program rdf implicit integer*4 (i-n) implicit real*8 (a-h,o-z) character*80 fn,s80,s81 character*1 ok logical lex,lex2,lex3,lex4,lcd,lmd,lex5,lall,lwater,lper,lsmooth, 1lex6,lsum common/ll/dmin,iq,lcd,lmd,ia common/c66/lex6,natu,nato real*8,allocatable::ro(:),ra(:),rh(:),rm(:),r(:),wso(:) integer*4,allocatable::nlist(:),iza(:),isu(:) c inquire(file='RDF.OPT',exist=lex) inquire(file='STRUCT.LST',exist=lex2) inquire(file='FILE.X',exist=lex3) inquire(file='TRACE.PRN',exist=lex4) inquire(file='SOLU.LST',exist=lex5) inquire(file='WEIGHT.LST',exist=lex6) c read nat only: if(lex2)then open(2,file='STRUCT.LST',status='old') read(2,200)fn 200 format(a80) close(2) open(4,file=fn) read(4,*)nat close(4) else if(lex3)then open(2,file='FILE.X',status='old') read(2,200)fn read(2,*)nat close(2) else call report('Geometry not found') endif endif write(6,*)nat,' atoms' allocate(r(3*nat),iza(nat),nlist(nat),isu(nat),wso(nat)) write(6,6000) 6000 format(//,' Radial distribution function calculation',/,/ 1' RDF.OPT ... option file',/, 1' STRUCT.LST ... list of MD files (Tinker)',/, 1' FILE.X ... list of geometries (alternative)',/, 1' SOLU.LST ... list of solute atoms ',/, 1' WEIGHT.LST ... list of solvent atomic weights ',/) ifa=0 px=0.0d0 py=0.0d0 pz=0.0d0 lall=.false. lsum=.false. lsmooth=.false. dmi=0.0d0 dm=10.0d0 nd=200 lwater=.true. lmd=.false. lcd=.false. lper=.false. nmax=0 nmin=0 if(lex)then open(9,file='RDF.OPT') 1000 read(9,900,end=99,err=99)s80 900 format(a80) c central atom: if(s80(1:4).eq.'ATOM')then read(9,*)nra if(nra.lt.0)then read(9,*)(nlist(i),i=1,-nra) write(6,*)'Atomic list read' endif endif c minimal distance: if(s80(1:4).eq.'RMIN')read(9,*)dmi c maximal distance: if(s80(1:4).eq.'RMAX')read(9,*)dm c number of poins: if(s80(1:2).eq.'NP')read(9,*)nd c water solvent: if(s80(1:5).eq.'WATER')read(9,*)lwater c follow atom: if(s80(1:6).eq.'FOLLOW')read(9,*)ifa c all solute atoms as central: if(s80(1:4).eq.'LALL')read(9,*)lall c one RDF as average of all: if(s80(1:4).eq.'LSUM')read(9,*)lsum c mass density: if(s80(1:4).eq.'MASS')read(9,*)lmd c charge density: if(s80(1:4).eq.'CHAR')read(9,*)lcd c smoothing: if(s80(1:4).eq.'SMOO')read(9,*)lsmooth c periodic boundary conditions: if(s80(1:4).eq.'LPER')read(9,*)lper c periodic box dimensions: if(s80(1:4).eq.'PERI')read(9,*)px,py,pz c minimal number of geometry taken if(s80(1:4).eq.'NMIN')read(9,*)nmin c maximal number of geometry taken if(s80(1:4).eq.'NMAX')read(9,*)nmax if(s80(1:3).eq.'END')goto 99 backspace 9 read(9,900)s81 if(s80.eq.s81)then write(6,900)s80 write(6,*)'unknown option' stop endif goto 1000 99 close(9) else write(6,*)' Select the atom number:' read(5,*)nra write(6,*) 1 ' Minimum and maximum distance (A) and number of intervals:' read(5,*)dmi,dm,nd write(6,*)' Periodic conditions (Y/N)?' read(5,'(a)')ok lper=ok.eq.'y'.or.ok.eq.'Y' if(lper)then write(6,*)' Periodic dimensions (X, Y, Z):' read(5,*)px,py,pz endif write(6,*)' Follow atom or h2o (#IZ/0)?' read(5,*)ifa endif if(lex5)then open(55,file='SOLU.LST') read(55,*)natu read(55,*)(isu(i),i=1,natu) close(55) write(6,*)natu,' atoms in SOLU.LST' endif if(lex6)then open(55,file='WEIGHT.LST') read(55,*)nato read(55,*)(wso(i),i=1,nato) close(55) write(6,*)nato,' Atomic solvent weights in WEIGHT.LST' endif dd=(dm-dmi)/dble(nd) write(6,6002)dd 6002 format(' Tick: ',f6.3,' A') c all solute atoms if(lall)then if(nra.lt.0)then call report('All solute and atom list not compatible options') else nstart=isu(1) nend=isu(natu) do 2001 i=nstart,nend 2001 nlist(i)=isu(i-nstart+1) endif else if(nra.gt.0)then nstart=nra nend=nra nlist(nstart)=nra else nstart=1 nend=-nra endif endif allocate(ra(nd),rh(nd),rm(nd),ro(nd)) call vz(ra,nd) call vz(rh,nd) call vz(rm,nd) call vz(ro,nd) write(6,60000)nra,dmi,dm,nd,ifa,lwater,lper,lall,lmd,lcd,lex5, 1lex2,lex3,lex6,lsmooth,nmax,nmin,px,py,pz 60000 format(' central atom',i10,/, 1 ' RMIN ',f10.3,' RMAX ',f10.5,/, 1 ' NPOINTS ',I10, ' FOLLOW ',I10,/, 1 ' WATER ',L10, ' LPER ',L10,/, 1 ' LALL ',L10, ' MASS ',L10,/, 1 ' CHARGE ',L10, ' SOLU.LST ',L10,/, 1 ' STRUCT.LST ',L10, ' FILE.X ',L10,/, 1 ' WEIGHT.LST ',L10, ' LSMOOTH ',L10,/, 1 ' NMAX ',I10, ' NMIN ',I10,/, 1 ' BOX ',3f10.3,/) write(6,*)nend-nstart+1,' atoms to watch:' write(6,60001)(nlist(i),i=nstart,nend) 60001 format(10i8) if(lsum)then write(6,*)'Average distribution requested' call wlist(nmin,nmax,nat,nlist,iza,nstart,nend,nd,isu,natu,ifa, 1 lsmooth,lwater,lex2,lex4,lex5,lex6,lper, 2 ra,rh,rm,ro,r,wso,px,py,pz,dd,dmi) else write(6,*)'separate distributions' call wsep(nmin,nmax,nat,nlist,iza,nstart,nend,nd,isu,natu,ifa, 1 lsmooth,lall,lwater,lex2,lex4,lex5,lex6,lper, 2 ra,rh,rm,ro,r,wso,px,py,pz,dd,dmi) endif end subroutine report(s) character*(*) s write(6,*)s write(6,*)'****** Program stopped. *********' stop end subroutine ppp(q1,q2,p,s) implicit real*8 (a-h,o-z) 1 if( s*(q1-q2).gt.1.5d0*p)then q1=q1-s*p*3.0d0 goto 1 endif return end subroutine cpc(xw,yw,zw,px,py,pz,xa,ya,za) c apply periodic box implicit real*8 (a-h,o-z) call ppp(xw,xa,px, 1.0d0) call ppp(yw,ya,py, 1.0d0) call ppp(zw,za,pz, 1.0d0) call ppp(xw,xa,px,-1.0d0) call ppp(yw,ya,py,-1.0d0) call ppp(zw,za,pz,-1.0d0) return end subroutine wrf(sn,nd,dmi,rh,dd) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) dimension rh(*) character*(*)sn open(3,file=sn) do 1 i=1,nd 1 write(3,4000)dmi+dd/2.0d0+dble(i-1)*dd,rh(i) 4000 format(f10.5,' ',g12.6) close(3) write(6,*)sn return end subroutine addx(nd,xw,yw,zw,xa,ya,za,dmi,dd,r,j,wso,isu) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) dimension r(*),amas(89),wso(*),isu(*) logical lcd,lmd common/ll/dmin,iq,lcd,lmd,ia data amas/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ logical lex6 common/c66/lex6,natu,nato pi=4.0d0*atan(1.0d0) d2=(xw-xa)**2+(yw-ya)**2+(zw-za)**2 dist=sqrt(d2) if(dist.lt.dmin)dmin=dist ni=int((dist-dmi)/dd)+1 if(ni.le.nd.and.ni.ge.1)then j=j+1 if(lcd)then r(ni)=r(ni)+dble(iq)/(4.0d0*pi*d2) else if(lmd)then c to have density in g/ml: r(ni)=r(ni)+dble(amas(iq))/(4.0d0*pi*d2*dd*6.022045d23*1.d-24) else if(lex6)then c multiply by atomic weight c suppose that solute atoms are isu(1)...isu(natu) c number of atoms in solute and solvent molecules :natu,nato c current atom ia, find its number in one solvent molecule: if(ia.lt.isu(1))then ns1=ia-(ia/nato)*nato else if(ia.gt.isu(natu))then ns1=ia-((ia-natu)/nato)*nato-natu else call report('cannot happen') endif endif r(ni)=r(ni)+1.0d0/(4.0d0*pi*d2)*wso(ns1) else r(ni)=r(ni)+1.0d0/(4.0d0*pi*d2) endif endif endif endif return end subroutine addxs(nd,xw,yw,zw,xa,ya,za,dmi,dd,r,j,wso,isu) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c version with smoothing dimension r(*),amas(89),isu(*),wso(*) logical lcd,lmd common/ll/dmin,iq,lcd,lmd,ia data amas/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ logical lex6 common/c66/lex6,natu,nato pi=4.0d0*atan(1.0d0) d2=(xw-xa)**2+(yw-ya)**2+(zw-za)**2 dist=sqrt(d2) if(dist.lt.dmin)dmin=dist ni=int((dist-dmi)/dd)+1 c c center of the grid cell: c=dmi+(dble(ni)-0.5d0)*dd v=dist-c if(v.gt.0)then np=ni+1 wp=v/dd else np=ni-1 wp=-v/dd endif w=1.0d0-wp if(ni.le.nd.and.ni.ge.1)then j=j+1 if(lcd)then r(ni)=r(ni)+dble(iq)/(4.0d0*pi*d2)*w else if(lmd)then c to have density in g/ml: r(ni)=r(ni)+dble(amas(iq)) 1 /(4.0d0*pi*d2*dd*6.022045d23*1.0d-24)*w else if(lex6)then c multiply by atomic weight c suppose that solute atoms are isu(1)...isu(natu) c number of atoms in solute and solvent molecules :natu,nato c current atom ia, find its number in one solvent molecule: if(ia.lt.isu(1))then ns1=ia-(ia/nato)*nato else if(ia.gt.isu(natu))then ns1=ia-((ia-natu)/nato)*nato-natu else call report('cannot happen') endif endif r(ni)=r(ni)+1.0d0/(4.0d0*pi*d2)*w*wso(ns1) else r(ni)=r(ni)+1.0d0/(4.0d0*pi*d2)*w endif endif endif endif if(np.le.nd.and.np.ge.1)then if(lcd)then r(np)=r(np)+dble(iq)/(4.0d0*pi*d2)*wp else if(lmd)then c to have density in g/ml: r(np)=r(np)+amas(iq)/(4.0d0*pi*d2*dd*6.022045d23*1.0d-24)*wp else if(lex6)then if(ia.lt.isu(1))then ns1=ia-(ia/nato)*nato else if(ia.gt.isu(natu))then ns1=ia-((ia-natu)/nato)*nato-natu else call report('cannot happen') endif endif r(np)=r(np)+1.0d0/(4.0d0*pi*d2)*wp*wso(ns1) else r(np)=r(np)+1.0d0/(4.0d0*pi*d2)*wp endif endif endif endif return end function symboltoatnumber(s3) implicit none character*3 s3,s integer symboltoatnumber,jstart,jend,iend,istart,i character*2 atsy(89) DATA ATSy/'H ','HE','LI','BE','B ','C ','N ','O ','F ','NE', 3'NA','MG','AL','SI','P ','S ','CL','AR', 4'K ','CA','SC','TI','V ','CR','MN','FE','CO','NI','CU','ZN', 4 'GA','GE','AS','SE','BR','KR', 5'RB','SR','Y ','ZR','NB','MO','TC','RU','RH','PD','AG','CD', 5 'IN','SN','SB','TE','I ','XE', 6'CS','BA','LA', 6 'CE','PR','ND','PM','SM','EU','GD','TB','DY','HO', 6 'ER','TM','YB','LU', 6'HF','TA','W ','RE','OS','IR','PT','AU','HG', 6 'TL','PB','BI','PO','AT','RN', 7'FR','RA','AC'/ s=s3 jstart=1 do 1 istart=1,len(s) 1 if(s(istart:istart).ne.' ')goto 2 2 do 3 iend=len(s),1,-1 3 if(s(iend:iend).ne.' ')goto 4 4 do 6 i=istart,iend 6 if(ichar(s(i:i)).gt.96)s(i:i)=char(ichar(s(i:i))-32) do 5 i=1,89 if(atsy(i)(2:2).eq.' ')then jend=1 else jend=2 endif if(atsy(i)(jstart:jend).eq.s(istart:iend))then symboltoatnumber=i return endif 5 continue write(6,*)s,' symbol not found' symboltoatnumber=0 return end subroutine vz(v,n) real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end subroutine dse(sn,is,ie) character*(*) sn integer*4 is,ie do 11 is=1,len(sn) 11 if(sn(is:is).ne.' ')goto 12 12 do 13 ie=len(sn),1,-1 13 if(sn(ie:ie).ne.' ')goto 14 14 return end subroutine wsep(nmin,nmax,nat,nlist,iza,nstart,nend,nd,isu,natu, 1ifa, 2lsmooth,lall,lwater,lex2,lex4,lex5,lex6,lper, 3ra,rh,rm,ro,r,wso,px,py,pz,dd,dmi) implicit none integer*4 nlist(*),nstart,nend,kk,ira,nd,isu(*),natu,nmc, 1nmct,nwmo,icvar,nat,nmin,ia,iza(*),symboltoatnumber,nwm, 1jo,jm,jh,ja,ix,iy,iz,ic,i,nmax,iq,jj,nh,ni,ifa,is,ie,is2, 1ie2,iend,istart,nh1,nh2,ir character*80 fn,sn,sn2,sn3 character*3 sy real*8 ra(*),rh(*),ro(*),rm(*),r(*),xa,ya,za,dmin,doh,x,y,z, 1xj,yj,zj,d2,xw,yw,zw,wso(*),xw2,yw2,zw2,xw1,yw1,zw1,xm,ym,zm, 2px,py,pz,sp,oa,om,dd,dmi logical lall,lex2,lex4,lex5,lex6,lper,lwater,lsmooth,lcd,lmd common/ll/dmin,iq,lcd,lmd,ia do 799 kk=nstart,nend ira=nlist(kk) if(lall)write(6,*)'looking at atom ',ira call vz(ra,nd) call vz(rh,nd) call vz(rm,nd) call vz(ro,nd) c open for averaging: if(lex2)then open(2,file='STRUCT.LST',status='old') write(6,*)'STRUCT.LST opened' else write(6,*)'FILE.X opened' open(2,file='FILE.X',status='old') endif write(6,*) c if file with minimal distances exists, rewrite the actual version if(lex4)open(44,file='TRACE.PRN') c check that solute atoms are together if(lex6)then if(isu(natu)-isu(1)+1.ne.natu)then call report('Forbidden ordering of solute atoms for lex6') endif endif nmc=0 nmct=0 nwmo=0 icvar=0 c if per. box, scan 3x3x3 superbox: if(lper)then ir=27 else ir=1 endif 1 continue if(lex2)then c Tinker files: read(2,200,end=500,err=500)fn 200 format(a80) open(4,file=fn) read(4,*)nat if(nmin.ne.0.or.nmct+1.lt.nmin)then c faster read for unanalyzed geoms: do 211 ia=1,nat 211 read(4,*) else c normal read: do 21 ia=1,nat read(4,400)sy,(r(3*(ia-1)+ix),ix=1,3) 400 format(7x,a3,1x,3f12.6) 21 iza(ia)=symboltoatnumber(sy) endif close(4) else c FILE.X: read(2,200,end=500,err=500)fn read(2,*,end=500,err=500)nat if(nmin.ne.0.and.nmct+1.lt.nmin)then c faster read for unanalyzed geoms: do 212 ia=1,nat 212 read(2,*) else c normal read: do 2 ia=1,nat 2 read(2,*)iza(ia),(r(3*(ia-1)+ix),ix=1,3) endif endif xa=r(3*(ira-1)+1) ya=r(3*(ira-1)+2) za=r(3*(ira-1)+3) nmct=nmct+1 dmin=1.0d9 c c maximum square of O-H distance: doh=1.85d0 c number of water molecules: nwm=0 c c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee if(nmin.eq.0.or.nmct.ge.nmin)then if(nmax.eq.0.or.nmct.le.nmax)then c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee jo=0 jm=0 jh=0 ja=0 nmc=nmc+1 c ix=-1 iy=-1 iz=-1 do 61 ic=1,ir do 6 ia=1,nat c exclude investigated atoms: do 612 i=nstart,nend 612 if(nlist(i).eq.ia)goto 6 c exclude solute atoms: if(lex5)then do 611 i=1,natu 611 if(ia.eq.isu(i))goto 6 endif iq=iza(ia) x=r(3*(ia-1)+1)+px*dble(ix) y=r(3*(ia-1)+2)+py*dble(iy) z=r(3*(ia-1)+3)+pz*dble(iz) c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww if(lwater)then c c decide, whether this is water O: if(iq.eq.8)then c ia is oxygen nh=0 do 62 jj=1,nat if(iza(jj).eq.1)then xj=r(3*(jj-1)+1)+px*dble(ix) yj=r(3*(jj-1)+2)+py*dble(iy) zj=r(3*(jj-1)+3)+pz*dble(iz) d2=(x-xj)**2+(y-yj)**2+(z-zj)**2 if(d2.lt.doh)nh=nh+1 if(d2.lt.doh.and.nh.eq.1)nh1=jj if(d2.lt.doh.and.nh.eq.2)nh2=jj endif 62 continue if(nh.eq.2)then c ia is bond to two hydrogens nh1 and nh2 nwm=nwm+1 c c RDF - O xw=x yw=y zw=z if(lper)call cpc(xw,yw,zw,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) else call addx(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) endif c c RDF - H xw1=r(3*(nh1-1)+1)+px*dble(ix) yw1=r(3*(nh1-1)+2)+py*dble(iy) zw1=r(3*(nh1-1)+3)+pz*dble(iz) if(lper)call cpc(xw1,yw1,zw1,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw1,yw1,zw1,xa,ya,za,dmi,dd,rh,jh,wso,isu) else call addx(nd,xw1,yw1,zw1,xa,ya,za,dmi,dd,rh,jh,wso,isu) endif xw2=r(3*(nh2-1)+1)+px*dble(ix) yw2=r(3*(nh2-1)+2)+py*dble(iy) zw2=r(3*(nh2-1)+3)+pz*dble(iz) if(lper)call cpc(xw2,yw2,zw2,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw2,yw2,zw2,xa,ya,za,dmi,dd,rh,jh,wso,isu) else call addx(nd,xw2,yw2,zw2,xa,ya,za,dmi,dd,rh,jh,wso,isu) endif c c RDF - mass center xm=(xw1+xw2+16.0d0*xw)/18.0d0 ym=(yw1+yw2+16.0d0*yw)/18.0d0 zm=(zw1+zw2+16.0d0*zw)/18.0d0 if(lper)call cpc(xm,ym,zm,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xm,ym,zm,xa,ya,za,dmi,dd,rm,jm,wso,isu) else call addx(nd,xm,ym,zm,xa,ya,za,dmi,dd,rm,jm,wso,isu) endif c c RDF - angle/water orientation c H c / c O ............ c .. \ c .. H c A sp= (xa-xw)*(xm-xw)+(ya-yw)*(ym-yw)+(za-zw)*(zm-zw) oa=sqrt((xa-xw)*(xa-xw)+(ya-yw)*(ya-yw)+(za-zw)*(za-zw)) om=sqrt((xm-xw)*(xm-xw)+(ym-yw)*(ym-yw)+(zm-zw)*(zm-zw)) d2=sp if(oa*om.gt.0.0d0)d2=acos(sp/oa/om)*180.0d0/4.0d0/atan(1.0d0) ni=int((nd*d2)/180.0d0)+1 if(ni.le.nd.and.ni.ge.1)then ra(ni)=ra(ni)+1.0d0 ja=ja+1 else write(6,*)d2,' angle outside limits' endif c endif c if(nh.eq.2) endif c if(iq.eq.8) c c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww else if(ifa.eq.iq.or.lcd.or.lmd.or.lex6)then c specific atom type nwm=nwm+1 xw=x yw=y zw=z if(lper)call cpc(xw,yw,zw,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) else call addx(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) endif endif endif c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww 6 continue c loop over atoms ix=ix+1 if(ix.gt.1)then ix=-1 iy=iy+1 if(iy.gt.1)then iy=-1 iz=iz+1 endif endif 61 continue c loop over periodic boxes if(nwm.ne.nwmo.and.nwmo.ne.0)icvar=icvar+1 nwmo=nwm if(lex4)write(44,440)nmc,dmin 440 format(i10,f12.5) c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee endif endif c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee c skip rest if limited number of geometries was required: if(nmct.gt.nmax.and.nmax.ne.0)goto 500 c goto 1 500 close(2) c if(icvar.gt.0)write(6,*)'Warning: Variable number of H2O', 1' molecules ',icvar,' times !' write(6,*)nmc,' geometries' write(6,*)nmct,' total geometries' write(6,*) write(6,*)'Last (super)box:' write(6,*)nwm,' waters' write(6,*)jo, ' oxygens' write(6,*)jh, ' hydrogens' write(6,*)jm, ' mass centers' write(6,*)ja, ' angles' c if(nmc.gt.1)then do 102 i=1,nd if(jh.gt.0)rh(i)=rh(i)/dble(nmc) if(ja.gt.0)ra(i)=ra(i)/dble(nmc) if(jm.gt.0)rm(i)=rm(i)/dble(nmc) 102 if(jo.gt.0)ro(i)=ro(i)/dble(nmc) endif write(sn,309)ira 309 format(i80) call dse(sn,is,ie) if(lwater)then call wrf(sn(is:ie)//'hrdf.prn',nd,dmi,rh,dd) call wrf(sn(is:ie)//'ordf.prn',nd,dmi,ro,dd) call wrf(sn(is:ie)//'mrdf.prn',nd,dmi,rm,dd) call wrf(sn(is:ie)//'ardf.prn',nd+1,-0.0d0,ra,180.0d0/dble(nd)) else sn3=' ' if(lall)write(sn3,*)ira call dse(sn3,istart,iend) if(iend.lt.istart)iend=istart if(lmd)then call wrf(sn3(istart:iend)//'md.prn',nd,dmi,ro,dd) else if(lcd)then call wrf(sn3(istart:iend)//'cd.prn',nd,dmi,ro,dd) else if(lex6)then call wrf(sn3(istart:iend)//'id.prn',nd,dmi,ro,dd) else write(sn2,309)ifa call dse(sn2,is2,ie2) call wrf( 1 sn(is:ie)//'_'//sn2(is2:ie2)//'rdf.prn', 2 nd,dmi,ro,dd) endif endif endif endif 799 continue if(lex4)close(44) return end subroutine wlist(nmin,nmax,nat,nlist,iza,nstart,nend,nd,isu,natu, 1ifa, 2lsmooth,lwater,lex2,lex4,lex5,lex6,lper, 3ra,rh,rm,ro,r,wso,px,py,pz,dd,dmi) implicit none integer*4 nlist(*),nstart,nend,kk,ira,nd,isu(*),natu,nmc, 1nmct,nwmo,icvar,nat,nmin,ia,iza(*),symboltoatnumber,nwm, 1jo,jm,jh,ja,ix,iy,iz,ic,i,nmax,iq,jj,nh,ni,ifa,is2, 1ie2,nh1,nh2,ir,nn character*80 fn,sn2 character*3 sy real*8 ra(*),rh(*),ro(*),rm(*),r(*),xa,ya,za,dmin,doh,x,y,z, 1xj,yj,zj,d2,xw,yw,zw,wso(*),xw2,yw2,zw2,xw1,yw1,zw1,xm,ym,zm, 2px,py,pz,sp,oa,om,dd,dmi logical lex2,lex4,lex5,lex6,lper,lwater,lsmooth,lcd,lmd common/ll/dmin,iq,lcd,lmd,ia call vz(ra,nd) call vz(rh,nd) call vz(rm,nd) call vz(ro,nd) nn=nend-nstart+1 c open for averaging: if(lex2)then open(2,file='STRUCT.LST',status='old') write(6,*)'STRUCT.LST opened' else write(6,*)'FILE.X opened' open(2,file='FILE.X',status='old') endif write(6,*) c if file with minimal distances exists, rewrite the actual version if(lex4)open(44,file='TRACE.PRN') c check that solute atoms are together if(lex6)then if(isu(natu)-isu(1)+1.ne.natu)then call report('Forbidden ordering of solute atoms for lex6') endif endif nmc=0 nmct=0 nwmo=0 icvar=0 c if per. box, scan 3x3x3 superbox: if(lper)then ir=27 else ir=1 endif 1 continue if(lex2)then c Tinker files: read(2,200,end=500,err=500)fn 200 format(a80) open(4,file=fn) read(4,*)nat if(nmin.ne.0.or.nmct+1.lt.nmin)then c faster read for unanalyzed geoms: do 211 ia=1,nat 211 read(4,*) else c normal read: do 21 ia=1,nat read(4,400)sy,(r(3*(ia-1)+ix),ix=1,3) 400 format(7x,a3,1x,3f12.6) 21 iza(ia)=symboltoatnumber(sy) endif close(4) else c FILE.X: read(2,200,end=500,err=500)fn read(2,*,end=500,err=500)nat if(nmin.ne.0.and.nmct+1.lt.nmin)then c faster read for unanalyzed geoms: do 212 ia=1,nat 212 read(2,*) else c normal read: do 2 ia=1,nat 2 read(2,*)iza(ia),(r(3*(ia-1)+ix),ix=1,3) endif endif nmct=nmct+1 dmin=1.0d9 c c maximum square of O-H distance: doh=1.85d0 c number of water molecules: nwm=0 c c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee if(nmin.eq.0.or.nmct.ge.nmin)then if(nmax.eq.0.or.nmct.le.nmax)then c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee jo=0 jm=0 jh=0 ja=0 nmc=nmc+1 c ix=-1 iy=-1 iz=-1 do 61 ic=1,ir do 6 ia=1,nat c exclude investigated atoms: do 612 i=nstart,nend 612 if(nlist(i).eq.ia)goto 6 c exclude solute atoms: if(lex5)then do 611 i=1,natu 611 if(ia.eq.isu(i))goto 6 endif iq=iza(ia) x=r(3*(ia-1)+1)+px*dble(ix) y=r(3*(ia-1)+2)+py*dble(iy) z=r(3*(ia-1)+3)+pz*dble(iz) c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww if(lwater)then c c decide, whether this is water O: if(iq.eq.8)then c ia is oxygen nh=0 do 62 jj=1,nat if(iza(jj).eq.1)then xj=r(3*(jj-1)+1)+px*dble(ix) yj=r(3*(jj-1)+2)+py*dble(iy) zj=r(3*(jj-1)+3)+pz*dble(iz) d2=(x-xj)**2+(y-yj)**2+(z-zj)**2 if(d2.lt.doh)nh=nh+1 if(d2.lt.doh.and.nh.eq.1)nh1=jj if(d2.lt.doh.and.nh.eq.2)nh2=jj endif 62 continue if(nh.eq.2)then c ia is bond to two hydrogens nh1 and nh2 nwm=nwm+1 c c RDF - O xw=x yw=y zw=z do 799 kk=nstart,nend ira=nlist(kk) xa=r(3*(ira-1)+1) ya=r(3*(ira-1)+2) za=r(3*(ira-1)+3) if(lper)call cpc(xw,yw,zw,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) else call addx(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) endif c c RDF - H xw1=r(3*(nh1-1)+1)+px*dble(ix) yw1=r(3*(nh1-1)+2)+py*dble(iy) zw1=r(3*(nh1-1)+3)+pz*dble(iz) if(lper)call cpc(xw1,yw1,zw1,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw1,yw1,zw1,xa,ya,za,dmi,dd,rh,jh,wso,isu) else call addx(nd,xw1,yw1,zw1,xa,ya,za,dmi,dd,rh,jh,wso,isu) endif xw2=r(3*(nh2-1)+1)+px*dble(ix) yw2=r(3*(nh2-1)+2)+py*dble(iy) zw2=r(3*(nh2-1)+3)+pz*dble(iz) if(lper)call cpc(xw2,yw2,zw2,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw2,yw2,zw2,xa,ya,za,dmi,dd,rh,jh,wso,isu) else call addx(nd,xw2,yw2,zw2,xa,ya,za,dmi,dd,rh,jh,wso,isu) endif c c RDF - mass center xm=(xw1+xw2+16.0d0*xw)/18.0d0 ym=(yw1+yw2+16.0d0*yw)/18.0d0 zm=(zw1+zw2+16.0d0*zw)/18.0d0 if(lper)call cpc(xm,ym,zm,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xm,ym,zm,xa,ya,za,dmi,dd,rm,jm,wso,isu) else call addx(nd,xm,ym,zm,xa,ya,za,dmi,dd,rm,jm,wso,isu) endif c c RDF - angle/water orientation c H c / c O ............ c .. \ c .. H c A sp= (xa-xw)*(xm-xw)+(ya-yw)*(ym-yw)+(za-zw)*(zm-zw) oa=sqrt((xa-xw)*(xa-xw)+(ya-yw)*(ya-yw)+(za-zw)*(za-zw)) om=sqrt((xm-xw)*(xm-xw)+(ym-yw)*(ym-yw)+(zm-zw)*(zm-zw)) d2=sp if(oa*om.gt.0.0d0)d2=acos(sp/oa/om)*180.0d0/4.0d0/atan(1.0d0) ni=int((nd*d2)/180.0d0)+1 if(ni.le.nd.and.ni.ge.1)then ra(ni)=ra(ni)+1.0d0 ja=ja+1 else write(6,*)d2,' angle outside limits' endif c 799 continue endif c if(nh.eq.2) endif c if(iq.eq.8) c c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww else if(ifa.eq.iq.or.lcd.or.lmd.or.lex6)then c specific atom type nwm=nwm+1 xw=x yw=y zw=z do 899 kk=nstart,nend ira=nlist(kk) xa=r(3*(ira-1)+1) ya=r(3*(ira-1)+2) za=r(3*(ira-1)+3) if(lper)call cpc(xw,yw,zw,px,py,pz,xa,ya,za) if(lsmooth)then call addxs(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) else call addx(nd,xw,yw,zw,xa,ya,za,dmi,dd,ro,jo,wso,isu) endif 899 continue endif endif c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww 6 continue c loop over atoms ix=ix+1 if(ix.gt.1)then ix=-1 iy=iy+1 if(iy.gt.1)then iy=-1 iz=iz+1 endif endif 61 continue c loop over periodic boxes if(nwm.ne.nwmo.and.nwmo.ne.0)icvar=icvar+1 nwmo=nwm if(lex4)write(44,440)nmc,dmin 440 format(i10,f12.5) c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee endif endif c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee c skip rest if limited number of geometries was required: if(nmct.gt.nmax.and.nmax.ne.0)goto 500 c goto 1 500 close(2) c if(icvar.gt.0)write(6,*)'Warning: Variable number of H2O', 1' molecules ',icvar,' times !' write(6,*)nmc,' geometries' write(6,*)nmct,' total geometries' write(6,*) write(6,*)'Last (super)box:' write(6,*)nwm,' waters' write(6,*)jo, ' oxygens' write(6,*)jh, ' hydrogens' write(6,*)jm, ' mass centers' write(6,*)ja, ' angles' do 103 i=1,nd rh(i)=rh(i)/dble(nn) ra(i)=ra(i)/dble(nn) rm(i)=rm(i)/dble(nn) 103 ro(i)=ro(i)/dble(nn) c if(nmc.gt.1)then do 102 i=1,nd if(jh.gt.0)rh(i)=rh(i)/dble(nmc) if(ja.gt.0)ra(i)=ra(i)/dble(nmc) if(jm.gt.0)rm(i)=rm(i)/dble(nmc) 102 if(jo.gt.0)ro(i)=ro(i)/dble(nmc) endif if(lwater)then call wrf('hrdf.prn',nd,dmi,rh,dd) call wrf('ordf.prn',nd,dmi,ro,dd) call wrf('mrdf.prn',nd,dmi,rm,dd) call wrf('ardf.prn',nd+1,-0.0d0,ra,180.0d0/dble(nd)) else if(lmd)then call wrf('md.prn',nd,dmi,ro,dd) else if(lcd)then call wrf('cd.prn',nd,dmi,ro,dd) else if(lex6)then call wrf('id.prn',nd,dmi,ro,dd) else write(sn2,309)ifa 309 format(i80) call dse(sn2,is2,ie2) call wrf(sn2(is2:ie2)//'rdf.prn',nd,dmi,ro,dd) endif endif endif endif if(lex4)close(44) return end