[gmx-users] Statistical uncertainty in gromacs

[Justin Lemkul]

On 2/21/13 6:58 PM, Juliette N. wrote:
> Hi Matthew,
>
> Thanks for your reply. I tried g_analyze as you suggested:
>
> 1) I am wondering why the average given by g_energy and g_analyze are not
> identical. I tried the following:
>
> g_energy -f  Potential. edr -o Potential.xvg and extracted the average then
> provided Potential.xvg as input to g_analyze ( I did this since g_analyze
> seems to read xvg files only). and issued: g_analyze –f.Potential.xvg–av.
>
> However the averages are not identical.
>
>
> 2) Also I tried g_energy -f  Potential. edr -o Potential.xvg and from
> g_energy I can see  Energy                      Average   *Err.Est.*
> RMSD  Tot-Drift
>
>
> are given by default. Then tried g_analyze -ee and noticed that the *
> Err.Est.*  reported by g_energy and the one from g_analyze -ee are not
> identical!

See the message below, printed by g_energy.

>
> Could you please explain why is this happening and how  *Err.Est.* is
> estimated? .
>
> g_energy says:
> An error estimate of the average is given based on a block averages over 5
> blocks using the full-precision averages. The error estimate can be
> performed over multiple block lengths with the options -nbmin and -nbmax. *

The following section explains the discrepancy in your results:

> Note* that in most cases the energy files contains averages over all MD
> steps, or over many more points than the number of frames in energy file.
> This makes the g_energy statistics output more accurate than the
> .xvg<http://manual.gromacs.org/online/xvg.html>output
>

Thus, the only way to achieve absolute agreement between .xvg and .edr analysis 
is to set nstenergy = 1, which will result in really big .edr files, and the 
difference in practice should be pretty small.

> I cant seem to understand how error is obtained.:( If I am reading the
> frames from say , -b 4000 -5000, what happens if I dont provide blocks
> etc...
>

Well, there are default values that are used for -nbmin and -nbmax, so if you 
don't specify any values manually, those are used.  The default is 5.  Without 
looking into the code, it would seem that block averaging of some sort is being 
done.  The error estimate for g_analyze -ee is documented more thoroughly and is 
described in the paper mentioned in the help description.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================