[gmx-developers] SHAKE / LINCS crash

["Pablo García Risueño"]

> 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
>

Dear prof. Hess

Thank you very much for the reply, now it is clear for me.

Best regards.


> --
> Gromacs Developers mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers or
> send a mail to gmx-developers-request at gromacs.org.
>


Beste Grüße.


--

Dr. Pablo García Risueño

Institut für Physik und IRIS Adlershof, Humboldt Universität zu Berlin,
Zum Großen Windkanal 6, 12489 Berlin, Germany

Tel. +49 030 209366369