[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.
Carl Christian Kjelgaard Mikkelsen, Ph.D
Department of Computing Science and HPC2N
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the gromacs.org_gmx-developers