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
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
