3.4 Molecular Dynamics in other Ensembles
NVT, NPH, NPT, ...
In the basic MD setting the macroscopic observables
$N,V,E,\vec{p}$ are held constant (where $\vec{p}$
is the total linear momentum.)
This defines the so-called "MD ensemble" which can be shown to be
thermodynamically equivalent to the other ensembles.
Other ensembles may be of advantage for specific problems:
$NVT$ is useful for comparisons with theory and MC simulations;
non-equilibrium simulations require a thermostat, thus $T=const$;
phase transitions are best studied in an $NPT$ ensemble, etc.
Subsections
NVT Molecular Dynamics
Initially (around 1980) stochastic procedures were used, in which the
particle velocities were regularly reset to random values picked from
the appropriate Maxwell-Boltzmann distribution. [Andersen 1980,
Heyes 1983]
Later, deterministic and even reversible equations of motion
were developed that nevertheless serve to roam phase space in a
more or less "canonical" manner. These are known as
Gaussian and Nose-Hoover dynamics, respectively.
Gauss dynamics
Gaussian, or isokinetic, dynamics was introduced by
Hoover and Evans in 1982. This is the idea:
Let
$
\sigma\{ \vec{v}\} \equiv \sum_{i=1}^{N}
\frac{\textstyle m}{\textstyle 2}v_{i}^{2}-E_{0}
$
be a (non-holonomic) constraint to the Newtonian motion of the
$N$
-particle system. This is only one constraint equation for
$3N$
variables, thus we have some freedom of choice.
A dynamical constraint may always be put in terms of
constraint forces. A particularly economical and "impartial"
prescription for the definition of these forces was given by
C. F. Gauss in 1829:
- From $\sigma\{ \vec{v}\}=0$ for all times we have
$
0=\dot{\sigma}
=\sum_{i}\left( \frac{\textstyle \partial}{\textstyle \partial \vec{v}_{i}}
\sigma \right) \cdot \dot{\vec{v}_{i}}
= m \sum_{i} \vec{v}_{i} \cdot \dot{\vec{v}_{i}}
\equiv m \vec{v} \cdot \dot{\vec{v}}
$
Since
$\dot{\vec{v}}=\left( \vec{K}+\vec{Z}\right)/m$
(with the Newtonian forces $\vec{K}$
and the yet unknown constraint forces
$\vec{Z}$) we have
$
\vec{v} \cdot \left( \vec{K}+\vec{Z}\right) = 0
$
- Among all possible sets of constraint forces let us choose the one with
the smallest norm. Requiring $\left\vert \vec{Z}^{0}\right\vert^{2}=min$
we find, upon variation, that
$\vec{Z}^{0} \cdot \delta \vec{Z}=0$. Since
for all physically allowed variations we have
$\delta \vec{Z} \cdot \vec{v}=0$
we conclude that
$\vec{Z}^{0}=-\lambda \vec{v}$. Inserting this in
$\vec{v} \cdot \left( \vec{K}+\vec{Z}^{0}\right) = 0$
we find
$
\vec{Z}_{i}^{0}=
- \frac{\textstyle \sum_{j}\vec{v}_{j}\cdot \vec{K}_{j}}{\textstyle \sum_{j}v_{j}^{2}} \vec{v}_{i}
= - \frac{\textstyle \vec{v}\cdot\vec{K}}{\textstyle |\vec{v}|^{2}} \vec{v}
$
- The equations of motion according to Gauss and Hoover are thus
$ \begin{eqnarray}
\dot{\vec{r}_{i}} &=& \vec{v}_{i}
\\
\dot{\vec{v}_{i}}&=& \frac{1}{m}
\left[ \vec{K}_{i}
- \frac{\sum_{j}\vec{v}_{j}\cdot \vec{K}_{j}}{\sum_{j}v_{j}^{2}} \vec{v}_{i}
\right]
\end{eqnarray}
$
These e.o.m. lead to a deterministic trajectory with
$E_{kin}=const$.
Note that $\vec{v}$ and $\vec{K}$
are needed at the same time step;
thus the usual Verlet algorithm is not applicable here. Instead, a
Predictor-Corrector algorithm or Runge-Kutta should be used.
Nosé-Hoover dynamics
Thermostats:
The constraint force used in isokinetic dynamics may be thought of as
the action of a particularly efficient thermostat, in which
not only the average kinetic energy (i. e. the temperature)
but even the momentary value of $E_{kin}$ is kept constant at all times.
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{\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. Nosé
had proposed a set of equations of motion that allow for a certain degree
of variation in the total energy of the system. He could show that with
these dynamics, and under a few very general conditions, the spread around
the given $E_{0}$ is such that a canonical distribution is approached.
Hoover rewrote Nosé's original equations of motion in the following
simple way:
$ \begin{eqnarray}
\dot{\vec{r}} &=& \frac{\vec{p}}{m}
\\
\dot{\vec{p}}&=&
\vec{K}_{i} - \zeta \vec{p}
\\
\dot{\zeta}&=& \frac{2}{Q}
\left[ E_{kin} - E_{0} \right]
\end{eqnarray} $
where we have omitted the particle index
$i$. The coupling parameter
$Q^{-1}$
represents the efficiency of the thermal bath holding the
temperature constant.
$\zeta$
is a kind of "viscosity" which,
however, may take on negative values.
For systems of many degrees of freedom (particles) this is sufficient
to produce a pseudo-canonical sequence of states.
(See also Posch, Hoover, Vesely Phys.Rev.A 33(1986)4253.)
For low-dimensional systems it may be necessary to use two NH thermostats
in tandem [MARTYNA 92].
For
$Q \rightarrow 0$ the spread of $E_{kin}$ around $E_{0}$ approaches
zero, and the Nosé-Hoover dynamics becomes identical to Gauss dynamics.
NPH Molecular Dynamics
Andersen (1980) introduced an additional ("synthetic") energy
that was coupled to changes of the system volume. Putting
$\Delta E_{pot}=P_{0}V$ and $\Delta E_{kin}=Q\dot{V}^{2}/2$
(with a generalized mass $Q$ which may be visualized as the mass of some
piston) he wrote the Hamiltonian as
$
H=E_{pot}+\Delta E_{pot}+E_{kin}+\Delta E_{kin}
$
Introducing scaled position vectors $\vec{s}_{i} \equiv \vec{r}_{i} V^{-1/3}$
he derived the equations of motion
$
\begin{eqnarray}
\ddot{\vec{s}_{i}}&=&\frac{\vec{K}_{i}}{m V^{1/3}}
- \frac{2}{3} \frac{\dot{V}}{V}\dot{\vec{s}_{i}}
\\
\ddot{V}&=&\left( P-P_{0}\right)/Q
\end{eqnarray}
$
These equations of motion conserve the enthalpy, as appropriate in
an isobaric ensemble.
Parrinello and Rahman (1980) extended the NPH method of Andersen
to allow for non-isotropic stretching and shrinking of box sides.
Important applications are structural phase transitions in solids.
Scaling of the position vectors now follows the equation
$\vec{r}_{i}= H \cdot \vec{s}_{i}$ instead of $=V^{1/3} \vec{s}_{i}$.
The matrix
$ H\equiv \left\{ \vec{h}_{1},\vec{h}_{2},\vec{h}_{3}\right\}$
describes the anisotropic transformation of the basic cell. The cell
volume is
$
V \equiv \left| H\right|
= \vec{h}_{1} \cdot \left[\vec{h}_{2} \times \vec{h}_{3}\right]
$
The additional terms in the Hamiltonian are
$
\Delta E_{pot} = P_{0} V\;\;\;\;{\rm and} \;\;\;\;
\Delta E_{kin} =
\frac{\textstyle Q}{\textstyle 2}\sum_{\alpha}\sum_{\beta}\dot{H}_{\alpha\beta}^{2}
$
and the equations of motion are
$ \begin{eqnarray}
\ddot{s}_{i}&=&\frac{1}{m} H^{-1} \cdot \vec{K}_{i}
- G^{-1}\cdot \dot{ G} \cdot \dot{\vec{s}}
\\
\ddot{ H}&=&\frac{1}{Q} \left[ P- I P_{0}\right]
V ( H^{-1})^{T}
\end{eqnarray} $
where
$ G= H^{T} \cdot H$ is a metric tensor,
$Q$ is a virtual mass (of dimension mass), and the
pressure/stress tensor is defined as
$
( P)_{\alpha \beta} = \frac{1}{V}
\left[ \sum_{i}m\left( H \cdot \vec{s}_{i}\right)_{\alpha}^{T}
\left( H \cdot \vec{s}_{i}\right)_{\beta}
+ \frac{1}{2}\sum_{i}\sum_{j}
\left( H\cdot \vec{s}_{ij}\right)_{\alpha}
\left( \vec{K}_{ij}\right)_{\beta}
\right]
$
Morriss and Evans (1983) devised another type of NPH dynamics.
They suggested to constrain the pressure not by an inert piston but
by a generalized constraint force in the spirit of Gaussian dynamics.
The same idea may be carried over to constrain
$E_{kin}$ in addition to the pressure. In this way one arrives at a
$NPT$ molecular dynamics procedure. (See Allen-Tildesley, Ch. 7.)
NPT Molecular Dynamics
The $NPT$ dynamics of Morriss and Evans (see end of previous
subsection) was actually a $NPE_{kin}$ dynamics, with zero fluctuation
of the kinetic energy. A genuine $NPT$ dynamics would constrain
$E_{kin}$ to stay near a desired average but with - hopefully
gaussian distributed - fluctuations around that value.
Such a method was conceived by Hoover (1985).
Extending the Hoover-Nosé idea to keep both the pressure and the
kinetic energy near a desired value $P_{0}$, he wrote the equations
of motion as
$ \begin{eqnarray}
\dot{\vec{s}_{i}} &=& \vec{v}_{i} V^{-1/3}
\\
\dot{\vec{v}_{i}} &=& \frac{1}{m}\vec{K}_{i}-(\chi+\zeta) \vec{v}_{i}
\\
\dot{\zeta} &=& \frac{2}{Q} \left[ E_{kin}-E_{0} \right]
\\
\chi &=& \frac{\dot{V}}{3V}
\\
\dot{\chi} &=& \frac{1}{t_{P}^{2} kT} \left[ P-P_{0} \right]
\end{eqnarray} $
where $t_{P}$ and $Q$ are relaxation times of the pressure and
temperature regulators, respectively.
Exercises
- Gauss Dynamics:
- Write a Gaussian MD program to simulate the following system:
3 particles without pairwise interaction are moving in a
one-dimensional harmonic potential. The sum of their kinetic energies
is kept constant. Draw the $x,\dot{x}$
trajectory of the first particle. Use Runge-Kutta.
- Nosé-Hoover thermostat:
- The equations of motion of the Nosé oscillator with
$k=m=E_{0}=1$ are
$ \begin{eqnarray}
\dot{x}&=&v
\\
\dot{v}&=&-x-\zeta v
\\
\dot{\zeta}&=&\frac{\textstyle 1}{\textstyle Q} \left[ v^{2}-2\right]
\end{eqnarray} $
Integrate this set of equations, using Runge-Kutta, and draw
the $x,v$ trajectory.
vesely may-06