[gmx-developers] Integrators

[Michael Shirts]

I wanted to e-mail you all a little bit of the analysis of integrators
I've been doing in the framework of Gromacs.  A very useful way of
writing reversible explicit integrators is the Trotter factorization,
where the integrator is broken up into a number of Liouville states
using the relation e^i(A+B)dt \approx (e^i[B/Pdt]e^i[Adt/P])^p (order
dt^P+1,) with P = 2.

In the Trotter decomposition, the basic outline is to break up the
pressure and temperature integrator updates into Liouville operator
half-steps.  In the Martyna-Tuckerman-Klein scheme, pressure and
temperature control steps are done outside of the position / velocity
update.

Velocity verlet can be schematically set up as:

1. Half Trotter P/T step.
2. v(t+dt/2) = v(t) + dt/2 *F(t)
3. r(t+dt) = r(t) + dt*v(t+dt/2)
4. v(t+dt) = v(t+dt/2) + dt/2*F(t+dt)
5. Half Trotter P/T step.

Note that the Trotter steps for two consecutive overall steps can be
performed sequentially, as they use the same temperature and pressure.
In GROMACS, with the current bookkeeping, this translates as:

1. Calculate force F(t)
2.  Do first update
    a. v(t-dt/2) -> v(t)
    b. constrain the velocities at the full time step
3. Write trajectories (with r(t), v(t), F(t))
4. Compute K(t) and P(t) from v(t) and virial V(t).
5. Do first 1/2 Trotter coupling step.
6. Do second 1/2 Trotter coupling step.
7. Do second update
   a. v(t) -> v(t+dt/2)
   b. r(t) -> r(t+dt)
   c. constrain positions and half step velocities.

One important note that it appears that if the trajectories are not
output, and the kinetic energy does not need to be computed then,
there does not need to be a second constraint step, as the velocities
become essentially a helper variable for Stoermer integration, and
only the position needs to be constrainted.  This becomes relevant
since using the Trotter decomposition framework, it's possible to
update the pressure and temperature coupling only every N time steps
(as described in the 1996 Martyna-Tuckerman-Tobias-Klein paper --
although some care needs to be taken in that case as simulations can
become unstable if they are done too infrequently.

It is also important to mote that it may be more accurate to compute
the kinetic energy and pressure via the averaged 1/2 step kinetic
energies (square, then average) than from the averaged 1/2 step
velocities (average, then square).  This is true whether the
simulations are done with velocity verlet or leapfrog.  If we want to
compute the kinetic energy using the averaged two half step
velocities, instead of the full step velocity, we would need to
iterate 4,5,6, since the temperature control would need to be altered,
requiring extra communication. This is obviously rather
undesireable. Using the K(t) temperature is trivial, and requires no
iteration, but may lead to temperatures and pressures that are too
low.

Alternately, we could approximate the half step kinetic energies by
computing v(t+dt/2) in the absence of the Trotter scaling step. It's
not immediately clear if this would be a better or worse
approximation.  Assuming relatively weak coupling, it seems like it
would be pretty good.  This is just used to adjust the effective
temperature, so it would not affect the energy conservation, just the
effective temperature of the simulation.  For example, with a box of
water, 2fs time steps can apparently lead to a on-time step K.E. off
by 8 K.  If we approximate the 0.5[K(-dt/2) + K(+dt/2)] energy, then
if scaling is 0.1% per time step, the KE's will be off by about 0.2%,
leading to a temperature difference of about 0.6 K instead of 8 K,
which would be much better.  If we make a decent guess of what the
scaling was going to be (from the Nose-Hoover variables of the last step),
we could do even better.

For leapfrog, things get more complicated. The leapfrog algorithm, in
Trotter decomposition form, is:

1. Half Trotter P/T step (requires virial and pressure at t).
2. r(t+dt/2) = r(t) + dt/2*v(t)
3. v(t+dt) = v(t) + dt*F(r(t+dt/2))
4. r(t+dt) = r(t+dt/2) + dt/2*v(t+dt)
5. Half Trotter P/T step (requires virial and pressure at t+dt).

Again, the Trotter steps can be put back to back without a problem.
But in this case, we can see that the force and virial are calculated
at different time steps, which doesn't really work for calculating the
virial efficiently.    Generally, this is done in leapfrog with
a dt/2 offset, leading to:

1. Calculate force F(t)
2. Write trajectories (with r(t), v(t-dt/2), F(t))
3. Do update
   b. v(t+dt/2) = v(t-dt/2) + F(r(t))
   a. r(t+dt/2) = r(t) + dt/2*v(t+dt/2)
   c. constrain positions and half step velocities.
4. Compute K(t+dt/2) and P(t+dt/2). Note we have K(t+dt/2) exactly,
though it's not clear if K(t+dt/2) is a replacement for
0.5*[K(t+dt/2)+K(t-dt/2)]
5. Half Trotter P/T step (requires virial and pressure at t+dt/2).
6. Half Trotter P/T step (requires virial and pressure at t+dt/2).
7. Do update:
   a. r(t+dt) = r(t+dt) + dt/2*v(t+dt/2)
   b. constrain positions.

The virial must be computed at the half step, but the force
calculation must be done at the full step, which is very hard to get
around in terms of time savings.  One could somehow calculate an
approximate virial to calculate the pressure, instead of the true
pressure, but it's not clear how good a calculation that would be.

So it appears that velocity verlet has a advantage over leapfrog in
the pressure computation, as it less straightfrward to construct the
correct coupling states in the leapfrog case.

So, in the near term, I'm going to set up only velocity verlet using
the trotter decomposition, keeping the leapfrog functionality the
same.  I'll try to keep things open enough in the Trotter
decomposition framework that it is relatively easy to add in only
computing pressure / kinetic energy every N time steps, or adding
multiple short range inner loops (for example, updated bonds every 0.5
fs, nonbonds every 2 fs, etc.)

Please let me know any thoughts!

Cheers,
Michael