Previous Thread
Next Thread
Print Thread
Different MSD from COOR ANALYSIS and CORREL
#34584 11/23/14 02:32 AM
Joined: Mar 2007
Posts: 148
J
jeffrey Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Mar 2007
Posts: 148
Hi all,

I performed a simulation of glucose solution (42 glucoses in a TIP3 water box) and want to compute diffusion coefficient of glucose molecules. For diffusion coefficient calculation, I followed the two methods proposed by Rick and Lennart in a previous relevant thread (http://www.charmm.org/ubbthreads/ubbthreads.php?ubb=showflat&Number=23584&page=1). However, I got different slope values from COOR ANALYSIS and CORREL. Here are the main commands in the input script:
-----------coor analysis
crystal define cubic @xlat @xlat @xlat 90.00 90.00 90.00
crystal build Noper 0 cutoff @xcut
image byres xcen 0.0 ycen 0.0 zcen 0.0 sele all end
traj firstu 51 nunit 1 begin 101000 skip 1000 stop 20000000
open unit 26 write form name O5-all.msd
coor anal select type o5 end - ! Glucose oxygens
firstu 51 nunit 1 nskip 1000 begin 101000 stop 20000000 -
imsd 26 -
rspin 0.0 rspout 999.9 -
ncors 500 -
xbox @xlat ybox @xlat zbox @xlat
-----------

----------correl
open unit 8 write unform name merge.dcd
OPEN READ UNIT 9 FILE NAME ../glucose.dcd
merge unfold xfluct firstu 9 nunit 1 outputu 8
close unit 8

open unit 8 write unform name merge.dcd

correl maxtime 20000000 maxatoms 180 maxseries 180
enter xy1 atom xyz mass sele segid aglc .and. resid 1 .and. (type o5) end
enter xy2 atom xyz mass sele segid aglc .and. resid 2 .and. (type o5) end
enter xy3 atom xyz mass sele segid aglc .and. resid 3 .and. (type o5) end
.
.
.
enter xy42 atom xyz mass sele segid aglc .and. resid 42 .and. (type o5) end
traj firstu 8 nunit 1 begin 101000 skip 1000 stop 20000000
corfun xy1 xy1 difference
write corr unit 11 dumb time
corfun xy2 xy2 difference
write corr unit 12 dumb time
corfun xy3 xy3 difference
write corr unit 13 dumb time
.
.
.
corfun xy42 xy42 difference
write corr unit 52 dumb time
----------average was computed from the 42 time series to get the correl-MSD.

Attached pls find the MSD-t plot (bold red: MSD-t from COOR ANAL; bold black: average MSD-t from CORREL; other lines: different time series from CORREL).
I suppose the two methods can get same results.Could you pls have a check of my input to see why I got very different slopes?

Thanks for your time.

Jeff

Attached Files MSD.png
Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34585 11/23/14 11:55 AM
Joined: Sep 2003
Posts: 4,771
Likes: 2
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 4,771
Likes: 2
correl may not be PBC-aware. You can test this by removing all PBCrelated stuff (CRYSTAL, IMAGE, XBOX ...) and run again. The correl results should be unchanged and the coor analysis should be similar (identical?) to the correl results.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34586 11/23/14 04:30 PM
Joined: Sep 2003
Posts: 8,447
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,447
You should remove the XFLUCT from the MERGE COOR UNFOLD (the UNFOLD is needed because CORREL is not PBC-aware).

The script appears incomplete; I don't see a command which opens UNIT 8 for reading.

To calculate diffusion properly, you should use the molecular center of mass, and not a single atom.

See also J. Phys. Chem. B, 114, 12501-12507 (2010).


Rick Venable
computational chemist

Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34587 11/23/14 11:00 PM
Joined: Mar 2007
Posts: 148
J
jeffrey Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Mar 2007
Posts: 148
Thanks for all comments and suggestions. Here are some new results from different tests.
1. trajectory merge without XFLUCT generates exactly the same results as COOR ANALYSIS;
2. without PBC related stuff, CHARMM crashed at the traj merge command with a message "***** CRYStal required for XFLUct, UNFOld".
3. the unit 8 for merged dcd should be opened by the command"
open read unit 8 file name merge.dcd
". Not the one I listed in the original post.

After read the paper recommended by Rick, I wonder if NPT ensemble can be used in calculation of diffusion coefficient using eq. (2) of the JPCb paper, D = D_pbc + K_b*T*xi/(6*PI*L*eta), with xi=2.837297, eta-the viscosity of the solution, L-the length of the cubic box in simulation. Under NPT, can L be averaged using all frames recorded for the calculation?

Many thanks.

Jeff

Attached Files msd-t-new.png
Last edited by jeffrey; 11/23/14 11:02 PM.
Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34588 11/24/14 03:30 AM
Joined: Sep 2003
Posts: 8,447
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,447
The problem is calculating eta from the simulation, which I believe requires NVT conditions.


Rick Venable
computational chemist

Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34589 11/24/14 10:12 AM
Joined: Sep 2003
Posts: 4,771
Likes: 2
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 4,771
Likes: 2
Some would argue that for transport properties you should use NVE.

Plotting D_pbc vs 1/L should give you both eta and D.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34591 11/24/14 05:10 PM
Joined: Sep 2003
Posts: 8,447
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,447
One can also compute eta directly from a single simulation via the off-diagonal elements of the pressure tensor, which is the approach we've used.


Rick Venable
computational chemist

Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34592 11/24/14 06:42 PM
Joined: Mar 2007
Posts: 148
J
jeffrey Offline OP
Forum Member
OP Offline
Forum Member
J
Joined: Mar 2007
Posts: 148
Thanks for the suggestions. I will try different schemes to see if the results are consistent. If using D_pbc vs 1/L to compute D and eta, is there any special reason for the NVE ensemble? Or can NVT be adopted?

Have a nice day.

Jeff

Re: Different MSD from COOR ANALYSIS and CORREL
jeffrey #34595 11/26/14 09:05 PM
Joined: Sep 2003
Posts: 4,771
Likes: 2
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 4,771
Likes: 2
Once your system is equilibrated there is little reason to use anything but NVE; this is system (and program) dependent, and e.g. membranes may be different. Using an integrator that interferes with the velocities can affect transport properties.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Moderated by  chmgr, John Legato, rmv 

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

PHP: 5.6.33-0+deb8u1 Page Time: 0.014s Queries: 34 (0.007s) Memory: 0.9640 MB (Peak: 1.0921 MB) Data Comp: Off Server Time: 2020-05-28 18:15:17 UTC
Valid HTML 5 and Valid CSS