[gmx-users] Water molecules cannot be settled: V-rescale is the cause

[Mark Abraham]

On 20/09/2012 7:10 AM, Ladasky wrote:
> After weeks of trying various conditions, I found my problem.  Here, from my
> position-restrained MDP file, are the two relevant lines:
>
> ; Temperature coupling is on
> tcoupl            = V-rescale ; Weak coupling for initial equilibration
> [snip]
> ; Pressure coupling is on
> pcoupl        = Berendsen     ; Pressure coupling on in NPT, also weak
> coupling
>
> This is a change that I forgot to mention in my list of changes, but in fact
> I did post the full MDP file.  I switched to using two DIFFERENT algorithms
> for temperature and pressure coupling when I performed the initial
> position-restrained equilibration.  This frequently sets up an INSTABILITY
> which manifests itself anywhere from a few picoseconds to a full nanosecond
> into the standard equilibrium MD run.
>
> I am giving this number in real time, rather than in compute cycles, because
> I tried cutting my step time to 1.0 fs from my usual 2.0 fs -- and I still
> had a water molecule dislocate abruptly, at nearly the same physical
> location and simulated time point.

This all sounds much like an issue with the topology or starting 
configuration. As you might know from the advice here 
http://www.gromacs.org/Documentation/How-tos/Steps_to_Perform_a_Simulation, 
having to use a smaller time step is a routine procedure in case of 
difficulties. I have seen systems where 0.5 fs without bond constraints 
was necessary to relax issues with the initial conformation, and using 
position restraints on (say) all heavy atoms would have made this 
impossible. Your protocol from July was also asking for trouble by 
generating velocities and moving straight into an NPT ensemble with 
position restraints. Starting with an NVT ensemble can be better idea, 
particularly if the volume is not quite right.

>
> I will share my full debugging process if anyone is interested.  It took me
> several runs, all of them over 24 hours long, to figure this out.
>
> Now, WHY would I have switched from a Berendsen temperature-coupling
> algorithm to the V-rescale algorithm?  Because of this cautionary note that
> I started receiving in GROMACS 4.5 when I started the position-restrained
> mdrun:
>
> NOTE 1 [file pr.mdp]:
>    The Berendsen thermostat does not generate the correct kinetic energy
>    distribution. You might want to consider using the V-rescale thermostat.
>
> I took this as a warning that a deeper study of the algorithms had revealed
> a flaw in Berendsen, and so I sought to use what GROMACS was advising me
> would be a better choice.  This is apparently MISLEADING advice when doing
> the initial equilibration.

Your observations on one system are not enough to reach this conclusion. 
v-rescale is normally quite appropriate for equilibration. The above 
hint is one that only matters when you wish to perform proper 
equilibrium sampling. Multiple phases of equilibration are normal, 
particularly in tricky cases, to gradually approach the conditions under 
which you wish to sample, via those that help deal with trouble with 
starting conditions.

>    If you are using the Berendsen pressure
> algorithm, you must also use the Berendsen temperature algorithm.

Have you controlled for the fact you were using position restraints 
before reaching this conclusion? Are the atoms you were restraining in 
fact in useful locations?

Merely switching to the Berendsen temperature algorithm can be lucky 
enough to lead to stable equilibration. Pathological starting conditions 
are sometimes quite sensitive. Numerical integration can be a tricky 
business, and it is difficult to draw general conclusions from limited data.

> I know that mdrun cannot be sophisticated enough to know which phase of my
> run I am performing, but I hope that the GROMACS documentation will be
> revised to discuss the appropriate use of the V-rescale thermostat.  I just
> misused it,and wasted a lot of time as a consequence.  I hope that some
> discussion of its proper use would be written and summarized on a page like,
> for example, this one:
>
> http://www.gromacs.org/Documentation/Terminology/Equilibration
>
> It already states: "Using, i.e., the Nosé-Hoover thermostat for initial NVT
> equilibration can lead to wild oscillations of the temperature, with the
> system ultimately blowing up."  A similar warning should be added about the
> V-rescale thermostat in NPT equilibration, ESPECIALLY if GROMACS is going to
> issue messages that recommend that you use it!
>
> There is a page for the Berendsen barostat/thermostat:
>
> http://www.gromacs.org/Documentation/Terminology/Berendsen
>
> There is apparently no corresponding page for the V-rescale thermostat.
> Perhaps one should be written.

Most of the documentation of GROMACS is a volunteer effort. I'm sorry 
you found it inadequate, but I do not think your conclusions are 
supported by your evidence.

Regards,

Mark