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

Carl Christian Kjelgaard Mikkelsen spock at cs.umu.se
Mon Apr 4 11:54:41 CEST 2016

```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>
```