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

[Mark Abraham]

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.