<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">
<div style="direction: ltr;font-family: Courier New;color: #000000;font-size: 10pt;">
Hi,<br>
<br>
I encountered a possible bug in simulations with the velocity-Verlet<br>
integrator. The cause is not yet clear to me and I hope that somebody<br>
on the list might have an idea.<br>
<br>
The simulated system is a membrane protein in a POPE<br>
bilayer with water and ions in an NpT ensemble. The symptom of the<br>
possible bug is a repeated large decrease and subsequent decrease of<br>
the membrane area from the normal value of ca 12.5 nm squared up to<br>
ca. 23 nm squared and back.The protein conformation is not visibly affected.<br>
The effect is also much less pronounced when using Stockholm-lipids<br>
and Amber99SB-ILDN as force fields. The same behavior can be observed<br>
for slightly different setups of the system (protonation, phosphorylation)<br>
and for a second, different system of similar size. I repeated the<br>
simulations with the leap-frog/md integrator which does not show this<br>
effect. The system also quickly recovers and restores the normal<br>
membrane area within 10s of ns if a simulation began with the<br>
velocity-Verlet integrator is continued with the leap-frog integrator.<br>
<br>
Plots of the box edge length in in the membrane plane for<br>
the two integrators for the same starting states can be found here:<br>
<br>
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/membrane_protein/Box-X.equi.md-vv.eps<br>
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/membrane_protein/Box-X.equi.md.eps<br>
<br>
The 15 different curves correspond to slightly different setups<br>
of the same system in terms of protonation and phosphorylation.<br>
<br>
This plot: <br>
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/membrane_protein/Epot.equi.md-vv_vs_md.eps<br>
shows that the potential energy is rather constant in both cases but<br>
at a markedly higher value in case of the velocity-Verlet integrator.<br>
Temperature and pressure are very stable / only show normal fluctuations<br>
around the preset values.<br>
<br>
<br>
I already made some initial attempts of debugging but have the uneasy<br>
feeling that the issue might not be simple.<br>
I started from the most recent commit of the git branch release-4-6 as<br>
of yesterday compiled in double precision. Tests were run in serial.<br>
As test system, I chose a very simple toy model which should be solvable<br>
analytically. I placed two oppositely charged ions in a distance of<br>
two nanometers in vacuum. The initial velocities are zero.<br>
I used a very small time step of 1 attosecond (.mdp option dt = 0.000001)<br>
to minimize differences in the integration accuracy of the two integrators.<br>
Each run comprises 10,000 steps. The initial velocities, coordinates and<br>
forces are verified to be identical for all runs.<br>
The two integrators yield the same results if no thermostat is used but<br>
markedly different results if the velocity-rescaling thermostat is added.<br>
The same random number seed for the thermostat is used in both runs<br>
(.mdp option ld_seed according to the manual). With this setup, I would<br>
have expected nearly identical trajectories for both integrators.<br>
Did I overlook something?<br>
<br>
Plots of the kinetic and potential energies for these runs can be found here:<br>
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/2_ions/2_ions.Ekin.eps<br>
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/2_ions/2_ions.Epot.eps<br>
<br>
The toy model with the raw output files can be found here:<br>
http://wwwuser.gwdg.de/~tullman/possible_md-vv_bug/2_ions/2_ions.tar.bz2<br>
<br>
<br>
For the moment I set up the simulations with the leap-frog integrator,<br>
but of course it would still be nice to understand the problem and<br>
have the bug fixed, if there is one. Redmine does not yet seem to have<br>
the exact same problem documented and most issues with the velocity-<br>
Verlet integrator are already older.<br>
<br>
Any ideas would be highly appreciated ...<br>
<br>
<div>Best,<br>
Thomas.<br>
<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"><br>
</div>
</div>
</body>
</html>