[gmx-users] g_msd uncertainty
Ivan Gladich
ivan.gladich at marge.uochb.cas.cz
Tue Jun 21 14:30:53 CEST 2011
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>
More information about the gromacs.org_gmx-users
mailing list