[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat

[Shirts, Michael R. (mrs5pt)]

> Obviously we agree on the dubious nature of the linear drift and that its origin from reduced precision round-off errors is doubtful.

I'm not sure I said this.  I thought you said that drift was zero in double precision, but present in single precision.

Since gromacs uses essentially the same code for single and double precision (the type 'real' is used, which is defined as single and double precision at compile time, except in some cases where double is always used), then it sounded like the reason was indeed round-off errors.

The question is how much NVT or NPT ensembles are affected by this lack of precision (obviously, you can run long timescale NVE without energy conservation).  That's an open question; for many problems, this source of error is much smaller than other sources of error, but there could be exceptions.

If I misunderstood and you saw drift in double precision, that's another story.

Best,
~~~~~~~~~~~~
Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434) 243-1821

From: Bernhard <b.reuter at uni-kassel.de<mailto:b.reuter at uni-kassel.de>>
Date: Tuesday, July 21, 2015 at 7:40 AM
To: "gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>" <gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>>, "michael.shirts at virginia.edu<mailto:michael.shirts at virginia.edu>" <michael.shirts at virginia.edu<mailto:michael.shirts at virginia.edu>>
Subject: Re: [gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat

Dear Michael,

thank you very much for your very helpfull answer.

In my opinion the occurence of a linear energy drift of this size could indicate a bug in the program.
So I startet a more rigorous investigation and would like to share some preliminary results:

Graph of Verlet buffer-size vs. energy-drift size for single and double precision:
https://www.dropbox.com/s/e56916inlm0ym48/buffersize-vs-total-energy-drift.png?dl=0

Graph of the linear total energy drift for a buffer-size of 0.02nm:
https://www.dropbox.com/s/ifq3v3jwfzn4goh/4_nve_100ps_rlist-1.02_Total-Energy.png?dl=0

Exemplary mdp and log files for a buffer-size of 0.02nm:
https://www.dropbox.com/s/peamt27d2exhclc/4_nve_100ps_rlist-1.02.mdp?dl=0
https://www.dropbox.com/s/70ictqb6nbs94wk/4_nve_100ps_rlist-1.02.log?dl=0

The investigated protein-water system consists of 22765 atoms, the AMBER99SB-ILDN force field with TIP3P water was used, all simulations were in the NVE-ensemble, GROMACS 4.6.7 was used.

I will repeat the tests with GROMACS 5.0 and for different test systems (i.e. the lysozyme system from the popular lysozyme tutorial of Justin Lemkul).

Best,
Bernhard

Am 20/07/15 um 00:50 schrieb Shirts, Michael R. (mrs5pt):
> Do I have to switch to double precision if I care about energy conservation, integrator symplecticity, phase space volume conservation and ergodicity?

This sounds a like a good idea.  If you are doing tests where this matters, use double precision. Sounds like the safest.

> Since bigger round-off errors by reduced precision shouldn't accumulate linearly but at worst with Sqrt(N): Shouldn't one be worried about the occurence of a linear systematic error by only changing the precision from double to single in a calculation?

Reduced precision errors only would be linear if the errors are uncorrelated, but it's not clear to me why roundoff errors would be uncorrelated.

> But if you have a constant downward drift of energy you must consider that there is less phase space volume at lower energies - so there is no volume conservation in phase space.

Correct, for NVE.  For NVT,  the conserved energy is a bookkeeping number, it has nothing to do with the current phase space of the system. The thermostat is pumping in more energy so that the kinetic energy remains consistent with the desired temperature.  We then actually have a steady state system, rather than an equilibrium system.  The question is, how different is this distribution from the true equilibrium distribution?

This is generally testable.  For thermodynamic  calculations (which is what one presumably is intrested in with a thermostat, rather than the dynamics ), what really matters is 1) whether the correct distribution is obtained within noise and 2) whether the sampling is ergodic.  2) is very hard to answer, but 1) can be checked by

https://github.com/shirtsgroup/checkensemble

With the theory described here:

http://dx.doi.org/10.1021/ct300688p

Gromacs in single precisions seems to behave fine statistically for systems of a few hundred atoms.

I suspect that there are subtle phenomena where the lack of exact symplecticness matters.  I also believe from my testing (no full paper on this) that there aren't very many that occur in highly chaotic systems with hundreds of particles at NIT.

I bet there are cases with just a few particles where the problems could become very obvious, however.

Best,
~~~~~~~~~~~~
Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu<mailto:michael.shirts at virginia.edu>
(434) 243-1821












-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20150721/f50ca8f4/attachment-0001.html>