time step must be taken considerably smaller than in ordinary molecular

dynamics.

For algorithmic details and applications of ab initio molecular dynamics

the reader is referred to Marx and Hutter (2000).

6.3.2 Simple molecular force ¬elds

The simplest force ¬elds, useful for large molecular systems, but not aiming

at detailed reproduction of vibrational spectroscopic properties, contain the

following elements:

• Atoms are the mass points that move in the force ¬eld. In united-atom

approaches some hydrogen atoms are “incorporated” into the atom to

which they are bound. In practice this is used for hydrogen atoms bound

to aliphatic carbon atoms; the resulting CH2 or CH3 groups are “united

atoms,” acting as a single moving mass.

• Atoms (or united atoms) are also the source points for the di¬erent terms

in the force ¬eld description. This means that the various contributions

to the forces are expressed as functions of the atomic positions.

• There are two types of interactions: bonded interactions between dedi-

cated groups of atoms, and non-bonded interactions between atoms, based

on their (changing) distance. These two types are computationally di¬er-

ent: bonded interactions concern atoms that are read from a ¬xed list, but

atoms involved in non-bonded interactions ¬‚uctuate and must be updated

regularly. Non-bonded interactions are assumed to be pairwise additive.

Bonded interactions are of the following types:

(i) A covalent bond between two atoms is described by a harmonic po-

154 Molecular dynamics

tential of the form:4

Vb (r i , r j ) = 1 kb (r ’ b)2 , (6.16)

2

where

r = |r i ’ r j |, (6.17)

and kb and b are parameters which di¬er for each bond type.

The harmonic potential may be replaced by the more realistic Morse

potential:

Vmorse (r i , r j ) = Dij [1 ’ exp(’βij (rij ’ b))]2 . (6.18)

Other forms contain harmonic plus cubic terms.

(ii) A covalent bond angle is described by a harmonic angular potential

of the form

Va (r i , r j , r k ) = 1 kθ (θ ’ θ0 )2 , (6.19)

2

where

r ij · r kj

θ = arccos , (6.20)

rij rjk

or by the simpler form

Va = 1 k (cos θ ’ cos θ0 )2 . (6.21)

2

(iii) Dihedral angles φ are de¬ned by the positions of four atoms i, j, k, l

as the angle between the normals n and m to the two planes i, j, k

and j, k, l:

n·m

φ = arccos (6.22)

nm

where

n = r ij — r kj

m = r jk — r lk . (6.23)

The dihedral potential is given by a periodic function

Vd (φ) = kφ (1 + cos(nφ ’ φ0 )). (6.24)

This makes all minima equal (e.g., the trans and the two gauche

states for a threefold periodic dihedral, as between two sp3 carbon

atoms). The actual di¬erence between the minima is caused by the

The GROMOS force ¬eld uses a quartic potential of the form V = (kb b’2 /8)(r 2 ’ b2 )2 , which

4

for small deviations is virtually equivalent to (6.16), but computationally much faster because

it avoids computation of a square root.

6.3 Force ¬eld descriptions 155

introduction of an extra interaction between atoms i and l, called the

1-4 interaction.

Instead of using a 1-4 interaction, one may also use a set of peri-

odic functions with di¬erent periodicity, or a set of powers of cosine

functions, as in the Ryckaert“Bellemans potential

5

Cn cosn φ.

VRB (φ) = (6.25)

n=0

(iv) In order to keep planar groups (as aromatic rings) planar and prevent

molecules from ¬‚ipping over to their mirror images, improper dihe-

drals are de¬ned, based on four atoms i, j, k, l and given a harmonic

restraining potential:

Vimproper (ξ) = 1 kξ (ξ ’ ξ0 )2 . (6.26)

2

Bonded interactions, if they are so sti¬ that they represent high-frequency

vibrations with frequency ν kB T /h, can be replaced by constraints. In

practice this can only be done for bond length constraints, and in some cases

for bond-angle constraints as well (van Gunsteren and Berendsen, 1977).

The implementation of constraints is described in Chapter 15, Section 15.8

on page 417.

Non-bonded interactions are pair-additive, and a function of the distance

rij = r between the two particles of each pair. Pairs that are already involved

in bonded interactions are excluded from the non-bonded interaction; this

concerns 1-2 and 1-3 interactions along a covalently-bonded chain. The 1-4

interactions are either excluded, or used in modi¬ed form, depending on

the dihedral functions that are used. Non-bonded interactions are usually

considered within a given cut-o¬ radius, unless they are computed as full

lattice sums over a periodic lattice (only in the case of periodic boundary

conditions). They are of the following types:

(i) Lennard“Jones interactions describe the short-range repulsion and

the longer-range dispersion interactions as

C12 C6

’ 6,

vLJ (r) = (6.27)

r12 r

which can be alternatively expressed as

σ σ

12 6

’

vLJ (r) = 4µ . (6.28)

r r

The treatment of the long-range part of the dispersion will be sepa-

rately considered below in Section 6.3.4. The r’12 repulsion term is

156 Molecular dynamics

of rather arbitrary shape, and can be replaced by the more realistic

exponential form

vrep = A exp(’Br), (6.29)

which, combined with the dispersion term, is usually referred to as

the Buckingham potential.

(ii) Coulomb interactions between charges or partial charges on atoms:

qi qj

VC (r) = fel (6.30)

µr r

Here fel = (4πµ0 )’1 and µr is a relative dielectric constant, usually

taken equal to 1, but in some force ¬elds taken to be a function of r

itself (e.g., equal to r measured in ˚) to mimic the e¬ect of dielectric

A

screening. The latter form must be considered as ad hoc without

physical justi¬cation. Special care is needed for the treatment of the

long-range Coulomb interaction, which is separately described below

in Section 6.3.5. The partial charges are often derived from empir-

ical dipole and quadrupole moments of (small) molecules, or from

quantum calculations. A simple Mulliken analysis of atomic charges

resulting from the occupation of atomic orbitals does not su¬ce; the

best partial charges are potential-derived charges: those that repro-

duce the electric potential in the environment of the molecule, with

the potential determined from a high-level ab initio quantum cal-

culation. Once the potential has been determined on a grid, and