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

[Michael Shirts]

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