[gmx-users] g_msd uncertainty

[Ivan Gladich]

Dear Gromacs community
          I would like to better understand how the uncertainty on the 
diffusion coefficient is calculated in GROMACS.
I am calculating water diffusivity for different water models at 
different temperatures.
I am using Gomacs 4.5.4


Reading the manual, I am a bit confused on how g_msd calculate the 
uncertainty on the average.
Indeed, the manual told me (pag 285-286)

1)

"g_msd computes the mean square displacement (MSD) of atoms from a set 
of initial positions. This
provides an easy way to compute the diffusion constant using the 
Einstein relation. The time between the
reference points for the MSD calculation is set with -trestart. The 
diffusion constant is calculated by
least squares fitting a straight line (D*t + c) through the MSD(t) from 
-beginfit to -endfit (note that t
is time from the reference positions, not simulation time). An error 
estimate given, which is the difference
of the diffusion coefficients obtained from fits over the two halves of 
the fit interval."


2) By the way, few lines below


"The diffusion coefficient is determined by linear regression of the MSD,
where, unlike for the normal output of D, the times are weighted 
according to
the number of reference points, i.e. short times have a higher weight. Also
when -beginfit=-1,fitting starts at 10% and when -endfit=-1, fitting goes to
90%.Using this option one also gets an accurate error estimate based on the
statistics between individual molecules. Note that this diffusion 
coefficient
and error estimate are only accurate when the MSD is completely linear
between -beginfit and -endfit."

Based on this last paragraph, I suppose that g_msd calculates for each 
molecule the diffusion coefficient and then it calculates the average 
and standard deviation from the daset of the diffusion coefficients of 
each molecules.

So, if I compute

g_msd -f traj.xtc -s topol.tpr -n index.ndx.ndx -o msd.xvg -b 0 -e 
100000 -beginfit 15000 -endfit 25000 -trestart 10 -type z

at the top of my output file msd.xvg I can read

# MSD gathered over 100000 ps with 10001 restarts
# Diffusion constants fitted from time 1500 to 25000 ps
# D[all_100000] = 0.1230 (+/- 0.0060) (1e-5 cm^2/s)


So, 0.0060 is one standard deviation based on the statistic between 
different molecules (second method) or is calculated using the 
difference of the diffusion coefficient on the two halves of the fitting 
region (first method)?

I was searching in the mail list without success...thanks for any help

Ivan

-- 
------
Ivan Gladich, Ph.D.
Postdoctoral Fellow
Academy of Sciences of the Czech Republic
Institute of Organic Chemistry and Biochemistry AS CR, v.v.i.
Flemingovo nám. 2.
166 10 Praha 6
Czech Republic

Tel: +420775504164
e-mail: ivan.gladich at uochb.cas.cz
web page:http://www.molecular.cz/~gladich/
-----

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110621/0036e96b/attachment.html>