[gmx-users] Statistical uncertainty in gromacs

[Justin Lemkul]

On 2/21/13 8:19 PM, Juliette N. wrote:
> Hi Justin,
>
> When one issues g_energy –f *.edr –nmol X –b XXX, g_energy reads the frames
> from -f edr. file which is already written each nstenergy = 1000 steps for
> instance. So g_energy is not reading all MD steps. and I guess the same
> steps are printed to .xvg as out of g_energy. Please help me realize why
> g_energy is more accurate. I am confused :( Where are these all MD steps
> written that is stated in g_energy - h
>
> Also :
>
> "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"
>
> By energy file we mean .edr file, right?
>
> Why does it say *energy file* contain averages over all MD ...than the
> number of frames in *energy file*?!
>
> I am not able to differentiate between energy file and energy file?
>

The .edr file will save specific energy terms at the interval specified by 
nstenergy, but it also saves energy history that is accumulated over all frames. 
  So in reality, .edr files still have knowledge of all MD frames, independent 
of what you set as nstenergy.  You can only output values every nstenergy frames 
for obvious reasons.

-Justin

> Please help!
>
> Sorry for naive questions!
>
> On 21 February 2013 19:23, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> 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<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<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
>>
>> ==============================**==========
>>
>> --
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users>
>> * Please search the archive at http://www.gromacs.org/**
>> Support/Mailing_Lists/Search<http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
>> * Please don't post (un)subscribe requests to the list. Use the www
>> interface or send it to gmx-users-request at gromacs.org.
>> * Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>
>>
>
>
>

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

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

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