<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
</head>
<body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; color: rgb(0, 0, 0); font-size: 14px; font-family: Calibri, sans-serif;">
<div>
<div>> Obviously we agree on the dubious nature of the linear drift and that its origin from reduced precision round-off errors is doubtful.</div>
<div><br>
</div>
<div>I'm not sure I said this. I thought you said that drift was zero in double precision, but present in single precision. </div>
<div><br>
</div>
<div>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. </div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>
<div>If I misunderstood and you saw drift in double precision, that's another story.</div>
<div><br>
</div>
<div>Best,</div>
~~~~~~~~~~~~<br>
Michael Shirts<br>
Associate Professor<br>
Department of Chemical Engineering<br>
University of Virginia<br>
<a href="michael.shirts@virginia.edu">michael.shirts@virginia.edu</a><br>
(434) 243-1821</div>
</div>
<div><br>
</div>
<span id="OLK_SRC_BODY_SECTION">
<div style="font-family:Calibri; font-size:14pt; text-align:left; color:black; BORDER-BOTTOM: medium none; BORDER-LEFT: medium none; PADDING-BOTTOM: 0in; PADDING-LEFT: 0in; PADDING-RIGHT: 0in; BORDER-TOP: #b5c4df 1pt solid; BORDER-RIGHT: medium none; PADDING-TOP: 3pt">
<span style="font-weight:bold">From: </span>Bernhard <<a href="mailto:b.reuter@uni-kassel.de">b.reuter@uni-kassel.de</a>><br>
<span style="font-weight:bold">Date: </span>Tuesday, July 21, 2015 at 7:40 AM<br>
<span style="font-weight:bold">To: </span>"<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>" <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>>, "<a href="mailto:michael.shirts@virginia.edu">michael.shirts@virginia.edu</a>"
<<a href="mailto:michael.shirts@virginia.edu">michael.shirts@virginia.edu</a>><br>
<span style="font-weight:bold">Subject: </span>Re: [gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat<br>
</div>
<div><br>
</div>
<div>
<div bgcolor="#FFFFFF" text="#000000">Dear Michael,<br>
<br>
thank you very much for your very helpfull answer.<br>
<br>
In my opinion the occurence of a linear energy drift of this size could indicate a bug in the program.<br>
So I startet a more rigorous investigation and would like to share some preliminary results:<br>
<br>
Graph of Verlet buffer-size vs. energy-drift size for single and double precision:<br>
<a class="moz-txt-link-freetext" href="https://www.dropbox.com/s/e56916inlm0ym48/buffersize-vs-total-energy-drift.png?dl=0">https://www.dropbox.com/s/e56916inlm0ym48/buffersize-vs-total-energy-drift.png?dl=0</a><br>
<br>
Graph of the linear total energy drift for a buffer-size of 0.02nm:<br>
<a class="moz-txt-link-freetext" href="https://www.dropbox.com/s/ifq3v3jwfzn4goh/4_nve_100ps_rlist-1.02_Total-Energy.png?dl=0">https://www.dropbox.com/s/ifq3v3jwfzn4goh/4_nve_100ps_rlist-1.02_Total-Energy.png?dl=0</a><br>
<br>
Exemplary mdp and log files for a buffer-size of 0.02nm:<br>
<a class="moz-txt-link-freetext" href="https://www.dropbox.com/s/peamt27d2exhclc/4_nve_100ps_rlist-1.02.mdp?dl=0">https://www.dropbox.com/s/peamt27d2exhclc/4_nve_100ps_rlist-1.02.mdp?dl=0</a><br>
<a class="moz-txt-link-freetext" href="https://www.dropbox.com/s/70ictqb6nbs94wk/4_nve_100ps_rlist-1.02.log?dl=0">https://www.dropbox.com/s/70ictqb6nbs94wk/4_nve_100ps_rlist-1.02.log?dl=0</a><br>
<br>
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.<br>
<br>
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).<br>
<br>
Best,<br>
Bernhard<br>
<br>
<div class="moz-cite-prefix">Am 20/07/15 um 00:50 schrieb Shirts, Michael R. (mrs5pt):<br>
</div>
<blockquote cite="mid:D1D19173.6F07F%25mrs5pt@eservices.virginia.edu" type="cite">
<div>
<div>
<div>> Do I have to switch to double precision if I care about energy conservation, integrator symplecticity, phase space volume conservation and ergodicity?</div>
</div>
<div><br>
</div>
<div>This sounds a like a good idea. If you are doing tests where this matters, use double precision. Sounds like the safest. </div>
<div><br>
</div>
<div>> 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?</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>> 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.</div>
<div><br>
</div>
<div>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? </div>
<div><br>
</div>
<div>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 </div>
</div>
<div><br>
</div>
<div><a moz-do-not-send="true" href="https://github.com/shirtsgroup/checkensemble">https://github.com/shirtsgroup/checkensemble</a></div>
<div><br>
</div>
<div>With the theory described here:</div>
<div><br>
</div>
<div><a moz-do-not-send="true" href="http://dx.doi.org/10.1021/ct300688p">http://dx.doi.org/10.1021/ct300688p</a></div>
<div><br>
</div>
<div>Gromacs in single precisions seems to behave fine statistically for systems of a few hundred atoms.</div>
<div><br>
</div>
<div>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. </div>
<div><br>
</div>
<div>I bet there are cases with just a few particles where the problems could become very obvious, however.</div>
<div><br>
</div>
<div>Best,</div>
<div>~~~~~~~~~~~~</div>
<div>Michael Shirts</div>
<div>Associate Professor</div>
<div>Department of Chemical Engineering</div>
<div>University of Virginia</div>
<div><a moz-do-not-send="true" href="mailto:michael.shirts@virginia.edu">michael.shirts@virginia.edu</a></div>
<div>(434) 243-1821</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset> <br>
</blockquote>
<br>
</div>
</div>
</span>
</body>
</html>