of a suitable parent model. For simplicity, the model should consist of

sites (atoms and possibly virtual sites) with partial charges only, without

higher multipoles. In order to represent real charge distributions faithfully,

smeared charge distributions should be allowed. Polarizability should be

realized by moving charge distributions attached to atoms, with special care

to correctly represent the induced charge distributions in close proximity

of external charges. The reason is not only simplicity, but it is also less

likely that induced dipoles or ¬‚uctuating charges will eventually prove to

be adequate, as they do not include polarization induced by exchange and

dispersion forces (see the discussion on page 178). Intermolecular repulsion

and dispersion interactions should be largely centered on the moving shells,

possibly re¬ned with short-range atom-based corrections. The model should

34 The three-body dispersion is a result of third-order perturbation and is inversely proportional

222

to r12 r23 r13 (Axilrod and Teller, 1948).

188 Molecular dynamics

be parameterized by ¬tting to high-level ab initio calculations, with possible

empirical ¬ne tuning.

6.3.10 QM/MM approaches

Force ¬elds of the type described above are not suitable for systems of

molecules that undergo chemical transformations. Con¬gurations along re-

action paths di¬er so much from the covalently bonded stable states that

they cannot be described by simple modi¬cations of the force ¬elds intended

for stable molecules. For such con¬gurations quantum calculations are re-

quired. However, in complex systems in which reactions take place, as in

enzymes and other catalysts, the part that cannot be described by force

¬elds is generally quite limited and it would be an overkill (if at all possible;

however, see the ab initio MD treated in Section 6.3.1) to treat the whole

system by quantum-chemical methods.

For such systems it is possible to use a hybrid method, combining quan-

tum calculations to obtain energies and forces in a limited fragment, em-

bedded in a larger system for which energies and forces are obtained from

the much simpler force ¬eld descriptions (see Fig. 6.11). These methods are

indicated as QM/MM (quantum mechanics/ molecular mechanics) meth-

ods. The principle was pioneered by Warshel and Levitt (1976), and has

been widely applied since the mid-1980s.35 Most simulation packages al-

low coupling with one or more quantum chemistry programs, which include

semi-empirical, as well as density functional and ab initio treatments of a

fragment. QM/MM can be used for optimizations, for example of transition

states and reaction paths, but also for dynamic trajectories. The forces and

energies can be generated every time step, producing a very costly dynamics

trajectory, or the QM calculations can be done at selected instances from

approximate trajectories. Note that the quantum part is always solved in

the Born“Oppenheimer approximation: QM/MM methods do not produce

quantum dynamics.

The coupling between the QM and the MM part must be carefully mod-

eled. In general, the total energy consists of three contributions, arising from

the QM part, the MM part and the interaction between the two. When the

35 For a survey see Gao and Thompson (1998). Chandra Singh and Kollman (1986) introduced the

cap atom, which is usually hydrogen and which is kept free to readjust. The coupling between

the QM and MM part is described by Field et al. (1990). A layered approach, allowing

layers of increasing quantum accuracy, has been introduced by Svensson et al. (1996). In the

“AddRemove” scheme of Swart (2003), the capping atom is ¬rst added on the line pointing to

the MM link atom, included in the QM, and then its interaction is removed from the sum of

QM and MM energies. Zhang and Yang (1999) use a pseudobond to a one-free-valence atom

with an e¬ective core potential instead of a cap atom.

6.4 Solving the equations of motion 189

MM

XQ XM

cap

QM

Figure 6.11 A quantum mechanical fraction embedded in a molecular mechanics

environment. Covalent bonds that cross the QM-MM boundary are replaced by a

cap atom in the QM part.

QM part consists of a separate molecule, the intermolecular interactions

couple the two systems. When the QM part consists of a molecular frag-

ment, the covalent bond(s) between the fragment and the MM environment

must be replaced by some other construct, most often a “cap” atom or

pseudo-atom.

In principle there is a discrepancy between the levels of treatment of the

QM part, which includes induction by the electric ¬elds of the environments,

and MM part in the case of non-polarizable MM force ¬elds. The reader is

referred to Bakowies and Thiel (1996) for a detailed evaluation of the various

ways the QM and MM systems can interact.

6.4 Solving the equations of motion

Given a conservative force ¬eld, we know the forces acting on atoms and on

virtual interaction sites. The forces on virtual sites are ¬rst redistributed

over the atoms from which the virtual sites are derived, so we end up with

(cartesian) forces on mass points. The description may contain constraints

of bond lengths and possibly bond angles.

While it is possible to write the equations of motion in internal coordi-

nates (see Chapter 15), these tend to become quite complicated and the

general recommendation for atomic systems is to stick to cartesian coordi-

nates, even in the presence of constraints. Modern constraint methods are

robust and e¬cient (see Section 15.8 on page 417). However, for completely

rigid molecules the use of quaternions may be considered (see page 413).

Consider a system of N particles with mass mi , coordinates r i and a

de¬ned recipe to compute the total potential energy Epot = V (r) and the

190 Molecular dynamics

forces F i (r) = ’∇i V (r), given the set of all coordinates. Let coordinates

and velocities v i be known at time t. We assume that there are no con-

straints. The equations of motion are simply Newton™s equations (which are

Hamiltonian as well):

™

ri = vi ,

™

v i = F i /mi . (6.79)

The total energy Etot = K + V will be conserved:

dEtot d 1 dV (r) ‚V

mi v i · v i + · vi

mi v 2 + ™

= =

i

dt dt 2 dt ‚r i

i i i

‚V

vi · F i + · vi

= = 0. (6.80)

‚r i

i

Properly solving these equations of motion will produce a microcanonical

or N, V, E ensemble. In practice there will be errors that cause deviations

from the ideal behavior: the ¬nite time step will cause integration errors

and the total energy will not be exactly conserved; errors in forces (e.g., due

to truncation) will produce pseudorandom disturbances that cause energies

to drift. Since the temperature is determined by the equipartition theorem

saying that K = 3 N kB T , the temperature may drift even when equilibrium

2

has been attained. Therefore there are always modi¬cations to the pure

Newtonian equations of motion needed to generate long stable trajectories.

The equations of motion are solved in time steps ”t. Three important

considerations in¬‚uence the choice of algorithm:

(i) Time reversibility, inherent in the Newtonian equations of motion,

should be conserved.

(ii) The generated trajectories should conserve volume in phase space,

and in fact also wedge products (area) dq§dp in general, i.e., the algo-

rithm should be symplectic (see Chapter 17, page 495). This is impor-

tant to conserve equilibrium distributions in phase space, because de-

viation from symplectic behavior will produce time-dependent weight

factors in phase space. This importance of the symplectic prop-

erty has been emphasized in the 1990s (Leimkuhler and Reich, 1994;

Leimkuhler and Skeel, 1994) and is now widely recognized.