CHARMM Development Project Forums CHARMM Community Script Archive self-diffusion-coeff.inp

*FILENAME: self-diffusion-coeff.inp
*PURPOSE: compute Mean Square Displacement vs time to get diffusion coeffcient D for water
* via the Einstein relation: 6*D is the slope of MSD vs t at long t.
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
*
!unix environment variable CHM_HOME points to CHARMM installation directory
read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp
read psf card name my_psf.psf

open unit 51 read unform name my_traj1.cor
open unit 52 read unform name my_traj2.cor

open unit 21 write form name msd.dat
write title unit 21
* TIME MSD
*

coor anal select type oh2 end - ! what atoms to look at, these are the TIP3 oxygens
firstu 51 nunit 2 skip 50 - ! trajectory specification
imsd 21 - ! flag to do the MSD analysis
rspin 0.0 rspout 999.9 - ! we are interested in ALL waters
ncors 200 - ! compute MSD out to NCORS*SKIP steps
xbox 40.0 ybox 45.0 zbox 45.0 ! and we did use PBC, simple rectangular box
! if trajectory is from constant pressure simulation
! actual box size is taken from trajectory

Hi Prof Nilsson,

Quote:

---

! if trajectory is from constant pressure simulation
! actual box size is taken from trajectory

---

Is it referring to @FILE-box-size.str? Or is there some way of getting the sizes instantaneously? from traj?

the statement "if trajectory is from constant pressure simulation
actual box size is taken from trajectory" means exactly what it says, namely that if you have a trajectory that was performed using constant pressure (in which case the box size may vary), the box size at each step (which in this case is stored in the trajectory) is taken from the trajectory. you do not have to do anything about it.

Can this be modified for computing the self-diffusion within RHDO pbc? If so, would you be so kind as to reveal the modifications?

Thanks a bunch!

The current implementation of PBC corrections in the COOR ANAL code is very simple, and only applies to orthorhombic lattices.

The problems with computing diffusion coefficients using the Einstein relation in PBC are: 1/ mean square displacements (MSD) at long times are underestimated due to the limited distance particles can travel since they are confined to the box by the recentering, and 2/ you can get spuriously high MSDs when molecules are moved across the box by the recentering. It is the second effect that is corrected by the PBC code in COOR ANAL (the first effect is usually not a problem because you can get a good estimate of the slope of MSD vs t before it sets in).

You can try approximately three different approaches:
1/ Implement PBC corrections that can also handle RHDO in subroutine corfuncb of source/solana.src (you can see how the same thing is done in source/manip/hbanal.src in the most recent version of CHARMM). When done, and satisfied that it works, please submit changes. I will try to find the time to put this in for the next release.

2/ Undo the recentering events using MERGE UNFOLD on your trajectory before the analysis.

3/ Restrict the MSD calculations to molecules that are not in the vicinity of the RHDO borders (to avoid problems of the second kind): Set RSPOut to a value well inside your RHDO. This is the simplest way, but you probably have enough data to get a reasonable estimate of the diffusion coefficient anyway; whether this is the case is highly system-dependent.

Thank you for the helpful response! I will consider implementation, if other research things go smoothly as to release my time constraints.

-Demian

The MERGE COOR UNFOLD operation is also currently restricted to orthonormal shapes (CUBIC, TETRagonal, ORTHorhombic).

Dear All,

In the script it shows that the trajectory was saved every 50 steps. My question is how often we should save the trajectory for diffusion coefficient calculation of TIP3 water. Is there any rule or something?

What is more, I understand that only the linear part of MSD should be used for diffusion coefficient calculation. However, is there any recommend value of ncors to calculate the diffusion of water base on your experience?

Thanks,
Leo

please read the relevant original literature!
for water you can get a good estimate of Dtransl for the 1-20 ps part of MSD(t)

Thanks, Lennart. I do appreciate your help.

-Leo

I would like to know, weather I can use above script for truncated octahedren box or not. If not what are the changes should I make in the script. My second question is " what does this "good estimate of Dtransl" means or what is the range of good estimate."

I will really appreciate if some one can help me.

Sorry, truncated octahedron is also not directly supported by this code. See previous posts in this thread for suggestions what to do.
"good" means to within a few percent of the value that would be obtained using much longer times for the fit of the straight line to the MSD(t) curve.

Dear all,

I am sorry that my question repeats something that was already wrote in a post above concerning the calculation of water diffusion coefficient based on MSD using coor anal command, but I am having trouble understanding the connection between the:
- number of frames of a trajectory,
- the frequency of saving coordinates in the trajectory,
- the timestep used in the dynamics simulation,
- the skip value used when reading the trajectory in coor anal command and
- ncors value in coor anal.
Are there any guidelines on selecting skip and ncors values depending on the analyzed trajectory features mentioned above (timestep, frequency of saving coordinates, number of frames of the analyzed trajectory)?
Thank you!

A trajectory of 1 ns simulated with a 2 fs time step, and with coordinates saved every 500 steps (NSAVC 500 in the dynamics command) will contain 1000 frames (1000*0.002ps*500=1000ps=1ns).
If you read this trajectory with NSKIP 2500 only every fifth frame (ie, every 5ps) will be read.
NCORS is the number of points in the MSD(t) curve.

For the above trajectory the NCORS 100 and NSKIP 1000 would compute MSD(t) from t=0 to t= 100 * 1000 * 0.002ps=200ps

You have to compute MSD(t) to a large enough t that you clearly see the linear part. Using a small NSKIP increases the computational time, but this is usually not a problem.

Dear Sir,

Thank you for your answer, I hope I understand it now.
I calculated the MSD for the same trajectory using different combinations of ncors and nskip, each time charmm gives a different estimate of diffusion coefficient, of the order E00 to E-06, and even negative values. I want to ask how reliable this estimation is, or maybe I should plot the data and do the linear fit myself.
I attached such a plot of MSD vs time and the plot is not that linear. Do you think it looks OK? Should I fit only the first 10 ps? In any way I fit it, I can not reproduce the estimated diffusion coefficient from charmm output.
I would be very grateful if you could help me on this matter.

You should of course do the fit yourself, especially with curve like the one in the graph you attached, which does look a bit strange.
Try to reproduce some published diffusion coefficient calculation for a box of water.