On 20/09/2012 5:08 PM, Ladasky wrote:
Mark Abraham wrote
This all sounds much like an issue with the topology or starting
configuration.
And it is for that reason that, during my debugging process, I switched back
from the chimeric protein structures that I was building myself to a
standard PDB structure for which I had previously completed a 5-ns stable
run. I wanted to eliminate any problems that the starting configuration
might have caused.
Fair enough, and good strategy to simplify while learning.
Mark Abraham wrote
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.
I'll keep in mind that even 1.0 fs is sometimes not short enough, when
working far from equilibrium. That being said, the fact that I saw almost
the same water molecule misbehave at almost the same simulated time point,
whether I used a 2.0 fs or a 1.0 fs time step, gave me the strong hint that
my simulation was both reproducible -- and wrong. And thus, that my time
step was probably not the critical factor.
Every time you change the ensemble (i.e. .mdp settings), which includes
beginning, the system is generally not in equilbrium and that can lead
to forces that physically jar it. One prefers to choose the longest
stable time step for a long simulation that remains in equilibrium, in
order to be efficient with computing resources. It doesn't follow that
the preparation should use that time step. Particularly if the forces
are large when not in equilibrium, subsequent displacements are large
and the simulation can fail immediately, or set up weird resonance
effects (think Tacoma Narrows bridge) that break things later.
The same water molecule having the problem suggests a local effect,
rather than a system-wide effect. Thus the starting position and not the
protocol might be to blame, although a gentler protocol can often
overcome a poor starting position. In particular, an isolated water
molecule placed by genbox in the middle of a protein can show these
kinds of problems. For this kind of reason visualizing the results of
computational procedures is always essential.
Mark Abraham wrote
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.
If you read my reply to Justin, you will know that I'm not deliberately
trying to skip any steps. I have made minimal modifications to the MDP
files that were part of the GROMACS 3.3 lysozyme tutorial shell script.
That protocol included a position-restrained solvent equilibration step, and
I used it as-is.
OK, but the authors probably knew that they'd set things up to have the
right density for their system, so that NPT with PR would succeed
without problems. Doing multiple equilibration phases for pedagogical
reasons might have meant the users of the tutorial spent more time doing
that and not focusing on whatever the tutorial's objective really was.
If the point of the tutorial was not to discuss equilibration issues, it
might not make a sound choice for a template.
I'm not married to the need to go straight to NPT, or even to use position
restraints during solvent equilibration. If the starting configuration of
my protein of interest shifts while I'm making adjustments to the solvent, I
don't mind. I'm trying to watch it fold, and the end point is far more
important to me than the starting point.
Then getting rid of PR entirely is plausible.
Did Berendsen himself ever post to this mailing list? I remember coming
across this signature line on some GROMACS-related post: "Why do YOU use
constraints?" In my case, the answer is, "because that's what my tutorial
recommended."
It's one of the quotes GROMACS tools print out as they exist, and
presumably an in-joke from years past.
Mark
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.
Mark Abraham wrote
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.
I understand that hint, in hindsight. It led me down a bit of a primrose
path.
Mark Abraham wrote
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.
Merely switching to the Berendsen temperature algorithm can be lucky
enough to lead to stable equilibration. Pathological starting conditions
are sometimes quite sensitive.
Well, I've been consistently lucky then, up to the point that I switched to
V-rescale. But I would rather be right than lucky. I don't care if my
simulation takes an extra day to run, as long as the darn thing doesn't
crash.
Thanks again for your help!
--
View this message in context:
http://gromacs.5086.n6.nabble.com/Re-Water-molecules-cannot-be-settled-why-tp4999302p5001134.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
--
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
