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

[Ullmann, Thomas]

Hi,

I encountered a possible bug in simulations with the velocity-Verlet
integrator. The cause is not yet clear to me and I hope that somebody
on the list might have an idea.

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.
The effect is also much less pronounced when using Stockholm-lipids
and Amber99SB-ILDN as force fields. The same behavior can be observed
for slightly different setups of the system (protonation, phosphorylation)
and for a second, different system of similar size. I repeated the
simulations with the leap-frog/md integrator which does not show this
effect. The system also quickly recovers and restores the normal
membrane area within 10s of ns if a simulation began with the
velocity-Verlet integrator is continued with the leap-frog integrator.

Plots of the box edge length in in the membrane plane for
the two integrators for the same starting states can be found here:

http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/membrane_protein/Box-X.equi.md-vv.eps
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/membrane_protein/Box-X.equi.md.eps

The 15 different curves correspond to slightly different setups
of the same system in terms of protonation and phosphorylation.

This plot:
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/membrane_protein/Epot.equi.md-vv_vs_md.eps
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.


I already made some initial attempts of debugging but have the uneasy
feeling that the issue might not be simple.
I started from the most recent commit of the git branch release-4-6 as
of yesterday compiled in double precision. Tests were run in serial.
As test system, I chose a very simple toy model which should be solvable
analytically. I placed two oppositely charged ions in a distance of
two nanometers in vacuum. 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.
The two integrators yield the same results if no thermostat is used but
markedly different results if the velocity-rescaling thermostat is added.
The same random number seed for the thermostat is used in both runs
(.mdp option ld_seed according to the manual). With this setup, I would
have expected nearly identical trajectories for both integrators.
Did I overlook something?

Plots of the kinetic and potential energies for these runs can be found here:
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/2_ions/2_ions.Ekin.eps
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/2_ions/2_ions.Epot.eps

The toy model with the raw output files can be found here:
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/2_ions/2_ions.tar.bz2


For the moment I set up the simulations with the leap-frog integrator,
but of course it would still be nice to understand the problem and
have the bug fixed, if there is one. Redmine does not yet seem to have
the exact same problem documented and most issues with the velocity-
Verlet integrator are already older.

Any ideas would be highly appreciated ...

Best,
Thomas.

--------------------------------------------------------------------------------------------------------------------------
R. Thomas Ullmann, PhD
Theoretical & Computational Biophysics
Max Planck Institute for Biophysical Chemistry
Am Fassberg 11
37077 Göttingen, Germany
thomas.ullmann at mpibpc.mpg.de
www.bisb.uni-bayreuth.de/People/ullmannt
--------------------------------------------------------------------------------------------------------------------------

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