molecules, and valid for a wide range of environments and conditions. It

is often the non-additivity of constituent terms, and the omission of im-

portant contributions, that renders the terms non-transferable. Since most

150 Molecular dynamics

force ¬elds contain parameters adjusted to empirical observations, an error

or omission in one term is compensated by changes in other terms, which

are then not accurate when used for other con¬gurations, environments or

conditions than those for which the parameters were adjusted.

Ideally, ab initio quantum calculations should provide a proper potential

energy surface for molecules and proper descriptions for the interaction be-

tween molecules. Density functional theory (DFT) has “ for small systems

“ advanced to the point that it is feasible to compute the energy and forces

by DFT at every time step of the molecular motion and thus evolve the sys-

tem dynamically. The “ab initio molecular dynamics” method of Car and

Parrinello (1985), described in Section 6.3.1, employs a clever method to

solve the electronic and nuclear equations simultaneously. Other quantum

approximations that scale more linearly with the number of particles, such

as the divide-and-conquer and the tight-binding approximations (see Section

4.9) are candidates for on-the-¬‚y quantum calculations of energies and force

during dynamic evolution. In general, however, for large systems such direct

methods are not e¬cient enough and simpler descriptions of force ¬elds are

required.

There are several reasons that quantum calculations on isolated molecules

do not su¬ce to produce reliable force ¬elds, and empirical adjustments are

still necessary:

• For condensed systems, interaction with the (in¬nite) environment must

be properly accounted for,

• The force ¬eld description must necessarily be simpli¬ed, if possible to

additive local terms, and this simpli¬cation involves approximation of the

full quantum potential energy surface.

• Even high-quality ab initio calculations are not accurate enough to pro-

duce overall accuracies better than kB T , as required to yield accurate

thermodynamic properties. Note that kB T = 2.5 kJ/mol for room tem-

perature, while an error of 6 kJ/mol in the free energy di¬erence of two

states corresponds to an error of a factor of 10 in concentrations of com-

ponents participating in an equilibrium.

• Small e¬ects that are not incorporated in the Born“Oppenheimer quan-

tum mechanics, such as nuclear quantum e¬ects, must be accounted for.

The choice must be made whether or not such corrections are applied to

the result of calculations, before parameter adjustments are made. If they

are not applied afterwards, the quantum e¬ects are mimicked by adjust-

ments in the force ¬eld contributions. In fact, this consideration does not

6.3 Force ¬eld descriptions 151

only apply to quantum corrections, but to all e¬ects that are not explicitly

accounted for in the force ¬eld description.

Still, it is through the study of quantum chemistry that insight is obtained

in the additivity of constituent terms and the shape of the potential terms

can be determined. Final empirical adjustments of parameters can then

optimize the force ¬eld.

6.3.1 Ab-Initio molecular dynamics

Considering nuclei as classical point particles and electrons as providing a

force ¬eld for the nuclear motion in the Born“Oppenheimer approximation,

one may try to solve the energies and forces for a given nuclear con¬guration

by quantum-chemical methods. The nuclear motion may then be advanced

in time steps with one of the standard molecular dynamics algorithms. For

e¬ciency reasons it is mandatory to employ the fact that nuclear con¬gu-

rations at successive time steps are very similar and therefore the solutions

for the electronic equations are similar as well.

In a seminal article, Car and Parrinello (1985) described the simultaneous

solution of the nuclear equations of motion and the evolution of the wave

function in a density-functional description. The electron density n(r) is

written in terms of occupied single-particle orthonormal Kohn“Sham orbitals

(K“S orbitals, see Section 4.7 for a more detailed description);

|ψi (r)|2 ,

n(r) = (6.10)

i

where each ψl (r) is a linear combination of well-chosen basis functions. Car

and Parrinello chose as basis functions a set of plane waves exp(ik · r) com-

patible with the

periodic boundary conditions. Thus every K“S orbital is a vector in re-

ciprocal space. A point of the Born“Oppenheimer potential energy surface

is given by the minimum with respect to the K“S orbitals of the energy

functional (4.59):

2

—

E=’ dr ψi (r)∇2 ψi (r) + U [n(r); R]. (6.11)

2me

i

Here the ¬rst term is the kinetic energy of the electrons and the second term

is a density functional, containing both the electron-nuclear and electron-

electron interaction. The latter consists of electronic Coulomb interactions,

and exchange and correlation contributions. The K“S wave functions ψi

152 Molecular dynamics

(i.e., the plane wave coe¬cients that describe each wave function) must be

varied to minimize E while preserving the orthonormality conditions

—

dr ψi (r)ψj (r) = δij . (6.12)

The nuclear coordinates are constant parameters in this procedure. Once

the minimum has been obtained, the forces on the nuclei follow from the

gradient of E with respect to the nuclear coordinates. With these forces the

nuclear dynamics can be advanced to the next time step.

The particular innovation introduced by Car and Parrinello lies in the

method they use to solve the minimization problem. They consider an ex-

tended dynamical system, consisting of the nuclei and the K“S wave funct-

ions. The wave functions are given a ¬ctitious mass μ and a Lagrangian

(see (15.2)) is constructed:

1 ™2

L= dr |ψi |2 + MI RI ’ E(ψ, R). (6.13)

2

I

This Lagrangian, together with the constraints (6.12), generate the following

equations of motion (see Section 15.8):

‚E

¨

μψi (r, t) = ’ + Λik ψk , (6.14)

‚ψi

k

¨

MI RI = ’∇RI E , (6.15)

where I numbers the nuclei, MI are the nuclear masses and RI the nu-

clear coordinates. The Λik are Lagrange multipliers introduced in order to

satisfy the constraints (6.12). The equations are integrated with a suitable

constraint algorithm (e.g., Shake, see Section 15.8).

When the ¬ctitious masses of the wave functions are chosen small enough,

the wave function dynamics is much faster than the nuclear dynamics and

the two types of motion are virtually uncoupled. Reducing the “velocities”

and consequently the kinetic energy or the “temperature” of the wave funct-

ions will cause the wave functions to move close to the B“O minimum energy

surface (in the limit of zero temperature the exact minimum is reached). In

fact, this dynamic cooling, reminiscent of the “simulated annealing” method

of Kirkpatrick et al. (1983), is an e¬ective multidimensional minimization

method.

In practice the temperature of the electronic degrees of freedom can be

kept low enough for the system to remain close to the B“O energy surface,

even when the temperature of the nuclei is considerably higher. Because

both systems are only weakly coupled, heat exchange between them is very

6.3 Force ¬eld descriptions 153

weak. Both systems can be coupled to separate thermostats (see Section

6.5) to stabilize their individual temperatures. There is a trade-o¬ between

computational e¬ciency and accuracy: when the “wave function mass” is

small, wave functions and nuclei are e¬ectively uncoupled and the system

can remain accurately on its B“O surface, but the electronic motions become

fast and a small time step must be used. For larger masses the motions of

the electronic and nuclear degrees of freedom will start to overlap and the

B“O surface is not accurately followed, but a larger time step can be taken.

In any case the Car“Parrinello method is time-consuming, both because a