[gmx-developers] Brainstorming ideas for minimally painless incorporation of an MC barostat (later MC moves in general)

[Shirts, Michael R. (mrs5pt)]

Hi, all-

> >       a. For a MC barostat, no communication between nodes is
>required, as
> >atoms are guaranteed stay on the same node -- everything only changes by
> >scaling

> This is not really true, since when constraints are present, these need
> to be applied after scaling, which leads to (minor) coordinate changes
> and requires communication. But we can probably do this without changing
> the decomposition.

Hmm.  Maybe this can be avoided.  If the centers of coupled constraint
system were scaled, rather than all atoms, then no iterative constraining
would be required. . . If these constraint data structures already contain
these sets of coupled constraints, then it would be easy to do this.  Note
that this is very close to molecular scaling.

If all atoms were scaled, then yes, any coordinate change would be minor
and could likely be handled by the same communication algorithm.



>I don't think any changes are needed for the forces.
>The constraining and calculating virial contributions are probably the
>most tricky parts.

One great thing about a MC barostat is that no virtual contribution is
needed for the accept/reject step.  One never calculates the pressure,
only the change in (E(T(x,p) + PV(T(x,p))), where T(x,p) is some
transformation of the coordinates and the momentum. where V(T(x,p)) is the
new volume after the transformation, and E(T(x,p)) is the total energy
after the transformation.  The pressure is the applied pressure, and is a
constant.   You also need to include a Jacobian factor calculated from
T(x,p) in the accept/reject criteria.

There are two main choices I've seen for T(x,p) -- one leaves the momentum
in place, and just scales the coordinates.  For isotropic, no constraints,
then the Jacobian factor (this may be off by a factor of V) is
(Vold_/Vnew)^N-1.  If the momentum are scaled in the opposite direction of
the coordinates, then the Jacobians cancel.  Parrinello-Rahman is
essentially a infinitesimal version of this second process.  For
constraints, then the exponent on V changes.

The virial would only need to be calculated in the normal place, to report
the average pressure, and would not be required to be calculated in the
accept-reject step.


The OpenMM developers have been playing around with choices for the MC
barostat, and I can get their opinions on which of these choices has
worked best.

There biggest reason I haven't dived into this yet is that I'm still
trying to get my mind around how the data structures for the force
calculation should be saved and destroyed when accepting and rejecting in
a efficient way. That's the biggest roadblock for me.  Once I understand
how to do that (or someone else makes any changes necessary for that),
then I can probably at least get a single node version up and running and
validated for correct statistical mechanics in a relatively short time.

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