[gmx-developers] re: Energy conservation

[Michael Shirts]

Hi, all-

> energy conservation seems to be a recurring topic. We are of course
> aware that different algorithms can degrade it. E.g. if you use
> constraints you will *never* get perfect energy conservation.

If the constraints are sufficiently tight, and implemented correctly,
then energy conservation is actually possible -- again, I need to
provide some more data here -- I agree that hearsay is not definitive
proof.  I was able to get very good conservation with TINKER -- I'll
try to get some data together.

> The leap-frog integration algorithm is fine though. We have done microsecond
> simulations of water clusters with almost perfect energy conservation using:
> - double precision
>- no cut-offs (vacuum simulations)
>- no thermostats
>- settle for constraints

This is great to nkow.  I think that 1) precision (double vs. single)
and 2) treatment of cutoffs are probably the biggest contributing
factors to lack of energy conservation.  The majority of simulations
are run under single precision / cutoff situations, so eventually, if
possible, energy conservation should be set as a goal for these
simulations as well.

> To the degree that energy conservation is not perfect (<1e-4) this is
> probably due to using a finite timestep.

Actually, finite timestep symplectic integrators conserve energy
exactly.  More precisely, a finite precision symplectic integrator
conserves the energy of a Hamiltonian that is a slight perterbation of
the original, infintely small timestep.  The energy error will be
constrainted to be a small distance away from the correct Hamiltonian,
for very long time scales.  So, lack of energy conservation is not
actually a result of finite precision timesteps.

Some refs: (http://www.unige.ch/~hairer/preprints/santander.pdf), and
I've put a few more preprints from Ernst Hairer here
http://www.columbia.edu/~mrs2144/verletX.pdf), where X is 1, 2 ,3.

> We have similar results for
> proteins with H-bonds constrained (using SHAKE) in water clusters, in
> vacuum, but much shorter simulations.

It might be good to have a validation data, test cases in the CVS.   I
think one of the problems is that we say "I've checked this" but if
the data and the mdp's used to generate the data is not easily
acessible, we end up talking past each other.  Might this be a feature
other people are interested in?  I think it would make things easier. 
  We can then work from the state of "This is exactly what it does
NOW" and "This is what a change could do" or "Here's what my data from
a similar system with another code are".  More productive than my
posts so far, certainly :)

> Michael, you mention "drifting off the energy surface". What does that
> mean? Are there two (or multiple) independent energy surfaces between
> which systems can not move if energy is conserved?

I'm referring to the manifold induced by Newtonian simulations using
the Hamiltonian.  It actually does evolve along a manifold in
H-dimensional space -- or with finite timsteps, a bounded distance
from this manifold.  For small systems (like the Henon-Heiles
problem), this can be easily visualized -- all-atom MD, not so much). 
If you conserve energy sufficiently, you remain bounded to near this
manifold.  The structure of such as space has lots of good
mathematical properties, like incompressibility of phase space, etc.

> Also, if people are convinced of the necessity of energy conservation, that's fine but not
> science. Is there any proof of these claims you mention? People tend to
> have lots of opinions on these things, but not many have the patience
> necessary for systematically testing and publishing them.

There is a signficant amount of research on this physics and applied
mathematics (some applied mathematics papers linked here).  But you're
right; it comes down to data, and I'll try to track more of that down
(in addition to the papers linked).

>I agree that not many published simulations are reproducible, but that
> usually is because people are sloppy in the methods section or in the
> description of results.

Complete agreement on this -- the vast majority of the lack of
reproducibility is because of sloppy methods and descriptions.  Errors
from improper intergration mostly show themselves in small biases, not
large discrepancies.  And the best value is still probably in
validation of code -- bugs just can't hide.

> By the way, how many people do actually try?

Do you want me to name names?  :) Try to get energy conservation? 
I've had discussions about this with Jay Ponder (TINKER), people in
Bruce Berne's (Sim package), Bill Swope and Jed Pitera (IBM's Blue
Matter code), people at D.E. Shaw Research (they've built a code
called DESMOND, you probably haven't heard of it yet, but it is going
to get some press relatively soon -- it will probably be speed
competive with GROMACS), and obviously in the Pande Lab.

> Sorry if this sounds like a rant,  but I would like to make clear that
> we do take these things seriously, and we're open to criticism that's
> based on facts, but not criticism that's based on opinions.

I agree that naming names is, at best, a shorthand for providing data,
and that good data is necessary.  I'd be happy to talk to these people
and either get them to post, or to get references for why energy
conservation is important, because, you're right, opinions alone are
only suggestive of  data, and not data themselves.  But I think
perceptions are important to address, as well.

I think that in the end, there won't be many changes in the results if
energy conservation is improved.  But I think there is real value in
making sure the underlying physics are rock-solid.  I'm not by a long
shot implying that GROMACS results are wrong, at all.  I'm just saying
that it would be nicer and more comforting to have some good
conserved-quantity properties underlying the simulation, to give one
less thing to worry about.  Having the energy conservation tests in
the CVS, including a .edr slot for the conserved quantity of the
ensemble, are both good things to put in eventually.    As long as we
agree that it's a good thing to work towards, then I'm fine putting it
lower on the priority list.

Best,
Michael