program flow implicit none character*80 fn character*5 so integer*4 nd(3),nper,nats,nc,nmc,nat,ia,ix,icu,ipr,i real*8 aaxis,baxis,caxis,tper,tsn,r0(3),dr(3),x,y,z,xr,yr,zr, 1dx,dy,dz,fac,aaxis2,baxis2,caxis2,acage,bcage,ccage,ave,ive,xi, 2weight,avogadro,xx,yy,zz,facd logical lsmooth,lstand real*8,allocatable::vx(:),vy(:),vz(:),ro(:),ra(:),r(:),vi(:) integer*4,allocatable::alist(:),izq(:),iza(:) real*4 amas(89) 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/ data avogadro/6.02214129D23/ c write(6,6000) 6000 format(/,' MD flow monitoring',/,/, 1 ' Input:',/, 1 ' FLOW.OPT... list of options',/, 1 ' FILE.X ... snapshot geometries',/,/, 1 ' Output:',/, 1 ' vx.cub ... vx density',/, 1 ' vy.cub ... vy density',/, 1 ' vz.cub ... vz density',/, 1 ' vi.cub ... vi density',/, 1 ' v.cub ... v density',/, 1 ' ro.cub ... density',/) c c just for allocation, do not read all options: open(9,file='FLOW.OPT') 98 read(9,900,err=99,end=99)so 900 format(a5) if(so.eq.'ATOMS')then read(9,*)nats allocate(alist(nats),ra(3*nats),izq(nats)) endif goto 98 99 close(9) c read rest of options: call readopt(fn,lsmooth,aaxis2,baxis2,caxis2, 1aaxis,baxis,caxis,acage,bcage,ccage,r0, 1tper,tsn,nd,nats,nper,alist,lstand) c c check that everything is initialize and write control message: call checkv(aaxis,baxis,caxis,acage,bcage,ccage,r0, 1tper,tsn,nd,nats,nper,nc,dr,lstand,lsmooth) allocate(vx(nc),vy(nc),vz(nc),ro(nc),vi(nc)) do 101 ix=1,nc vi(ix)=0.0d0 vx(ix)=0.0d0 vy(ix)=0.0d0 vz(ix)=0.0d0 101 ro(ix)=0.0d0 c average velocity: ave=0.0d0 ive=0.0d0 nmc=0 c icu - current geometry: icu=0 c ipr - reference geometry i ntime - tsn * nper: ipr=icu-nper c start reading the snapshots: open(4,file=fn) 1 read(4,*,end=400,err=400) read(4,*,end=400,err=400)nat nmc=nmc+1 c initialize coordinates for nper geometries: if(nmc.eq.1)allocate(r(3*nat*(nper+1)),iza(nat)) icu=icu+1 if(icu.gt.nper+1)icu=1 ipr=ipr+1 if(ipr.gt.nper+1)ipr=1 if(mod(nmc,10).eq.0)write(6,6003)nmc,icu,ipr,ave/ive 6003 format(i12,2i3,f10.4) do 4 ia=1,nat 4 read(4,*)iza(ia),(r(3*nat*(icu-1)+3*(ia-1)+ix),ix=1,3) if(nmc.eq.1)then do 3 ia=1,nats 3 izq(ia)=iza(alist(ia)) do 5 ia=1,nats do 5 ix=1,3 5 ra(ix+3*(ia-1))=r(ix+3*(alist(ia)-1)) write(6,*)nats,' atoms selected from the first geometry' endif if(nmc.gt.nper)then do 6 ia=1,nat c exclude listed atoms: if(lstand)then do 7 i=1,nats 7 if(alist(i).eq.ia)goto 6 endif x =r(3*nat*(icu-1)+3*(ia-1)+1) y =r(3*nat*(icu-1)+3*(ia-1)+2) z =r(3*nat*(icu-1)+3*(ia-1)+3) xr=r(3*nat*(ipr-1)+3*(ia-1)+1) yr=r(3*nat*(ipr-1)+3*(ia-1)+2) zr=r(3*nat*(ipr-1)+3*(ia-1)+3) call fixp(x,xr,aaxis2,aaxis) call fixp(y,yr,baxis2,baxis) call fixp(z,zr,caxis2,caxis) c velocity in A/fs: dx=x-xr dy=y-yr dz=z-zr xi=dsqrt(dx**2+dy**2+dz**2) ave=ave+xi ive=ive+1.0d0 xx=(x+xr)/2.0d0 yy=(y+yr)/2.0d0 zz=(z+zr)/2.0d0 weight=dble(amas(iza(ia))) if(lsmooth)then call addxsmooth(nd,xx,yy,zz,r0,dr,vi,xi,aaxis,baxis,caxis) call addxsmooth(nd,xx,yy,zz,r0,dr,vx,dx,aaxis,baxis,caxis) call addxsmooth(nd,xx,yy,zz,r0,dr,vy,dy,aaxis,baxis,caxis) call addxsmooth(nd,xx,yy,zz,r0,dr,vz,dz,aaxis,baxis,caxis) call addxsmooth(nd,x,y,z,r0,dr,ro,weight,aaxis,baxis,caxis) else call addx(nd,xx,yy,zz,r0,dr,vi,xi, 1 aaxis,baxis,caxis) call addx(nd,xx,yy,zz,r0,dr,vx,dx, 1 aaxis,baxis,caxis) call addx(nd,xx,yy,zz,r0,dr,vy,dy, 1 aaxis,baxis,caxis) call addx(nd,xx,yy,zz,r0,dr,vz,dz, 1 aaxis,baxis,caxis) call addx(nd,x,y,z,r0,dr,ro,weight, 1 aaxis,baxis,caxis) endif 6 continue endif goto 1 400 close(4) fac=1.0d0/dble(nmc) c density from g/mol/A^3 to g/ml: facd=1.0d24/dble(nmc)/avogadro/dr(1)/dr(2)/dr(3) do 103 ix=1,nc vi(ix)=vi(ix)*fac/tper vx(ix)=vx(ix)*fac/tper vy(ix)=vy(ix)*fac/tper vz(ix)=vz(ix)*fac/tper 103 ro(ix)=ro(ix)*facd write(6,6002)ave/ive,tper 6002 format('average shift ',f10.3,' A per ',f5.1,' fs') call wrf('vx.cub',nats,izq,ra,dr,r0,nd,vx) call wrf('vi.cub',nats,izq,ra,dr,r0,nd,vi) call wrf('vy.cub',nats,izq,ra,dr,r0,nd,vy) call wrf('vz.cub',nats,izq,ra,dr,r0,nd,vz) call wrf('ro.cub',nats,izq,ra,dr,r0,nd,ro) do 102 ix=1,nc 102 ro(ix)=dsqrt(vx(ix)**2+vy(ix)**2+vz(ix)**2) call wrf('v.cub',nats,izq,ra,dr,r0,nd,ro) stop end subroutine report(s) character*(*) s write(6,*)s write(6,*)'****** Program stopped. *********' stop end subroutine wrf(sn,nats,izq,ra,dr,r0,nd,rh) implicit none integer*4 nats,izq(*),nd(3),i,ix,iy,iz real*8 ra(*),rh(*),r0(3),dr(3),bohr character*(*)sn data bohr/0.529177d0/ c open(3,file=sn) write(3,3000)nats,(r0(ix)/bohr,ix=1,3),nd(1),dr(1)/bohr,0.0,0.0, 1nd(2),0.0,dr(2)/bohr,0.0,nd(3),0.0,0.0,dr(3)/bohr 3000 format(' Nuclear density',/,' SCF Total Density',/, 13(i5,3f12.6,/),i5,3f12.6) do 2 i=1,nats 2 write(3,3001)izq(i),dble(izq(i)),(ra(ix+3*(i-1))/bohr,ix=1,3) 3001 format(i5,4f12.6) do 1 ix=1,nd(1) do 1 iy=1,nd(2) 1 write(3,4000)(rh(iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2)),iz=1,nd(3)) 4000 format(6E13.5) close(3) write(6,*)sn return end subroutine addxsmooth(nd,xw,yw,zw,r0,dr,r,weight, 1aaxis,baxis,caxis) implicit none real*8 r(*),r0(3),dr(3),xw,yw,zw,vx,vy,vz,wx,wy,wz, 1w(27),a14,dr1,dr2,dr3,an,p,weight,aaxis,baxis,caxis,down integer*4 j,nd(3),ix,i,k,iy,iz,ni,ii,jj,kk,nd1,nd2,nd3,i9,j3 c dr1=dr(1) dr2=dr(2) dr3=dr(3) nd1=nd(1) nd2=nd(2) nd3=nd(3) if(xw.lt.r0(1))xw=xw+aaxis if(yw.lt.r0(2))yw=yw+baxis if(zw.lt.r0(3))zw=zw+caxis if(xw.gt.r0(1)+dble(nd1-1)*dr1)xw=xw-aaxis if(yw.gt.r0(2)+dble(nd2-1)*dr2)yw=yw-baxis if(zw.gt.r0(3)+dble(nd3-1)*dr3)zw=zw-caxis ix=int((xw-r0(1))/dr1)+1 wx=xw-dr1*(0.50d0+dble(ix-1))-r0(1) iy=int((yw-r0(2))/dr2)+1 wy=yw-dr2*(0.50d0+dble(iy-1))-r0(2) iz=int((zw-r0(3))/dr3)+1 wz=zw-dr3*(0.50d0+dble(iz-1))-r0(3) c occasional error check - should not happen: if(ix.le.0.or.ix.gt.nd1)return if(iy.le.0.or.iy.gt.nd2)return if(iz.le.0.or.iz.gt.nd3)return c c zero cube contributions do 2 i=1,27 2 w(i)=0.d0 a14=0.0d0 c loop over the all 27 cubes vx=-dr1-dr1 do 1 i=-1,1 vx=vx+dr1 vy=-dr2-dr2 i9=(i+1)*9 do 1 j=-1,1 vy=vy+dr2 vz=-dr3-dr3 j3=(j+1)*3+i9 do 1 k=-1,1 vz=vz+dr3 if(k.ne.0.or.j.ne.0.or.i.ne.0)then down=(vx*vx+vy*vy+vz*vz) if(down.ne.0.0d0)then p=(wx*vx+wy*vy+wz*vz)/down else p=0.0d0 endif if(dabs(p).gt.1.d0)then write(6,*)xw,yw,zw,ix,iy,iz,dr1,dr2,dr3,r0,vx,vy,vz iy=int((yw-r0(2))/dr2)+1 write(6,*)ix,'xori' call fixper(ix,nd2) write(6,*)ix,'xfixed' write(6,*)iy,'yori' call fixper(iy,nd2) write(6,*)iy,'yfixed' write(6,*)iz,'zori' call fixper(iz,nd2) write(6,*)iz,'zfixed' write(6,*)p write(6,*)wx,wy,wz call report('error') endif if(p.gt.0.d0)then w(k+2+j3)=p w(14)=w(14)+(1.0d0-p) a14=a14+1.0d0 endif endif 1 continue c c adjust central cube: w(14)=w(14)/a14 c c normalize total contribution to the weight: an=0.0d0 do 3 i=1,27 3 an=an+w(i) if(an.gt.0.0d0)then do 4 i=1,27 4 w(i)=w(i)/an*weight endif c c increment RDF: do 5 i=-1,1 ii=ix+i call fixper(ii,nd1) i9=(i+1)*9 do 5 j=-1,1 jj=iy+j call fixper(jj,nd2) j3=(j+1)*3+i9 do 5 k=-1,1 kk=iz+k call fixper(kk,nd3) ni=kk+(jj-1+(ii-1)*nd2)*nd3 5 r(ni)=r(ni)+w(k+2+j3) return end subroutine fixper(i,n) implicit none integer*4 i,n if(i.lt.1)then i=i+n-1 else if(i.gt.n)then i=i-n+1 endif endif return end subroutine addx(nd,xw,yw,zw,r0,dr,r,weight, 1aaxis,baxis,caxis) implicit none real*8 r(*),r0(3),dr(3),xw,yw,zw,weight, 1dr1,dr2,dr3,aaxis,baxis,caxis integer*4 nd(3),ix,iy,iz,ni,nd1,nd2,nd3 dr1=dr(1) dr2=dr(2) dr3=dr(3) nd1=nd(1) nd2=nd(2) nd3=nd(3) if(xw.lt.r0(1))xw=xw+aaxis if(yw.lt.r0(2))yw=yw+baxis if(zw.lt.r0(3))zw=zw+caxis if(xw.gt.r0(1)+dble(nd1-1)*dr1)xw=xw-aaxis if(yw.gt.r0(2)+dble(nd2-1)*dr2)yw=yw-baxis if(zw.gt.r0(3)+dble(nd3-1)*dr3)zw=zw-caxis ix=int((xw-r0(1))/dr(1))+1 iy=int((yw-r0(2))/dr(2))+1 iz=int((zw-r0(3))/dr(3))+1 c occasional error check - should not happen: if(ix.le.0.or.ix.gt.nd(1))return if(iy.le.0.or.iy.gt.nd(2))return if(iz.le.0.or.iz.gt.nd(3))return ni=iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2) r(ni)=r(ni)+weight return end subroutine fixp(z,zr,caxis2,caxis) implicit none real*8 z,zr,caxis2,caxis if(z.gt.zr+caxis2)z=z-caxis if(z.lt.zr-caxis2)z=z+caxis return end subroutine readopt(fn,lsmooth,aaxis2,baxis2,caxis2, 1aaxis,baxis,caxis,acage,bcage,ccage,r0, 1tper,tsn,nd,nats,nper,alist,lstand) implicit none character*(*)fn character*5 so real*8 aaxis,baxis,caxis,acage,bcage,ccage,r0(*),z0, 1tper,tsn,aaxis2,baxis2,caxis2 integer*4 nats,nd(*),nper,alist(*),ia logical lsmooth,lstand z0=0.0d0 fn='FILE.X' lsmooth=.true. lstand=.false. tper=z0 tsn=z0 nper=0 aaxis=z0 baxis=z0 caxis=z0 acage=z0 bcage=z0 ccage=z0 nd(1)=0 nd(2)=0 nd(3)=0 c origin of the cage: r0(1)=-99.99d0 c number of atoms of the standing molecules nats=0 c lsmooth - smoothing of the distributions c lstand - exclude listed atoms from the statistics c tsn [fs] - timem interval between the snapshots c nper - number in snapshots in integration period c tper = nper * tsn, v = (x(t+tper)-x(t))/tper write(6,*)'FLOW.OPT:' open(9,file='FLOW.OPT') 98 read(9,900,err=99,end=99)so write(6,900)so 900 format(a5) if(so.eq.'NPERI')read(9,*)nper if(so.eq.'A-AXI')read(9,*)aaxis if(so.eq.'B-AXI')read(9,*)baxis if(so.eq.'C-AXI')read(9,*)caxis if(so.eq.'A-CAG')read(9,*)acage if(so.eq.'B-CAG')read(9,*)bcage if(so.eq.'C-CAG')read(9,*)ccage if(so.eq.'X-ORI')read(9,*)r0(1) if(so.eq.'Y-ORI')read(9,*)r0(2) if(so.eq.'Z-ORI')read(9,*)r0(3) if(so.eq.'X-POI')read(9,*)nd(1) if(so.eq.'Y-POI')read(9,*)nd(2) if(so.eq.'Z-POI')read(9,*)nd(3) if(so.eq.'TSNAP')read(9,*)tsn if(so.eq.'SMOOT')read(9,*)lsmooth if(so.eq.'STAND')read(9,*)lstand if(so.eq.'FILEN')read(9,9001)fn 9001 format(a80) if(so.eq.'ATOMS')then read(9,*)nats,(alist(ia),ia=1,nats) endif goto 98 99 close(9) if(acage.eq.0.0d0)then acage=aaxis bcage=baxis acage=caxis write(6,*)'cage set according to the periodic box' endif c origin of the cage: if(r0(1).eq.-99.99d0)then r0(1) = -aaxis*0.5d0 r0(2) = -baxis*0.5d0 r0(3) = -caxis*0.5d0 write(6,6001) 6001 format(/,' cage origin set to -MDside/2',/) endif tper=dble(nper)*tsn aaxis2=aaxis/2.0d0 baxis2=baxis/2.0d0 caxis2=caxis/2.0d0 return end subroutine checkv(aaxis,baxis,caxis,acage,bcage,ccage,r0, 1tper,tsn,nd,nats,nper,nc,dr,lstand,lsmooth) implicit none real*8 z0,aaxis,baxis,caxis,acage,bcage,ccage,r0(*), 1tper,tsn,dr(*) integer*4 nd(*),nats,nper,nc logical lstand,lsmooth z0=0.0d0 if(aaxis.eq.z0)call report('aaxis not set') if(baxis.eq.z0)call report('baxis not set') if(caxis.eq.z0)call report('caxis not set') if(acage.eq.z0)call report('aacage not set') if(bcage.eq.z0)call report('bacage not set') if(ccage.eq.z0)call report('cacage not set') if(tper.eq.z0)call report('tper not set') if(nper.eq.0)call report('nper not set') if(tsn.eq.z0)call report('tsn not set') write(6,*) write(6,4031)'A-AXIS MD box [A ]:',aaxis write(6,4031)'B-AXIS [A ]:',baxis write(6,4031)'C-AXIS [A ]:',caxis write(6,4031)'A-CAGE cub dimension [A ]:',acage write(6,4031)'B-CAGE [A ]:',bcage write(6,4031)'C-CAGE [A ]:',ccage write(6,4031)'X-ORIGIN cage origin [A ]:',r0(1) write(6,4031)'Y-ORIGIN cage origin [A ]:',r0(2) write(6,4031)'Z-ORIGIN cage origin [A ]:',r0(3) write(6,4031)'TPERIODIC integration period [fs]:',tper write(6,4031)'TSNAP time between snapshots [fs]:',tsn 4031 format(a39,f8.1) write(6,4003)'SMOOTH smoothing of dist:',lsmooth write(6,4003)'STAND exclude listed ats:',lstand 4003 format(a25,l10) write(6,4004)'X-POINTS number :',nd(1) write(6,4004)'Y-POINTS number :',nd(2) write(6,4004)'Z-POINTS number :',nd(3) write(6,4004)'ATOMS number :',nats write(6,4004)'NPERI snapshots in per. :',nper 4004 format(a25,i10) dr(1)=acage/dble(nd(1)-1) dr(2)=bcage/dble(nd(2)-1) dr(3)=ccage/dble(nd(3)-1) nc=nd(1)*nd(2)*nd(3) return end