[gmx-developers] Velocity verlet version of Gromacs with trotter-decomposition NPT

[Michael Shirts]

Hi, all-

After much excruciating debugging, and 10 months or so after I thought
it would be done, I have a velocity Verlet Trotter decomposition based
implementation including both Nose-Hoover chains and isotropic
pressure control, all of which work with constraints.  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.  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.

It 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. It is
synchronized with the main git distribution as of September 25, 2009.
I have not played around with multiprocessor jobs yet, however.
Rerunning does not quite work for the ekin and pressures in the
trotter-decomposition NPT version because of the integration of the
Nose-hoover chains; that might take some extra work to sort out, but
potential energies are fine, as is rerunning without NPT.

The code is in my public repository at:
ssh://git.gromacs.org/~mrshirts/public_git/gromacs_mrs.git

It would be very good to have Berk look over it and see how much can
be ported over as is (hopefully most of it!)  There's a few code
architecture decisions that he would likely have some better ideas on;
and I'm certain there are a few bugs that will pop up.  I'm personally
willing to do whatever needs to be done to get these features in and
working now so that I don't have to spend a day or two synchronizing
and debugging every time a big change is made to the code base!

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.

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.  These
can be done over a day or two as needed, certainly in time 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.

-Semiisotropic scaling: This would be a bit harder:  I'd have to sit
down and do the new derivations for the trotter decomposition.  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.

-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.  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 pretty easy to code in
without constraints.  With constraints, however, it requires solving a
self consistent equation, which can be done -- it just takes some
extra effort.  Might take some more mucking around.  I'm not exactly
certain how important this is in the near term, could definitely be
considered in the longer term.

-Nose-Hoover with chains with the leapfrog Nose-Hoover algorithm.
Leapfrog uses a linearized version, whereas my implementation follows
Martyna-Kline-Tuckerman and is not; so it would take a bit of work,
but wouldn't be too bad.

-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=50 fs tau_t=50 fs 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.

-Two major structural components were splitting update() in to several
different components, since I need to iterate over several different
components of the updating algorithm, and grouping all of the "gather
information from multiple states" into a single compute_globals call,
since this needed to be done multiple times per loop -- there's a
bunch of messy flag passing that could probably be simplified here.

-There's a number of smaller changes here and there -- these include
including a few more variables in structures to avoid passing things
back and forth as much, but these are the biggest ones.

-The number of Nose-Hoover chains is hard coded for simplicity.  This
should probably be switched eventually, just on principle.  Will
require changing code in a number of places, but should be
straightforward.

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.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821