[gmx-developers] Various integrator questions . . .

[Michael Shirts]

Hi, all-

> Yes, being trained by Berendsen, please allow me to be biased. First I
> think the algorithms should be there because of choice, and users should
> be encouraged to take an informed decision.

As a Platonic ideal -- as time goes on, code accumulates a large
number of different methods, complicating things immensely.  Methods
that are discovered to not be particularly useful should be removed in
order to keep the code base workable and flexible, and even to prevent
people who don't understand the consequences of various choices from
making bad black-box decisions.  If discovered afterwards that such a
method is really needed, then can be added back in.

In this -specific- case, however, I agree that the overhead of keeping
Berendsen is very low, and it can easily be included at this time. 
Erik Lindhal pointed out in a separate e-mail that it is very stable
and quick in relaxing to near-equilibrium ensembles.

> Second, the Nose-Hoover thermostat gives the correct ensemble for your
> system plus "something else". In fact it therefore does not give you the
> correct ensemble for what you are interested in.

I think Michel addressed this before -- for the physical variables,
you get the correct canonical
distribution.

> Third it is very useful to strive for proper algorithms, in particular
> for reasons of reproducibility of simulations. However, the accuracy of
> the results is dependent on many variables, and in particular the force
> field introduces systematic errors that cannot be compensated for by
> other means. Moreover, small errors in e.g. the temperature coupling
> algorithm are negligible compared to systematic bias due to force field,
> finite system size etc.

I'm really not certain that this is the case.  Michel presented some
examples.   Certainly, for the infinite N limit, all ensembles
converge to the right distribution (if, of course, the method results
in a well-defined ensemble!)  But it's better to eliminate as many
sources of error as possible.  We'll have more time to do more science
if we can eliminate various forms of error that need to be accounted
for.

> In summary I'm all for better algorithms, but I vote for keeping the old
> ones as well. This means that we'll have to look into ways to make the
> integration part more modular, such that we, as an old friend used to
> say, can separate the Math from the Physics.

Different temperature controls are usually pretty easy to add on in a
very modular way.  Different pressure controls can be more complicated
-- especially with leapfrog, because the velocity and pressure are not
defined at the same time.


Michel said:
> Like you say, in MD, there are so many parameters, and so few checks to
> constrain these parameters. If we have the chance to detect that
> something might be a source of error, we should correct it. Good
> integrators are available
> already in the literature (see references in a previous mail).

Yes.

> A little cautious pointer: The Nose-Poincare thermostat (Bond et al, J
> Comp Phys 151 : 114, 1999; Leimkuhler and Sweet, JCP 121(1) : 108, 2004)
> might become the method of choice.  It is hamiltonian like the Nose
> thermostat, but cures the problem of the uneven time sampling. Due to
> it's hamiltonian structure, it allows to build formally smplectic
> integrators, which have all the nice properties one might dream of  (eg.
> mathematcally ;-) bounded energy drift).

I'll have to check this out.  Thanks for pointing it out.  Symplectic
integrators are always very good, and should be used whenever
possible,


Berk said:

> I changed the kinetic energy calculation in the CVS from being based on
> the average of the n-1/2 and n+1/2 velocities to n-1/2 and n+1/2 Ekin.
> With the old method the discretization error in Ekin (always causing
> Ekin to be too small) was 4 times larger than with the new code.
> The error causes the temperature to appear as too low, so the temperature
> coupling scales it up too much, which effectively results in an ensemble
> at a too high temperature.

I'll put in a plug again for Velocity Verlet, which avoids this whole
issue, as the velocities are determined at the same time as the
coordinates.  I'll start working on implementing my local version on
the 3.3 codebase, so that people can check it out.

Best,
Michael