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

[Only registered users see links. ] (monika) wrote in message news:<246572e3.0311191948.3dad217b@posting.google. com>...
[snip...]

Looks very professional, but there's
no author's name, or company, no disclaimer,
or even a version number.

However...GOOGLE: tinker molecular modeling ~3090 hits

There's a copyright on the script. Why not
contact the people who distribute it?


Mark