[gmx-developers] single precision, double precision hybrid GROMACS
lindahl at sbc.su.se
Fri Apr 7 18:38:21 CEST 2006
> 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.
NVE for simple liquids and similar systems can certainly be useful,
but for large molecules you might be problems due to a slowly
decreasing potential energy. In NVE that will lead to a temperature
> 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?
Currently we use a few double precision _scalar_ variables in places
where we've found and isolated accuracy errors. Virial summation and
the correction force due to constraints comes to mind, but there are
Doing the entire update and constraints in double is something
completely different. It includes double precision arrays for both
coordinates and forces, and accounts for several percent of the
flops. I could still accept having a separate option for that, but
before doing so I want to know where the error comes from, and that
it is indeed necessary to do that entire process in double -
otherwise we will eventually end up doing everything in double
instead of fixing the actual problem!
What needs to be done is this:
1. Create a simple test system (water can be used to test LINCS/SHAKE
too if the SETTLE section is removed from the itp file), and make
sure you conserve energy in double precision.
2. Run it in single too, and compare the accuracy of coordinates,
velocities, and forces in three places: before integration, after
integration but before constraints, and finally after constraints.
3. It isn't magic - when you find the first difference it must be
possible to track the propagation and determine exactly which parts
need to be done in double, or if it can be rewritten to work in single.
> 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 agree completely - that was actually the main reason I wanted to
make sure we can conserve perfectly in double (and single w/o
constraints) a couple of years ago.
> 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'm not saying you're wrong, but this is a classical difference of
opinion between "physicists" and "chemists", or rather depending on
the size of systems you simulate :-) Just as truly ergodic
thermostats or accurate ensemble barostats it seems to be much less
important for big systems.
However, energy conservation _per_se_ doesn't say much. Without
cutoffs we would get perfect conservation using cos(r) as a coulomb
potential, and even standard switch functions introduce an unphysical
bulge in the force.
The PME accuracy is actually related to this. We would conserve
energy better if you applied a (small) switch/shift function on top
of the Ewald-adjusted interaction, but doing so the total error in
the force compared to the exact result would increase compared to the
current sharp cutoff. There isn't a given answer - we might need to
allow both options in the future.
Anyway, it's far from an impossible problem, and it would be very
desirable to conserve energy perfectly even with constraints in
single. Doing the whole integration+constraints in double was just my
first attempt to isolate the problem, and then I haven't found time
to continue. I'd be happy to assist with comments and ideas, though.
More information about the gromacs.org_gmx-developers