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

[Berk Hess]

I understand that the virial is not needed for the pressure scaling 
itself, but it will usually be calculated for reporting at a scaling 
step, so we need to get it right, also with constraints.

If we don't want to calculate the force 3x when a MC step is rejected, 
we can directly reuse the storing mechanism implemented in minimize.c.

Cheers,

Berk

On 07/17/2014 04:20 PM, Shirts, Michael R. (mrs5pt) wrote:
> 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
>
>
>