Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
#403 10/23/03 10:51 AM
Joined: Sep 2003
Posts: 4,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,881
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
C
Forum Member
Offline
Forum Member
C
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,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,881
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
dmr Offline
Forum Member
Offline
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,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,881
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
dmr Offline
Forum Member
Offline
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,650
Likes: 26
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,650
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
L
leo Offline
Forum Member
Offline
Forum Member
L
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,881
Likes: 12
lennart Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 4,881
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
L
leo Offline
Forum Member
Offline
Forum Member
L
Joined: Oct 2004
Posts: 22
Thanks, Lennart. I do appreciate your help.

-Leo

Page 1 of 2 1 2

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u3 Page Time: 0.009s Queries: 35 (0.006s) Memory: 0.7821 MB (Peak: 0.8763 MB) Data Comp: Off Server Time: 2023-06-05 19:32:32 UTC
Valid HTML 5 and Valid CSS