[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


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.


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


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)


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


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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20160404/e11e7f8d/attachment.html>

More information about the gromacs.org_gmx-developers mailing list