face is required to enable the simulation of motion on that surface. This

is the ¬rst task: design a suitable force ¬eld from which the forces on each

atom, as well as the total energy, can be e¬ciently computed, given the

positions of each atom.1 Section 6.3 describes the principles behind force

¬elds, and emphasizes the di¬culties and insu¬ciencies of simple force ¬eld

descriptions. But before considering force ¬elds, we must de¬ne the system

with its boundary conditions (Section 6.2). The way the interactions over

the boundary of the simulated systems are treated is in fact part of the total

potential energy description.

The force ¬eld descriptions take the covalent structure of molecules into

account. They are not valid when chemical reactions may take place, chang-

ing the covalent structure, or the redox state, or even the protonation state

of a molecule. In such cases at least the reactive part of the molecular sys-

tem should be treated di¬erently, e.g., by quantum-chemical methods. The

1 The name “atom” is used interchangeably with “nucleus;” as the electronic motion is not

separately considered, the di¬erence is immaterial.

139

140 Molecular dynamics

“ab initio molecular dynamics” method (Section 6.3.1) solves the electronic

and nuclear equation simultaneously. Methods that use quantum-chemical

approaches for a subsystem, embedded in a larger system described by a

“standard” force ¬eld, are called QM/MM methods, and are the subject of

Section 6.3.10.

The methods to solve the equations of motion are described in Section 6.4.

We focus on methods that retain cartesian coordinates for the atomic po-

sitions, as those are by far the easiest to implement, and generally also

the fastest to execute. Internal constraints are then implemented by spe-

cial algorithms. A review of the principles of classical mechanics and a more

complete description of the special techniques for rigid-body and constrained

dynamics is given in Chapter 15.

In Section 6.5 coupling of the simulated system to external temperature

and pressure baths is described. This includes extended system methods

that extend the system with extra degrees of freedom allowing the control

of certain variables within the system. Such controls can also be used to

invoke non-equilibrium molecular dynamics, driving the system away from

thermal equilibrium in a speci¬ed way, and allowing the direct study of

transport properties.

Some straightforward applications are discussed in Section 6.7. The com-

putation of macroscopic quantities that depend on the extent of accessible

phase space rather than on microscopic averages, viz. entropy and free en-

ergy, is left for the next chapter.

6.2 Boundary conditions of the system

The ¬rst item to consider is the overall shape of the system and the boundary

conditions that are applied. Long-range interactions, notably of electrostatic

origin, are dependent on such conditions. Only isolated molecules in the

dilute gas phase are an exception, but these have a very limited interest in

practice. In general we are concerned with systems consisting of molecules

in the condensed phase, with a size much larger than the system we can

a¬ord to simulate. Thus the simulated system interacts over its boundaries

with an environment that is (very) di¬erent from vacuum. The important

consideration is that the environment must respond to changes in the system;

such a response is not only static, involving an average interaction energy,

but is also dynamic, involving time-dependent forces reacting to changes in

the system.

6.2 Boundary conditions of the system 141

y

T

¡ ¡ ¡ ¡

¡ ¡ ¡ ¡

e e e

¡ ¡ ¡ ¡

¡ ¡ ¡ ¡

¡ ¡ ¡ ¡

jb ¡ b b¡

¡ ¨¡

¨¡

¡ ¡ ¡

d

d ¨

b ¡ du ¨

¡ ¡

! ¡

ey e

¨i

¡ ¡

¡ e

¡ ¡ ¡

!

·b¡ ¡

¡ ¡ ¡ ¡

e

¡ ¡

¡ ¡ ¡ ¡

er

b¡ b¡

¡ ¡

¡ j¡

¡E

¡ ¡ Ex

E¡ ¡

¡

ξa e x ¡a e

e

¡ ¡ ¡

¡ ¡ ¡ ¡

¡ ¡ ¡ ¡

¡ b¡ b¡ b¡

¡ ¡ ¡ ¡

¡ ¡ ¡ ¡

Figure 6.1 Periodic box (2D, for simplicity) determined by base vectors a and b.

Particles i and j in the unit cell are given with their images; the nearest image to

i is not j, but j in the NW neighboring cell.

6.2.1 Periodic boundary conditions

The simplest, and most often applied, conditions are periodic boundary con-

ditions (Fig. 6.1). The system is exactly replicated in three dimensions, thus

providing a periodic lattice consisting of unit cells. Each unit cell can have

an arbitrary triclinic shape, de¬ned by three basis vectors a, b, c, with arbi-

trary angles ±, β, γ (± is the angle between b and c, etc.) between the basis

vectors. If there is only one angle di¬erent from 90—¦ , the cell is monoclinic; if

all angles are 90—¦ , the cell is rectangular, and if in addition all basis vectors

have equal length, the cell is cubic. Note that the unit cell of a periodic

system is not uniquely de¬ned: the origin can be placed arbitrarily, the unit

cell can be arbitrarily rotated, and to each base vector a linear combination

of integer multiples of other base vectors may be added. The volume of the

unit cell does not change with any of these operations. For example, if c is

changed into c + na a + nb b (na , nb = 0, ±1, ±2, . . .), the volume (see below)

V = (a — b) · c does not change since a and b are both perpendicular to

a — b. One can always choose a unit cell for which the projection of b on a

is smaller than 1 a, and make sure that γ ≥ 60—¦ .

2

Several other cell shapes have been devised, such as the truncated octahe-

dron and the rhombic dodecahedron, both of which pack in three dimensions

and have a smaller volume than a cubic cell for the same minimal distance

142 Molecular dynamics

(a) (b)

(c)