next up previous

4.3.1 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
$\displaystyle \ddot{\vec{r}_{1}}$ $\textstyle =$ $\displaystyle \vec{b}_{1}+\frac{\lambda}{m_{1}}   \vec{r_{12}}$ (4.6)
$\displaystyle \ddot{\vec{r}_{2}}$ $\textstyle =$ $\displaystyle \vec{b}_{2}-\frac{\lambda}{m_{2}}   \vec{r_{12}}$ (4.7)

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
    $\displaystyle \vec{r}_{1}^{ n+1}$ $\textstyle =$ $\displaystyle \vec{r}_{1}{'}+\frac{\lambda}{m_{1}}\vec{r}_{12}^{ n}$ (4.8)
    $\displaystyle \vec{r}_{2}^{ n+1}$ $\textstyle =$ $\displaystyle \vec{r}_{2}{'}-\frac{\lambda}{m_{2}}\vec{r}_{12}^{ n}$ (4.9)

    with undetermined $\lambda$, requiring that the corrected positions fulfill the constraint equations we find
- 2 \frac{\lambda}{\mu} \left(\vec{r}_{12}^{ n}\cdot\vec{r}_{12}{'}\right)
+ \frac{\lambda^{2}}{\mu^{2}} l^{2}=0
\end{displaymath} (4.10)

    with $1/\mu_{12} \equiv 1/m_{1}+1/m_{2}$.
This equation for $\lambda/\mu$ may be solved exactly:
\frac{\lambda}{\mu} =
\frac{\left(\vec{r}_{12}^{ n}\cdot\v...\vec{r}_{12}{'}\right)^{2}}{l^{4}}
- \frac{\sigma'}{l^{2}}}
\end{displaymath} (4.11)

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) = \mbox{$\left( \begin{array}{cc}c&-s\...
...} 0.049917 \ \vspace{-9 pt}\ -0.497502 \end{array} \right)$}
\end{displaymath} (4.12)

\vec{r}_{2}(\Delta t) = \mbox{$\left( \begin{array}{cc}c&-s\...
...} -0.049917 \ \vspace{-9 pt}\ 0.497502 \end{array} \right)$}
\end{displaymath} (4.13)

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}{'} =\mbox{$\left( \begin{array}{r} 0.049917 \ \...
...} -0.049917 \ \vspace{-9 pt}\ 0.502498 \end{array} \right)$}
\end{displaymath} (4.14)

with $\vec{r}_{12}{'}=\mbox{$\left( \begin{array}{r} -0.099834 \ \vspace{-9 pt}\ 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{\lambda}{\mu} = 1.00996 - \sqrt{1.004996^{2}-0.01998} = 0.009992
\end{displaymath} (4.15)

Since $\mu=1/2$ we have $\lambda=0.004996$. Insert this in the correction formula 4.8-4.9 to find
\vec{r}_{1}(\Delta t) =\mbox{$\left( \begin{array}{r} 0.0499...
...} -0.049917 \ \vspace{-9 pt}\ 0.497502 \end{array} \right)$}
\end{displaymath} (4.16)

which is identical to the exact values.

Thermostats / Nosé-Hoover method:
In a nonequilibrium process work must be invested which systematically heats up the system. In lab experiments a thermostat takes care of this; in simulations the problem of thermostating is non-trivial. The temperature of a molecular dynamics sample is not an input parameter to be manipulated at will; rather, it is a quantity to be ``measured'' in terms of an average of the kinetic energy of the particles,

\langle E_{kin}\rangle \equiv \langle \sum_{i}\frac{mv_{i}^{2}}{2}\rangle
= d  \frac{NkT}{2}

($d \dots $ dimension). Suggestions as to how one can maintain a desired temperature in a dynamical simulation range from repeatedly rescaling all velocities (``brute force thermostat'') to adding a suitable stochastic force acting on the molecules. Such recipes are unphysical and introduce an artificial trait of irreversibility and/or indeterminacy into the microscopic dynamics.

A interesting deterministic method of thermostating a simulation sample was introduced by Shuichi Nosé and William Hoover. Under very mild conditions the following augmented equations of motion will lead to a correct sampling of the canonical phase space at a given temperature $T_{0} $:
$\displaystyle \dot{\mbox{$\bf v$}_{i}}$ $\textstyle =$ $\displaystyle \frac{1}{m}\mbox{$\bf K$}_{i}-\xi   \mbox{$\bf v$}_{i}$  
$\displaystyle \dot{\xi}$ $\textstyle =$ $\displaystyle \frac{2}{Q}\left[ E_{kin}-3NkT_{0}/2 \right]$  

The coupling parameter $Q$ describes the inertia of the thermostat. The generalized viscosity $\xi (t)$ is temporally varying and may assume negative values as well.

For systems of many degrees of freedom (particles) this is sufficient to produce a pseudo-canonical sequence of states. For low-dimensional systems it may be necessary to use two NH thermostats in tandem [MARTYNA 92].

next up previous
F. J. Vesely / University of Vienna