Franz J. Vesely > MolSim Tutorial > Geometrical Constraints / SHAKE;
 
 





 
< >
Vesely MolSim Tutorial




3.3 Geometrical constraints / SHAKE method:

Internal geometrical constraints, e.g. rigid bond lengths (see table 1.2):
  • Stiff harmonic bonds - uneconomical
  • SHAKE method by Ryckaert et al. [ RYCKAERT 77] - better


SHAKE: Consider the smallest non-trivial Kramers chain consisting of three atoms connected by massless rigid bonds of lengths $l_{1,2}$. Constraint equations:
$ \sigma_{1}( r_{12})=\left| r_{12}\right|^{2}-l_{1}^{2}=0 \;\;\;{\rm and}\;\;\; \sigma_{2}( r_{23})=\left| r_{23}\right|^{2}-l_{2}^{2}=0 $

with $ r_{12} \equiv r_{2}- r_{1}$ etc.

Lagrange constraint forces: keep the bond lengths constant. The constraint force on atom $1$ is parallel to $ r_{12}$. Atom $2$ is subject to two constraint forces along $- r_{12}$ and $ r_{23}$. Atom $3$ is kept in line by a force along $- r_{23}$:

$ \begin{eqnarray} \ddot{ r}_{1}&=& b_{1}+\frac{a_{1}}{m_{1}} r_{12} \\ \ddot{ r}_{2}&=& b_{2}+\frac{1}{m_{2}} \left[-a_{1} r_{12}+a_{2} r_{23}\right] \\ \ddot{ r}_{3}&=& b_{3}-\frac{a_{2}}{m_{3}} r_{23} \end{eqnarray} $

where $ b_{1..3}$ are the physical accelerations due to Lennard-Jones or other pair potentials.

Procedure:
  1. Let the positions $r_{i}^{n}$ be given at time $t_{n}$. Integrate the equations of motion for one time step with all a_{k}$ set to zero, i. e. without considering the constraint forces; denote the resulting preliminary positions (at time $n+1$) as $ r_{i}'$. These will not yet fulfill the constraint equations; instead, the values of $\sigma_{1}( r_{12}')$ and $\sigma_{2}( r_{23}')$ will have some nonzero values $\sigma_{1}'$, $\sigma_{2}'$.

  2. Now we make the correction ansatz

    $ \begin{eqnarray} r_{1}^{n+1}&=& r_{1}'+\frac{a_{1}}{m_{1}} r_{12}^{n}\\ r_{2}^{n+1}&=& r_{2}'+\frac{1}{m_{2}} \left[-a_{1} r_{12}^{n}+a_{2} r_{23}^{n}\right]\\ r_{3}^{n+1}&=& r_{3}'-\frac{a_{2}}{m_{3}} r_{23}^{n} \end{eqnarray} $ 3.1


    3.2


    3.3


    with undetermined $a_{k}$, requiring that the corrected positions fulfill the constraint equations:

    $ \begin{eqnarray} \sigma_{1}' - 2 \frac{a_{1}}{\mu_{12}} \left( r_{12}^{n}\cdot r_{12}'\right) + 2 \frac{a_{2}}{m_{2}} \left( r_{23}^{n}\cdot r_{12}'\right) +\left[ \dots \right]^{2} &=& 0\\ \sigma_{2}' + 2 \frac{a_{1}}{m_{2}} \left( r_{12}^{n}\cdot r_{23}'\right) - 2 \frac{a_{2}}{\mu_{23}} \left( r_{23}^{n}\cdot r_{23}'\right) +\left[ \dots \right]^{2} &=& 0 \end{eqnarray} $


    with $1/\mu_{12} \equiv 1/m_{1}+1/m_{2}$; the terms $\left[ \dots \right]^{2}$ are quadratic in $a_{1,2}$.

  3. Instead of solving these two quadratic equations for the unknowns $a_{1,2}$ we ignore, for the time being, the small quadratic terms. The remaining linear equations are solved iteratively, meaning that this system of linear equation is solved to arrive at an improved estimate for $a_{1,2}$ which is again inserted in 3.1-3.3 leading to a new set of linearized equations etc., until the absolute values of $a_{1,2}$ are negligible; generally, this will occur after a very few iterations.

    Solving the linearized equations involves a matrix inversion. To avoid this we introduce one more simplification:
    In passing through the chain from one end to the other we consider only one constraint per atom:
    • First the bond $ r_{12}$ is repaired by displacing $1$ and $2$.
    • By repairing the next bond $r_{23}$ we disrupt the first bond again; this is accepted in view of further iterations.
    • By going through the chain several times we reduce both the error introduced by neglecting the quadratic terms and the error due to considering only one constraint at a time.


    In our case the procedure is:
    $ a_{1}=\frac{\textstyle \mu_{12}}{\textstyle 2} \frac{\textstyle \sigma_{1}'}{\textstyle r_{12}'\cdot r_{12}^{n}} \;\;\;\;\; a_{2}=\frac{\textstyle \mu_{23}}{\textstyle 2} \frac{\textstyle \sigma_{2}'}{\textstyle r_{23}'\cdot r_{23}^{n}} $

    insert this in 3.1-3.3 and iterate until $a_{1,2}$ are negligible.
The technique is easily generalized to longer chains.

Molecules and robots:
Robot arms made up of several successive links and joints bear some resemblance to chain molecules. A standard problem of robotics, the inverse kinematic problem, may therefore be solved a la simulation.[KASTENMEIER 86]

The principle of the so-called constrained dynamics/stochastic optimization technique is as follows:
  • Define a required "world trajectory"for the final element in the chain.
  • If the robot is redundant, meaning that it has more degrees of freedom (links and joints) than necessary, the problem is underdetermined: different combinations of movements by the individual joints produce identical paths of the final member.
  • This redundancy may be exploited to fulfill additional requirements, such as an overall minimum of angular accelerations (i. e. mechanical wear) in the joints, avoidance of obstacles, etc.
  • The last element is now moved, one time step at a time, along its requested trajectory, and SHAKE is invoked to have the preceding joints follow.
  • The additional requirements are combined into a cost function which is minimized, at each time step, using a stochastic search (simulated annealing or simple random search).



Constraint dynamics with one bond length constraint

Molecular Dynamics for Linear Diatomics: Constraint method

The classical method to simulate diatomics is due to Singer (1977) and Fincham (1984), as also described in Allen-Tildesley, Computer Simulation of Liquids (Clarendon, Oxford 1987).

However, we prefer to use the method of "Constraint Dynamics" - which is particularly simple when applied to one single constraint.

Assume a diatomic with potential forces $\vec{F}_{1,2}$ acting on the two atoms, respectively, and write $\vec{b}_{1} \equiv \vec{F}_{1}/m_{1}$ etc. The only constraint is $\sigma(\vec{r}_{12}) = \vert\vec{r}_{12}\vert^{2}-l^{2}=0$. The equations of motion are then

$ \begin{eqnarray} \ddot{\vec{r}_{1}}&=&\vec{b}_{1}+\frac{\lambda}{m_{1}} \vec{r_{12}} \\ \ddot{\vec{r}_{2}}&=&\vec{b}_{2}-\frac{\lambda}{m_{2}} \vec{r_{12}} \end{eqnarray} $


True to the idea of constraint dynamics, proceed as follows:
  1. Given the positions $\vec{r}_{i}^{n}$ at time $t_{n}$, integrate the equations of motion for one time step only considering the potential forces. The resulting positions are denoted $\vec{r}_{i}{'}$. These preliminary position vectors will not fulfill the constraint equation; rather, the value of $\sigma' \equiv \sigma (\vec{r}_{i}{'})$ will be non-zero.

  2. Making the correction ansatz

    $ \begin{eqnarray} \vec{r}_{1}^{n+1}&=&\vec{r}_{1}{'}+\frac{\lambda}{m_{1}}\vec{r}_{12}^{n}\\ \vec{r}_{2}^{n+1}&=&\vec{r}_{2}{'}-\frac{\lambda}{m_{2}}\vec{r}_{12}^{n} \end{eqnarray} $ 3.4


    3.5


    with undetermined $\lambda$, requiring that the corrected positions fulfill the constraint equations we find

    $ \sigma' - 2 \frac{\textstyle \lambda}{\textstyle \mu} \left(\vec{r}_{12}^{n}\cdot\vec{r}_{12}{'}\right) + \frac{\textstyle \lambda^{2}}{\textstyle \mu^{2}} l^{2}=0 $

    with $1/\mu_{12} \equiv 1/m_{1}+1/m_{2}$.


  3. This equation for $\lambda/\mu$ may be solved exactly:

    $ \frac{\textstyle \lambda}{\textstyle \mu} = \frac{\textstyle \left(\vec{r}_{12}^{n}\cdot\vec{r}_{12}{'}\right)}{\textstyle l^{2}} \pm \sqrt{\frac{\textstyle \left(\vec{r}_{12}^{n}\cdot\vec{r}_{12}{'}\right)^{2}}{\textstyle l^{4}} - \frac{\textstyle \sigma'}{\textstyle l^{2}}} $
    where the negative sign prevails.



Simple example: Free rotation of a 2D dumbbell

For this problem the constraint method should be exact, since neither the Verlet propagator nor the constraint reconstruction introduce errors.

At $t_{n}=0$ let $\vec{r}_{1}=(0 / -0.5)$, $\vec{r}_{2}=(0 / 0.5)$; the masses are $m_{1}=m_{2}=1$, the angular velocity is $\omega=1$, and a time step of $\Delta t=0.1$ is assumed. Then the exact positions at time $\Delta t$ are

$ \vec{r}_{1}(\Delta t) = \left( \begin{array}{cc}c&-s\\s&c\end{array} \right) \cdot \left( \begin{array}{r} 0 \\ -0.5 \end{array} \right) =\left( \begin{array}{r} 0.049917 \\-0.497502 \end{array} \right) $

and

$ \vec{r}_{2}(\Delta t) = \left( \begin{array}{cc}c&-s\\s&c\end{array} \right) \cdot \left( \begin{array}{r} 0 \\0.5 \end{array} \right) =\left( \begin{array}{r} -0.049917 \\0.497502 \end{array} \right) $

where we have written $c \equiv \cos \omega \Delta t$ and $s \equiv \sin \omega \Delta t$.

Now let us test the constraint dynamics procedure against this exact result. For the Verlet integrator we need $\vec{r}_{i}(- \Delta t)$, for which we take the exact values $\vec{r}_{1}(-\Delta t) =\{-0.049917 /-0.497502\}$ and $\vec{r}_{2}(-\Delta t) =\{0.049917 /0.497502\}$. After the free flight the new positions are

$ \vec{r}_{1}{'} =\left( \begin{array}{r} 0.049917 \\-0.502498 \end{array} \right) \;\;\;\; \vec{r}_{2}{'} =\left( \begin{array}{r} -0.049917 \\0.502498 \end{array} \right) $

with $\vec{r}_{12}{'}=\left( \begin{array}{r} -0.099834 \\1.05000 \end{array} \right)$ and $\sigma'=0.019983$. Using $\vec{r}_{12}{'} \cdot \vec{r}_{12}(0) = 1.004996$ and $l=1$ we find for $\lambda/\mu$

$ \frac{\textstyle \lambda}{\textstyle \mu} = 1.00996 - \sqrt{1.004996^{2}-0.01998} = 0.009992 $

Since $\mu=1/2$ we have $\lambda=0.004996$. Insert this in the correction formula 3.4-3.5 to find

$ \vec{r}_{1}(\Delta t) =\left( \begin{array}{r} 0.049917 \\-0.497502 \end{array} \right) \;\;\;\; \vec{r}_{2}(\Delta t) =\left( \begin{array}{r} -0.049917 \\0.497502 \end{array} \right) $

which is identical to the exact values.


vesely may-06

< >