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

[Michael Shirts]

In the version checked in last night, the Nose-Hoover chains had a
fixed length defined at compile time.  Now, the number of particles in
the Nose-Hoover chain can be set in the mdp file, using the keyword
"nhchains".  The default is currently 5, but can be changed by a
#define in include/types/state.h

I noticed a slight problem with the first frame of the checkpoint file
for velocity verlet -- the trajectory was correct, but the pressure
and virial in the first frame after a restart was incorrect.

Also, because the output for Nose-Hoover chains can be significant if
a large number are selected, the positions and velocities of the
chains will now only be written to the edr file if
GMX_NOSEHOOVER_CHAINS is defined in the environment.



On Sat, Dec 19, 2009 at 2:04 AM, Michael Shirts
<michael.shirts at virginia.edu> 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
>