program rdfa implicit none integer*4 nd,jh,ig,i,natu,ir,ii,jj,kk,ia,ib,nat,ni,it,na,nb, 1s1,s2,ic real*8 px,py,pz,dmi,dm,dd,xa,ya,za,xb,yb,tol,pi,d2,dist,so logical lper,lsol,ls1,ls2 character*80 fn real*8,allocatable::rh(:),r(:) integer*4,allocatable::iza(:),isu(:),bl(:),al(:) c pi=4.0d0*atan(1.0d0) tol=1.0d-4 so=1.0d9 write(6,600) 600 format(/,' RDFA .. radial distribution function',/,/, 1 ' g(R_A1..A2)',/,/, 1 ' Input: FILE.X geometry',/, 1 ' RDFA.OPT options',/, 1 ' Output: g.prn RDF',/,/, 1' Options variable/default:',/, 1' -------------------------',/, 1' RMIN dmi/0 minimal distance / A',/, 1' RMAX dm /10 maximal distance / A',/, 1' NP nd/200 number of points',/, 1' SOLU lsol/f define solute atoms in file SOLU.LST',/, 1' 1SOL ls1/t include solute in first atom',/, 1' 2SOL ls2/t include solute in second atom',/, 1' LPER lper/f perodic boundary conditions',/, 1' PERI px py pz / 0 0 0 periodic dimensions',/, 1' FIRS na/0,n1,.... number of first atoms and their list',/, 1' SECO nb/0,n1,.... number of second atoms and their list',/, 1' END end of options',/,/, 1' By default, makes RDF between all pairs of atoms.',/, 1' If na <>0, first atom is limited:',/, 1' For na>0, only listed atomic numbers (one for hydrogen,', 1' etc.) are taken.',/, 1' For na<0, only listed atoms are taken', 1' (from 1..nat in FILE.X).',/, 1' (Similarly for nb and the second atom.)',/) open(2,file='FILE.X',status='old') read(2,*) read(2,*)nat write(6,*)nat,' atoms' close(2) allocate(r(3*nat),iza(nat),al(nat),bl(nat)) call readopt(px,py,pz,dmi,dm,nd,lper,lsol,dd, 1na,nb,al,bl,nat,ls1,ls2) allocate(rh(nd)) if(lsol)then open(55,file='SOLU.LST') read(55,*)natu allocate(isu(natu)) read(55,*)(isu(i),i=1,natu) close(55) write(6,655)natu,ls1,ls2 655 format(i6,' atoms in SOLU.LST, in first, second:',2l2) else write(6,656) 656 format(' No solute') endif rh=0.0d0 jh=0 ig=0 if(lper)then ir=1 else ir=0 endif open(2,file='FILE.X') 1 continue read(2,200,end=500,err=500)fn 200 format(a80) read(2,*,end=500,err=500)nat ig=ig+1 do 2 ia=1,nat 2 read(2,*)iza(ia),(r(3*(ia-1)+i),i=1,3) do 101 ia=1,nat if(lsol)then ic=0 do 108 i=1,natu 108 if(isu(i).eq.ia)ic=1 if(ls1.and.ic.eq.1)goto 1021 if(.not.ls1.and.ic.eq.0)goto 1021 goto 101 endif if(na.ne.0)then if(na.gt.0)then do 102 i=1,na 102 if(al(i).eq.iza(ia))goto 1021 goto 101 else do 107 i=1,-na 107 if(al(i).eq.ia)goto 1021 goto 101 endif endif 1021 xa=r(3*(ia-1)+1) ya=r(3*(ia-1)+2) za=r(3*(ia-1)+3) do 104 ib=1,nat if(lsol)then ic=0 do 109 i=1,natu 109 if(isu(i).eq.ib)ic=1 if(ls2.and.ic.eq.1)goto 1031 if(.not.ls2.and.ic.eq.0)goto 1031 goto 104 endif if(nb.ne.0)then if(nb.gt.0)then do 103 i=1,nb 103 if(bl(i).eq.iza(ib))goto 1031 goto 104 else do 106 i=1,-nb 106 if(bl(i).eq.ib)goto 1031 goto 104 endif endif 1031 it=3*(ib-1) do 105 ii=-ir,ir xb=(r(it+1)+px*dble(ii)-xa)**2 do 105 jj=-ir,ir yb=(r(it+2)+py*dble(jj)-ya)**2+xb do 105 kk=-ir,ir d2=(r(it+3)+py*dble(kk)-za)**2+yb dist=dsqrt(d2) if(dist.lt.so)then so=dist s1=ia s2=ib endif ni=int((dist-dmi)/dd)+1 if(ni.le.nd.and.ni.ge.1)then jh=jh+1 rh(ni)=rh(ni)+1.0d0/(4.0d0*pi*d2) endif 105 continue 104 continue 101 continue goto 1 500 close(2) c if(jh.gt.0)then rh=rh/dble(jh) call wrf('g.prn',nd,dmi,rh,dd) endif write(6,6001)ig,' geometries' write(6,6001)jh,' pairs ' 6001 format(i14,a11) write(6,6002)s1,iza(s1),s2,iza(s2),so 6002 format(' Shortest distance between ',i5,' (',i2, 1') -',i6,' (',i2,'):',f10.3,' A') 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 readopt(px,py,pz,dmi,dm,nd,lper,lsol,dd, 1na,nb,al,bl,nat,ls1,ls2) implicit none real*8 px,py,pz,dmi,dm,dd integer*4 nd,nat,al(nat),bl(nat),na,nb,i logical lper,lex,lsol,ls1,ls2 character*80 s80,s81 al=0 bl=0 na=0 nb=0 px=0.0d0 py=0.0d0 pz=0.0d0 dmi=0.0d0 dm=10.0d0 nd=200 lper=.false. lsol=.false. ls1=.true. ls2=.true. inquire(file='RDFA.OPT',exist=lex) if(lex)then open(9,file='RDFA.OPT') 1000 read(9,900,end=99,err=99)s80 900 format(a80) c solute defined if(s80(1:4).eq.'SOLU')read(9,*)lsol 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 include solute in first: if(s80(1:4).eq.'1SOL')read(9,*)ls1 c include solute in second: if(s80(1:4).eq.'2SOL')read(9,*)ls2 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 atom list, first: if(s80(1:4).eq.'FIRS')read(9,*)na,(al(i),i=1,abs(na)) c atom list, second: if(s80(1:4).eq.'SECO')read(9,*)nb,(bl(i),i=1,abs(nb)) 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) endif dd=(dm-dmi)/dble(nd) write(6,60000)dmi,dm,dd,nd,lsol,lper 60000 format(' RMIN ',f10.3,' RMAX ',f10.3,/, 1 ' TICK ',f10.3,' NPOINTS ',I10 ,/, 1 ' SOLU.LST ',L10 ,' LPER ',L10 ,/) if(lper)write(6,60001)px,py,pz 60001 format(' BOX ',3f10.3,/) if(na.gt.0)write(6,607)'First atom types :' if(na.lt.0)write(6,607)'First atom numbers :' 607 format(1x,a20) if(na.ne.0)write(6,608)(al(i),i=1,abs(na)) 608 format(10i6) if(nb.gt.0)write(6,607)'Second atom types :' if(nb.lt.0)write(6,607)'Second atom numbers:' if(nb.ne.0)write(6,608)(bl(i),i=1,abs(nb)) return end