|
Joined: Sep 2003
Posts: 4,883 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,883 Likes: 12 |
*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
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: May 2004
Posts: 222
Forum Member
|
Forum Member
Joined: May 2004
Posts: 222 |
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?
CHARMM 30b1 driven by
1/ Xeon (32 bits)
2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel
3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum.
4/ Gracious support from the forum.
|
|
|
|
Joined: Sep 2003
Posts: 4,883 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,883 Likes: 12 |
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.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Apr 2004
Posts: 9
Forum Member
|
Forum Member
Joined: Apr 2004
Posts: 9 |
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!
|
|
|
|
Joined: Sep 2003
Posts: 4,883 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,883 Likes: 12 |
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.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Apr 2004
Posts: 9
Forum Member
|
Forum Member
Joined: Apr 2004
Posts: 9 |
Thank you for the helpful response! I will consider implementation, if other research things go smoothly as to release my time constraints.
-Demian
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
The MERGE COOR UNFOLD operation is also currently restricted to orthonormal shapes (CUBIC, TETRagonal, ORTHorhombic).
Rick Venable computational chemist
|
|
|
|
Joined: Oct 2004
Posts: 22
Forum Member
|
Forum Member
Joined: Oct 2004
Posts: 22 |
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
|
|
|
|
Joined: Sep 2003
Posts: 4,883 Likes: 12
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 4,883 Likes: 12 |
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)
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
|
|
|
Joined: Oct 2004
Posts: 22
Forum Member
|
Forum Member
Joined: Oct 2004
Posts: 22 |
Thanks, Lennart. I do appreciate your help.
-Leo
|
|
|
|
|