[gmx-developers] Re: Integrators

[Michael Shirts]

Well, unfortunately, it looks like the case with constraints is more
complicated than I had initially thought.  The issue is that there's
coupling between the Trotter steps and the constraint steps, which
means there needs to be iteration over the constraints.  Fortunately,
this is discussed in Appendix D of the 1996
Martyna-Tuckerman-Tobias-Klein paper.

The procedure presented there goes like this (using a = F/m, and h = \delta t).

1)  Trotter half-step [using P(t-1) and K(t-1)], using
    the Lagrangian multipliers from RATTLE of the previous step.
2)  r(t+h)  = r(t) + h*v(t) + h^2/2*a(t)
3)  v(t+h/2) = v(t) + h/2*a(t)
4)  Constrain r(t+h), and recalculate the Lagrangian multipliers
5)  Iterate 1-4
6)  Calculate a(t+h)
7)  v(t+h)  = v(t+h/2) + h/2*a(t)
8)  Second Trotter half-step, using the P(t) with the forces from the Lagrangian
    multiplier obtained after step 5.
9)  Constrain v(t+h) using RATTLE, to determine the full step velocities.
10) Repeat steps 7-9 until the constraints are satisfied, and the box
    speed is consistent.

Note that the even though the constraints can be satisfied with one
step, the Trotter step is no longer correct in this case, because it
requires the velocities -- which were not correctly constrainted when
we started.  So we must iterate.

Translating this into the current GROMACS bookkeeping, we get:

1)  Calculate a(t)
2)  v(t)   = v(t-h/2) + h/2*a(t)
3)  A Trotter half-step, using the P with the forces from the Lagrangian
    multiplier obtained after step 5.
4)  Constrain v(t+h) using RATTLE, to determine v(t).
5)  Repeat steps 2-4 until the constraints are satisfied, and the box
speed is consistent.
6)  Output r(t),v(t),a(t)
7)  Another Trotter half-step, using the Lagrangian multipliers from
RATTLE of the previous step.
8)  r(t+h)  = r(t) + h*v(t) + h^2/2*a(t)
9)  v(t+h/2) = v(t) + h/2*a(t)
10) Constrain r(t+h), and recalculate the Lagrangian multipliers.
11) Iterate 7-10.
12) Output P(t), K(t), and U(t).

This is obviously a bit more complicated.  We must do several
iterations over the constraints, and several summations of the kinetic
energy each time we do the constraints, leading to a lot of communication.

I also still need to translate this into the constraint formalism used
in Gromacs currently, where the Lagrangian multipliers are not
specifically stored.

One thing that might not totally kill the scaling is that it -looks-
like we still don't have to pressure / temperature couple after every
step.  Perfectly legitimate dynamics can be constructed by only
pressure coupling every N steps, though the bounds to the
approximation needs to be tested pretty strictly -- it could be that
it starts having issues if it's more than 2-3 steps at a time between
pressure "bumps".  Kinetic energy and multiple constraint iterations
don't need to be done in between.

Any thoughts as to what the right order of putting things might be?  I
think it's worth getting this in, even if it's not going to be the
default for large scale simulations.   But it will take a bit of time.
 I can probably go ahead and get everything except the iteration over
the pressure coupling pretty easily -- it's pretty much in, I just
need to tie together the loose ends I left by trying to get
constraints working and debug a bit.

Anything beyond that will take a bit longer. . .

Cheers,
Michael