1c 2c 3c ############################################################ 4c ## COPYRIGHT (C) 1995 by Yong Kong & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################ 7c 8c ################################################################ 9c ## ## 10c ## program diffuse -- find liquid self-diffusion constant ## 11c ## ## 12c ################################################################ 13c 14c 15c "diffuse" finds the self-diffusion constant for a homogeneous 16c liquid via the Einstein relation from a set of stored molecular 17c dynamics frames; molecular centers of mass are unfolded and mean 18c squared displacements are computed versus time separation 19c 20c the estimate for the self-diffusion constant in 10-5 cm**2/sec 21c is printed in the far right column of output and can be checked 22c by plotting mean squared displacements as a function of the time 23c separation; values for very large time separation are inaccurate 24c due to the small amount of data 25c 26c 27 program diffuse 28 use atomid 29 use atoms 30 use bound 31 use inform 32 use iounit 33 use molcul 34 use usage 35 implicit none 36 integer i,j,k,m 37 integer nframe,iframe 38 integer iarc,start,stop 39 integer step,skip,size 40 integer, allocatable :: list(:) 41 integer, allocatable :: ntime(:) 42 real*8 xmid,ymid,zmid 43 real*8 xold,yold,zold 44 real*8 xdiff,ydiff,zdiff 45 real*8 xr,yr,zr,weigh 46 real*8 tstep,dunits,delta 47 real*8 xvalue,yvalue,zvalue 48 real*8 rvalue,dvalue,counts 49 real*8, allocatable :: xmsd(:) 50 real*8, allocatable :: ymsd(:) 51 real*8, allocatable :: zmsd(:) 52 real*8, allocatable :: xcm(:,:) 53 real*8, allocatable :: ycm(:,:) 54 real*8, allocatable :: zcm(:,:) 55 logical exist,query 56 character*240 record 57 character*240 string 58c 59c 60c perform the standard initialization functions 61c 62 call initial 63c 64c open the trajectory archive and read the initial frame 65c 66 call getarc (iarc) 67c 68c get numbers of the coordinate frames to be processed 69c 70 start = 1 71 stop = 100000 72 step = 1 73 query = .true. 74 call nextarg (string,exist) 75 if (exist) then 76 read (string,*,err=10,end=10) start 77 query = .false. 78 end if 79 call nextarg (string,exist) 80 if (exist) read (string,*,err=10,end=10) stop 81 call nextarg (string,exist) 82 if (exist) read (string,*,err=10,end=10) step 83 10 continue 84 if (query) then 85 write (iout,20) 86 20 format (/,' Numbers of First & Last Frame and Step', 87 & ' Increment : ',$) 88 read (input,30) record 89 30 format (a240) 90 read (record,*,err=40,end=40) start,stop,step 91 40 continue 92 end if 93c 94c get the time increment between frames in picoseconds 95c 96 tstep = -1.0d0 97 call nextarg (string,exist) 98 if (exist) read (string,*,err=50,end=50) tstep 99 50 continue 100 if (tstep .le. 0.0d0) then 101 write (iout,60) 102 60 format (/,' Enter the Time Increment in Picoseconds', 103 & ' [1.0] : ',$) 104 read (input,70) tstep 105 70 format (f20.0) 106 end if 107 if (tstep .le. 0.0d0) tstep = 1.0d0 108c 109c get the atom parameters, lattice type and molecule count 110c 111 call field 112 call unitcell 113 call lattice 114 call katom 115 call molecule 116c 117c perform dynamic allocation of some local arrays 118c 119 size = 40 120 allocate (list(size)) 121c 122c find atoms and molecules to be excluded from consideration 123c 124 call active 125 if (nuse .eq. n) then 126 do i = 1, size 127 list(i) = 0 128 end do 129 i = 0 130 do while (exist) 131 call nextarg (string,exist) 132 if (exist) then 133 read (string,*,err=80,end=80) list(i+1) 134 i = i + 1 135 end if 136 end do 137 80 continue 138 if (i .eq. 0) then 139 write (iout,90) 140 90 format (/,' Numbers of any Atoms to be Excluded : ',$) 141 read (input,100) record 142 100 format (a240) 143 read (record,*,err=110,end=110) (list(i),i=1,size) 144 110 continue 145 end if 146 i = 1 147 do while (list(i) .ne. 0) 148 list(i) = max(-n,min(n,list(i))) 149 if (list(i) .gt. 0) then 150 k = list(i) 151 if (use(k)) then 152 use(k) = .false. 153 nuse = nuse - 1 154 end if 155 i = i + 1 156 else 157 list(i+1) = max(-n,min(n,list(i+1))) 158 do k = abs(list(i)), abs(list(i+1)) 159 if (use(k)) then 160 use(k) = .false. 161 nuse = nuse - 1 162 end if 163 end do 164 i = i + 2 165 end if 166 end do 167 end if 168c 169c perform deallocation of some local arrays 170c 171 deallocate (list) 172c 173c alter the molecule list to include only active molecules 174c 175 do i = 1, nmol 176 do j = imol(1,i), imol(2,i) 177 k = kmol(j) 178 if (.not. use(k)) imol(1,i) = 0 179 end do 180 end do 181 k = 0 182 do i = 1, nmol 183 if (imol(1,i) .ne. 0) then 184 k = k + 1 185 imol(1,k) = imol(1,i) 186 imol(2,k) = imol(2,i) 187 molmass(k) = molmass(i) 188 end if 189 end do 190 nmol = k 191 write (iout,120) nmol 192 120 format (/,' Total Number of Molecules :',i16) 193c 194c count the number of coordinate frames in the archive file 195c 196 abort = .false. 197 rewind (unit=iarc) 198 nframe = 0 199 do while (.not. abort) 200 call readxyz (iarc) 201 nframe = nframe + 1 202 end do 203 nframe = nframe - 1 204 rewind (unit=iarc) 205 stop = min(nframe,stop) 206 nframe = (stop-start)/step + 1 207 write (iout,130) nframe 208 130 format (/,' Number of Coordinate Frames :',i14) 209c 210c perform dynamic allocation of some local arrays 211c 212 allocate (ntime(nframe)) 213 allocate (xmsd(nframe)) 214 allocate (ymsd(nframe)) 215 allocate (zmsd(nframe)) 216 allocate (xcm(nmol,nframe)) 217 allocate (ycm(nmol,nframe)) 218 allocate (zcm(nmol,nframe)) 219c 220c get the archived coordinates for each frame in turn 221c 222 write (iout,140) 223 140 format (/,' Reading the Coordinates Archive File :',/) 224 nframe = 0 225 iframe = start 226 skip = start 227 do while (iframe.ge.start .and. iframe.le.stop) 228 do j = 1, skip-1 229 call readxyz (iarc) 230 end do 231 iframe = iframe + step 232 skip = step 233 call readxyz (iarc) 234 if (n .eq. 0) goto 160 235 nframe = nframe + 1 236 if (mod(nframe,100) .eq. 0) then 237 write (iout,150) nframe 238 150 format (4x,'Processing Coordinate Frame',i13) 239 end if 240c 241c unfold each molecule to get its corrected center of mass 242c 243 do i = 1, nmol 244 xmid = 0.0d0 245 ymid = 0.0d0 246 zmid = 0.0d0 247 do j = imol(1,i), imol(2,i) 248 k = kmol(j) 249 weigh = mass(k) 250 xmid = xmid + x(k)*weigh 251 ymid = ymid + y(k)*weigh 252 zmid = zmid + z(k)*weigh 253 end do 254 weigh = molmass(i) 255 xmid = xmid / weigh 256 ymid = ymid / weigh 257 zmid = zmid / weigh 258 if (nframe .eq. 1) then 259 xold = xmid 260 yold = ymid 261 zold = zmid 262 else 263 xold = xcm(i,nframe-1) 264 yold = ycm(i,nframe-1) 265 zold = zcm(i,nframe-1) 266 end if 267 xr = xmid - xold 268 yr = ymid - yold 269 zr = zmid - zold 270 if (use_bounds) call image (xr,yr,zr) 271 xcm(i,nframe) = xold + xr 272 ycm(i,nframe) = yold + yr 273 zcm(i,nframe) = zold + zr 274 end do 275 end do 276 160 continue 277 close (unit=iarc) 278 if (mod(nframe,100) .ne. 0) then 279 write (iout,170) nframe 280 170 format (4x,'Processing Coordinate Frame',i13) 281 end if 282c 283c increment the squared displacements for each frame pair 284c 285 do i = 1, nframe 286 ntime(i) = 0 287 xmsd(i) = 0.0d0 288 ymsd(i) = 0.0d0 289 zmsd(i) = 0.0d0 290 end do 291 do i = 1, nframe-1 292 do j = i+1, nframe 293 m = j - i 294 ntime(m) = ntime(m) + 1 295 do k = 1, nmol 296 xdiff = xcm(k,j) - xcm(k,i) 297 ydiff = ycm(k,j) - ycm(k,i) 298 zdiff = zcm(k,j) - zcm(k,i) 299 xmsd(m) = xmsd(m) + xdiff*xdiff 300 ymsd(m) = ymsd(m) + ydiff*ydiff 301 zmsd(m) = zmsd(m) + zdiff*zdiff 302 end do 303 end do 304 end do 305c 306c perform deallocation of some local arrays 307c 308 deallocate (xcm) 309 deallocate (ycm) 310 deallocate (zcm) 311c 312c get mean squared displacements and convert units; 313c conversion is from sq. Ang/ps to 10-5 sq. cm/sec 314c 315 dunits = 10.0d0 316 do i = 1, nframe-1 317 counts = dble(nmol) * dble(ntime(i)) 318 xmsd(i) = xmsd(i) * (dunits/counts) 319 ymsd(i) = ymsd(i) * (dunits/counts) 320 zmsd(i) = zmsd(i) * (dunits/counts) 321 end do 322c 323c estimate the diffusion constant via the Einstein relation 324c 325 write (iout,180) 326 180 format (/,' Mean Squared Displacements and Self-Diffusion', 327 & ' Constant :', 328 & //,5x,'Time Gap',6x,'X MSD',7x,'Y MSD',7x,'Z MSD', 329 & 7x,'R MSD',4x,'Diff Const', 330 & /,7x,'(ps)',9x,'(/2)',8x,'(/2)',8x,'(/2)',8x,'(/6)', 331 & 5x,'(x 10^5)',/) 332 do i = 1, nframe-1 333 delta = tstep * dble(i) 334 xvalue = xmsd(i) / 2.0d0 335 yvalue = ymsd(i) / 2.0d0 336 zvalue = zmsd(i) / 2.0d0 337 rvalue = (xmsd(i) + ymsd(i) + zmsd(i)) / 6.0d0 338 dvalue = rvalue / delta 339 write (iout,190) delta,xvalue,yvalue,zvalue,rvalue,dvalue 340 190 format (f12.2,4f12.2,f12.4) 341 end do 342c 343c perform deallocation of some local arrays 344c 345 deallocate (ntime) 346 deallocate (xmsd) 347 deallocate (ymsd) 348 deallocate (zmsd) 349c 350c perform any final tasks before program exit 351c 352 call final 353 end 354