describing a corresponding lattice. Each unit cell contains parts of four hexagons

and each hexagon is spread over four cells. (b) The 3D rhombic octahedron with 12

faces. (c) A stereo pair depicting a packed array of rhombic dodecahedra and the

triclinic unit cell of the corresponding lattice (parallel view: left picture for left eye).

Figures reproduced from Wassenaar (2006) by permission of Tsjerk Wassenaar,

Groningen.

to a neighboring cell. However, they can all be de¬ned in a triclinic periodic

lattice, and there is no need to invoke such special unit cells (Bekker, 1997).

Figure 6.2 shows how a periodic, but more near-spherical shape (as the

hexagon in two dimensions or the rhombic octahedron in three dimensions)

packs in a monoclinic or triclinic box.

Coordinates of particles can always be expressed in a cartesian coordinate

system x, y, z, but expression in relative coordinates ξ, ·, ζ in the periodic

unit cell can be convenient for some purposes, e.g., for Fourier transforms

6.2 Boundary conditions of the system 143

(see page 331). These latter are the contravariant components2 of the posi-

tion vector in the oblique coordinate system with unit vectors a, b, c.

r = xi + yj + zk (6.1)

r = ξa + ·b + ζc. (6.2)

Here i, j, k are cartesian unit vectors, and a, b, c are the base vectors of

the triclinic unit cell. The origins of both coordinate systems coincide. For

transformation purposes the matrix T made up of the cartesian components

of the unit cell base vectors is useful:3

⎛ ⎞

ax bx cx

T = ⎝ ay by cy ⎠ . (6.3)

az bz cz

This matrix can be inverted when the three base vectors are linearly inde-

pendent, i.e., when the volume of the cell is not zero. It is easily veri¬ed

that

r = Tρ (6.4)

ρ = T’1 r, (6.5)

where r and ρ denote the column matrices (x, y, z)T and (ξ, ·, ζ)T , respec-

tively. These equations represent the transformations between cartesian and

oblique contravariant components of a vector (see also page 331). Another

characteristic of the oblique coordinate system is the metric tensor g, which

de¬nes the length of a vector in terms of its contravariant components:

(dr)2 = (dx)2 + (dy)2 + (dz)2 = gij dξi dξj , (6.6)

i,j

where ξi stands for ξ, ·, ζ. The metric tensor is given by

g = TT T. (6.7)

Finally, the volume V of the triclinic cell is given by the determinant of the

transformation vector:

V = (a — b) · c = det T. (6.8)

This determinant is also the Jacobian of the transformation, signifying the

2 One may also de¬ne covariant components, which are given by the projections of the vector

onto the three unit vectors, but we shall not need those in this context. We also do not

follow the convention to write contra(co)variant components as super(sub)scripts. In cartesian

coordinates there is no di¬erence between contra- and covariant components.

In many texts this transformation matrix is denoted with the symbol h.

3

144 Molecular dynamics

Table 6.1 Unit cell de¬nitions and volumes for the cubic box, the rhombic

dodecahedron and the truncated octahedron (image distance d)

Box Box Box vectors Box angles

±βγ

type volume abc

d 0 0

90—¦ 90—¦ 90—¦

3 0 d 0

cubic d

0 0 d

rhombic ⎛ ⎞

d 0 d/2

⎝0 ⎠ 60—¦ 60—¦ 90—¦

0.707d3 d √d/2

dodecahedron

0 0 2d/2

truncated ⎛ ⎞

’d/3

d d/3

√ √

⎝0 ⎠ 70.5—¦ 70.5—¦ 70.5—¦

0.770d3 2 2d/3 √2d/3

octahedron

0 0 6d/3

modi¬cation of the volume element dx dy dz:

‚(x, y, z)

dx dy dz = dξ d· dζ = det T dξ d· dζ. (6.9)

‚(ξ, ·, ζ)

In practice, it is always easier to express forces and energies in cartesian

rather than oblique coordinates. Oblique coordinates are useful for manip-

ulation with images.

In many applications, notably proteins in solvent, the optimal unit cell has

a minimal volume under the condition that there is a prescribed minimum

distance between any atom of the protein and any atom of any neighboring

image. This condition assures that the interaction between images (which

is an artefact of the periodicity) is small, while the minimal volume min-

imizes the computational time spent on the less interesting solvent. For

approximately spherical molecules the rhombic dodecahedron is the best,

and the truncated octahedron the second-best choice (see Table 6.1). The

easy, and therefore often used, but suboptimal choice for arbitrary shapes is

a properly chosen rectangular box. As Bekker et al. (2004) have shown, it

is also possible to automatically construct an optimal molecular-shaped box,

that minimizes the volume while the distances between atoms of images re-

main larger than a speci¬ed value. An example is given in Fig. 6.3. For the

particular protein molecule shown, and with the same minimum distance

between atoms of images, the volumes of the cube, truncated octahedron,

rhombic dodecahedron and molecular-shaped boxes were respectively 817,

6.2 Boundary conditions of the system 145

629, 578 and 119 nm3 ; even more dramatically the numbers of solvent (wa-

ter) molecules were respectively 26 505, 20 319, 18 610 and 3 474, the latter

reducing the computational time for MD simulation from 9h41 for the cubic

box to 1h25 for the molecular-shaped box. Using a molecular-shaped box

in MD simulations, it is mandatory to prevent overall rotation of the cen-

tral molecule in the box in order not to destroy the advantages of the box

shape when time proceeds. An algorithm to constrain overall rotation and

translation due to Amadei et al. (2000) can be easily implemented and it

has been shown that rotational constraints do not to in¬‚uence the internal

dynamics of the macromolecule (Wassenaar and Mark, 2006).

Artifacts of periodic boundary conditions

Periodic boundary conditions avoid the perturbing in¬‚uence of an arti¬cial

boundary like a vacuum or a re¬‚ecting wall, but add the additional artefact

of periodicity. Only if one wishes to study a crystal, periodic boundary con-

ditions are natural, but even then they suppress all motions in neighboring

unit cells that are di¬erent from the motions in the central cell. In fact,

periodic boundary conditions suppress all phenomena that “ in reciprocal

space “ concern k-vectors that do not ¬t multiples of the reciprocal basis

vectors in (see Section 12.9 on page 331 for a description of the reciprocal

lattice) 2π/a— , 2π/b— , 2π/c— .

Periodic boundary conditions imply that potential functions are also peri-

odic. For example, the Coulomb energy between two charges q1 (r 1 ), q2 (r 2 )

contains the interaction with all images, including the self-images of each

particle. Algorithms to implement lattice sums for Coulomb interactions are

described in Chapter 13. This leads to the requirement of overall electrical

neutrality of the unit cell to avoid diverging electrostatic energies, to cor-

rection terms for the self-energy and to considerable artifacts if the system