6.3 Molecular Dynamics Simulation
Fundamentally different methods for
(a) Hard bodies with impulsive collisons and
(b) Particles interacting via smooth potential forces.
Subsections
6.3.1 Hard Spheres / Hard disks
Set up the hard spheres/disks on a lattice; then assign
random initial velocities to the $N$ particles.
Adjust these velocities such that
- the total kinetic energy is consistent with some desired
temperature: $E_{k}=3NkT/2$
- the total momentum (conserved in the simulation) equals zero.
Note that the adjustment of the temperature is only temporary; it will
have to be repeated several times before thermal equilibrium is reached;
even then, $T$ will continue to fluctuate.
Now calculate, for each pair of particles
$(i,j)$ in the system, the time $t_{ij}$ it would take that pair to meet:
$
t_{ij} = \frac{\textstyle -b - \sqrt{b^{2}-v^{2}(r^{2}-d^{2})}}{\textstyle v^{2}}
$
where $d$ is the sphere diameter, $r$ is the distance between the centers of
$i$ and $j$, and
$
\begin{eqnarray}
b & = & ( r_{j}- r_{i}) \cdot ( v_{j}- v_{i})
\\
v & = & \left| ( v_{j}- v_{i})\right|
\end{eqnarray}
$
Note that the above formula for $t_{ij}$ derives from a straightforward
solution of the quadratic equation $r_{ij}^{2}(t)-d^{2}=0$.
As discussed in
the appendix, this solution may give rise to numerical errors. A more
secure alternative is, with the same meaning for
$b$ and $v$,
$
q = - \left[ b +sgn(b) \sqrt{b^{2}-v^{2}(r^{2}-d^{2})} \right],
\;\;\;\;t_{ij} =
min \left\{
\frac{\textstyle q}{\textstyle v^{2}},
\frac{\textstyle r^{2}-d^{2}}{\textstyle q} \right\}
$
(Since $b$ must be negative for a pair to meet in the future, we
may restrict the calculation to the case $sgn(b)<0$.)
Set up two arrays that contain, for each particle
$i$, the smallest positive collision time
$t(i)=min(t_{ij})$ and the next collision partner $j(i)$.
If a particle has no collision partner at positive times
we simply set $j(i)$ and $t(i)=[\infty]$, i.e. the largest
representable number.
This double loop over all $N(N-1)/2$ pairs need be performed only once,
at the start of the simulation.
Now find the smallest element in the table of free flight times,
calling it $t(i_{0})$. This is the time until the very next encounter
between two particles in the system. The indices of those two particles
are named $i_{0}$ and $j_{0}$.
During the time $t(i_{0})$
all particles perform a free flight, thus:
$
r_{i}\longrightarrow r_{i}+ v_{i} \cdot t(i_{0})
$
and all $t(i) \longrightarrow t(i)-t(i_{0})$.
Now an elastic collision between $i=i_{0}$ and
$j=j_{0}$ occurs, resulting in the new velocities
$
v_{i}'= v_{i} + \Delta v , \;\;\;
v_{j}'= v_{j} - \Delta v
$
where
$
\Delta v =b \frac{\textstyle r_{ij}}{\textstyle d^{2}}
$
Since $i_{0}$ and $j_{0}$ have new flight directions and speeds,
all pairwise collision times $t_{ij}$ involving these two must be
recalculated. This means that $2N-3$ pairs have to be scanned.
This completes the basic hard sphere/disc MD step. Now all
$t_{i}$ are once more searched for the smallest element, etc.
Here is the procedure in succinct form:
Molecular dynamics simulation of hard spheres:
Immediately after a collision, for each particle $i$ in the system the time
$t(i)$ to its next collision and the partner $j(i)$ at that collision is
assumed to be known.
-
Determine the smallest positive element $t(i_{0})$ among the $t(i)$,
identify the corresponding particle $i_{0}$ and its collision partner
$j_{0}\equiv j(i_{0})$.
-
Let all particles follow their free flight paths for a period $\Delta t \equiv t(i_{0})$;
subtract $\Delta t$ from each $t(i)$.
-
Perform the elastic collision between $i_{0}$ and $j_{0}$; after the
collision these spheres have the new velocities
$
v' = v \pm \Delta v , \;\;\;{\rm with} \;\;
\Delta v = b \frac{\textstyle r_{ij}}{\textstyle d^{2}}
$
-
Recalculate all times $t(i)$ that involve either $i_{0}$ or $j_{0}$,
i.e. for $i=i_{0}$, $i=j(i_{0})$, $i=j_{0}$, and $i=j(j_{0})$.
-
Go to (1).
At low densities the large free path may create problems with the periodic
boundary conditions, some particle suddenly appearing where it overlaps
another. One therefore limits the time allowed for free flight such that
for each particle and each coordinate $\alpha$ the free flight displacement
fulfills $\Delta x_{\alpha} \equiv v_{\alpha} \Delta t \leq L/2-d$.
|
|
Figure 6.5: Molecular dynamics of hard spheres
|
Project MD (Hard Disks):
(See also
here)
For a two-dimensional system of hard disks, write subroutines to
a) set up an initial configuration (simplest, though not best: square
lattice;) b) calculate $t(i)$ and $j(i)$; c) perform a pair collision.
Combine these subroutines into an MD code. To avoid the difficulty
mentioned at the end of the preceding figure, one might simply use
reflecting boundary conditions, doing a "billiard dynamics" simulation.
6.3.2 Continuous Potentials
For continously varying pair potentials we have for a particle
$i$ at any time $t$
$
\ddot{ r}_{i}(t) = \frac{\textstyle 1}{\textstyle m}
\sum \limits_{j\neq i} K_{ij}(t)\;\;\;\;{\rm with} \;\;\;\;
K_{ij}\equiv -\nabla_{i}u(r_{ij})
$
Considering the Lennard-Jones interaction, we find for the pair
force
$
K_{ij}= -24 \frac{\textstyle \epsilon}{\textstyle \sigma^{2}}
\left[ 2 \left(
\frac{\textstyle r_{ij}}{\textstyle \sigma}\right)^{-14}\right.
-\left. \left(\frac{\textstyle r_{ij}}{\textstyle \sigma}\right)^{-8}\right]
r_{ij}
$
where
$ r_{ij} \equiv r_{j}- r_{i}$.
The above-mentioned nearest image convention (NIC) is used in the
evaluation of the force acting on a particle.
Having determined this total force, the equation of motion
for particle $i$ may be numerically integrated.
A widely used technique is Verlet's algorithm
$
r_{i}(t_{n+1})= 2 r_{i}(t_{n})- r_{i}(t_{n-1})
+ b_{i}(t_{n})(\Delta t)^{2}
$
(with $ b_{i}\equiv \sum \limits_{j\neq i} K_{ij}/m$.)
PROJECT MD (LENNARD-JONES):
(See also
here)
Augment the subroutine module ENERGY such that it computes, for each
Lennard-Jones particle $i$ in the system, the total force exerted on it
by all other particles $j$:
$ K_{i} \equiv \sum \limits_{j \neq i} K_{ij}$, with $ K_{ij}$
as given above; remember to apply the nearest image convention.
Write a subroutine MOVE to integrate the equations of motion by a suitable
algorithm such as Verlet's. Having advanced each particle for one
time step, do not forget to apply periodic boundary conditions
to retain them all in the simulation box.
Write a main routine that puts the subroutines STARTCONF, ENERGY and
MOVE to work. Test your first MD code by monitoring the mechanically
conserved quantities.
Do a number of MD steps - say, $50-100$ - and average the
quantity $| v^{*}|^{2}$ to estimate the actual temperature.
To adjust the temperature to a desired value, scale all
velocity components of all particles in a suitable way. Repeat
this procedure up to 10 times. After $500 - 100$ steps the
fluid will normally be well randomized in space, and the temperature
will be steady - though fluctuating slightly.
6.3.3 Beyond Basic Molecular Dynamics
- orientation dependent potentials
- ionic or multipolar interactions
- polymer chains (geometrically constrained molecules)
- nonequilibrium dynamics (e.g. laminary flow); see
[EVANS 86],
[HOOVER 99]
- thermostatted dynamics (Nosé-Hoover)
Geometrical constraints / SHAKE method:
Internal geometrical constraints, e.g. rigid bond lengths
(see table
6.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}
\;\;\;{\rm and}\;\;\;
\sigma_{2}( r_{23})=\left| r_{23}\right|^{2}-l_{2}^{2}
$
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{\textstyle a_{1}}{\textstyle m_{1}} r_{12}
\\
\ddot{ r}_{2}&=& b_{2}+\frac{\textstyle 1}{\textstyle m_{2}}
\left[-a_{1} r_{12}+a_{2} r_{23}\right]
\\
\ddot{ r}_{3}&=& b_{3}-\frac{\textstyle a_{2}}{\textstyle m_{3}} r_{23}
\end{eqnarray}
$
where $ b_{1..3}$ are the physical accelerations due to Lennard-Jones
or other pair potentials.
Procedure:
- 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}'$.
- 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}
$
(6.4-6.6)
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}$.
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:
- 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 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:
- 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).
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 \limits_{i}
\frac{\textstyle mv_{i}^{2}}{\textstyle 2}\rangle
= d \frac{\textstyle NkT}{\textstyle 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}$:
$
\begin{eqnarray}
\dot{ v_{i}}&=&\frac{\textstyle 1}{\textstyle m} K_{i}-\xi v_{i}
\\
\dot{\xi} &=&\frac{\textstyle 2}{\textstyle Q}\left[ E_{kin}-3NkT_{0}/2 \right]
\end{eqnarray}
$
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].
vesely 2006
|