[gmx-developers] SHAKE / LINCS crash

[Berk Hess]

On 2015-01-30 11:59, "Pablo García Risueño" wrote:
>> On 2015-01-29 08:43, David van der Spoel wrote:
>>> On 2015-01-28 18:13, "Pablo García Risueño" wrote:
>>>>> On Tue, Jan 27, 2015 at 11:26 AM, "Pablo García Risueño" <
>>>>> Risueno at physik.hu-berlin.de> wrote:
>>>>>
>>>>>> Dear Gromacs developers
>>>>>>
>>>>>> My name is Pablo G. Risueño, I am a postdoc researcher at the
>>>>>> Humboldt
>>>>>> University of Berlin. I am currently working with some colleagues
>>>>>> in the
>>>>>> development of new constrained MD schemes.
>>>>>>
>>>>>> A few years ago I had a conversation with prof. B. Hess, where he
>>>>>> told
>>>>>> me
>>>>>> that both SHAKE and LINCS can make a MD simulation to crash
>>>>>> (specially
>>>>>> for
>>>>>> large systems with many time steps). We would like to analyse the
>>>>>> problem,
>>>>>> so we would like to have some example(s) of simulation where such
>>>>>> crash
>>>>>> happens. We would need the Gromacs input file, and the
>>>>>> coordinates/velocities files at the step before the crashing. May
>>>>>> anybody
>>>>>> help us?
>>>>>>
>>>>> Hi,
>>>>>
>>>>> This is one of the easier problems in molecular dynamics. :-) You can
>>>>> probably take anybody's MD tutorial, multiply the time step by 10 and
>>>>> observe the simulation blow up from numerical instability. The
>>>>> constraint
>>>>> algorithm is often the first thing to fail, but it is not itself the
>>>>> problem. It seems to me that 99+% of the time the user has not
>>>>> prepared
>>>>> the
>>>>> system or model physics well enough, and most of the rest of the
>>>>> cases are
>>>>> code bugs.
>>>>>
>>>>> It is not clear to me that even if a constraint algorithm could be
>>>>> designed
>>>>> to satisfy such cases (of poor user input), that having it would be
>>>>> a good
>>>>> thing. It is irritating to have a simulation blow up because you
>>>>> were too
>>>>> rough with it, but it could well be worse to have your simulation
>>>>> succeed
>>>>> with only subtle evidence of the original problem.
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Mark
>>>> Dear Mark
>>>>
>>>> Thank you very much for your reply. However, if I understand it
>>>> correctly,
>>>> the crash of SHAKE I am speaking about is not a consequence of poor
>>>> definition of the input, but of nearly-singularity in the constraint
>>>> system of equations to solve. I guess that the (statistical) fails in
>>>> LINCS may be due to the fact that the correction that LINCS does may
>>>> not
>>>> suffice to satisfy the constraints for some given configurations.
>>>>
>>> Both shake and lincs only correct bond lengths along the direction of
>>> the original bond. However, if due to an unphysical situation, e.g.
>>> large forces or too long time steps, the bond rotates too much, let us
>>> say 90 degrees, the problem becomes unsolvable. Mathematically that
>>> could correspond to a singularity in the equations, but the equations
>>> become unstable due to other problems as Mark indicated.
>> There is one, unimportant, regime where LINCS (not SHAKE) will fail.
>> LINCS has the restriction that the eigenvalues of the constraint
>> coupling matrix should be smaller than 1, otherwise the expansion
>> doesn't converge (see the P-LINCS paper). For normal simulations where
>> only bonds are constrained all eigenvalues are smaller than 1. But when
>> you would additionally constrain all angles in a medium or large
>> molecules, eigenvalues larger than 1 will appear. However, constraining
>> all angles in larger molecules is not very useful.
>>
>> Cheers,
>>
>> Berk
> Dear prof. Hess
>
> Thank you very much for your reply. We had a conversation in Sept 2011 in
> Zaragoza, where you told me about some errors both in SHAKE and LINCS that
> happened specially for large systems after many simulation time steps.
> After the explanations of prof. van der Spoel and prof. Feenstra, I find
> myself a bit confused, so if you don't mind, I would like to ask you for
> further clarification.
>
> I understand the problem of "nearly orthogonality of bonds" pointed by
> prof. van der Spoel. However, prof. Feenstra told me that when his
> simulation crashes, he solves it by restarting the simulation in the
> latest non-crashed step and increasing the number of figures of the atomic
> positions. If this can solve the problem, I guess that it cannot be due to
> orthogonality of bonds, but rather to constraint matrix singularity (I
> expect the constraint matrix to be somewhat chaotic, so the increase of
> the number of figures could solve the nearly-singularity). However, I am
> not sure on what is actually the origin of the SHAKE and LINCS errors that
> you told me about in the past. May you explain me a bit more on it?
>
This is exactly the case David is referring to. It can happen, 
especially with large time steps and/or with not full equilibrated 
systems, that a bonds rotates close to 90 degrees. Then it might be 
problematic to get the bonds to the correct length, since the component 
orthogonal to the original length is nearly the whole bond length. This 
is not an issue of the matrix. Resetting to some steps before avoids the 
issue in most cases, since MD is chaotic and you will thus not end up in 
exactly the same state (unless you have binary reproducibility).

Berk