[gmx-developers] Shake question?

[Shirts, Michael (mrs5pt)]

Hi, all-

A couple questions about the shake implementation.

Right now, it appears that shake tolerance (set as shake_tol) is the
tolerance on the relative distance from the constraint squared -- from the
manual, one might guess that shake_tol is the maximum relative deviation of
a constrained bond.  I think expressing an the parameter in terms of
distances instead of distances**2 would make more sense.

Also, currently in the cshake code, the constraint is considered satisfied
if iconv = 0, where:

iconv   = fabs(diff)*tt[ll];

Diff is the difference between the square bond length and the square desired
bond length, and tt[ll] currently is 1/(2*tol) -- rather than the tol*tol
which I think would be a bit more clear where tol is equal to shake_tol.

If tt is too big, then its possible that there's a 1 in about 4 billion
chance that iconv will be zero even if the tolerance isn't satisfied.
Actually, this appears to be the case even with the current code - if
ir->shake_tol is less than 10^-13, then you'd start seeing the looping
around the integers in a long int.

I'd propose changing the definition of 1/(2*tol) to 1/(tol*tol), and to
avoid the problem of integer wrapping by capping shake_tol internally in
readir.c , so that  the hardware limit for "long int" isn't violated.

I can go ahead and make these changes, but I wanted to ask around and verify
there wasn't something I was missing!

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