<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=utf-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<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>
</body>
</html>