<html>
<head>
<meta content="text/html; charset=windows-1252"
http-equiv="Content-Type">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<div class="moz-cite-prefix">Hi,<br>
<br>
You have identified the correct function.<br>
The exact Lagrange multiplier is described in a comment in the
code:<br>
scaled_lagrange_multiplier Scaled Lagrange multiplier for each
constraint (-2 * eta from p. 336<br>
* of the paper,
divided by the constraint distance)<br>
<br>
But SHAKE is almost never used in GROMACS, because LINCS is faster
and especially because SHAKE is ill suited for parallelization
over MPI. The only other constraint method I know of that would be
suitable is CG MSHAKE, for a comparison with SHAKE and LINCS see:<br>
<a class="moz-txt-link-freetext" href="http://link.springer.com/article/10.1140%2Fepjst%2Fe2011-01525-9">http://link.springer.com/article/10.1140%2Fepjst%2Fe2011-01525-9</a><br>
But there would not be any significant advantage over LINCS.<br>
<br>
Cheers,<br>
<br>
Berk<br>
<br>
On 2016-04-04 11:54, Carl Christian Kjelgaard Mikkelsen wrote:<br>
</div>
<blockquote cite="mid:570239DE.8090704@cs.umu.se" type="cite">
<meta http-equiv="Content-Type" content="text/html;
charset=windows-1252">
<font face="Courier New, Courier, monospace">Hello,<br>
<br>
I wish to formulate a contribution to GROMACS. I am writing code
for applying Newton's method to solving nonlinear constraint
equations in the context of molecular dynamics using specialized
sparse direct solvers for the relevant linear systems.<br>
<br>
<b>NOTATION:</b><br>
<br>
I would like to know *exactly* which nonlinear equation is being
solved by the GROMACS subroutine cshake(). This routine accepts
unconstrained atomic coordinates (rprime) and returns
constrained atomic coordinates (rnew) along with a so-called
scaled Lagrange multiplier.<br>
<br>
My analysis of your code cshake() indicates that the exact
equation being solved is the vector equation<br>
<br>
g(rprime - inv(M)*Dg(r)'*lambda) = 0<br>
<br>
where <br>
<br>
g is the constraint function, see below<br>
rprime is the unconstrained positions at the next time step<br>
r is the constrained positions at the current time step<br>
lambda is the vector of the scaled Lagrange multipliers, and<br>
Dg(r)' is the transpose of the Jacobian of the constraint
function g at the current time step.<br>
M is the diagonal mass matrix <br>
<br>
If m is the number of atoms and n is the number of constraints,
then the Jacobian Dg is an n by 3m sparse matrix of partial
derivatives. The diagonal mass matrix M has dimension 3m.<br>
<br>
If bond k is between atoms a and b, with square bond length
sigma2_k, then one of many possibilities is <br>
<br>
g_k(r) = 0.5*(sigma2_k - |r_a - r_b|^2)<br>
<br>
where <br>
<br>
r_a is the 3d coordinates of atom a<br>
r_b is the 3d coordinates of atom b<br>
<br>
The corresponding constraint equation is simply the scalar
equation g_k(r) = 0.<br>
<br>
If I apply a nonlinear Gauss-Seidel iteration to this system of
n equations using one step of Newton's method for the individual
scalar equations, then I (appear) to arrive at the code which is
implemented in cshake(). The norm being used is the Euclidian
norm. <br>
<br>
The new constrained coordinates are <br>
<br>
rnew:= rprime - Dg(r)'*lambda<br>
<br>
<b>REQUEST:</b><br>
<br>
Please confirm or deny if I have correctly identified both<br>
<br>
1) the exact constraint function used by cshake<br>
2) the exact equation being solved by cshake<br>
<br>
I have to stress the importance of even the smallest detail such
as the scaling factor of 0.5 or the order of the vectors r_a and
r_b. The new constrained positions are unaffected by our
choices, but the *value* of the scaled Lagrange multiplier is
not. I understand that GROMACS uses the scaled Lagrange
multiplier outside of cshake, so any alternative to cshake()
must necessarily return exactly the same value.<br>
<br>
Kind regards<br>
<br>
Carl Christian<br>
<br>
-----------------------------------------<br>
Carl Christian Kjelgaard Mikkelsen, Ph.D<br>
Department of Computing Science and HPC2N<br>
Umeaa University<br>
Sweden</font> <br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<br>
</blockquote>
<br>
</body>
</html>