et al. (1999), and also by Paricaud et al. (2005), both authors using polariz-

abilities on all three atoms. However, shell models with charge distributions

for the shells are far more natural and e¬cient. In a revealing study, ¬tting

electrostatic models to ab initio potentials, Tanizaki et al. (1999) showed

that point charge models tend to become counterintuitive (as positive charge

on lone-pair positions) unless shielding distributed charges are used. As the

MCDHO model of Saint-Martin et al. (2000) shows, a single shell for the

molecule su¬ces, neglecting the polarizabilities on hydrogen. This model

does not only reproduce liquid behavior, but has excellent second virial co-

e¬cients (see Fig. 6.9) and excellent liquid“vapor coexistence and critical

behavior (Hern´ndez-Cobos et al., 2005). But also with the MCDHO model

a

the polarization is too strong in the condensed phase, leading to a high

dielectric constant,30 and some re¬nement seems necessary.

6.3.8 Energies and forces for polarizable models

Consider a collection of ¬xed source terms, which may be charges qi only or

include higher multipoles as well. The sources may be charge or multipole

density distributions. They produce an electric ¬eld E 0 (r) at a position

r, which diverges at a point source itself. In addition there will be induced

dipoles μs or shells with charge qk at sites r s (these may be atoms or virtual

s

k k

sites), with the shell connected by a harmonic spring with spring constant

s

kk to an atom or virtual site r k . The latter may or may not carry a ¬xed

charge as well. These dipoles or shells (which may be points or distributions)

29 See Jedlovszky and Richardi (1999) for a comparison of water models in the critical region.

30 Saint-Martin et al. (2005).

182 Molecular dynamics

produce an extra ¬eld E ind (r) at a position r. The total energy can best

be viewed as the sum of the total electrostatic interaction plus a positive

polarization energy Vpol needed to create the induced dipoles μ:

(μs )2

k

Vpol = (6.67)

2±k

k

or to stretch the springs

1s s

kk |r k ’ r k |2 .

Vpol = (6.68)

2

k

The total electrostatic energy consists of the source“source interaction Vqq

and the source“dipole plus dipole“dipole interaction Vqμ + Vμμ , or for shells

the source“shell plus shell“shell interactions Vqs + Vss . In these sums every

pair should be counted only once and, depending on the model used, some

neighbor interactions must be excluded (minimally a shell has no Coulomb

interaction with the site(s) to which it is bound and dipoles do not interact

with charges on the same site). The form of the potential function for a pair

interaction depends on the shape of the charge or dipole distribution. The

polarizations (dipoles or shell positions) will adjust themselves such that the

total potential energy is minimized. For dipoles this means:

‚Vtot

= 0, Vtot = Vqq + Vqμ + Vμμ + Vpol , (6.69)

‚μsk

or, for shells:

‚Vtot

= 0, Vtot = Vqq + Vqs + Vss + Vpol . (6.70)

‚r s

k

It is easily shown that this minimization leads to

μs = ±k E tot (r s ), (6.71)

k k k

or, for shells:

s

qk tot s

rs = r k + s E k (r k ). (6.72)

k

kk

Here E tot is minus the gradient of Vtot . The shell displacement corresponds

to an induced dipole moment equal to q s times the displacement and hence

corresponds to a polarizability

(qk )2

s

±k = s. (6.73)

kk

6.3 Force ¬eld descriptions 183

The polarization energy appears to compensate just half of the electrostatic

interaction due to induced charges or shells.

This, in fact, completes the equations one needs. The forces on site i are

determined as minus the gradient to the position r i of the total potential en-

ergy.31 The electrostatic interactions between arbitrary charge distributions

are not straightforward. The interaction energy between a point charge q at

a distance r from the origin of a spherical charge distribution q s w(r) (with

∞ 2

0 4πr w(r) dr = 1) is given (from Gauss™ law) by integrating the work

done to bring the charge from in¬nity to r:

∞ r

qq s qq s s

1

s 2

Vqq = 4πr w(r ) dr = φ (r). (6.74)

r2

4πµ0 4πµ0

r 0

In general the function φs (r) and its derivative (to compute ¬elds) can be

best tabulated, using cubic splines for interpolation (see Chapter 19). For a

Slater-type exponential decay of the density:32

1 ’2r/»

w(r) = e , (6.75)

π»3

the potential function is given by

r ’2r/»

φ(r) = 1 ’ 1 +

e . (6.76)

»

For two interacting charge distributions the two-center integral needs to be

tabulated, although Carrillo-Trip et al. (2003) give an approximate expres-

sion that gives good accuracy at all relevant distances (note that shells are

never very close to each other):

ss »1 »2 + (»1 + »2 )2

q1 q2

[1 ’ (1 + z)e’2z ), z =

Vss ≈ r12 . (6.77)

(»1 + »2 )3

r12

6.3.9 Towards the ideal force ¬eld

General force ¬elds for molecular systems are not of su¬cient quality for

universal use. They are often adapted to speci¬c interactions and the terms

are not su¬ciently transferable. The latter is clear from the ever increas-

ing number of atom types deemed necessary when force ¬elds are re¬ned.

A large amount of work has been done on models for water, being impor-

tant, small, polar, polarizable, hydrogen-bonding, and endowed with a huge

31 In the dipolar case the induced dipoles are kept constant in taking the derivatives. This is

correct because the partial derivatives to the dipole are zero, provided the dipoles have been

completely relaxed before the force is determined.

32 » is the decay length of the corresponding wave function; 32% of the charge density is within

a distance » from the origin.

184 Molecular dynamics

knowledge body of experimental properties. The excessive number of pub-

lished water models indicate a state of confusion and a lack of accepted

guidelines for the development of force ¬elds. Still, it seems that an ad-