[gmx-developers] Velocity verlet option in Gromacs with trotter-decomposition NPT

[David Mobley]

Congrats on the accomplishments! Let me know when it's starting to look like the kinks are worked out and we may try to start running with it and doing some testing on this end. 

Thanks!
David

On Dec 19, 2009, at 1:04 AM, Michael Shirts wrote:

> After much excruciating debugging, and a year or so after I thought it
> would be done, I have merged in my velocity Verlet Trotter
> decomposition based
> implementation including both Nose-Hoover chains and isotropic
> pressure control.  It was the pressure control + constraints that held
> back everything for most of that time, but it appears to be working
> perfectly now.   By working, it conserves the conserved quantity (NVE,
> NVT, NPH, NPT), with quadratic dependence on error in the timestep as
> well good as the NVT simulation.  These are not symplectic algorithms,
> so there is some possibility of drift in the conserved quantity,
> although it must average to zero.  This works best in double, of
> course -- in single, it conserves energy as well as current methods,
> though starts to have
> problems with constraints for small timesteps.
> 
> This code also preserves previous MD methods without any differences
> in the output caused by more than order of operations changes in the
> numerical precision (at least, all that I've checked -- there may be
> some combinations of options that still have bugs).  It has also been
> verified with checkpoints, and passes all the regression tests.
> Rerunning is not quite working for ekin and pressure for leapfrog with
> constraints, but it was not working before; for velocty verlet, ekin
> is correct, but pressure and virial are not.
> 
> Implementing velocity verlet required a myriad small changes in
> bookkeeping, and it's almost certain I've missed some combination of
> options that will mess things up.  Please let me know ASAP if you
> encounter any problems so I can get them resolved.
> 
> Basic documentation.  The new integrator is "md-vv".  There is another
> integrator "md-vv-avek" which I will describe later in this email. The
> pressure control option for both vv integrators is 'MTTK', the
> Martyna-Tuckerman-Tobias-Klein, with some later modifications.  This
> is all the new keywords that are added.
> 
> Additional features in this code:
> 
> I have worked out velocity verlet temperature control using  KE(t) =
> 0.5*[(KE(t+dt/2) + KE(t-dt/2)] instead of KE(t) = \sum ( v(t-dt/2) +
> v(t+dt/2))^2. Using
> average KE vs. average velocities gives worse fluctuations in the
> conserved quantity but an average temperature that depends much less
> on the timestep size.  Without temperature control, it doesn't matter
> -- you're getting the same trajectory whichever definition you use,
> but
> with temperature control it can mean that for large steps, you are not
> keeping the system at the temperature you think you are.  For small
> timesteps, this difference is not large, but if you're taking 4 fs
> timsteps, it could cause problems in terms of the temperature that you
> -think- your simulation is running at.   This is done using the
> 'md-vv-avek' integrator.  Note that without temperature coupling, this
> is numerically identical to leapfrog, other than starting at a
> different place in the timestep.
> 
> -Virial Temperature: There's a question as to exactly how much the ave
> ekin vs ave velocity matters.  As a comparison, I've implemented the
> temperature calculated from the virial, which is a purely structural
> measure of temperature.  However, it converges much more slowly than
> the ave_ekin definition of temperature, so I haven't been able to do
> full comparisons yet.  This is turned on by setting the environment
> variable "GMX_VIRIAL_TEMPERATURE".
> 
> -Nose-Hoover with chains with the the velocity verlet integrators.
> Note that right now, the number of elements in the chain is in a
> #define, so can only be changed at compile time.  It wouldn't be that
> much more difficult. Note that leapfrog does not have NH chains;
> Leapfrog uses a linearized version, whereas my implementation follows
> Martyna-Kline-Tuckerman and is not; so it would take a bit of work to
> adapt this with chains.
> 
> There are a number of additional features to implement that have not
> yet been added, but will be relatively straightforward to implement in
> a unified code framework instead of working on a code branch, and can
> be done in the order that people require them.   These include:
> 
> -Documentation of the methods -- I've been working off of a Mark
> Tuckerman preprint; I did the iteration a bit differently than he did,
> and I have not entirely finished the .tex documentation files.  I will
> get these prepared gradually; certainly for 4.1
> 
> -Nonisotropic scaling: This will be pretty trivial to put in.  It's
> mostly a matter of replacing scalars with matrices that are simple to
> compute and set up.  Will take a couple of days, but is
> straightforward, as the theory is worked out by Tuckerman et al
> already.
> 
> -Semiisotropic scaling: This would be a bit harder:  I'd have to sit
> down and do the new derivations for the trotter decomposition such
> that two dimensions are fixed, but one is floating.  It's essentially
> a remapping of a 2x2 matrix, so it can be done.  After the exact
> equations are worked out, it can be implemented in the nonisotropic
> scaling framework, so that won't require any additional coding at that
> point.  Will come after nonisotropic scaling, might take some thought.
> 
> -Changing pressure / thermostat coupling to only every N steps, thus
> requiring virial and ekin calculations, by rewriting the order of the
> Trotter decomposition. Not yet implemented, but -should- be (unless
> I'm missing something) pretty simple to do.  This should mean greatly
> improved scaling with pressure and temperature control.  You would
> have to have very long tau_t and tau_p for this to work, even with
> something fairly low like N=10; but for long simulations, long tau_t
> and tau_p are fine in any case.
> 
> -Converging the iterative equations; Right now these are controlled
> with some heuristic guidelines which may occasionally lead to lack of
> convergence -- these will likely need to be adjusted slightly.
> 
> Some more detailed coding issues:
> 
> -I had to significantly rearrange the code of the do_md loop and the
> update() function to get this to work.  I had to have iterative loops
> to solve some of the implicit equations for constraint + NPT, so I
> introduced a number of helping functions to do this.  I used a
> modified secant search instead of self-consistent iteration ion, which
> greatly improves the convergence -- in most cases, it converges in
> just 6-8 iterations for even short temperature and pressure time
> constants (for example, a tau_p=0.05 ps tau_t=0.05 ps combination).
> For 1
> ps tau_t/tau_p, it takes 4 iterations pretty much every time; with
> some fiddling, this could probably be reduced to 3 or so.
> 
> Theoretical note:
> 
> You can write either velocity verlet or leapfrog as a trotter
> decomposition; it's just a matter of where you start from.  However,
> in most cases leapfrog requires things before you calculate them, so
> doing pressure control is very difficult.  Velocity verlet requires
> them after you calculate them, which usually makes it cleaner.  I can
> explore various other combinations as time permits and people are
> interested.
> 
> Next: there are also a number of free energy features that I'd like to
> add in as soon as these integrator changes are incorporated.  These
> should be much simpler to add since I've got them going in my 3.1.4
> version already and they are pretty portable --- unlike this batch of
> changes, which is pretty heavily woven throughout a number of
> different parts of the code.
> 
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu
> (434)-243-1821
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.

David Mobley
dmobley at gmail.com
504-383-3662