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

[Shirts, Michael R. (mrs5pt)]

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.

Sure, but in what way?  For example, is the right ensemble sampled in NVT?  That would allow one to separate out issues with just NPT.

How about the second, simpler test described at the end of the mail (two particles only)?

That's almost certainly not a good test of anything integrator related.  You can't test ensembles when looking at the approach to equilibrium.  I look at those two graphs and see perfectly valid approaches to equilibrium, especially if vv is using a different random number seed in each. Half step velocities will not be interpreted the same way by both problems, so looking at the first picosecond or so will likely not yield any useful information; it will be entirely equilibration artifacts from the different starting points. There's also a real ergodicity problem, in that if COM motion is being removed, there are only a few degrees of freedom for a symmetric particle system, none of which can exchange kinetic energy with each other, making all thermostats other than langevin suspect.  A nice little box of water or LJ fluid would be much preferable a test, especially since it would allow separate tests of NVT and NPT (in fact, I suspect an issue with non isotropic pressure control and + md-vv)

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).

I'm not sure how being dominated by different initial conditions really helps anything. To me, it looks like several valid approaches to equilibrium, with different velocity verlet seeds.

Overall:
My guess is that you are using md-vv with Parinello-Rahnman in semiisotropic mode, which I think is supposed to have a warning, but perhaps it got missed.   Perhaps if we saw the mdp's of the system that was clearly having area problems, it would be possible to tell what what combination of options might be causing problems.  If you think it was md vs md-vv, test a box of water or LJ particles over 50 or so times tau_t, and see if they have the same thermodynamics.

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