Computational Techniques in Theoretical Physics
Section 18
Molecular dynamics in canonical and grand canonical ensembles
Key Point
• Canonical ensmble
• Grand canonical ensmble
• Nose-Hoover method.

Canonical Ensemble:

Closed systems:
A closed system can exchange heat with its surroundings.

Chacateristics:

Physical description:

• Total number of particles fixed.
• Volume of system is fixed.
• Total energy fluctuating.
• Average energy (temperature) fixed at equilibrium.

Mathematical description:

sum P(x) = 1

< E > = sum E(x) P(x)
= sum H(x) P(x)

where H(x) is the Hamiltonian of the system.

Canonical ensemble is a statistical description of a closed system in equilibrium:

Entropy is maximum with the condition that the number of particles and average energy is constant:

d  int dx [ a P(x) + b H(x) P(x) - kB P(x) ln cN P(x) ] = 0

From this equation, we can obtain:

P(x) = exp{-H(x)/kBT} / cN Z

Z = int dx exp{ -H(x)/kBT }  / cN

Z is the partition function. All physical quantities can be derived from it:

Free energy:

F = - kBT ln Z
Entropy:

S = - dF / dT  (with V and N fixed).
Internal energy:

< E > = F + S T
Pressure:

P = dF / dV  (with T and N fixed)
Heat Capacity:

CV = d< E > / dT = T dS / dT
(with V and N fixed)

Grand canonical ensemble:

Open systems:

An open system can exchange both heat and matter with its surroundings

Characteristics:

Physical description:

• Total energy and particle number fluctuating.
• Average energy (temperature), pressure and everage particle number fixed at equilibrium.

Mathematical description:

< N > = sum N P(x

< E > = sum E(x) P(x)
= sum H(x) P(x)

where H(x) is the Hamiltonian of the system.

Grand canonical ensemble is a statistical description of a closed system in equilibrium:

Entropy is maximum with the condition of fixed pressure, average number of particles and average energy:

d  int dx [ a P(x) + b H(x) P(x)  + c N P(x)
- kB P(x) ln cN P(x) ] = 0

From this equation, we can obtain:

P(x) = exp{-[H(x) - mu N]/kBT} / cNZ'

Z' = int dx exp{ -[H(x) - mu N]/kBT }  / cN

Z' is the grand partition function, and mu is the chemical potential. All physical quantities can be derived from Z':

Free energy:

G = - kBT ln Z'
Entropy:

S = - dG / dT  (with P and mu fixed).
Internal energy:

< E > = F + S T
Volume:

P =  -dG / dP (with T and N fixed)
Heat Capacity:

Cp =  - T d2G /dT2
(with P and N fixed)

Nose-Hoover Molecular Dynamics for canonical ensembles

General principles:
In statistical mechanics, we deal most conveniently with canonical ensemble, which has a fixed number of particles N, fixed volume V, and fixed temperature T.

It is also possible to simulate such systems by molecular dynamics. Such algorithms introduce addition forces to the Newton's equations. The forces are introduced purely for the sake of producing correct ensemble, they are not real forces. Thus these systems are somewhat artificial and less fundamental.

There are a number of methods which realizes canonical system, one of them is to rescale velocities at each time step. Here we'll just introduce the elegant method of Nose-Hoover

Algorithm for isothermal molecular dynamics (canonical ensemble)
The new dynamical equation proposed by Nose and Hoover in 1984 is the following:

dqi/dt = pi/m

dpi/dt = F(q) - z pi

dz/dt = sum i [ pi2/mkBT - 1 ] / N t'd

The notations p and q stand for the vector (p1, ..., pN) and (q1, ..., qN). The frictional coefficient z is a scalar variable introduced to realize the canonical ensemble. The parameter t' is an arbitrary constant, d is dimension, and N is the number of particles.

Note that if we set z = 0, we reduce the equations to Hamilton's equations of motion. The temperature T appears in the equation, which sets the equilibrium temperature of the system.

It can be shown that the canonical distribution augmented by a Gaussian distribution in the frictional coefficient:

P(p, q, z) ~ exp[ - H(p,q)/kBT - 1/2 z2 t'2 ]

is a steady state distribution of the Nose-Hoover dynamics. This guarantees that the time average over the trajectories of the system is equal to the average over the canonical distribution.

We can rewrite the Nose-Hoover dynamics in a more familiar Newtonian form:

m dvi/dt = Fi - z m vi

dz/dt = [K/K0 - 1] / t'2

K = sum i 1/2 mvi2

K0 = 1/2 NkBTd

Here K is total kinetic energy, and K0 is average kinetic energy. This set of equations can be solved with an algorithm similar to that of Verlet as

ri(t+h) = 2ri(t) - ri(t-h) + h2 Fi(t) / m
- z(t)/2h [ri(t+h) - ri(t-h)]

z(t+h) = z(t) + h/t'2 [ K(t+h/2)/K0 - 1 ]

K(t+h/2) = m/2h2 sum i ( ri(t+h) - ri(t) )2

The difference scheme here is also exactly time-reversible:

If you start from ri(-h), ri(0), and z(0), the system develops to ri((n-1)h), ri(nh), and z(nh), after n steps with positive h.

Now if you run the system backwards, namely, starting from ri(nh) as the previous position, ri((n-1)h) as the current position and z((n-1)h), applying the recursion relation after n steps with h negative, you get back exactly the starting values.

This cannot be true if a fourth-order Runge-Kutta method were used. This time-reversal symmetry is the key ingredient for long-time stability of the algorithm. In the usual microcanonical molecular dynamics, total energy is a constant. This is no longer so in the canonical dynamics. Nevertheless, there exists an analogous quantity which is a constant of the motion:

C = K + V + K0(t'z)2 + 2K0 int z(t'') dt''

This equation can be used to check the stability of integration schemes.

Grand canonical ensemble
(Isothermal-isobaric)  molecular dynamics

An isothermal-isobaric ensemble is a one with fixed number of particles N, temperature T, and pressure P0.

Most experimental conditions are of this type.

In an isothermal-isobaric ensemble volume is a variable. The canonical distribution is replaced by:

P(p,q,V) ~ exp{ -[ H(p,q) + P0 V ]/ kBT }

where H=K+V is total energy, P0 is the constant pressure. In this ensemble, the logarithm of the (isothermal-isobaric) grand partition function is -G/kB T, G is Gibbs free energy.

The generalization of the Newton's equations of motion is

d2 xi /dt2 = Fi/mL - (z+2e) dxi/dt

de/dt = V(P - P0)/kBTaV

e = 1/L dL/dt

dz/dt = [K/K0 - 1] / aT

K0 = 1/2 N kB T d

where xi = ri/L is the coordinate vector normalized by the length of the box; the volume of the box is V = Ld; aV and aT are two parameters; P and K are the instantaneous pressure and kinetic energy, defined by

PV = 2/d sum i 1/2 m vi2 + 1/d sum i<j rij * Fij

K = sum i 1/2 m vi2

Homework