Intramolecular bonds, if
$kT$ is comparable to the bond energy

Born-Huggins-Mayer

$u(r)=\frac{\textstyle{q_{1}q_{2}}}{\textstyle{4\pi \epsilon_{0}r}}
+ B e^{\textstyle{-\alpha r}} -
\frac{\textstyle{C}}{\textstyle{r^{6}}}-
\frac{\textstyle{D}}{\textstyle{r^{8}}}$

Ionic melts;
$q_{i}$ are the ion charges

Table 6.2:
Anisotropic model potentials in statistical-mechanical simulation:
$u(12)=u( r_{12}, e_{1}, e_{2})$

Hard dumbbells,
hard spherocylinders, etc.

$u(12)= \infty \;\;$
if overlap
$= 0 \;\;\;$
otherwise

First approximation to rigid molecules

Interaction site models, rigid

$u(12)=$
sum of isotropic pair energies
$\textstyle{u(r_{i(1),j(2)})}$,
where several interaction sites $i$ and $j$
are in fixed positions on molecules $1$ and $2$, respectively

Rigid molecules

Interaction site models with non-rigid bonds

$u(12)=$ sum of isotropic pair energies, both intra- and intermolecular

Non-rigid molecules

Kramers-type

$u(12)=$ sum of isotropic pair energies, exclusively between sites on different
molecules; certain intramolecular distances (bonds) and/or angles
are fixed

Flexible molecules, from ethane to biopolymers

Stockmayer

$u(12)=$ Lennard-Jones + point dipoles

First approximation to small polar molecules

Anisotropic Kihara

$u(12)=4 \epsilon \left[ \left(
\frac{\textstyle{\rho_{12}}}{\textstyle{\sigma}}\right)^{-12}
- \left(
\frac{\textstyle{\rho_{12}}}{\textstyle{\sigma}}\right)^{-6} \right]$
where
$\rho_{12}$ is the shortest distance between two linear rods

Rigid linear molecules
with distributed
Lennard-Jones interaction

Gay-Berne

$u(12)=$
$4 \epsilon(12)
\left[ \left(
\frac{\textstyle{r_{12}-\sigma(12)+\sigma_{0}}}{\textstyle{\sigma_{0}}}\right)^{-12}
\right.$
$- \left. \left(
\frac{\textstyle{r_{12}-\sigma(12)+\sigma_{0}}}{\textstyle{\sigma_{0}}}\right)^{-6} \right]$
where
$\sigma(12)$ and $\epsilon(12)$ depend on
$ r_{12}, e_{1},e_{2}$ and certain substance-specific
shape parameters

Liquid crystal molecules of (roughly) ellipsoidal shape.

At time $t$, combine the position vectors of the
$N$ atoms into a vector $\Gamma_{c}\equiv$
$\{ r_{1} \dots r_{N}\}$.
The set of all possible such vectors spans the
$3N$ -dimensional "configuration space"
$\Gamma_{c}$.

Let $a(\Gamma_{c})$ be a property of the $N$
-body system, depending on the positions of all particles
(i.e. of the microstate $\Gamma_{c}$).
The thermodynamic average of $a$ is given by

Problem: the probability density
$p(\Gamma_{c})$ is known only up to an indetermined normalizing factor.
Thus, $p_{can}(\Gamma_{c})$ $\propto exp[-E(\Gamma_{c})/kT]$
but the normalizing denominator $Q$ (the Configurational Partition
Function) in

Basic model of ferromagnetic solids:
atoms are fixed on the vertices of a lattice. Their
dipole vectors (spins) may have varying directions: either
up/down (Ising model) or any direction (Heisenberg model).

Microscopic configuration $\Gamma_{c}$: given by
the $N$ spins on the lattice.

Example:
Two-dimensional square Ising lattice; only the four nearest spins
contribute to the energy of spin $\sigma_{i}$
($= \pm 1$). (Three dimensions: six nearest neighbors).
The total energy is

$
E=-\frac{A}{2}\sum \limits_{i}^{N} \sum \limits_{j(i)}^{4\; or \;6} \sigma_{i}
\sigma_{j(i)}
$

($A$ ... coupling constant).

$\Longrightarrow$
Magnetic polarization

$M \equiv \sum \limits_{i}^{N} \sigma_{i}$

may be determined as a function of temperature.
If a magnetic field $H$ is applied, the additional potential energy is
$E_{H}=-H \sum \limits_{i} \sigma_{i}$.

6.1.2 Technical Matters

What units?

What about boundaries?

What initial confiuration?

How to adjust density and temperature

Units:
To avoid the handling of small numbers, choose units appropriate to
the model.

Lennard-Jones:
Energy unit $E_{0}=\epsilon$; length
$l_{0}=\sigma$.
The pair energy is then

where
$u^{*} \equiv u/\epsilon$ and $r^{*} \equiv r/\sigma$.

As the third mechanical unit, choose the atomic mass
$m_{0}=1AMU= 1.6606 \cdot 10^{-27}kg$.

The time unit is now the combination
$t_{0}=\sqrt{m_{0} \sigma^{2}/\epsilon}$.

Electrical charge: best measured in multiples of the
electron charge, $q_{0}= 1.602 \cdot 10^{-19} As$.

Number density:
$\rho=N/V$ is a large number; therefore we
reduce it by a suitable standard density,
$\rho_{0}/\sigma^{3}$:
$\rho^{*} \equiv N \sigma^{3}/V$.

Temperature:
$T_{0}=\epsilon/k$

Hard spheres:
no "natural" unit of energy; therefore choose
self-consistent time unit $t_{0}=\sqrt{m_{0}d^{2}/kT}$.

For hard spheres of diameter $d_{0} =2 r_{0}$
the customary standard density is
$\rho_{0}=\sqrt{2}/d_{0}^{3}$; thus:
$\rho^{*}=N d_{0}^{3}/V\sqrt{2}$.

Periodic boundary conditions (PBC) and nearest image
convention (NIC):
To avoid "wall effects" we surround the basic cell containing
the $N$ particles by periodic images of itself.

Spin lattices:
In each coordinate direction the last ("rightmost") spin in a row
interacts also with a right neighbor, which is taken to be identical to
the first (leftmost) spin in that row, and vice versa (see Figure).

Figure:
Periodic boundary conditions on a spin lattice

Fluid: Apply the following rule:

For each particle $i$ store, instead of any coordinate
$x_{i}$, the quantity

$(x_{i}+2L) \mbox{mod} L $

(with $L$ the side length of the cell).
When a particle leaves the cell to the right,
it is automatically replaced by a particle entering from the left, etc.

Nearest Image Convention:
In computing pair vectors between two particles
$i$ and $j$ one needs only differences of coordinate values.
If this coordinate difference
$\Delta x_{ij} \equiv x_{j}-x_{i}$
is larger than $L/2$, then the particle
$j$ will not be regarded as an interaction partner of $i$;
instead, its left periodic image with coordinate $x_{j}-L$
interacts with $i$.

The NIC rule may be implemented by a sequence of IF commands or
by the more compact code line

$
\Delta x = \Delta x - L \cdot \mbox{nint}\left( \frac{\textstyle \Delta x}{\textstyle L}\right)
$

where nint($a$) denotes the rounded value of
$a$, i.e. the integer nearest to $a$. This formulation of NIC is also better suited for vectorising
compilers.

Starting configuration:
On an Ising lattice, draw $N$ spin values with equal probabilities
for $+ 1$ and $-1$.

Molecules in disordered media:
Overlaps must be avoided.
$\Longrightarrow$
Place the molecules on a lattice, then "melt" this crystal before
the actual simulation run: Thermalization.

Population number in a cubic cell with face-centered cubic arrangement:
$4m^{3}$, with $m=1,,2,\dots$. Therefore typical particle numbers in
simulations are $N = 32, 108, 256, 500$ etc.

Adjusting density and temperature:
Given $N$, the density is adjusted to a desired value by shrinking or expanding
the volume: scale all coordinates by a suitable factor.

The temperature is a constant parameter in an MC simulation.
In molecular dynamics it must be adjusted in the following manner.
Since

$
T=m \langle| v|^{2} \rangle/3k
$

(with $k=1.3804 \cdot 10^{-23} J/deg$) or, in reduced units,

$
T^{*}=m^{*} \langle | v^{*}|^{2} \rangle/3
$

we first take the average of $| v^{*}|^{2}$
over a number of MD steps to determine the actual temperature of the simulated
system. Then we scale each velocity component according to

Since $T^{*}$ is a fluctuating quantity it can be adjusted
only approximately.

EXERCISE:
(See also
here)
Reduced units:
Consider a pair of Lennard-Jones particles with
$\epsilon=1.6537\cdot 10^{-21} J$
and $\sigma=3.405 \cdot 10^{-10} m$
(typical for Argon).
Let the two molecules be situated at a distance of
$3.2 \cdot10^{-10}m$ from each other.

- Calculate the potential energy of this arrangement.

- Do the same calculation using
$\epsilon$ and $\sigma$
as units of energy and length, respectively.
These parameters then vanish from the expression for the pair energy,
and the calculation is done with quantities of order $1$.

- With the above units for energy and length, together with the
atomic mass unit, compute the metric value of the self-consistent unit
of time? Let one of the particles have a metric speed
$v=500 m/s$, typical
of the thermal velocities of atoms or small molecules. What is the value
of $v$ in self-consistent units?

PROJECT MC/MD:
(See also
here)
As a first reusable module for a simulation program, write a code to
set up a cubic box of side length $L$
inhabited by up to $N =4 m^{3}$ particles
in a face-centered cubic arrangement. Use your favourite programming language
and make the code flexible enough to allow for easy change of volume
(i.e. density). Make sure that the lengths are measured in units
of $\sigma_{LJ}$. For later reference, let us call this subroutine
STARTCONF.

Advice: It is convenient to count the lower, left and front
face of the cube as belonging to the basic cell, while the three other
faces belong to the next periodic cells.

To test for correct arragement of the particles, compute the diagnostic

$S_{0} \equiv \frac{1}{3} \sum \limits_{i}^{N} \left[
\cos (4 m \pi x_{i}/L)+\cos (4 m \pi y_{i}/L)+
\cos (4 m \pi z_{i}/L) \right]
$

which is sometimes called "melting factor". For a fcc
configuration it should be equal to $N$ (why?).

By scaling all lengths, adjust the
volume such that the reduced number density becomes
$\rho^{*}=0.6$.

Augment the subroutine STARTCONF by a procedure that assigns random
velocities to the particles, making sure that the
sum total of each velocity component is zero.

PROJECT MC/MD:
(See also
here)
The second subroutine will serve to compute the total potential
energy in the system, assuming a Lennard-Jones interaction and applying
the nearest image convention: