tinker molecular modelling - code

Hi  I'm using a molecular modelling program called tinker.  One
execultable is called diffuse.x but it only computes the
self-diffusion constant for a homogeneous liquid using Einsteins
equation.  I don't really know any programing (rather haven't done so
in many, many years) to modify the code in order to compute the
diffusion of individual molecules.  Has anyone done this?

Here is a description of diffuse.x.  Sorry if this is excessive.
    ################################################################
c     ##                                                            ##
c     ##  program diffuse  --  find liquid self-diffusion constant  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "diffuse" finds the self-diffusion constant for a homogeneous
c     liquid via the Einstein relation from a set of stored molecular
c     dynamics frames; molecular centers of mass are unfolded and mean
c     squared displacements are computed versus time separation
c
c     the estimate for the self-diffusion constant in 10-5 cm**2/sec
c     is printed in the far right column of output and can be checked
c     by plotting mean square displacements as a function of the time
c     separation
c
c     diffusion values for very large time separation are inaccurate
c     due to the small amount of data; the current version requires
c     an orthogonal unit cell
c
c
      program diffuse
      implicit none
      include 'sizes.i'
      include 'atmtyp.i'
      include 'atoms.i'
      include 'boxes.i'
      include 'iounit.i'
      include 'molcul.i'
      integer maxmol,maxframe
      parameter (maxmol=800)
      parameter (maxframe=2000)
      integer i,j,k,m
      integer iarc,nframe
      integer next,freeunit
      integer ntime(maxframe)
      real*8 xmid,ymid,zmid
      real*8 xold,yold,zold
      real*8 xdiff,ydiff,zdiff
      real*8 dx,dy,dz,weight
      real*8 tstep,dunits,delta
      real*8 xvalue,yvalue,zvalue
      real*8 rvalue,dvalue,counts
      real*8 xmsd(maxframe)
      real*8 ymsd(maxframe)
      real*8 zmsd(maxframe)
      real*8 xcm(maxmol,maxframe)
      real*8 ycm(maxmol,maxframe)
      real*8 zcm(maxmol,maxframe)
      logical exist
     character*60 arcfile
      character*120 record
      character*120 string
c
c
c     perform the standard initialization functions
c
      call initial
c
c     try to get a filename from the command line arguments
c
      call nextarg (arcfile,exist)
      if (exist) then
         call basefile (arcfile)
         call suffix (arcfile,'arc')
         call version (arcfile,'old')
         inquire (file=arcfile,exist=exist)
      end if
c
c     ask for the user specified input structure filename
c
      dowhile (.not. exist)
         write (iout,10)
   10    format (/,' Enter Name of the Coordinate Archive File :  ',$)
         read (input,20)  arcfile
   20    format (a60)
         call basefile (arcfile)
         call suffix (arcfile,'arc')
         call version (arcfile,'old')
         inquire (file=arcfile,exist=exist)
      end do
c
c     read the first coordinate set in the archive
c
      iarc = freeunit ()
      open (unit=iarc,file=arcfile,status='old')
      call readxyz (iarc)
      nframe = 1
c
c     get the time increment between frames in picoseconds
c
      tstep = -1.0d0
      call nextarg (string,exist)
      if (exist)  read (string,*,err=30,end=30)  tstep
   30 continue
      if (tstep .le. 0.0d0) then
         write (iout,40)
   40    format (/,' Enter the Time Increment in Picoseconds',
     &              ' [1.0] :  ',$)
         read (input,50)  tstep
   50    format (f20.0)
      end if
      if (tstep .le. 0.0d0)  tstep = 1.0d0
c
c     get the unit cell axis lengths from keyfile or user
c
      call unitcell
      dowhile (xbox .eq. 0.0d0)
         write (iout,60)
   60    format (/,' Enter Unit Cell Axis Lengths :  ',$)
         read (input,70)  record
   70    format (a80)
         read (record,*,err=80,end=80)  xbox,ybox,zbox
   80    continue
         if (ybox .eq. 0.0d0)  ybox = xbox
         if (zbox .eq. 0.0d0)  zbox = xbox
      end do
c
c     set the half width values for the periodic box
c
      xbox2 = 0.5d0 * xbox
      ybox2 = 0.5d0 * ybox
      zbox2 = 0.5d0 * zbox
c
c     assign the atom parameters and count the molecules
c
      call field
      call katom
      call molecule
c
c     check for too many iudividual molecules in the system
c
      if (nmol .gt. maxmol) then
         write (iout,90)  maxmol
   90    format (/,' DIFFUSE  --  The Maximum of',i6,' Molecules',
     &              ' has been Exceeded')
         call fatal
      end if
c
c     store each molecular center of mass for the first frame
c
      do i = 1, nmol
         xmid = 0.0d0
        ymid = 0.0d0
         zmid = 0.0d0
         do j = imol(1,i), imol(2,i)
            k = kmol(j)
            weight = mass(k)
            xmid = xmid + x(k)*weight
            ymid = ymid + y(k)*weight
            zmid = zmid + z(k)*weight
         end do
         weight = molmass(i)
         xmid = xmid / weight
         ymid = ymid / weight
         zmid = zmid / weight
         xcm(i,nframe) = xmid
         ycm(i,nframe) = ymid
         zcm(i,nframe) = zmid
      end do
c
c     get the archived coordinates for each frame in turn
c
      write (iout,100)
  100 format (/,' Reading the Dynamics Trajectory Archive File :',/)
      dowhile (.true.)
         read (iarc,110,err=150,end=150)
  110    format ()
         do i = 1, n
            next = 1
            read (iarc,120)  record
  120       format (a120)
            read (record,*)  tag(i)
            call getword (record,name(i),next)
            string = record(next:120)
            read (string,*)  x(i),y(i),z(i)
         end do
         nframe = nframe + 1
         if (nframe .gt. maxframe) then
            write (iout,130)  maxframe
  130       format (/,' DIFFUSE  --  The Maximum of',i6,' Frames',
     &                 ' has been Exceeded')
            call fatal
         end if
c
c     unfold each molecule to get its corrected center of mass
c
         do i = 1, nmol
            xold = xcm(i,nframe-1)
            yold = ycm(i,nframe-1)
            zold = zcm(i,nframe-1)
            xmid = 0.0d0
            ymid = 0.0d0
            zmid = 0.0d0
            do j = imol(1,i), imol(2,i)
               k = kmol(j)
               weight = mass(k)
               xmid = xmid + x(k)*weight
               ymid = ymid + y(k)*weight
               zmid = zmid + z(k)*weight
            end do
            weight = molmass(i)
            xmid = xmid / weight
            ymid = ymid / weight
            zmid = zmid / weight
            dx = xmid - xold
            dy = ymid - yold
            dz = zmid - zold
            dowhile (dx .gt. xbox2)
               dx = dx - xbox
            end do
            dowhile (dx .lt. -xbox2)
               dx = dx + xbox
            end do
            dowhile (dy .gt. ybox2)
               dy = dy - ybox
            end do
            dowhile (dy .lt. -ybox2)
               dy = dy + ybox
            end do
            dowhile (dz .gt. zbox2)
               dz = dz - zbox
            end do
            dowhile (dz .lt. -zbox2)
               dz = dz + zbox
            end do
            xcm(i,nframe) = xold + dx
            ycm(i,nframe) = yold + dy
            zcm(i,nframe) = zold + dz
         end do
c
c     print a message after each group of frames is processed
c
         if (mod(nframe,100) .eq. 0) then
            write (iout,140)  nframe
  140       format (4x,'Processing Trajectory Frame',i11)
         end if
      end do
  150 continue
      close (unit=iarc)
      write (iout,160)  nframe
  160 format (/,' Total Number of Trajectory Frames',i8)
c
c     increment the squared displacements for each frame pair
c
      do i = 1, nframe
         ntime(i) = 0
         xmsd(i) = 0.0d0
         ymsd(i) = 0.0d0
         zmsd(i) = 0.0d0
      end do
      do i = 1, nframe-1
         do j = i+1, nframe
            m = j - i
            ntime(m) = ntime(m) + 1
            do k = 1, nmol
               xdiff = xcm(k,j) - xcm(k,i)
               ydiff = ycm(k,j) - ycm(k,i)
               zdiff = zcm(k,j) - zcm(k,i)
               xmsd(m) = xmsd(m) + xdiff*xdiff
               ymsd(m) = ymsd(m) + ydiff*ydiff
               zmsd(m) = zmsd(m) + zdiff*zdiff
            end do
         end do
      end do
c
c     get mean squared displacements and convert units;
c     conversion is from sq. Ang/ps to 10-5 sq. cm/sec
c
      dunits = 10.0d0
      do i = 1, nframe-1
         counts = dble(nmol) * dble(ntime(i))
         xmsd(i) = xmsd(i) * (dunits/counts)
         ymsd(i) = ymsd(i) * (dunits/counts)
         zmsd(i) = zmsd(i) * (dunits/counts)
      end do
c
c     estimate the diffusion constant via the Einstein relation
c
      write (iout,170)
  170 format (/,' Mean Squared Diffusion Distance and Self-Diffusion',
     &           ' Constant :',
     &        //,5x,'Time Step',5x,'X MSD',7x,'Y MSD',7x,'Z MSD',
     &           7x,'R MSD',4x,'Diff Const',/)
      do i = 1, nframe-1
         delta = tstep * dble(i)
         xvalue = xmsd(i) / 2.0d0
         yvalue = ymsd(i) / 2.0d0
         zvalue = zmsd(i) / 2.0d0
         rvalue = (xmsd(i) + ymsd(i) + zmsd(i)) / 6.0d0
         dvalue = rvalue / delta
         write (iout,180)  delta,xvalue,yvalue,zvalue,rvalue,dvalue
  180    format (f12.2,4f12.2,f12.4)
      end do

c     perform any final tasks before program exit
c
      call final
      end