Constraint (computational chemistry)
In computational chemistry, a constraint algorithm is a method for satisfying the Newtonian motion of a rigid body which consists of mass points.
The general steps involved are: (i) choose novel unconstrained coordinates (internal coordinates), (ii) introduce explicit constraint forces, (iii) minimize constraint forces implicitly by the technique of Lagrange multipliers or projection methods.
However, explicit constraint forces give rise to inefficiency; more computational power is required to get a trajectory of a given length.
Constraint algorithms achieve computational efficiency by neglecting motion along some degrees of freedom.
For instance, in atomistic molecular dynamics, typically the length of covalent bonds to hydrogen are constrained; however, constraint algorithms should not be used if vibrations along these degrees of freedom are important for the phenomenon being studied.
The motion of a set of N particles can be described by a set of second-order ordinary differential equations, Newton's second law, which can be written in matrix form where M is a mass matrix and q is the vector of generalized coordinates that describe the particles' positions.
If M constraints are present, the coordinates must also satisfy M time-independent algebraic equations where the index j runs from 1 to M. For brevity, these functions gi are grouped into an M-dimensional vector g below.
This problem was studied in detail by Joseph Louis Lagrange, who laid out most of the methods for solving it.
The drawback of this approach is that the equations may become unwieldy and complex; for example, the mass matrix M may become non-diagonal and depend on the generalized coordinates.
The two difficulties of this approach are that the constraints are not satisfied exactly, and the strong forces may require very short time-steps, making simulations inefficient computationally.
For example, the dihedral angles of a protein are an independent set of coordinates that specify the positions of all the atoms without requiring any constraints.
The difficulty of such internal-coordinate approaches is twofold: the Newtonian equations of motion become much more complex and the internal coordinates may be difficult to define for cyclic systems of constraints, e.g., in ring puckering or when a protein has a disulfide bond.
The original methods for efficient recursive energy minimization in internal coordinates were developed by Gō and coworkers.
in the next timestep, the Lagrange multipliers should be determined as the following equation, This implies solving a system of
Although there are a number of algorithms to compute the Lagrange multipliers, these difference is rely only on the methods to solve the system of equations.
The SHAKE algorithm was first developed for satisfying a bond geometry constraint during molecular dynamics simulations.
[10] The method was then generalised to handle any holonomic constraint, such as those required to maintain constant bond angles, or molecular rigidity.
The original SHAKE algorithm is capable of constraining both rigid and flexible molecules (eg.
water, benzene and biphenyl) and introduces negligible error or energy drift into a molecular dynamics simulation.
[13] One issue with SHAKE is that the number of iterations required to reach a certain level of convergence does rise as molecular geometry becomes more complex.
[13] Hence the CPU requirements of the SHAKE algorithm can become significant, particularly if a molecular model has a high degree of rigidity.
A later extension of the method, QSHAKE (Quaternion SHAKE) was developed as a faster alternative for molecules composed of rigid units, but it is not as general purpose.
It is worth mentioning that MSHAKE computes corrections on the constraint forces, achieving better convergence.
P-SHAKE computes and updates a pre-conditioner which is applied to the constraint gradients before the SHAKE iteration, causing the Jacobian
The M-SHAKE algorithm[21] solves the non-linear system of equations using Newton's method directly.
This improvement comes at a cost though, since the Jacobian is no longer updated, convergence is only linear, albeit at a much faster rate than for the SHAKE algorithm.
Several variants of this approach based on sparse matrix techniques were studied by Barth et al..[22] The SHAPE algorithm[23] is a multicenter analog of SHAKE for constraining rigid bodies of three or more centers.
It also allows rigid bodies to be linked with one or two common centers (e.g. peptide planes) by solving rigid body constraints iteratively in the same basic manner that SHAKE is used for atoms involving more than one SHAKE constraint.
An alternative constraint method, LINCS (Linear Constraint Solver) was developed in 1997 by Hess, Bekker, Berendsen and Fraaije,[24] and was based on the 1986 method of Edberg, Evans and Morriss (EEM),[25] and a modification thereof by Baranyai and Evans (BE).
This approximation only works for matrices with eigenvalues smaller than 1, making the LINCS algorithm suitable only for molecules with low connectivity.
