the present confusion. Let us attempt to clarify the principles involved and

make a series of rational choices.

We list a few considerations, not necessarily in order of importance:

• Simplicity: as such very simple models as SPC/E are already so successful

within their realm of applicability, it seems an overkill to devise very

sophisticated models with many multipoles and polarizabilities, and many

intermolecular interaction terms. However excellent such models may be,

they will not become adopted by the simulation community.

• Robustness, meaning validity in varying environments. Correct phase be-

havior can only be expected if the same model represents all phases well,

over a wide range of temperature and pressure. The same model should

also be valid in very di¬erent molecular environments.

• Transferability, i.e., the principle of constructing the model, and as much

as possible also the parameters, should be applicable to other molecules

with similar atoms.

• Accuracy, meaning the precision reached in reproducing experimental

properties. The nature of such properties can be thermodynamic (e.g.,

phase boundaries, free energies), static structural (e.g., radial distribution

functions), static response (e.g., dielectric constant, viscosity, di¬usion

constant) or dynamic response (e.g., dielectric dispersion, spectroscopic

relaxation times). The choice of properties and the required accuracy

depend on the purpose for which the force ¬eld is to be used. The good-

for-everything force ¬eld will be too complex and computationally inten-

sive to be useful for the simulation of large systems with limited accuracy

requirements. Thus there is not one ideal force ¬eld, but a hierarchy

depending on the application.

• Ab initio computability, meaning that the model parameters should in

principle be obtained from ¬tting to high-level quantum chemistry cal-

culations. This opens the way to construct reliable terms for unusual

molecules for which su¬cient experimental data are not readily available.

If ab initio parametrization does not yield acceptable results and empir-

ical re¬nement is necessary, the model designer should ¬rst consider the

question whether some interaction type or inherent approximation has

been overlooked.

Ideally, a very accurate parent force ¬eld, derivable from quantum calcu-

6.3 Force ¬eld descriptions 185

lations, should be constructed to form the top of the hierarchy. Then, by

constraining terms or properties that do ¬‚uctuate in the parent model, to

their average values in a given set of conditions, simpler child models can

be derived with more limited applicability. This process may be repeated

to produce even simpler and more restricted grandchildren. In this way

simple and e¬cient force ¬elds for a limited range of applications may be

derived without the need for extensive reparametrization of each new force

¬eld. This strategy has been advocated by Saint-Martin et al. (2005) and

shown to be successful for the MCDHO model.33 The full model has internal

¬‚exibility and is polarizable; if the ¬‚exibility is constrained at the average

value obtained in a given environment and under given conditions, a sim-

pler “child” model emerges that is valid for a range of similar environments

and conditions. The same applies to constraining induced dipoles to yield

a “grandchild”: a simple four-site e¬ective pair potential, valid in a more

limited range of conditions. In a similar fashion child models with simpler

force truncation schemes could be constructed, with corrections obtained

from the average long-range contributions within the parent model.

The parent model should explicitly express separate aspects: any omis-

sion will lead to e¬ective incorporation of the omitted aspect into terms of

another physical nature at the expense of robustness and transferability of

the model. We list a number of these e¬ects:

• Quantum character of nuclear motions These can be included in a ther-

modynamically correct (but dynamically questionable) way by replacing

each nucleus by a path integral in imaginary time, approximated by a

string of beads as explained in Chapter 5. “Grandparent” models with

complete path integral implementations are considerably more complex.

More approximately, but su¬cient for most applications, these quantum

e¬ects may be estimated by quantum corrections to second order in , as

detailed in Section 3.5, or they may be incorporated in a re¬nement of the

model that includes Feynman“Hibbs quantum widths of the nuclei. Such

second-order corrections are not su¬cient for oscillators with frequencies

much above kB T /h.

• Quantum character of high-frequency vibrations There are no good meth-

ods to include quantum vibrational states dynamically into a classical

system. The best one can do is to make an adiabatic approximation and

include the equilibrium distribution over quantum states as a potential of

mean force acting on the nuclei. Consider an single oscillating bond (the

theory is readily extended to coupled vibrations) between two masses m1

33 See also Hern´ndez-Cobos et al. (2005).

a

186 Molecular dynamics

and m2 , with the bond length as a general coordinate q. Let q0 be the

bond length for which there is no net force (potential force plus centrifugal

force); q0 depends on the environment and the velocities. The deviation

ξ = q ’ q0 can be separated as a quantum degree of freedom acting in a

quadratic potential with its minimum at ξ = 0. It will have an oscillator

(angular) frequency ω = k/μ, where k is the force constant and μ the

reduced mass m1 m2 /(m1 + m2 ). We now wish to treat the total system as

a reduced dynamic system with q constrained to q0 , and omit the quan-

tum degree of freedom ξ. This we call a ¬‚exible constraint because the

position q0 ¬‚uctuates with external forces and angular velocities. In order

to preserve the correctness of thermodynamic quantities, the potential of

the reduced (i.e., ¬‚exibly constrained) system must be replaced by the po-

tential of mean force V mf with respect to the omitted degree of freedom

(this is the adiabatic approximation):

1

ω + kB T ln(1 ’ e’

V mf = Vclass + ω/kB T

). (6.78)

2

Since in principle ω is a function of the classical coordinates through its

dependence on the force constant, the potential of mean force will lead

to forces on the classical system and energy exchange with the quan-

tum oscillator. However, this dependence is generally small enough to be

negligible. In that case it su¬ces to make posteriori corrections to the

total energy. The important point is to impose ¬‚exible constraints on

the system. As Hess et al. (2002), who give a more rigorous statistical-

mechanical treatment, have shown, the iterations necessary to implement

¬‚exible constraints do not impose a large computational burden when

polarization iterations are needed anyway.

It should be noted that fully ¬‚exible classical force ¬elds will add an

incorrect kB T per oscillator to the total energy. Appropriate corrections

are then mandatory.

• Intramolecular structural response to external interactions Both full ¬‚ex-

ibility and ¬‚exible constraints take care of intramolecular structural re-

sponse to external interactions. Generally the structural deviations are

small, but they represent a sizeable energetic e¬ect, as the restoring in-

tramolecular potentials are generally quite sti¬. Structural ¬‚uctuations

also react back on the environment and cause a coupling between mechan-

ical forces and polarization. Not only the structural response to external

forces should be taken into account, but also the structural response to

external electric ¬elds.

• Intramolecular electronic response to external interactions This, in fact,

6.3 Force ¬eld descriptions 187

is the polarizability response, through which the intramolecular charge

distribution responds to electrical (and in the case of shell models also

mechanical) forces from the environment. Their incorporation in parent

models is mandatory.

• Long-range dispersion forces With reference to the discussion in Section

6.3.4 on page 159, we may state that accurate parent models should in-

clude very long range dispersion interactions, preferably by evaluating the

corresponding lattice sum for periodic systems.

• Long-range electrical interactions From the discussion in Section 6.3.5 on

page 164 it is clear that parent models must evaluate long-range Coulomb

interactions as lattice sums for periodic systems, or else (e.g., for clusters)

use no truncation at all.

• Non-additivity of repulsion and attraction Repulsion between electronic

distributions based on exchange is not strictly pair-additive. It is expected

that the nonadditivity is at least partly taken care of by repulsions between

moving charge distributions, but this remains to be investigated. A parent

model should at least have an estimate of the e¬ect and preferably contain

an appropriate three-body term. Also the dispersion interaction is non-

additive, but this so-called Axilrod“Teller e¬ect34 is of order r’9 in the

distance and probably negligible.

• E¬ects of periodicity Parent models should at least evaluate the e¬ects of

periodicity and preferably determine the in¬nite box-size limit for all eval-

uated properties. Especially for electrolyte solutions e¬ects of periodicity

are a matter of concern.