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

[Michael Shirts]

Hi,

> last month (on 3rd of April) you wrote some thoughts about Parrinello-Rahman
> pressure control to this list. Do the problems you described mean that I
> cannot expect to get the right density after a long MD with PR pressure
> coupling? ( I calculated the system density with the g_energy tool.)

The current PR pressure control is not implemented correctly in the current
version of gromacs.  Two problems are the isotropic volume acceleration is
computed incorrectly, and there is no scaling of positions.  I am working with
Berk Hess and Erik Lindahl to try to get this corrected iin 4.0.

> Do I end up with a different density when using i.e. Berendsen pressure
> coupling?

You should check :).  It doesn't take too long.  Right now, since 1) PR is not
implemented correctly and 2) Berendsen pressure control does not actually
yield isobaric ensembles, I would expect to get a slightly different density
for most systems.  How significant the difference is depends on the system and
how precisely you need the answer.

> I observed some deviation in the density of simple PE systems after
> equilibration with GROMACS - isotropic PR pressure coupling and all bonds
> constrained - in comparison to some collegues using other simulation
> packages (both of us used identical force fields and simulation parameter as
> far as possible).

There are a number of other possible reasons the results would not be the
same, as many "standard" methods are implemented in different ways in
different codes.  Usually, all use reasonable (i.e., physical)
implementations, but that doesn't mean they are compatible.  Also, there are
still a number of bugs in virtually all MD packages.

In comparing the results of different codes, I go through this hierarchy.
Sorry if this is a bit of a review -- I'm sure you probably already do
something like this.

1) Does a single configuration gives the same energy with both simulation
codes?  You get the same energy to a reasonable number of digits for single
conformations (in my opinion, at least to single precision, so agreement to
one part in 10^-7 to 10^-8).

2) Do you get the same results (for a variety of ensemble variables) for two
NVE simuations with the different codes?  Because of differences machine
rounding, it is virtually certain that two simulations with different codes
starting from the same simulation give different trajectories.  However, they
should yield ensemble averages that agree to within the uncertainty.

3) If the NVE simulations agree, then try NVT.

4) If NVT agrees, then try NPT.

That's how I'd approach it.  In my experience, you'll start seeing differences
between any two codes in either step 1) or step 2), let alone step 3) and 4).

Cheers,
Michael Shirts
Pande Group / Folding at Home
Stanford University