[gmx-developers] Question on SHAKE equations

["Pablo García Risueño"]

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?

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