<html dir="ltr">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style id="owaParaStyle" type="text/css">P {margin-top:0;margin-bottom:0;}</style>
</head>
<body ocsi="0" fpstyle="1" style="word-wrap:break-word; font-family:Calibri,sans-serif; font-size:14px; color:rgb(0,0,0)">
<div style="direction: ltr;font-family: Courier New;color: #000000;font-size: 10pt;">
Hi Michael,<br>
<br>
I fully agree that the actual system is not suitable for testing, but it lead me to think that there might be a bug. Thanks for the paper, I already saw it, but I think that it is obvious that the large systems do not sample the correct ensemble.<br>
<br>
How about the second, simpler test described at the end of the mail (two particles only)?<br>
<br>
While it certainly is not physically too sound to thermostat the single moving particle alone, I thought it should be usable to test the technical function. The small time step should also make the results independent of the slightly different interpretation
of the input velocities, as verified without thermostat).<br>
<br>
Best,<br>
Thomas.<br>
<div><br>
<div style="font-family:Tahoma; font-size:13px">
<div style="font-family:Tahoma; font-size:13px">--------------------------------------------------------------------------------------------------------------------------<br>
R. Thomas Ullmann, PhD<br>
Theoretical & Computational Biophysics<br>
Max Planck Institute for Biophysical Chemistry<br>
Am Fassberg 11<br>
37077 Göttingen, Germany<br>
thomas.ullmann@mpibpc.mpg.de<br>
www.bisb.uni-bayreuth.de/People/ullmannt<br>
--------------------------------------------------------------------------------------------------------------------------<br>
</div>
</div>
</div>
<div style="font-family: Times New Roman; color: #000000; font-size: 16px">
<hr tabindex="-1">
<div style="direction: ltr;" id="divRpF818830"><font color="#000000" face="Tahoma" size="2"><b>From:</b> gromacs.org_gmx-developers-bounces@maillist.sys.kth.se [gromacs.org_gmx-developers-bounces@maillist.sys.kth.se] on behalf of Shirts, Michael R. (mrs5pt)
[mrs5pt@eservices.virginia.edu]<br>
<b>Sent:</b> Tuesday, July 01, 2014 6:24 AM<br>
<b>To:</b> gmx-developers@gromacs.org<br>
<b>Subject:</b> Re: [gmx-developers] possible bug connected to the velocity-Verlet integrator<br>
</font><br>
</div>
<div></div>
<div>
<div>Hi, Thomas-</div>
<div><br>
</div>
<div>One gold standard test Is the test published here: <a href="http://pubs.acs.org/doi/abs/10.1021/ct300688p" target="_blank">http://pubs.acs.org/doi/abs/10.1021/ct300688p</a>.</div>
<div><br>
</div>
<div>For 4.6.3, both leapfrog and velocity verlet worked for md/parrinello-rahman and md-vv/Martina-Tuckerman-Klein integrators, but there is certainly a possibility that there's some combination of options that could easily cause a problem that would not have
the right ensemble.</div>
<div><br>
</div>
<blockquote id="MAC_OUTLOOK_ATTRIBUTION_BLOCKQUOTE" style="border-left:#b5c4df 5 solid; padding:0 0 0 5; margin:0 0 0 5">
<div>shows that the potential energy is rather constant in both cases but at a markedly higher value in case of the velocity-Verlet integrator. Temperature and pressure are very stable / only show normal fluctuations around the preset values.</div>
</blockquote>
<div><br>
</div>
<div>For the test case, the difference is almost certainly because the same input velocities are interpreted differently depending on whether they are velocity verlet or leapfrog. The same input velocities will be interpreted in different ways. either as full
step or half-step velocities. If you run at with any temperature control, or even with generated velocities, the average potential energies should be much closer. So I don't think you are testing whatever issue is occurring in the membrane simulation by this
test. </div>
<div><br>
</div>
<div>
<div>> The initial velocities are zero.</div>
<div>> I used a very small time step of 1 attosecond (.mdp option dt = 0.000001)</div>
<div>> to minimize differences in the integration accuracy of the two integrators.</div>
<div>> Each run comprises 10,000 steps. The initial velocities, coordinates and</div>
<div>> forces are verified to be identical for all runs.</div>
</div>
<div><br>
</div>
<div>Running < 1 time period for any thermostat means that one is not really testing the thermostat; thermostats only really give the correct kinetic energy distribution if averaged over multiple time periods </div>
<div><br>
</div>
<div>> The simulated system is a membrane protein in a POPE</div>
<div>> bilayer with water and ions in an NpT ensemble. The symptom of the</div>
<div>> possible bug is a repeated large decrease and subsequent decrease of</div>
<div>> the membrane area from the normal value of ca 12.5 nm squared up to</div>
<div>
<div>> ca. 23 nm squared and back.The protein conformation is not visibly affected.</div>
</div>
<div><br>
</div>
<div>What pressure control algorithm is being used? Were any warnings printed in those logs? If I were to hazard a guess, it would be the interaction with pressure control and the integrator. </div>
<div>
<div><br>
</div>
<div>> The effect is also much less pronounced when using Stockholm-lipids</div>
<div>> and Amber99SB-ILDN as force fields.</div>
</div>
<div><br>
</div>
<div>Hmm. That doesn't make that much sense. It sounds like there's a lot of statistical non-reproducibility.</div>
<div><br>
</div>
<div>Best,</div>
<div>~~~~~~~~~~~~</div>
<div>Michael Shirts</div>
<div>Assistant Professor</div>
<div>Department of Chemical Engineering</div>
<div>University of Virginia</div>
<div><a href="mailto:michael.shirts@virginia.edu" target="_blank">michael.shirts@virginia.edu</a></div>
<div>(434)-243-1821</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
</div>
</div>
</div>
</body>
</html>