[gmx-developers] Question on SHAKE equations

["Pablo García Risueño"]

Dear Berk

I am a bit concerned about this point. Let us restrict ourselves to the
case of all H bond angles constrained. I was assuming that every H had one
single (arbitrary) constrained bond angle. This is, in the scheme below,
e.g. the angle C1-C2-H would be constrained, but not the H-C2-C3 angle.
Isn't this true?

            H
            |
        C1--C2--C3

In some literature I read that the vibration period of H bond angles is
the same as that of heavy atoms bond lengths (the Gromacs option h-angles
for constraints seems to assume this). Isn't this Gromacs option
proceeding as stated above?

Can you give me some reference where the procedure of 'virtual sites' to
constrain H bond angles is explained?

Thank you very much. Best regards


> On 6/10/13 12:15 , David van der Spoel wrote:
>> On 2013-06-10 10:57, "Pablo García Risueño" wrote:
>>> Thank you very much for the answer. It is said that:
>>>
>>> "A practical issue is that not many people use shake, since LINCS
>>> works in
>>> parallel too."
>>>
>>> But I thought that SHAKE was used when bond angles are constrained,
>>> this
>>> is when bond lengths of heavy atoms are constrained (since bond
>>> lengths of
>>> heavy atoms and bond angles of hydrogen have the same period). Isn't
>>> this
>>> correct?
>>>
>> That's up to the user, but LINCS can be used too for angles when the
>> number of iterations is increased a bit (although not too many people
>> use angle constraints either). Nevertheless gromacs supports the shake
>> algorithm and hence it should work.
> But note that constraining all-angles will lead to incompatible
> constraints, unless the molecule is very small.
> Constraining all angles involving hydrogens has similar issues. We use
> virtual sites to remove H-angle vibrations.
> So in practice the only molecule for which angle constraints are used is
> water, but there SETTLE is algorithm of choice.
>
> Cheers,
>
> Berk
>>> Thank you very much
>>>
>>>
>>>
>>>> On 2013-06-07 19:26, "Pablo García Risueño" wrote:
>>>>> Dear Gromacs developers
>>>>>
>>>>> I am studying the Gromacs implementation of the SHAKE algorithm. As
>>>>> far
>>>>> as
>>>>> I know, it corresponds to the cshake procedure, which lies in the
>>>>> shakef.c
>>>>> file, and the equation (5.6) of the SHAKE paper
>>>>> (see
>>>>> http://www.sciencedirect.com/science/article/pii/0021999177900985)
>>>>> is
>>>>> essentially in
>>>>>
>>>>> acor = omega*diff*m2[ll]/rrpr
>>>>> lagr[ll] += acor;
>>>>>
>>>>> where the relationship between variables is:
>>>>>
>>>>> SHAKE var               Gromacs var
>>>>>                             omega=1
>>>>> d^2 - (r')^2               diff
>>>>> 1/(2(1/m_i+1/mj))          m2[ll]
>>>>> \vec{r}_{ij}\cdot\vec{r}'_{ij}  rrpr
>>>>> g                         lagr[ll]
>>>>>
>>>>> being ll the constraint index. However, I think that the term
>>>>> proportional
>>>>> to g^2 in the equation (5.6) of SHAKE paper is not in the Gromacs
>>>>> code.
>>>>> This term is relatively small, but I was assuming that it is not
>>>>> negligible. Could some developer explain me whether this term is
>>>>> neglected, or it is included in the calculations elsewhere?
>>>>>
>>>>
>>>> I don't have the complete overview, so I will assume your
>>>> observation is
>>>> correct.
>>>> Since the algorithm is iterative the second term may not be needed,
>>>> although it might be converge in fewer iterations with it.
>>>>
>>>> A practical issue is that not many people use shake, since LINCS works
>>>> in parallel too.
>>>>
>>>>> Thank you very much. Best regards
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>>
>>>>> Dr. Pablo García Risueño
>>>>>
>>>>> Institut für Physik und IRIS Adlershof, Humboldt Universität zu
>>>>> Berlin,
>>>>> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>>>>>
>>>>> Tel. +49 030 209366369
>>>>>
>>>>
>>>>
>>>> --
>>>> David van der Spoel, Ph.D., Professor of Biology
>>>> Dept. of Cell & Molec. Biol., Uppsala University.
>>>> Box 596, 75124 Uppsala, Sweden. Phone:    +46184714205.
>>>> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
>>>> --
>>>> gmx-developers mailing list
>>>> gmx-developers at gromacs.org
>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>> Please don't post (un)subscribe requests to the list. Use the
>>>> www interface or send it to gmx-developers-request at gromacs.org.
>>>>
>>>
>>>
>>> --
>>>
>>> Dr. Pablo García Risueño
>>>
>>> Institut für Physik und IRIS Adlershof, Humboldt Universität zu Berlin,
>>> Zum Grossen Windkanal 6, 12489 Berlin, Germany
>>>
>>> Tel. +49 030 209366369
>>>
>>
>>
>
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
>


--

Dr. Pablo García Risueño

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

Tel. +49 030 209366369