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

[Berk Hess]

Hi,

You have identified the correct function.
The exact Lagrange multiplier is described in a comment in the code:
scaled_lagrange_multiplier    Scaled Lagrange multiplier for each 
constraint (-2 * eta from p. 336
  *                                             of the paper, divided by 
the constraint distance)

But SHAKE is almost never used in GROMACS, because LINCS is faster and 
especially because SHAKE is ill suited for parallelization over MPI. The 
only other constraint method I know of that would be suitable is CG 
MSHAKE, for a comparison with SHAKE and LINCS see:
http://link.springer.com/article/10.1140%2Fepjst%2Fe2011-01525-9
But there would not be any significant advantage over LINCS.

Cheers,

Berk

On 2016-04-04 11:54, Carl Christian Kjelgaard Mikkelsen wrote:
> 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/9bd2e447/attachment.html>