Constrained molecular dynamics

Constrained molecular dynamics is performed using the SHAKE algorithm.. In this algorithm, the Lagrangian for the system $\mathcal{L}$ is extended as follows:

$\mathcal{L}^*(\mathbf{q,\dot{q}}) = \mathcal{L}(\mathbf{q,\dot{q}}) + \sum_{i=1}^{r} \lambda_i \sigma_i(q),$

where the summation is over r geometric constraints, $\mathcal{L}^*$ is the Lagrangian for the extended system, and λ_i is a Lagrange multiplier associated with a geometric constraint σ_i:

$\sigma_i(q) = \xi_i({q})-\xi_i \;$

with ξ_i(q) being a geometric parameter and ξ_i is the value of ξ_i(q) fixed during the simulation.

In the SHAKE algorithm, the Lagrange multipliers λ_i are determined in the iterative procedure:

1. Perform a standard MD step (leap-frog algorithm):

   $v^{t+{\Delta}t/2}_i = v^{t-{\Delta}t/2}_i + \frac{a^{t}_i}{m_i} {\Delta}t$

   $q^{t+{\Delta}t}_i = q^{t}_i + v^{t+{\Delta}t/2}_i{\Delta}t$

2. Use the new positions q(t+Δt) to compute Lagrange multipliers for all constraints:

   ${\lambda}_k= \frac{1}{{\Delta}t^2} \frac{\sigma_k(q^{t+{\Delta}t})}{\sum_{i=1}^N m_i^{-1} \bigtriangledown_i{\sigma}_k(q^{t}) \bigtriangledown_i{\sigma}_k(q^{t+{\Delta}t})}$

3. Update the velocities and positions by adding a contribution due to restoring forces (proportional to λ_k):

   $v^{t+{\Delta}t/2}_i = v^{t-{\Delta}t/2}_i + \left( a^{t}_i-\sum_k \frac{{\lambda}_k}{m_i} \bigtriangledown_i{\sigma}_k(q^{t}) \right ) {\Delta}t$

   $q^{t+{\Delta}t}_i = q^{t}_i + v^{t+{\Delta}t/2}_i{\Delta}t$

4. repeat steps 2-4 until either |σ_i(q)| are smaller than a predefined tolerance (determined by SHAKETOL), or the number of iterations exceeds SHAKEMAXITER.

Ho to

Geometric constraints are introduced by defining one or more entries with the STATUS parameter set to 0d in the ICONST-file. Constraints can be used within a standard NVT or NpT MD setting introduced by MDALGO=1|2|3. Note that fixing geometric parameters related to lattice vectors is not allowed within an NVT simulation. Constraints can be combined with restraints, time-dependent bias potentials (Metadynamics)

---