c c ---------------------------------------------------------------------- c c in core map extension - Kevin Cowtan c C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C c ---------------------------------------------------------------------- c program ncsmask include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c integer size1,size2,sizes(2) character types(2)*1 c integer lunsto,idum external ncsmaskmain,lunsto,XYZINIT c data types /'R','R'/ c call ccpfyp call ccprcs (6, 'ncsmask', '$Date: 2002/04/12 09:01:49 $') call XYZINIT c idum=0 lmap = 1 lpt=lunsto(idum) c call rcards(size1,size2) c sizes(1)=size1 sizes(2)=size2 call ccpalc(ncsmaskmain,2,types,sizes) c call ccperr(0,' (ncsmask) - normal termination ') end c c ---------------------------------------------------------------------- c subroutine ncsmaskmain(nmap,map,nunit,unit) include 'ncsmask.fh' include 'params.fh' c mask array and optional unit cell array integer nmap,nunit real map(nmap),unit(nunit) c integer i character*240 title external trim c title=' ' c do 100 i=1,nmap map(i)=0.0 100 continue do 110 i=1,nunit unit(i)=0.0 110 continue c c read mask if (mskin.eq.1) call ccpmskin(map,title) c or coords if (xyzin.eq.1) call pdbmask(map) c peak counting if (numpks.gt.0) call peak(map) c expand/contract if (exprad.ne.0.0) call expand(map) c smooth if (smtrad.ne.0.0) call smooth(map) c isolate monomer if (monomer.ne.0.0) call makemono(map) c overlap removal if (overlap.ne.0) call over(map,unit) c work out limits of set region if (edgtrim.ne.0) call trim(map) c output call ccpmskout(map,title) c and comment c call comment(map) c return end c c ---------------------------------------------------------------------- c subroutine rcards(size1,size2) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c rcards - read input cards c integer size1,size2 c character namspg*10,nampg*10,namlau*10,title*80, + line*80,xyz*4,caxis(10)*10 c real iuvw(2,3),ouvw(2,3),symmat(4,4) real r,omega,phi,kappa,alpha,beta,gamma real xspace,yspace,zspace integer i,j,k,l,icsym logical lover,lntrm,lmono,linvm,leule,lpola,ltran,lomat,cont, + lspacing c c spacegroups with yxz=1 and zxy=2 axis ordering integer axis(230) data axis/2,2,2,2,1,1,1,1,1,2,1,1,1,1,1,2,2,2,1,2,2,1,2,207*1/ c data lover,lntrm,lmono,linvm/4*.false./ c xyz='XYZ' c c set defaults c iuvw(1,1)=-1.0e6 ouvw(1,1)=-1.0e6 cell(1)=0.0 cell(2)=0.0 cell(3)=0.0 cell(4)=90.0 cell(5)=90.0 cell(6)=90.0 caxis(1)='?' caxis(2)='?' caxis(3)='?' nspgrp=1 namspg=' ' nu=0 nv=0 nw=0 ncsym=0 icsym=1 monomer=0 overlap=0 edgtrim=1 numpks=-1 atmrad=3.0 exprad=0.0 smtrad=0.0 xspace=0.0 yspace=0.0 zspace=0.0 c do 910 i=1,4 do 900 j=1,4 r=0.0 if (i.eq.j) r=1.0 do 880 k=1,maxncs rsymncs(i,j,k)=r 880 continue do 890 k=1,96 rsym(i,j,k)=r 890 continue 900 continue 910 continue c c now parse input c 1000 call memoparse(.true.) c leule=.false. lpola=.false. ltran=.false. lomat=.false. call parsesubkey('ROTA','EULE',leule) call parsesubkey('ROTA','POLA',lpola) call parsesubkey('TRAN',' ',ltran) call parsesubkey('OMAT',' ',lomat) c c CELL call parsecell(cell) c SYMM call parsesymm(namspg,nspgrp,nampg,nsym,nsymp,rsym) c GRID lspacing = .false. call parsesubkey('GRID','SPAC',lspacing) if (.not. lspacing) then call parsesubint('GRID',' ',1,.true.,nu) call parsesubint('GRID',' ',2,.true.,nv) call parsesubint('GRID',' ',3,.true.,nw) else call parsesubreal('GRID','SPAC',1,.true.,xspace) call parsesubreal('GRID','SPAC',2,.false.,yspace) if (yspace.gt.0.0) + call parsesubreal('GRID','SPAC',3,.true.,zspace) end if c XYZLIM call parsesubreal('XYZL',' ',1,.true.,ouvw(1,1)) call parsesubreal('XYZL',' ',2,.true.,ouvw(2,1)) call parsesubreal('XYZL',' ',3,.true.,ouvw(1,2)) call parsesubreal('XYZL',' ',4,.true.,ouvw(2,2)) call parsesubreal('XYZL',' ',5,.true.,ouvw(1,3)) call parsesubreal('XYZL',' ',6,.true.,ouvw(2,3)) c AXIS call parsesubchar('AXIS',' ',1,.true.,caxis(1)) call parsesubchar('AXIS',' ',2,.true.,caxis(2)) call parsesubchar('AXIS',' ',3,.true.,caxis(3)) c OVER call parsesubkey('OVER',' ',lover) call parsesubint('OVER',' ',1,.false.,overlap) c PEAK call parseint('PEAK',numpks) c RADIUS call parsereal('RADI',atmrad) c EXPAND call parsereal('EXPA',exprad) c SMOOTH call parsereal('SMOO',smtrad) c NOTRIM call parsekey('NOTR',lntrm) c MONOMER call parsekey('MONO',lmono) C AVER call parseint('AVER',ncsym) call parsesubkey('AVER','INVE',linvm) c ROTA call parsesubreal('ROTA','EULE',1,.true.,alpha) call parsesubreal('ROTA','EULE',2,.true.,beta) call parsesubreal('ROTA','EULE',3,.true.,gamma) call parsesubreal('ROTA','POLA',1,.true.,omega) call parsesubreal('ROTA','POLA',2,.true.,phi) call parsesubreal('ROTA','POLA',3,.true.,kappa) call parsesubreal('ROTA','MATR',1,.true.,rsymncs(1,1,icsym)) call parsesubreal('ROTA','MATR',2,.true.,rsymncs(1,2,icsym)) call parsesubreal('ROTA','MATR',3,.true.,rsymncs(1,3,icsym)) call parsesubreal('ROTA','MATR',4,.true.,rsymncs(2,1,icsym)) call parsesubreal('ROTA','MATR',5,.true.,rsymncs(2,2,icsym)) call parsesubreal('ROTA','MATR',6,.true.,rsymncs(2,3,icsym)) call parsesubreal('ROTA','MATR',7,.true.,rsymncs(3,1,icsym)) call parsesubreal('ROTA','MATR',8,.true.,rsymncs(3,2,icsym)) call parsesubreal('ROTA','MATR',9,.true.,rsymncs(3,3,icsym)) c TRAN call parsesubreal('TRAN',' ',1,.true.,rsymncs(1,4,icsym)) call parsesubreal('TRAN',' ',2,.true.,rsymncs(2,4,icsym)) call parsesubreal('TRAN',' ',3,.true.,rsymncs(3,4,icsym)) c deal with special input keywords if (lomat) read (*,*)((rsymncs(i,j,icsym),i=1,3),j=1,4) if (leule) call eulertomatrix(alpha,beta,gamma,rsymncs(1,1,icsym)) if (lpola) call polartomatrix(omega,phi ,kappa,rsymncs(1,1,icsym)) if (ltran.or.lomat) icsym=icsym+1 c call parsediagnose(cont) if (cont) goto 1000 c find out where the input is coming form c problems may occur if xyzin or mskin has been previously assigned mskin=0 xyzin=0 call ugtenv('MSKIN',line) if (line.ne.' ') mskin=1 call ugtenv('XYZIN',line) if (line.ne.' ') xyzin=1 c put in a check if (mskin.eq.1 .and. xyzin.eq.1) then call ccperr(1, 'ERR: XYZIN and MSKIN both assigned') elseif (mskin.eq.1) then write (lpt, *) 'INFO: MSKIN assigned' elseif (xyzin.eq.1) then write (lpt, *) 'INFO: XYZIN assigned' else call ccperr(1, 'ERR: XYZIN or MSKIN must be assigned') endif c c Now process the input c c c try and get cell/symm from file if (mskin.ne.0) then call ccpmskhead(iuvw) jfmso(1)=jfmsi(1) jfmso(2)=jfmsi(2) jfmso(3)=jfmsi(3) elseif (xyzin.ne.0) then call pdbhead(iuvw) else goto 9060 endif c call calcmetric() c if (ncsym.ne.icsym-1) goto 9030 ncsym=max(ncsym,1) if (lover) overlap=max(overlap,1) if (lntrm) edgtrim=0 if (lmono) monomer=1 do 290 k=1,ncsym if (linvm) then call symcpy(rsymncs(1,1,k),symmat) call invsym(symmat,rsymncs(1,1,k)) endif 290 continue c c spacegroup info call msymlb(lmap,nspgrp,namspg,nampg,nsymp,nsym,rsym) c axis order do 150 i=1,3 jfmso(i)=index(xyz,caxis(i)(1:1)) 150 continue if (jfmso(1).le.0) then if (axis(mod(nspgrp,1000)).eq.1) then jfmso(1)=2 jfmso(2)=1 jfmso(3)=3 else jfmso(1)=3 jfmso(2)=1 jfmso(3)=2 endif endif c set number of grid points if (xspace.gt.0.0) then if (yspace*zspace.eq.0.0) then yspace=xspace zspace=xspace elseif (yspace.lt.0.0 .or. zspace.lt.0.0) then goto 9070 endif nu=nint(cell(1)/xspace) nv=nint(cell(2)/yspace) nw=nint(cell(3)/zspace) if (nu*nv*nw.eq.0) then write (lpt, *) 'INFO: grid spacing too large, defaulting to 1A' endif elseif (xspace.lt.0.0) then goto 9070 endif c if (nu*nv*nw.eq.0) then nu=nint(cell(1)) nv=nint(cell(2)) nw=nint(cell(3)) endif c if (xyzin.ne.0) then iuvw(1,1)=(iuvw(1,1)-(atmrad+2.0)/q(1,1))*nu iuvw(2,1)=(iuvw(2,1)+(atmrad+2.0)/q(1,1))*nu iuvw(1,2)=(iuvw(1,2)-(atmrad+2.0)/q(2,2))*nv iuvw(2,2)=(iuvw(2,2)+(atmrad+2.0)/q(2,2))*nv iuvw(1,3)=(iuvw(1,3)-(atmrad+2.0)/q(3,3))*nw iuvw(2,3)=(iuvw(2,3)+(atmrad+2.0)/q(3,3))*nw endif if (exprad.gt.0.0) then iuvw(1,1)=iuvw(1,1)-exprad/q(1,1)*nu iuvw(2,1)=iuvw(2,1)+exprad/q(1,1)*nu iuvw(1,2)=iuvw(1,2)-exprad/q(2,2)*nv iuvw(2,2)=iuvw(2,2)+exprad/q(2,2)*nv iuvw(1,3)=iuvw(1,3)-exprad/q(3,3)*nw iuvw(2,3)=iuvw(2,3)+exprad/q(3,3)*nw endif if (ouvw(1,1).lt.-0.9e6) then ouvw(1,1)=iuvw(1,1) ouvw(2,1)=iuvw(2,1) ouvw(1,2)=iuvw(1,2) ouvw(2,2)=iuvw(2,2) ouvw(1,3)=iuvw(1,3) ouvw(2,3)=iuvw(2,3) endif if ((ouvw(2,1)-ouvw(1,1))*(ouvw(2,2)-ouvw(1,2))* + (ouvw(2,3)-ouvw(1,3)).lt.8.5) then ouvw(1,1)=ouvw(1,1)*nu ouvw(2,1)=ouvw(2,1)*nu ouvw(1,2)=ouvw(1,2)*nv ouvw(2,2)=ouvw(2,2)*nv ouvw(1,3)=ouvw(1,3)*nw ouvw(2,3)=ouvw(2,3)*nw endif c iuvwlim11=int(iuvw(1,1)) iuvwlim21=int(iuvw(2,1)) iuvwlim12=int(iuvw(1,2)) iuvwlim22=int(iuvw(2,2)) iuvwlim13=int(iuvw(1,3)) iuvwlim23=int(iuvw(2,3)) ouvwlim11=int(ouvw(1,1)) ouvwlim21=int(ouvw(2,1)) ouvwlim12=int(ouvw(1,2)) ouvwlim22=int(ouvw(2,2)) ouvwlim13=int(ouvw(1,3)) ouvwlim23=int(ouvw(2,3)) c lu=iuvwlim21-iuvwlim11+1 lv=iuvwlim22-iuvwlim12+1 lw=iuvwlim23-iuvwlim13+1 c size1=lu*lv*lw size2=256 if (lover) size2=nu*nv*nw c return 9010 call ccperr(1,' rcards - error in AXIS card') 9020 call ccperr(1,' rcards - error in FROM card') 9030 call ccperr(1,' rcards - error in AVER card') 9040 call ccperr(1,' rcards - conflict between card and file input') 9060 call ccperr(1,' rcards - no input file on MSKIN or XYZIN') 9070 call ccperr(1,' rcards - GRID SPACING must be greater than zero') end c c ---------------------------------------------------------------------- c subroutine ccpmskhead(uvwlim) include 'ncsmask.fh' include 'io.fh' include 'params.fh' include 'crystal.fh' c CCPMSKhead - get limits of mask from file c real uvwlim(2,3) c integer mxyz(3),juvwi(3),m1(3),m2(3),msec,mode integer i,ifail,iprint real rmin,rmax,rmean,rrms character title*80 c ifail=0 iprint=1 call mrdhds(lmap,'MSKIN',title,msec,jfmsi,mxyz,m1(3),m1(1),m2(1), + m1(2),m2(2),cell,nspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint) m2(3)=m1(3)+msec-1 c do 110 i=1,3 juvwi(jfmsi(i))=i 110 continue c nu=mxyz(1) nv=mxyz(2) nw=mxyz(3) uvwlim(1,1)=m1(juvwi(1)) uvwlim(2,1)=m2(juvwi(1)) uvwlim(1,2)=m1(juvwi(2)) uvwlim(2,2)=m2(juvwi(2)) uvwlim(1,3)=m1(juvwi(3)) uvwlim(2,3)=m2(juvwi(3)) c return end c c ---------------------------------------------------------------------- c subroutine ccpmskin(mask,title) include 'ncsmask.fh' include 'io.fh' include 'params.fh' include 'crystal.fh' c CCPMSKIN - read a ccp4 logical mask and store in xyz order c real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) character title*(*) real lsec(0:maxsec) c integer ifast,imedm,islow,lfast,lmedm,lslow,ierr,ifail,iprint integer m1(3),m2(3),ifms(3),juvwi(3),mxyz(3) integer i,m,iu,iv,iw,dum(3) integer nsec,mode,nsymexp,lmod,mspgrp real rmin,rmax,rmean,rrms,u,v,w,rho c external lmod c c Note: juvw convert from fast/med/slow to u/v/w c jfms convert from u/v/w to fast/med/slow c c now open map header and read map, re-ordering it as necessary c ifail=0 iprint=0 call mrdhds(lmap,'MSKIN',title,nsec,dum,mxyz,m1(3),m1(1),m2(1), + m1(2),m2(2),cell,mspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint) m2(3)=m1(3)+nsec-1 do 110 i=1,3 juvwi(jfmsi(i))=i 110 continue c find out the mask grid dimensions lfast=m2(1)-m1(1)+1 lmedm=m2(2)-m1(2)+1 lslow=m2(3)-m1(3)+1 c if (lfast*lmedm.gt.maxsec) + call ccperr(1,' ccpmskin - Mask section > lsec: recompile') c now read the map in the order it is on file do 200 islow=0,lslow-1 c get a section ierr=0 call mgulpr(lmap,lsec,ierr) if (ierr.ne.0) call ccperr(1,' ccpmskin - ccp4 read error') c now sort into uvw mask ifms(3)=islow+m1(3) do 190 imedm=0,lmedm-1 ifms(2)=imedm+m1(2) do 190 ifast=0,lfast-1 ifms(1)=ifast+m1(1) iu=ifms(juvwi(1)) iv=ifms(juvwi(2)) iw=ifms(juvwi(3)) if (iu.ge.iuvwlim11.and.iu.le.iuvwlim21.and. + iv.ge.iuvwlim12.and.iv.le.iuvwlim22.and. + iw.ge.iuvwlim13.and.iw.le.iuvwlim23) then mask(iu,iv,iw)=lsec(ifast+lfast*imedm) endif 190 continue 200 continue c close the map file call mrclos(lmap) c write (lpt,910) 910 format (/' MASK READ SUCCESSFUL'//) c return end c c ---------------------------------------------------------------------- c subroutine ccpmskout(mask,title) include 'ncsmask.fh' include 'io.fh' include 'params.fh' include 'crystal.fh' c CCPMSKOUT - write a ccp4 logical mask and store in xyz order c real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) character title*(*) c real lsec(0:maxsec) c integer ifast,imedm,islow,lfast,lmedm,lslow,ierr,ifail,iprint integer m1(3),m2(3),ifms(3),juvwo(3),mxyz(3) integer i,m,iu,iv,iw,si,sn integer nsec,mode,nsymexp,lmod real rmin,rmax,rmean,rrms,u,v,w,u1,v1,w1,rho c external lmod c c Note: juvw convert from fast/med/slow to u/v/w c jfms convert from u/v/w to fast/med/slow c c now open map header and read map, re-ordering it as necessary c do 110 i=1,3 juvwo(jfmso(i))=i 110 continue c mxyz(1)=nu mxyz(2)=nv mxyz(3)=nw m1(juvwo(1))=ouvwlim11 m2(juvwo(1))=ouvwlim21 m1(juvwo(2))=ouvwlim12 m2(juvwo(2))=ouvwlim22 m1(juvwo(3))=ouvwlim13 m2(juvwo(3))=ouvwlim23 nsec=m2(3)-m1(3)+1 mode=0 c call mwrhdl(lmap,'MSKOUT',title,nsec,jfmso,mxyz,m1(3),m1(1),m2(1), + m1(2),m2(2),cell,nspgrp,mode) call msywrt(lmap,nsym,rsym) c c find out the mask grid dimensions lfast=m2(1)-m1(1)+1 lmedm=m2(2)-m1(2)+1 lslow=m2(3)-m1(3)+1 c if (lfast*lmedm.gt.maxsec) + call ccperr(1,' ccpmapin - Mask section > lsec: recompile') c now write the map in the new order sn=0 si=0 do 200 islow=0,lslow-1 c sort onto section ifms(3)=islow+m1(3) do 190 imedm=0,lmedm-1 ifms(2)=imedm+m1(2) do 190 ifast=0,lfast-1 ifms(1)=ifast+m1(1) iu=ifms(juvwo(1)) iv=ifms(juvwo(2)) iw=ifms(juvwo(3)) i=0 if (iu.ge.iuvwlim11.and.iv.ge.iuvwlim12.and.iw.ge.iuvwlim13.and. + iu.le.iuvwlim21.and.iv.le.iuvwlim22.and.iw.le.iuvwlim23) then if (mask(iu,iv,iw).gt.0.5) i=1 endif call ccpstb(i,lsec,ifast+lfast*imedm+1) sn=sn+1 si=si+i 190 continue c write the section call mspew(lmap,lsec) 200 continue c close the map file call mwclose(lmap) c write (lpt,550)sn,(100.0*sn)/(nu*nv*nw),si,(100.0*si)/(nu*nv*nw), + (100.0*si)/sn 550 format (' OUTPUT MASK STATISTICS:',/, + ' Mask contains ',i9,' gridpoints, or ',f6.1,' % of the cell',/, + 6x,'of which ',i9,' are set filling ',f6.1,' % of the cell',/, + 38x,'or ',f6.1,' % of the mask',/) c if (sn.gt.nu*nv*nw) write (lpt,560) 560 format ( + ' Warning: The mask box covers more than one unit cell',/) c return end c c ---------------------------------------------------------------------- c subroutine pdbhead(uvwlim) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c real b(6),uvwlim(2,3) c integer i,ifail,idum,iresn,iser,ixyzin,iz,nread,wwarn real biso,occ,x,y,z,u,v,w,umin,vmin,wmin,umax,vmax,wmax,border character atnam*4,resnam*4,resno*4,id*4,chnnam*1,altcod*1, + inscod*1,segid*4 c c open coordinate file ifail=0 IXYZIN = 0 call XYZOPEN('XYZIN','INPUT',' ',IXYZIN,ifail) if (ifail.ne.0) goto 9000 c call rbcell(cell,vol) if (cell(1).eq.0.0) goto 9010 call calcmetric() c uvwlim(1,1)=1.0e6 uvwlim(1,2)=1.0e6 uvwlim(1,3)=1.0e6 uvwlim(2,1)=-1.0e6 uvwlim(2,2)=-1.0e6 uvwlim(2,3)=-1.0e6 nread=0 wwarn=0 c Read atomic coordinates, convert to fractional 300 call XYZADVANCE(IXYZIN,0,0,*500,*500) call XYZATOM(IXYZIN,iser,atnam,resnam,chnnam,iresn,resno, + inscod,altcod,segid,iz,id) call XYZCOORD(IXYZIN,'O','U',x,y,z,occ,biso,b) c c Skip ANISOU cards c if (b(2).ne.0.0 .or. b(3).ne.0.0) then if (x.eq.0.0 .and. y.eq.0.0 .and. z.eq.0.0) goto 300 endif nread=nread+1 call orth2frac(x,y,z,u,v,w) uvwlim(1,1)=min(u,uvwlim(1,1)) uvwlim(1,2)=min(v,uvwlim(1,2)) uvwlim(1,3)=min(w,uvwlim(1,3)) uvwlim(2,1)=max(u,uvwlim(2,1)) uvwlim(2,2)=max(v,uvwlim(2,2)) uvwlim(2,3)=max(w,uvwlim(2,3)) if (resnam(1:3).eq.'HOH' .or. resnam(1:3).eq.'DOD' .or. + resnam(1:3).eq.'H20' .or. resnam(1:3).eq.'WAT') wwarn=1 goto 300 c done 500 call XYZCLOSE(IXYZIN) write (lpt,510)nread 510 format (' Number of atoms read: ',i7) c if (wwarn.eq.1) write (lpt,520) 520 format ( + /,' WARNING: this PDB file appears to contain solvent atoms.' + /,' DO you really want to use these in your mask?',/) c return 9000 call ccperr(1,' pdbhead - XYZIN error') 9010 call ccperr(1,' pdbhead - XYZIN has no CRYSTL card') end c c ---------------------------------------------------------------------- c subroutine pdbmask(mask) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c PDBMASK - make a protein mask from coordinates in a .pdb file c real b(6),mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c real biso,occ,x,y,z,u,v,w,du,dv,dw,ratom,r2cut,r2 integer i,ixyzin,j,k,i1,j1,k1,ru,rv,rw,nread,nused integer ifail,idum c r2cut=atmrad**2 ru=int(atmrad*nu/q(1,1))+1 rv=int(atmrad*nv/q(2,2))+1 rw=int(atmrad*nw/q(3,3))+1 c c open coordinate file ifail=0 IXYZIN = 0 call XYZOPEN('XYZIN','INPUT',' ',IXYZIN,ifail) if (ifail.ne.0) goto 9000 c nread=0 nused=0 c Read atomic coordinates, convert to fractional 300 call XYZADVANCE(IXYZIN,0,0,*500,*500) call XYZCOORD(IXYZIN,'O','U',x,y,z,occ,biso,b) c c Skip ANISOU cards c if (b(2).ne.0.0 .or. b(3).ne.0.0) then if (x.eq.0.0 .and. y.eq.0.0 .and. z.eq.0.0) goto 300 endif nread=nread+1 call orth2frac(x,y,z,u,v,w) c fill in mask within 2A of atom i=nint(u*nu) j=nint(v*nv) k=nint(w*nw) do 400 k1=max(k-rw,iuvwlim13),min(k+rw,iuvwlim23) do 400 j1=max(j-rv,iuvwlim12),min(j+rv,iuvwlim22) do 400 i1=max(i-ru,iuvwlim11),min(i+ru,iuvwlim21) du=real(i1)/nu-u dv=real(j1)/nv-v dw=real(k1)/nw-w r2=rlm11*du**2+rlm22*dv**2+rlm33*dw**2+ + rlm12*du*dv+rlm13*du*dw+rlm23*dv*dw if (r2.lt.r2cut) mask(i1,j1,k1)=1.0 400 continue c next atom nused=nused+1 go to 300 c c end of .pdb file c 500 write (lpt,510)nread,nused 510 format (' Number of atoms read: ',i6/ + ' Number of atoms used: ',i6//) c end of .pdb file call XYZCLOSE(IXYZIN) return c 9000 call ccperr(1,' (PDBMASK) Error opening file XYZIN') end c c ---------------------------------------------------------------------- c subroutine smooth(mask) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c SMOOTH - smooth the mask real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c real u,v,w,dist,wt integer maxnn parameter (maxnn=11*11*11) integer nn,du(maxnn),dv(maxnn),dw(maxnn),hist(0:maxnn) integer i,j,n,m,dn,dm,iu,iv,iw,ju,jv,jw c c do 50 i=0,maxnn hist(i)=0 50 continue c c make list of neighbours orthogonal pixel distances to nearby pixels c dn=int(smtrad*(nu+nv+nw)/(cell(1)+cell(2)+cell(3))+1.0) dn=min(dn,5) nn=0 do 100 jw=-dn,dn w=real(jw)/nw do 100 jv=-dn,dn v=real(jv)/nv do 100 ju=-dn,dn u=real(ju)/nu dist=sqrt(rlm11*u**2+rlm22*v**2+rlm33*w**2+ + rlm12*u*v +rlm13*u*w +rlm23*v*w ) if (dist.lt.smtrad) then nn=nn+1 du(nn)=ju dv(nn)=jv dw(nn)=jw endif 100 continue cxx write (*,*)dn,nn c if (nn.eq.0) then write (lpt,110) 110 format (' smooth - smoothing radius too small') return endif c c now smooth the mask by testing for each point if more than half are set c n=0 m=0 dn=0 dm=0 c c first set all points inside mask to 0.5, outside to -0.5 do 150 iw=iuvwlim13,iuvwlim23 do 150 iv=iuvwlim12,iuvwlim22 do 150 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.0.5) then mask(iu,iv,iw)=0.5 n=n+1 else mask(iu,iv,iw)=-0.5 endif 150 continue c c now replace each point by the count of set neighbours, c +ve if inside, -ve if outside do 200 iw=iuvwlim13,iuvwlim23 do 200 iv=iuvwlim12,iuvwlim22 do 200 iu=iuvwlim11,iuvwlim21 c only treat points inside mask (to speed things up) wt=0.0 do 190 j=1,nn ju=iu+du(j) jv=iv+dv(j) jw=iw+dw(j) if (ju.ge.iuvwlim11.and.ju.le.iuvwlim21.and. + jv.ge.iuvwlim12.and.jv.le.iuvwlim22.and. + jw.ge.iuvwlim13.and.jw.le.iuvwlim23) then if (mask(ju,jv,jw).gt.0.0) wt=wt+1.0 endif 190 continue mask(iu,iv,iw)=sign(wt,mask(iu,iv,iw)) 200 continue c c now make a histogram of values do 250 iw=iuvwlim13,iuvwlim23 do 250 iv=iuvwlim12,iuvwlim22 do 250 iu=iuvwlim11,iuvwlim21 i=nint(abs(mask(iu,iv,iw))) hist(i)=hist(i)+1 250 continue c c search the histogram do 300 i=maxnn-1,0,-1 hist(i)=hist(i)+hist(i+1) 300 continue j=0 do 350 i=1,maxnn if (abs(hist(i)-n).lt.abs(hist(j)-n)) j=i 350 continue cxx write (*,*)j c c now mark the map and calculate the changes do 450 iw=iuvwlim13,iuvwlim23 do 450 iv=iuvwlim12,iuvwlim22 do 450 iu=iuvwlim11,iuvwlim21 if (nint(abs(mask(iu,iv,iw))).ge.j) then if (mask(iu,iv,iw).le.0.0) dm=dm+1 mask(iu,iv,iw)=1.0 m=m+1 else if (mask(iu,iv,iw).gt.0.0) dn=dn+1 mask(iu,iv,iw)=0.0 endif 450 continue c write (lpt,550)n,m,dn,dm 550 format(' MASK SMOOTHING:'/ + ' Points in mask before smoothing : ',i9/ + ' Points in mask after smoothing : ',i9/ + ' Points removed from mask during smoothing : ',i9/ + ' Points added to mask during smoothing : ',i9/) c return end c c ---------------------------------------------------------------------- c subroutine expand(mask) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c EXPAND - expand/contract the mask real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c real step,sn,sm integer iu,iv,iw c c work out how may pixels to expand the mask step=exprad*(nu+nv+nw)/(cell(1)+cell(2)+cell(3)) c c set the mask range -0.5..0.5 c sn=0.0 do 100 iw=iuvwlim13,iuvwlim23 do 100 iv=iuvwlim12,iuvwlim22 do 100 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.0.5) then sn=sn+1.0 mask(iu,iv,iw)=0.5 else mask(iu,iv,iw)=-0.5 endif 100 continue c c mark up the mask with distance from boundary, c +ve values inside and -ve outside c call markdist(mask) c c now build the expanded mask c sm=0.0 do 300 iw=iuvwlim13,iuvwlim23 do 300 iv=iuvwlim12,iuvwlim22 do 300 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.-step) then sm=sm+1.0 mask(iu,iv,iw)=1.0 else mask(iu,iv,iw)=0.0 endif 300 continue c write (lpt,910)nint(sn),100.0*sn/(nu*nv*nw), + nint(sm),100.0*sm/(nu*nv*nw) 910 format (' MASK EXPAND:' + /' Before expand mask contained ',i8, + ' points or ',f6.1,' % of unit cell' + /' After expand mask contained ',i8, + ' points or ',f6.1,' % of unit cell'/) return end c c ---------------------------------------------------------------------- c subroutine peak(mask) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c integer maxpeak parameter (maxpeak=10000) integer i,j,k,iu,iv,iw,ju,jv,jw integer npeaks,ipeak(0:maxpeak),npeak(0:maxpeak),index(0:maxpeak) c do 10 i=0,maxpeak ipeak(i)=0 npeak(i)=1 index(i)=i 10 continue c k=1 do 100 iw=iuvwlim13,iuvwlim23 do 100 iv=iuvwlim12,iuvwlim22 do 100 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.0.5) then mask(iu,iv,iw)=k k=k+1 endif 100 continue c do 300 i=1,500 k=0 do 200 iw=iuvwlim23,iuvwlim13,-1 do 200 iv=iuvwlim22,iuvwlim12,-1 do 200 iu=iuvwlim21,iuvwlim11,-1 if (mask(iu,iv,iw).gt.0.5) then j=0 do 150 jw=max(iw-1,iuvwlim13),min(iw+1,iuvwlim23) do 150 jv=max(iv-1,iuvwlim12),min(iv+1,iuvwlim22) do 150 ju=max(iu-1,iuvwlim11),min(iu+1,iuvwlim21) j=max(j,nint(mask(ju,jv,jw))) 150 continue if (j.gt.nint(mask(iu,iv,iw))) then mask(iu,iv,iw)=j k=k+1 endif endif 200 continue if (k.eq.0) goto 400 300 continue call ccperr(1,' peak - unable to peak-search mask') c 400 npeaks=0 do 500 iw=iuvwlim13,iuvwlim23 do 500 iv=iuvwlim12,iuvwlim22 do 500 iu=iuvwlim11,iuvwlim21 i=nint(mask(iu,iv,iw)) do 450 k=0,npeaks if (i.eq.ipeak(k)) then npeak(k)=npeak(k)+1 goto 500 endif 450 continue npeaks=npeaks+1 if (npeaks.gt.maxpeak) call ccperr + (1,' peak - mask is too complex for peak counting') ipeak(npeaks)=i 500 continue write (lpt,510)npeaks 510 format (' Found',i6,' peaks') c do 600 i=1,npeaks-1 do 600 j=i+1,npeaks if (npeak(index(j)).gt.npeak(index(i))) then k=index(j) index(j)=index(i) index(i)=k endif 600 continue c write (lpt,710)npeak(0) do 700 i=1,min(npeaks,max(2*numpks,10)) write (lpt,720)i,npeak(index(i)) 700 continue 710 format (' Points outside mask ',i8) 720 format (' Peak ',i6,' size ',i8) c do 800 iw=iuvwlim13,iuvwlim23 do 800 iv=iuvwlim12,iuvwlim22 do 800 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.0.5) then j=nint(mask(iu,iv,iw)) mask(iu,iv,iw)=1.0 do 750 i=1,numpks if (j.eq.ipeak(index(i))) goto 800 750 continue mask(iu,iv,iw)=0.0 endif 800 continue c return end c c ---------------------------------------------------------------------- c subroutine over(mask,unit) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c TRIM - trim the mask to the set volume c real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) real unit(0:nu-1,0:nv-1,0:nw-1) c integer nsymtot,nsymmax parameter (nsymmax=720) real mat(4,4,nsymmax),minv(4,4),temp(4,4) integer nuvwlim(6,nsymmax) c integer i,j,k,iu,iv,iw,ju,jv,jw,ku,kv,kw real u,v,w,u1,v1,w1,su,sv,sw,sn,sm c integer lmod real rfracp external rfracp,lmod c c start nsymtot=nsym*ncsym if (nsymtot.gt.nsymmax) goto 9000 c c now convert ncs ops to fractional coords do 10 i=1,ncsym call symmatmul(rsymncs(1,1,i),q,temp) call symmatmul(qt,temp,rsymncs(1,1,i)) 10 continue c c make all the combined symmetry operators do 21 j=1,ncsym do 20 i=1,nsym k=i+(j-1)*nsym call symmatmul(rsym(1,1,i),rsymncs(1,1,j),mat(1,1,k)) nuvwlim(1,k)=1000000 nuvwlim(2,k)=1000000 nuvwlim(3,k)=1000000 nuvwlim(4,k)=-1000000 nuvwlim(5,k)=-1000000 nuvwlim(6,k)=-1000000 20 continue 21 continue c c first find the bounding box of each molecule c do 200 i=1,nsymtot do 100 iw=iuvwlim13,iuvwlim23 w=real(iw)/nw do 100 iv=iuvwlim12,iuvwlim22 v=real(iv)/nv do 100 iu=iuvwlim11,iuvwlim21 u=real(iu)/nu if (mask(iu,iv,iw).gt.0.5) then call symuvw(mat(1,1,i),u,v,w,u1,v1,w1) ju=nint(u1*nu) jv=nint(v1*nv) jw=nint(w1*nw) nuvwlim(1,i)=min(nuvwlim(1,i),ju) nuvwlim(2,i)=min(nuvwlim(2,i),jv) nuvwlim(3,i)=min(nuvwlim(3,i),jw) nuvwlim(4,i)=max(nuvwlim(4,i),ju) nuvwlim(5,i)=max(nuvwlim(5,i),jv) nuvwlim(6,i)=max(nuvwlim(6,i),jw) endif 100 continue 200 continue c c NOW START THE OVERLAP REMOVAL LOOP c do 800 k=1,overlap c c clear unit cell mask, set the input mask range -0.5..0.5 c do 250 iw=0,nw-1 do 250 iv=0,nv-1 do 250 iu=0,nu-1 unit(iu,iv,iw)=-0.5 250 continue c do 300 iw=iuvwlim13,iuvwlim23 do 300 iv=iuvwlim12,iuvwlim22 do 300 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.0.5) then mask(iu,iv,iw)=0.5 else mask(iu,iv,iw)=-0.5 endif 300 continue c c mark up the mask with distance from boundary, c +ve values inside and -ve outside c call markdist(mask) c c now search map, generating equivalents c do 500 i=1,nsymtot call invsym(mat(1,1,i),minv) c loop through each point in this mask domain do 400 iw=nuvwlim(3,i),nuvwlim(6,i) w=real(iw)/nw kw=lmod(iw,nw) do 400 iv=nuvwlim(2,i),nuvwlim(5,i) v=real(iv)/nv kv=lmod(iv,nv) do 400 iu=nuvwlim(1,i),nuvwlim(4,i) u=real(iu)/nu ku=lmod(iu,nu) c rotate back to mask frame call symuvw(minv,u,v,w,u1,v1,w1) ju=nint(nu*u1) jv=nint(nv*v1) jw=nint(nw*w1) if (ju.ge.iuvwlim11.and.ju.le.iuvwlim21.and. + jv.ge.iuvwlim12.and.jv.le.iuvwlim22.and. + jw.ge.iuvwlim13.and.jw.le.iuvwlim23) then c mark this cell position with the corresponding mask value if (mask(ju,jv,jw).gt.unit(ku,kv,kw)) then unit(ku,kv,kw)=mask(ju,jv,jw) else if (mask(ju,jv,jw).eq.unit(ku,kv,kw)) then unit(ku,kv,kw)=mask(ju,jv,jw)+0.1 endif endif 400 continue 500 continue c c now go back to the mask, take out any overlap c and restore to 0..+1 range c sn=0.0 sm=0.0 do 600 iw=iuvwlim13,iuvwlim23 w=real(iw)/nw kw=lmod(iw,nw) do 600 iv=iuvwlim12,iuvwlim22 v=real(iv)/nv kv=lmod(iv,nv) do 600 iu=iuvwlim11,iuvwlim21 u=real(iu)/nu ku=lmod(iu,nu) if (mask(iu,iv,iw).gt.0.0) sn=sn+1.0 if (unit(ku,kv,kw).gt.mask(iu,iv,iw)) mask(iu,iv,iw)=-0.5 if (mask(iu,iv,iw).gt.0.0) sm=sm+1.0 if (mask(iu,iv,iw).gt.0.0) then mask(iu,iv,iw)=1.0 else mask(iu,iv,iw)=0.0 endif 600 continue c write (lpt,710)nsymtot,nint(sn),100.0*sn/(nu*nv*nw), + nint(sm),100.0*sm/(nu*nv*nw) 710 format (' OVERLAP REMOVAL:' + /' Total number of copies of mask in unit cell: ',i3 + /' Before overlap mask contained ',i8, + ' points or ',f6.1,' % of unit cell' + /' After overlap mask contained ',i8, + ' points or ',f6.1,' % of unit cell'/) c c NOW GO BACK FOR NEXT CYCLE c (why does it take more than one?) c 800 continue c return 9000 call ccperr(1,'over - You have a ridiculous number of symm ops!') end c c ---------------------------------------------------------------------- c subroutine markdist(mask) include 'ncsmask.fh' include 'crystal.fh' include 'params.fh' c markdist - mark up the mask, marking each point with the (approximate) c distance to the nearest mask boundary. Max error ~ 15% (cubic) c Distances are in approximate grid points. real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c real olddif,newdif,maxht,minht,uphill,dnhill,dist,cut real u,v,w,distab(-1:1,-1:1,-1:1) integer iu,iv,iw,ju,jv,jw,ku,kv,kw c c make table of orthogonal pixel distances to nearby pixels do 100 jw=-1,1 w=jw/cell(3) do 100 jv=-1,1 v=jv/cell(2) do 100 ju=-1,1 u=ju/cell(1) distab(ju,jv,jw)=sqrt( + rlm11*u**2+rlm22*v**2+rlm33*w**2+ + rlm12*u*v +rlm13*u*w +rlm23*v*w ) 100 continue distab(0,0,0)=1.0e12 cxx write (*,'(3(3f8.3/)/)')distab c c and mark up - iterative refinement cut=0.0 newdif=0.0 150 olddif=newdif+0.001 maxht=0.0 minht=0.0 c loop to increase stuff inside mask do 200 iw=iuvwlim13,iuvwlim23 do 200 iv=iuvwlim12,iuvwlim22 do 200 iu=iuvwlim11,iuvwlim21 c only treat points inside mask (to speed things up) if (abs(mask(iu,iv,iw)).gt.cut) then uphill=cut+1.5 dnhill=-cut-1.5 do 180 jw=-1,1 do 180 jv=-1,1 do 180 ju=-1,1 ku=iu+ju kv=iv+jv kw=iw+jw dist=distab(ju,jv,jw) c mark with highest/lowest neighbour+distance if (ku.ge.iuvwlim11.and.ku.le.iuvwlim21.and. + kv.ge.iuvwlim12.and.kv.le.iuvwlim22.and. + kw.ge.iuvwlim13.and.kw.le.iuvwlim23) then uphill=min(uphill,mask(ku,kv,kw)+dist) dnhill=max(dnhill,mask(ku,kv,kw)-dist) else c edges of box count as outside mask boundaries as well uphill=min(uphill,-0.5+dist) endif 180 continue if (mask(iu,iv,iw).gt.0.0) then mask(iu,iv,iw)=uphill else mask(iu,iv,iw)=dnhill endif endif maxht=max(maxht,mask(iu,iv,iw)) minht=min(minht,mask(iu,iv,iw)) 200 continue c we can increase the cutoff to speed things up cut=cut+1.0 newdif=maxht-minht cxx write (6,*)' maxht,minht ',maxht,minht if (newdif.gt.olddif) goto 150 return end c c ---------------------------------------------------------------------- c subroutine trim(mask) include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' include 'params.fh' c TRIM - trim the mask to the set volume c real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c integer i,iu,iv,iw c ouvwlim11=iuvwlim21 ouvwlim12=iuvwlim22 ouvwlim13=iuvwlim23 ouvwlim21=iuvwlim11 ouvwlim22=iuvwlim12 ouvwlim23=iuvwlim13 c do 100 iw=iuvwlim13,iuvwlim23 do 100 iv=iuvwlim12,iuvwlim22 do 100 iu=iuvwlim11,iuvwlim21 if (mask(iu,iv,iw).gt.0.5) then ouvwlim11=min(ouvwlim11,iu) ouvwlim21=max(ouvwlim21,iu) ouvwlim12=min(ouvwlim12,iv) ouvwlim22=max(ouvwlim22,iv) ouvwlim13=min(ouvwlim13,iw) ouvwlim23=max(ouvwlim23,iw) endif 100 continue c write (lpt,150) + iuvwlim11,iuvwlim21,iuvwlim12,iuvwlim22,iuvwlim13,iuvwlim23, + ouvwlim11,ouvwlim21,ouvwlim12,ouvwlim22,ouvwlim13,ouvwlim23 150 format(' MASK TRIM:' + /' Mask dimensions before trim:' + /' U ',i4,' :',i4/' V ',i4,' :',i4/' W ',i4,' :',i4 + /' After trim:' + /' U ',i4,' :',i4/' V ',i4,' :',i4/' W ',i4,' :',i4/) c return end c c ---------------------------------------------------------------------- c subroutine makemono(mask) include 'ncsmask.fh' include 'crystal.fh' include 'params.fh' include 'io.fh' c makemono - make a monomer from a multimer for simple proper sym real mask(iuvwlim11:iuvwlim21,iuvwlim12:iuvwlim22, + iuvwlim13:iuvwlim23) c integer i,j,k,iu,iv,iw real org(3),axis(3),axis2(3),r(3),s(3),dum(3) real u,v,w,x,y,z,xr,xs,sn,sm,sumstd,phi real a11,a12,a13,a21,a22,a23,a31,a32,a33,twopi real rmod c external rmod c c twopi=2.0*pi sumstd=0.0 do 50 i=1,3 org(i)=0.0 axis(i)=0.0 axis2(i)=0.0 50 continue c c first check the syms and calc the symmetry origin and axis c .. origin do 100 k=1,ncsym call symuvw(rsymncs(1,1,k),0.0,0.0,0.0,x,y,z) org(1)=org(1)+x org(2)=org(2)+y org(3)=org(3)+z 100 continue c .. axis do 200 k=2,ncsym a11=rsymncs(1,1,k)-1.0 a12=rsymncs(1,2,k) a13=rsymncs(1,3,k) a21=rsymncs(2,1,k) a22=rsymncs(2,2,k)-1.0 a23=rsymncs(2,3,k) a31=rsymncs(3,1,k) a32=rsymncs(3,2,k) a33=rsymncs(3,3,k)-1.0 if (abs(a12*a21-a11*a22).gt.0.1) then dum(1)=-(a12*a23-a13*a22)/(a12*a21-a11*a22) dum(2)=-(a13*a21-a11*a23)/(a12*a21-a11*a22) dum(3)=1.0 else if (abs(a23*a32-a22*a33).gt.0.1) then dum(1)=1.0 dum(2)=-(a23*a31-a21*a33)/(a23*a32-a22*a33) dum(3)=-(a21*a32-a22*a31)/(a23*a32-a22*a33) else call ccperr(1,' monomer - get Kevin to write some more code') endif do 190 i=1,3 axis(i)=axis(i)+dum(i) axis2(i)=axis2(i)+dum(i)**2 190 continue 200 continue c do 210 i=1,3 org(i)=org(i)/ncsym axis(i)=axis(i)/(ncsym-1) axis2(i)=axis2(i)/(ncsym-1) sumstd=sumstd+sqrt(max(axis2(i)-axis(i)**2,0.0)) 210 continue cxx write (*,*)org,axis c if (sumstd.gt.0.05) call ccperr + (1,' monomer - ncs ops do not appear to be a single proper set') c c now make a vector dum not parallel to axis (fudge for now) dum(1)=1.0 dum(2)=2.0 dum(3)=3.0 c make r perpendicular to dum and axis call crossprod(axis,dum,r) c make s perpendicular to r and axis call crossprod(axis,r,s) c c now search the map removing points outside angle 0...2pi/nncs sn=0.0 sm=0.0 do 400 iw=iuvwlim13,iuvwlim23 w=real(iw)/nw do 400 iv=iuvwlim12,iuvwlim22 v=real(iv)/nv do 400 iu=iuvwlim11,iuvwlim21 u=real(iu)/nu call frac2orth(u,v,w,x,y,z) x=x-org(1) y=y-org(2) z=z-org(3) xr=x*r(1)+y*r(2)+z*r(3) xs=x*s(1)+y*s(2)+z*s(3) phi=rmod(atan2(xs,xr),twopi) if (mask(iu,iv,iw).gt.0.5) sn=sn+1.0 if (phi.gt.twopi/ncsym) mask(iu,iv,iw)=0.0 if (mask(iu,iv,iw).gt.0.5) sm=sm+1.0 400 continue c write (lpt,910)ncsym,nint(sn),100.0*sn/(nu*nv*nw), + nint(sm),100.0*sm/(nu*nv*nw) 910 format (' MONOMER GENERATION:' + /' Number of non-crystallographic syms: ',i3 + /' Before monomer-ing mask contained ',i8, + ' points or ',f6.1,' % of unit cell' + /' After monomer-ing mask contained ',i8, + ' points or ',f6.1,' % of unit cell'/) c return end c c ---------------------------------------------------------------------- c subroutine crossprod(a,b,c) include 'ncsmask.fh' real a(3),b(3),c(3) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end c c ---------------------------------------------------------------------- c subroutine frac2orth(u,v,w,x,y,z) include 'ncsmask.fh' include 'crystal.fh' real u,v,w,x,y,z x=q(1,1)*u+q(1,2)*v+q(1,3)*w+q(1,4) y=q(2,1)*u+q(2,2)*v+q(2,3)*w+q(2,4) z=q(3,1)*u+q(3,2)*v+q(3,3)*w+q(3,4) return end c c ---------------------------------------------------------------------- c subroutine orth2frac(x,y,z,u,v,w) include 'ncsmask.fh' include 'crystal.fh' real u,v,w,x,y,z u=qt(1,1)*x+qt(1,2)*y+qt(1,3)*z+qt(1,4) v=qt(2,1)*x+qt(2,2)*y+qt(2,3)*z+qt(2,4) w=qt(3,1)*x+qt(3,2)*y+qt(3,3)*z+qt(3,4) return end c c ---------------------------------------------------------------------- c subroutine calcmetric() include 'ncsmask.fh' include 'io.fh' include 'crystal.fh' c Calculate the orthonormal - crystallographic transform c integer i,j real a,b,c,alph,beta,gamm,sum,cosas,q123 c do 10 j=1,4 do 10 i=1,4 q(i,j)=0.0 qt(i,j)=0.0 10 continue c a=cell(1) b=cell(2) c=cell(3) alph=dtor*cell(4) beta=dtor*cell(5) gamm=dtor*cell(6) sum=0.5*(alph+beta+gamm) vol=2.0*a*b*c* + sqrt(sin(sum)*sin(sum-alph)*sin(sum-beta)*sin(sum-gamm)) cosas=(cos(beta)*cos(gamm)-cos(alph))/(sin(beta)*sin(gamm)) c orthogonalisation matrix q(1,1)=a q(1,2)=b*cos(gamm) q(1,3)=c*cos(beta) q(2,2)=b*sin(gamm) q(2,3)=-c*sin(beta)*cosas q(3,3)=c*sin(beta)*sqrt(1.0-cosas**2) q(4,4)=1.0 c and its inverse q123=q(1,1)*q(2,2)*q(3,3) qt(1,1)=1.0/q(1,1) qt(1,2)=-q(1,2)*q(3,3)/q123 qt(1,3)=(q(1,2)*q(2,3)-q(2,2)*q(1,3))/q123 qt(2,2)=1.0/q(2,2) qt(2,3)=-q(1,1)*q(2,3)/q123 qt(3,3)=1.0/q(3,3) qt(4,4)=1.0 c write (lpt,4)vol 4 format (' Unit cell volume:',f10.1) write (lpt,1)((q(i,j),j=1,3),(qt(i,j),j=1,3),i=1,3) 1 format (' Orthogonalisation matrices:',3(/1x,3f8.3,8x,3f8.3)) c Calculate real space metric tensor rlm11=q(1,1)*q(1,1)+q(2,1)*q(2,1)+q(3,1)*q(3,1) rlm12=2.0*(q(1,1)*q(1,2)+q(2,1)*q(2,2)+q(3,1)*q(3,2)) rlm13=2.0*(q(1,1)*q(1,3)+q(2,1)*q(2,3)+q(3,1)*q(3,3)) rlm22=q(1,2)*q(1,2)+q(2,2)*q(2,2)+q(3,2)*q(3,2) rlm23=2.0*(q(1,2)*q(1,3)+q(2,2)*q(2,3)+q(3,2)*q(3,3)) rlm33=q(1,3)*q(1,3)+q(2,3)*q(2,3)+q(3,3)*q(3,3) write (lpt,2)rlm11,rlm22,rlm33,rlm12,rlm13,rlm23 2 format (' Real space metric tensor: '/1x,6f9.1) c Calculate reciprocal space metric tensor rcm11=qt(1,1)*qt(1,1)+qt(1,2)*qt(1,2)+qt(1,3)*qt(1,3) rcm12=2.0*(qt(1,1)*qt(2,1)+qt(1,2)*qt(2,2)+qt(1,3)*qt(2,3)) rcm13=2.0*(qt(1,1)*qt(3,1)+qt(1,2)*qt(3,2)+qt(1,3)*qt(3,3)) rcm22=qt(2,1)*qt(2,1)+qt(2,2)*qt(2,2)+qt(2,3)*qt(2,3) rcm23=2.0*(qt(2,1)*qt(3,1)+qt(2,2)*qt(3,2)+qt(2,3)*qt(3,3)) rcm33=qt(3,1)*qt(3,1)+qt(3,2)*qt(3,2)+qt(3,3)*qt(3,3) write (lpt,3)rcm11,rcm22,rcm33,rcm12,rcm13,rcm23 3 format (' Reciprocal space metric tensor:'/1x,6f9.6) c return end c c ---------------------------------------------------------------------- c subroutine eulertomatrix(alpha,beta,gamma,matrix) include 'ncsmask.fh' c calculate rotation matrix from Euler angles real alpha,beta,gamma,matrix(4,4) c real ca,cb,cg,sa,sb,sg c ca = cos(dtor*alpha) sa = sin(dtor*alpha) cb = cos(dtor*beta) sb = sin(dtor*beta) cg = cos(dtor*gamma) sg = sin(dtor*gamma) matrix(1,1) = cg*cb*ca - sg*sa matrix(1,2) = -sg*cb*ca - cg*sa matrix(1,3) = sb*ca matrix(2,1) = cg*cb*sa + sg*ca matrix(2,2) = -sg*cb*sa + cg*ca matrix(2,3) = sb*sa matrix(3,1) = -cg*sb matrix(3,2) = sg*sb matrix(3,3) = cb return end c c ---------------------------------------------------------------------- c subroutine polartomatrix(omega,phi,kappa,matrix) include 'ncsmask.fh' c calculate rotation matrix from Euler angles real omega,phi,kappa,matrix(4,4) c real someg,comeg,sinph,cosph,sinchi,coschi,cosch1,dc1,dc2,dc3 c someg = sin(omega*dtor) comeg = cos(omega*dtor) sinph = sin(phi*dtor) cosph = cos(phi*dtor) dc1 = someg*cosph dc2 = someg*sinph dc3 = comeg sinchi = sin(kappa*dtor) coschi = cos(kappa*dtor) cosch1 = 1.0 - coschi matrix(1,1) = dc1*dc1*cosch1 + coschi matrix(1,2) = dc1*dc2*cosch1 - dc3*sinchi matrix(1,3) = dc3*dc1*cosch1 + dc2*sinchi matrix(2,1) = dc1*dc2*cosch1 + dc3*sinchi matrix(2,2) = dc2*dc2*cosch1 + coschi matrix(2,3) = -dc1*sinchi + dc2*dc3*cosch1 matrix(3,1) = -dc2*sinchi + dc3*dc1*cosch1 matrix(3,2) = dc2*dc3*cosch1 + dc1*sinchi matrix(3,3) = dc3*dc3*cosch1 + coschi return end c c ---------------------------------------------------------------------- c subroutine symmatmul(a,b,c) include 'ncsmask.fh' c symmetmaul - 4x4 matrix multiply, c=a*b real a(4,4),b(4,4),c(4,4) integer i,j,k c do 20 k=1,4 do 20 i=1,4 c(i,k)=0.0 do 20 j=1,4 c(i,k)=c(i,k)+a(i,j)*b(j,k) 20 continue c return end c c ---------------------------------------------------------------------- c subroutine symuvw(mat,u,v,w,u1,v1,w1) include 'ncsmask.fh' c symuvw - multiply a 4-matrix by a 3-vector real mat(4,4),u,v,w,u1,v1,w1 u1=mat(1,1)*u+mat(1,2)*v+mat(1,3)*w+mat(1,4) v1=mat(2,1)*u+mat(2,2)*v+mat(2,3)*w+mat(2,4) w1=mat(3,1)*u+mat(3,2)*v+mat(3,3)*w+mat(3,4) return end c c ---------------------------------------------------------------------- c subroutine symcpy(a,b) include 'ncsmask.fh' c symcpy - copy a symmetry matrix real a(16),b(16) integer i do 100 i=1,16 b(i)=a(i) 100 continue return end c c ---------------------------------------------------------------------- c real function rfracp(x) include 'ncsmask.fh' c fractional part -0.5:0.5 real x rfracp=x-real(nint(x)) return end c c ---------------------------------------------------------------------- c real function rmod(r,s) include 'ncsmask.fh' c fixed real mod real r,s,t c t=mod(r,s) if (t.ge.0.0) then rmod=t else rmod=t+s endif return end c c ---------------------------------------------------------------------- c integer function lmod(i,j) include 'ncsmask.fh' c fixed integer mod (optimised for 0<=i