[gmx-developers] Which nonlinear constraint equation is being solved in GROMACS subroutine cshake()?

[Carl Christian Kjelgaard Mikkelsen]

Hello,

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.

*NOTATION:*

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.

My analysis of your code cshake() indicates that the exact equation 
being solved is the vector equation

g(rprime - inv(M)*Dg(r)'*lambda) = 0

where

g      is the constraint function, see below
rprime is the unconstrained positions at the next time step
r      is the constrained positions at the current time step
lambda is the vector of the scaled Lagrange multipliers, and
Dg(r)' is the transpose of the Jacobian of the constraint function g at 
the current time step.
M      is the diagonal mass matrix

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.

If bond k is between atoms a and b, with square bond length sigma2_k, 
then one of many possibilities is

g_k(r) = 0.5*(sigma2_k - |r_a - r_b|^2)

where

r_a    is the 3d coordinates of atom a
r_b    is the 3d coordinates of atom b

The corresponding constraint equation is simply the scalar equation 
g_k(r) = 0.

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.

The new constrained coordinates are

rnew:= rprime - Dg(r)'*lambda

*REQUEST:*

Please confirm or deny if I have correctly identified both

1) the exact constraint function used by cshake
2) the exact equation being solved by cshake

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.

Kind regards

Carl Christian

-----------------------------------------
Carl Christian Kjelgaard Mikkelsen, Ph.D
Department of Computing Science and HPC2N
Umeaa University
Sweden
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20160404/e11e7f8d/attachment.html>