[gmx-users] Energy Drift

[Ben Reynwar]

I posted to the list a few days ago with an energy drift problem.
Mark Abraham helpfully suggested using all-bonds rather than h-bonds
which solved the problem.  I'm now trying to understand quite why that
helped so much.

The simulation is a protein of about 5000 atoms using GBSA, a time
step of 2 fs, and a cut-off of 1.6 for VdW, coulomb and GB.

I've run energy drift simulations using the md integrator, with no
thermostat, and using the ld integrator with a tau_t of 1000000.
Simulations were run changing various parameters.  The changes were:
 - dt reduced to 1.5
 - dt reduced to 1.0
 - cutoff increased to 2.0
 - constraint changed from h-bonds to all-bonds

Plot of the energy drifts can be seen at
http://www.reynwar.net/ben/gromacs/energy_drift.png.

Noticeable features are:
 - changing the time step makes no difference (in my last post I
claimed it did, which is why you should make plots rather than
eyeballing log files).
 - increasing the cut-off helps a lot.
 - changing constraint to all-bonds make a dramatic difference
 - using ld there is a downwards drift in the energy when using
all-bonds constraint.  The temperature is roughly 300 K and the set
point is 400 K so this downwards drift seems unlikely to be due to
coupling to the langevin thermostat.

My questions are:
 - why does the all-bonds constraint help so much?
 - why doesn't moving to a smaller time step help with this?
 - what is the cause of the downwards drift when using ld with all-bonds?

Cheers,
Ben