[gmx-developers] Question on SHAKE equations

[Berk Hess]

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