# Constrained molecular dynamics

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

${\displaystyle {\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, ${\displaystyle {\mathcal {L}}^{*}}$ is the Lagrangian for the extended system, and λi is a Lagrange multiplier associated with a geometric constraint σi:

${\displaystyle \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):
${\displaystyle v_{i}^{t+{\Delta }t/2}=v_{i}^{t-{\Delta }t/2}+{\frac {a_{i}^{t}}{m_{i}}}{\Delta }t}$
${\displaystyle q_{i}^{t+{\Delta }t}=q_{i}^{t}+v_{i}^{t+{\Delta }t/2}{\Delta }t}$
2. Use the new positions q(tt) to compute Lagrange multipliers for all constraints:
${\displaystyle {\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):
${\displaystyle v_{i}^{t+{\Delta }t/2}=v_{i}^{t-{\Delta }t/2}+\left(a_{i}^{t}-\sum _{k}{\frac {{\lambda }_{k}}{m_{i}}}\bigtriangledown _{i}{\sigma }_{k}(q^{t})\right){\Delta }t}$
${\displaystyle q_{i}^{t+{\Delta }t}=q_{i}^{t}+v_{i}^{t+{\Delta }t/2}{\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.

## How to

Geometric constraints are introduced by defining one or more entries with the STATUS parameter set to 0 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 (VASP would terminate with an error message). Constraints can be combined with restraints, time-dependent bias potentials (Metadynamics), monitored coordinates and other elements available within the context of MD.