[gmx-developers] Re: Energy conservation

[Michael Shirts]

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

Because the original system wasn't equilibrated?  If you start in an
equilibrated system, and have no/little energy drift, then there
shouldn't be an potential energy decrease. . .

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

Ah, excellent.  Thanks for the clarification.

> 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!

Very reasonable.

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

This is very true.  There are very likely fewer problems the larger
the system gets.   But it is, at the core, physics, and problems that
pop up in small systems have a tendency to appear in pathological
cases in large systems.   So, if one can avoid them without too much
cost, one less thing to worry about.

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

Well, yes, there is also the question of accuracy of the force field. 
But I see that as a qualitively different, even orthogonal question,
and definitely the biggest one.

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

Yes, I think this gets to my previous e-mail on the list -- I though
there already was one, but if there isn't then rlist>rcoul doesn't
make much sense.  This is something I'd be happy to look at -- though,
it will be further on my priority list as well ;)

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

I'm very happy with this.  As long as we recognize it's a good design
goal, then I think it will eventually happen, and that's good enough
for me.  I'll be happy to help out  -- though it is also a little
lower on my priority list for now as well! :)

Best,
Michael