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