[gmx-developers] Question on SHAKE equations

[Shirts, Michael (mrs5pt)]

Hi, all-

I asked this before, and I believe (though Berk would have to chime in, or
I'd have to take the time looking through which I don't have) that there's a
numerical shortcut in the calculation that is true in the limit of small
delta, which of course delta is.  If you actually open up the debugger and
look at the deviations from the constraints after, they are indeed at the
proper tolerance specified by the user.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821


> From: David van der Spoel <spoel at xray.bmc.uu.se>
> Reply-To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Date: Fri, 7 Jun 2013 22:28:35 +0200
> To: <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] Question on SHAKE equations
> 
> 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.