CHARMM Development Project Forums User Discussion & Questions Other User Discussion and Questions Different MSD from COOR ANALYSIS and CORREL

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

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.

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).

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

The problem is calculating eta from the simulation, which I believe requires NVT conditions.

Some would argue that for transport properties you should use NVE.

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

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.

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

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.