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

[Erik Lindahl]

Hi,

>
> 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  
increase...

>
> 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  
others too.

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.

Cheers,

Erik