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