[gmx-developers] PR pressure control (ANick at t-online.de)

[Michael Shirts]

> sorry for bothering you again, but I have used PR for many, many simulations
> during the last year (the final results seemed not to be too bad) and I
> still hope to be able to use some of them.

> Does this on the other hand mean that the ANISOTROPIC pressure coupling
> is working fine?

In isotropic pressure control, the box acceleration is computed incorrectly,
and so the average pressure will be incorrect.  With anisotropic pressure
coupling, the box acceleration is correct, and the average pressure will be
correct.  What is incorrect in both cases is that the particle positions are
not scaled, so instead of scaling in the positions when the box volume
changes, so rthta the density is always equal throughout the box, the box
sides "peel away" from the atoms -- which isn't idea, since it's supposed to
be isotropic with the periodic boundary conditions.  So the average pressure
will be correct, but the ensemble will not be NPT.

How bad is this in practice?  With a large coupling constant (longer time
scale), it won't be that bad, since it will always be very close to
equilibrium. My results seem pretty reasonable -- if you aren't doing high
precision measurements, it's probably OK, but I can't guarantee it.  You could
also use Berendsen pressure control, which also gives the correct average
pressure but does not generate the correct ensemble, but is still very common
-- so it's no worse that what most people are doing.

If you use anisotropic temperature control, the other thing to watch out for
is the box side lengths changing -- since there is no restoring force, the box
side lengths can drift, causing eventual failure if one dimension drops below
the cutoff length.  Anisotropic temperature control really isn't designed for
fluids.

> Today, I have reanalysed many systems and what I actually see is a
> systematic, rapid decrease of the system density when changing from
> isotropic to anisotropic PR pressure coupling in all cases where I used a
> small pressure coupling constant (tau_p=0.1). My systems are simple united
> atom PE systems.

Yeah, I'm not suprised.  Larger constants should solve much of the problem
(maybe tau_p=1.0?  I can't be sure) -- it will get you much closer to the
correct behavior.

It is very likely that this will all be solved in Gromacs 4.0 -- but it will
probably be too complex to patch until then :(

Cheers,
Michael