[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>
More information about the gromacs.org_gmx-developers
mailing list