program rtime c looks at FILE.X snapshots and analyzes interaction c of solvent with a solute implicit none integer*4 n1,n2,m,io,iesc,ia,ix,iargc,nat,ng,nl,nlo, 1imax,imin,smax,smin,mmax character*10 s real*8 d,dt,an,an2,time,atime,px,py,pz,atime2,tmax logical lex,notinlist,lwr real*8,allocatable::r(:),t(:) integer*4,allocatable::l(:),lo(:) c extracts columns from multiple files in F.LST if(iargc().ne.9)then write(6,*)' usage: rtime N1 N2 M d dt px py pz lwr' write(6,*) write(6,*)' N1: first atom of the solute molecule' write(6,*)' N2: last atom of the solute molecule' write(6,*)' M: number of atoms in a solvent molecule' write(6,*)' d: solvent-solute distance cutoff' write(6,*)' dt: time(fs) between snapshots' write(6,*)' pi: periodic box dimension' write(6,*)' lwr: (t/f) extensive output' write(6,*) else call getarg(1,s) read(s,*)n1 call getarg(2,s) read(s,*)n2 call getarg(3,s) read(s,*)m call getarg(4,s) read(s,*)d call getarg(5,s) read(s,*)dt call getarg(6,s) read(s,*)px call getarg(7,s) read(s,*)py call getarg(8,s) read(s,*)pz call getarg(9,s) read(s,*)lwr an=0.0d0 an2=0.0d0 time=0.0d0 atime=0.0d0 atime2=0.0d0 iesc=0 nlo=0 mmax=0 tmax=0.0d0 nl=0 inquire(file='FILE.X',exist=lex) if(lex)then c number of geometries: ng=0 c number of atoms: nat=0 open(9,file='FILE.X') 2 read(9,*,end=99,err=99) read(9,*,end=99,err=99)nat ng=ng+1 time=time+dt if(ng.eq.1)then write(6,*)nat,' atoms' allocate(r(3*nat),l(nat),lo(nat),t(nat)) do 5 io=1,nat t(io)=0.0d0 l(io)=0 5 lo(io)=0 endif do 1 ia=1,nat 1 read(9,*)r(3*(ia-1)+1),(r(3*(ia-1)+ix),ix=1,3) call makelist(nl,l,n1,n2,m,d,nat,r,px,py,pz) if(ng.eq.1)then smax=1 smin=1 imax=nl imin=nl else if(nl.gt.imax)then imax=nl smax=ng endif if(nl.lt.imin)then imin=nl smin=ng endif endif an=an+dble(nl) an2=an2+dble(nl**2) do 31 io=1,nl if(notinlist(l(io),lo,nlo))then if(lwr)write(6,*)l(io),' joined at ',time t(l(io))=time endif 31 continue do 3 io=1,nlo if(notinlist(lo(io),l,nl))then c molecule lo(io) escaped atime=atime+time-t(lo(io)) atime2=atime2+(time-t(lo(io)))**2 if(lwr)write(6,*)lo(io),' escaped after ',time-t(lo(io)) if(time-t(lo(io)).gt.tmax)then tmax=time-t(lo(io)) mmax=lo(io) endif iesc=iesc+1 endif 3 continue if(lwr)write(6,900)ng,nl,(l(io),io=1,nl) 900 format(20i4) if(lwr)write(6,901)(t(lo(io)),io=1,nl) 901 format(8x,20f9.1) c remember list: do 4 io=1,nl 4 lo(io)=l(io) nlo=nl goto 2 99 close(9) an=an/dble(ng) an2=an2/dble(ng) an2=dsqrt(an2-an**2) atime=atime/dble(iesc) atime2=atime2/dble(iesc) atime2=dsqrt(atime2-atime**2) write(6,*)ng,' geometries' write(6,6001)imax,smax,imin,smin,an,an2,atime,atime2, 1 tmax,mmax 6001 format(' maximum of ',i2,' solvents in snapshot ',i5,/, 1 ' minimum of ',i2,' solvents in snapshot ',i5,/, 1 ' average number of molecules around:',f12.3,/, 1 ' RMS variation :',f12.3,/, 1 ' average rest time (fs) :',f12.3,/, 1 ' RMS variation (fs) :',f12.3,/, 1 ' largest resting time (fs) :',f12.3,/, 1 ' for molecule :',i12,/) else write(6,*)'FILE.X does not exist' endif endif end c ============================================================ subroutine makelist(nl,l,n1,n2,m,d,nat,r,px,py,pz) c list of molecules that are around solute:l implicit none integer*4 nl,ia,n1,n2,l(*),m,nat,iap,imol,i,j,k logical notinlist real*8 r(*),d,x,y,z,xp,yp,zp,dist,px,py,pz, 1tx,ty,tz nl=0 do 1 ia=n1,n2 x=r(3*(ia-1)+1) y=r(3*(ia-1)+2) z=r(3*(ia-1)+3) do 1 iap=1,nat if(iap.lt.n1.or.iap.gt.n2)then tx=-2.0d0*px do 11 i=1,3 tx=tx+px xp=r(3*(iap-1)+1)+tx ty=-2.0d0*py do 11 j=1,3 ty=ty+py yp=r(3*(iap-1)+2)+ty tz=-2.0d0*pz do 11 k=1,3 tz=tz+pz zp=r(3*(iap-1)+3)+tz dist=dsqrt((x-xp)**2+(y-yp)**2+(z-zp)**2) if(dist.le.d)then if(iap.lt.n1)then imol=(iap-1)/m+1 else imol=(n1-1-1)/m+1+(iap-n2-1)/m+1 endif if(notinlist(imol,l,nl))then nl=nl+1 l(nl)=imol endif endif 11 continue endif 1 continue return end c ============================================================ function notinlist(imol,l,nl) implicit none logical notinlist,lo integer*4 imol,l(*),nl,i lo=.true. do 1 i=1,nl if(l(i).eq.imol)then lo=.false. goto 2 endif 1 continue 2 notinlist=lo return end c ============================================================