[gmx-developers] single precision, double precision hybrid GROMACS

[Michael Shirts]

Hi, all-

A couple clarifications, questions, and comments:  Dan Ensign works
for Vijay Pande, so I believe that his question is actually following
after conversations that Erik, Vijay, and I had back at Stanford.  So
I think Dan is talking about what Erik is talking about here.

There are two "accuracy" issues here.  One is the accuracy of the
model.  Clearly, molecular mechanics models are not accurate.  The
question with energy conservation is -not- about accuracy of the
simulation -- it's is about the physical fidelity of the simulation. 
If energy is not conserved, the simulations drifts off the manifold of
phase space, and no longer samples from the correct distribution.  How
far off does it get? Well, really, that's an uncontrollable factor. 
This doesn't just hold in NVE, it holds in any other ensemble.  If you
can use a thermostat, then there's less of a problem -- but there are
some kinetic variables (diffusion constants, for example) that can't
be properly measured with a thermostat on.   And it would be good to
run NVE for 10's of ns without it heating up!

But even if you run with temperature control (preferably a temperature
control with an ensemble defined!) you still drift away from the
ensemble.  You can -say- that it doesn't matter, and in many cases it
probably doesn't, but physically, something incorrect for the model is
happening, and any averages are no longer strictly ensemble averages. 
  From previous talks with Dan, I believe that he is is planning on
running NVE simulations of large proteins, so the need the speed of
single precision with the energy conservation of double-precision is
important if it is obtainable -- which it seems to be.

Mark and Erik -- you seem to contradict each other.   Erik talks about
some mixed single-double experiments that were performed on key parts
(i.e., integration and constraints).  Mark seems to think they are in
the code, Erik doesn't.  Could you clarify, Erik?   Are the two of you
talking about the same code?  Assuming it's not in yet, would it be
possible to put such a version into the code as a separate compiler
option? Perhaps Dan can help Erik in implementing?

I suspect implementing Erik's test cases would go a long way towards
fixing this (at least for tapered cutoffs -- I suspect there's some
more tweaking that would need to be done for PME based on some of the
energy conservation tests I've been doing) -- I doubt that there would
be a need for anything other than integrator/constraints to be in
double.

A more general comment, about energy conservation:  Another important
reason for making sure energy converges is for bug-checking.  Lack of
energy conservation can show very subtle bugs that tend to hide
themselves otherwise.  A thermostat covers a multitude of sins!

I've talked to a number of people (Bill Swope and Jed Pitera at IBM,
Hans Andersen at Stanford), that are convinced that good energy
conservation is vital for replicable, statistically valid simulations,
and they've convinced me.   If you aren't convinced, please let me
know, and I'll dig up some more details.   At the very least, there
are a decent number of fairly smart people who feel uncomfortable
about GROMACS unless it can show good energy conservation.  I
understand that changing things will take some time, but I really do
think that ability to conserve energy is a design goal that really
needs to be met, and should be included in tests.  I'm talking about
very strict design goals, such as 1) how long can it run before
statistically significant energy drift is introduced and 2) does the
RMSD of the energy scale as t^2 as t->0?  Are tests like this done
right now as part of the validation suite?  Even if energy
conservation doesn't exist, it would be nice to document exactly how
far off the code is.

Clearly, I need to do my part -- A few months back, I committed to
putting velocity verlet in (which is required for correct
implementation of pressure control, since you need velocities and
positions at the same time to get the correct pressure, and therefore
conserve the conserved quantity for the ensemble), and haven't been
able to find time to merge codes yet.  I'll try to clear some more
time up in the next little while.

Thanks,
Michael Shirts
Research Fellow
Columbia University