Constraint algorithms#
Constraints can be imposed in GROMACS using LINCS (default) or the traditional SHAKE method.
SHAKE#
The SHAKE 46 algorithm changes a
set of unconstrained coordinates
This action is consistent with solving a set of Lagrange multipliers in the constrained equations of motion. SHAKE needs a relative tolerance; it will continue until all constraints are satisfied within that relative tolerance. An error message is given if SHAKE cannot reset the coordinates because the deviation is too large, or if a given number of iterations is surpassed.
Assume the equations of motion must fulfill
For example,
where
The displacement due to the constraint forces in the leap-frog or
Verlet algorithm is equal to
SETTLE#
For the special case of rigid water molecules, that often make up more than 80% of the simulation system we have implemented the SETTLE algorithm 47 (sec. Constraint algorithms). The implementation of SETTLE in GROMACS is a slight modification of the original algorithm, in that it completely avoids the calculation of the center of mass of the water molecule. Apart from saving a few operations, the main gain of this is a reduction in rounding errors. For large coordinates, the floating pointing precision of constrained distances is reduced, which leads to an energy drift which usually depends quadratically on the coordinate. For SETTLE this dependence is now linear, which enables accurate integration of systems in single precision up to 1000 nm in size. But note that the drift due to SHAKE and LINCS still has a quadratic dependence, which limits the size of systems with normal constraints in single precision to 100 to 200 nm.
For velocity Verlet, an additional round of constraining must be done, to constrain the velocities of the second velocity half step, removing any component of the velocity parallel to the bond vector. This step is called RATTLE, and is covered in more detail in the original Andersen paper 48.
LINCS#
The LINCS algorithm#
LINCS is an algorithm that resets bonds to their correct lengths after an unconstrained update 49. The method is non-iterative, as it always uses two steps. Although LINCS is based on matrices, no matrix-matrix multiplications are needed. The method is more stable and faster than SHAKE, but it can only be used with bond constraints and isolated angle constraints, such as the proton angle in OH. Because of its stability, LINCS is especially useful for Brownian dynamics. LINCS has two parameters, which are explained in the subsection parameters. The parallel version of LINCS, P-LINCS, is described in subsection Constraints in parallel.
The LINCS formulas#
We consider a system of
where
In a numerical integration scheme, LINCS is applied after an unconstrained update, just like SHAKE. The algorithm works in two steps (see figure Fig. 8). In the first step, the projections of the new bonds on the old bonds are set to zero. In the second step, a correction is applied for the lengthening of the bonds due to rotation. The numerics for the first step and the second step are very similar. A complete derivation of the algorithm can be found in 49. Only a short description of the first step is given here.
A new notation is introduced for the gradient matrix of the constraint equations which appears on the right hand side of this equation:
Notice that
where
The derivation of this equation from (104) and (105) can be found in 49.
This first step does not set the real bond lengths to the prescribed
lengths, but the projection of the new bonds onto the old directions of
the bonds. To correct for the rotation of bond
where
This correction for rotational effects is actually an iterative process, but during MD only one iteration is applied. The relative constraint deviation after this procedure will be less than 0.0001 for every constraint. In energy minimization, this might not be accurate enough, so the number of iterations is equal to the order of the expansion (see below).
Half of the CPU time goes to inverting the constraint coupling matrix
The matrix
The matrix
This inversion method is only valid if the absolute values of all the
eigenvalues of
For molecules with all bonds constrained the eigenvalues of
where
The LINCS Parameters#
The accuracy of LINCS depends on the number of matrices used in the expansion (112). For MD calculations a fourth order expansion is enough. For Brownian dynamics with large time steps an eighth order expansion may be necessary. The order is a parameter in the mdp file. The implementation of LINCS is done in such a way that the algorithm will never crash. Even when it is impossible to to reset the constraints LINCS will generate a conformation which fulfills the constraints as well as possible. However, LINCS will generate a warning when in one step a bond rotates over more than a predefined angle. This angle is set by the user in the mdp file.