[gmx-developers] Thoughts on pressure control / constraints / leapfrog

[Michael Shirts]

Hi, all-

I've been working to fix Parinello-Rahman pressure control and
Andersen temperature control, and have encountered some issues in the
GROMACS code.

The first subset consists in pressure control algorithms. The currrent
Parrinello-Rahman algorithm in GROMACS has some problems.  First, in
the isotropic case of Parrinello-Rahman pressure coupling the second
derivative of the box volume is incorrectly calculated in line of
update.c, as the cross-terms are absent (coupling.c:173-184, gmx
3.2). This is fairly easy to fix, and appears to solve the problem
I've documented on the mailing list of the actual pressure being
different from the control temperature.

There is also no scaling of the positions as the box size changes (for
example, line 123 of update.c), which is what is in the original
definition of P-R calls for. I've been working on this, but have yet
to get it working correctly to ensure enthalpy conservation (NPH
ensemble).  As it is currently in gmx 3.2, the box is expanding, but
the molecules are not properly spreading apart, so there is a gap
around the edge of the box.

The second set is problems that I've encountered because of the
leapfrog integrator.  Without constraints, velocity verlet and
leapfrog verlet are completely equivalent.  However, WITH constraints,
the velocities in the two algorithms are no longer the same.  Also,
the equations appear to not be equivalent with extendend ensemble
methods of pressure control.

The difference occurs is that because with velocity verlet, to
correctly implement constraints, at time t, it is possible (necessary,
actually) to correct both the velocities and the positions such that
they obey constraints (i.e., SHAKE and RATTLE). In GROMACS currently,
this constrained velocity at t+0.5*dt is determined by taking the
difference between the constrained positions r(t+dt) - r(t) / dt.
This is not equivalent, and it doesn't appear that v(t+0.5*dt) is
actually obeying any constraints with respect to the bonds, but I
haven't finished with the math yet.  Does anybody have any references
that leapfrog and verlet are equivalent with constraints?

Additionally, pressure control is difficult becase the new coordinates
at time t+dt depend on the velocities at time t, as well as the
velocities from t+0.5*dt.  You can't just solve the equations in the
same way as you do to produce leapfrog, because the equations mix the
p and q and their derivitives in a way that newton's original
equations don't.  Velocity verlet solves this, as you have everything
at the same time.

I realize that one advantage of leapfrog is that lincs can be applied, which
is more easily paralellizable for large molecules that are split between
processors.  But I would argue that it may be useful to enable rattle as well,
and move to velocity verlet.  Perhaps lincs could be generalized to include
velocity constraints as well?  Also, rattle can be done iteratively per
molecule, which means that most of the waters don't need to be communicte to
run it.  Can't lincs be used with velocity verlet instead of leapfrog, with
velocities from position difference?  It should be no worse than the current
situation.

Thanks,
Michael