next up previous
Next: 6.4 Evaluation of Simulation Up: 6.3 Molecular Dynamics Simulation Previous: 6.3.2 Continuous Potentials


6.3.3 Beyond Basic Molecular Dynamics



Geometrical constraints / SHAKE method:
Internal geometrical constraints, e.g. rigid bond lengths (see table 6.2):

SHAKE: Consider the smallest non-trivial Kramers chain consisting of three atoms connected by massless rigid bonds of lengths $l_{1,2}$. Constraint equations:

\begin{displaymath}
\sigma_{1}(\mbox{$\bf r$}_{12})=\left\vert \mbox{$\bf r$}_{1...
...23})=\left\vert \mbox{$\bf r$}_{23}\right\vert^{2}-l_{2}^{2}=0
\end{displaymath}

with $\mbox{$\bf r$}_{12} \equiv \mbox{$\bf r$}_{2}-\mbox{$\bf r$}_{1}$ etc.

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

$\displaystyle \ddot{\mbox{$\bf r$}}_{1}$ $\textstyle =$ $\displaystyle \mbox{$\bf b$}_{1}+\frac{a_{1}}{m_{1}}\mbox{$\bf r$}_{12}$  
$\displaystyle \ddot{\mbox{$\bf r$}}_{2}$ $\textstyle =$ $\displaystyle \mbox{$\bf b$}_{2}+\frac{1}{m_{2}}
\left[-a_{1}\mbox{$\bf r$}_{12}+a_{2}\mbox{$\bf r$}_{23}\right]$  
$\displaystyle \ddot{\mbox{$\bf r$}}_{3}$ $\textstyle =$ $\displaystyle \mbox{$\bf b$}_{3}-\frac{a_{2}}{m_{3}}\mbox{$\bf r$}_{23}$  

where $\mbox{$\bf b$}_{1..3}$ are the physical accelerations due to Lennard-Jones or other pair potentials.

Procedure:
  1. Let the positions $\mbox{$\bf 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 $\mbox{$\bf r$}_{i}'$. These will not yet fulfill the constraint equations; instead, the values of $\sigma_{1}(\mbox{$\bf r$}_{12}')$ and $\sigma_{2}(\mbox{$\bf r$}_{23}')$ will have some nonzero values $\sigma_{1}'$, $\sigma_{2}'$.
  2. Now we make the correction ansatz
    $\displaystyle \mbox{$\bf r$}_{1}^{n+1}$ $\textstyle =$ $\displaystyle \mbox{$\bf r$}_{1}'+\frac{a_{1}}{m_{1}}\mbox{$\bf r$}_{12}^{n}$ (6.4)
    $\displaystyle \mbox{$\bf r$}_{2}^{n+1}$ $\textstyle =$ $\displaystyle \mbox{$\bf r$}_{2}'+\frac{1}{m_{2}}
\left[-a_{1}\mbox{$\bf r$}_{12}^{n}+a_{2}\mbox{$\bf r$}_{23}^{n}\right]$ (6.5)
    $\displaystyle \mbox{$\bf r$}_{3}^{n+1}$ $\textstyle =$ $\displaystyle \mbox{$\bf r$}_{3}'-\frac{a_{2}}{m_{3}}\mbox{$\bf r$}_{23}^{n}$ (6.6)

    with undetermined $a_{k}$, requiring that the corrected positions fulfill the constraint equations:
    $\displaystyle \sigma_{1}'
- 2 \frac{a_{1}}{\mu_{12}} \left(\mbox{$\bf r$}_{12}^...
...box{$\bf r$}_{23}^{n}\cdot\mbox{$\bf r$}_{12}'\right)
+\left[ \dots \right]^{2}$ $\textstyle =$ $\displaystyle 0$ (6.7)
    $\displaystyle \sigma_{2}'
+ 2 \frac{a_{1}}{m_{2}} \left(\mbox{$\bf r$}_{12}^{n}...
...box{$\bf r$}_{23}^{n}\cdot\mbox{$\bf r$}_{23}'\right)
+\left[ \dots \right]^{2}$ $\textstyle =$ $\displaystyle 0$ (6.8)

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

    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 6.4-6.6 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:



    In our case the procedure is:

    \begin{displaymath}
a_{1}=\frac{\mu_{12}}{2}\,
\frac{\sigma_{1}'}{\mbox{$\bf r$...
...\sigma_{2}'}{\mbox{$\bf r$}_{23}'\cdot\mbox{$\bf r$}_{23}^{n}}
\end{displaymath}

    insert this in 6.4-6.6 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:

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,

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

($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
Next: 6.4 Evaluation of Simulation Up: 6.3 Molecular Dynamics Simulation Previous: 6.3.2 Continuous Potentials
Franz J. Vesely Oct 2005
See also:
"Computational Physics - An Introduction," Kluwer-Plenum 2001