[gmx-developers] Re: New implementations of Parrinello-Rahman

[Michael Shirts]

> You said that the code that you implemented is not completely
> reversible. I think that that it would be worth to look in the
> bibliography for a truly reversible method.

I've implemented the original isotropic coupling scheme of Andersen
(J. Chem. Phys. 72 (1980), no. 4, 2384-2393). There is one
approximation that needs to be iterated in order to be completely
reversible, in the velocity of the box volume at the t+dt step, but
it's pretty minor, and you get excellent energy (or rather, enthalpy)
conservation even without it.

> Could you place your code in a CVS branch?

To do this correctly, it requires knowing the velocity and positions
at the same time.  This requires the use of velocity verlet, which
requires a complete reorganization of the iteration (i.e. update())
function.  I have this working on my personal code, which is from
3.1.4.  There's some nonsignificant debugging to do to merge them, and
since it's been decided it won't be in the 3.3 release, and I'm being
paid to work on another project, it's a bit lower on my priority list
right now.  I am pretty sure I'll be able to get it in in the next
couple of months.  It's working fine now locally, it just needs to be
integrated, which might take a bit of time.

> On the other hand, how do you handle constraints? Adjusting
> velocities so that they do not modify constraints is not difficult,
> and does not require iteration. But it leads to drifts. One must
> constraint positions, and doing so would require either dropping
> reversibility or iterating. Perhaps there is a method in the bibligraphy
> to apply constraints in a reversible way.

The constraints are applied to the positions after the rescaling.   
It's most efficient to just rescale the center of mass positions of
the molecules, and to use the molecular virial instead of the atomic
virial, but that's not really supported in GROMACS right now since
there is no separate list of molecules set up.  It might be possible
to use the molecular virial for the water molecules (or whatever the
solvent is), since they are already defined separately, and the atomic
for the other molecules.   I've heard that there's a very small issue
with this (heard through the grapevine from Martyna and Tuckerman - I
think it's for any set of constraints, not just constraints updated
this way, if anyone has the reference, let me know), but in my tests
of energy conservation, it's extremely small, and there's no drift in
the conserved quantity for many nanoseconds.

Best,
Michael