[gmx-developers] possible bug connected to the velocity-Verlet integrator

[Shirts, Michael R. (mrs5pt)]

Hi, Thomas-

One gold standard test Is the test published here: http://pubs.acs.org/doi/abs/10.1021/ct300688p.

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.

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.

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.

> The initial velocities are zero.
> I used a very small time step of 1 attosecond (.mdp option dt = 0.000001)
> to minimize differences in the integration accuracy of the two integrators.
> Each run comprises 10,000 steps. The initial velocities, coordinates and
> forces are verified to be identical for all runs.

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

> The simulated system is a membrane protein in a POPE
> bilayer with water and ions in an NpT ensemble. The symptom of the
> possible bug is a repeated large decrease and subsequent decrease of
> the membrane area from the normal value of ca 12.5 nm squared up to
> ca. 23 nm squared and back.The protein conformation is not visibly affected.

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.

> The effect is also much less pronounced when using Stockholm-lipids
> and Amber99SB-ILDN as force fields.

Hmm. That doesn't make that much sense.  It sounds like there's a lot of statistical non-reproducibility.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant 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/20140701/95cbe562/attachment-0001.html>