moment within the sphere gives a ¬eld in the origin; the total quadrupole

moment gives a ¬eld gradient, etc. When the system contains charges only,

and no explicit higher multipoles, we need only the reaction ¬eld to compute

the forces on each particle.13 Including the direct interactions, the force on

i is given by

⎡ ¤

⎢ r ij ⎥

⎢ ⎥

qi r ij

qj 3 ’

⎢ qj f (µr , κ) 3 ⎥ .

Fi = (6.50)

⎣ rc ¦

4πµ0 rij

j; R ij ¤r c j

Rij ¤rc

(i,j)∈exclst

/

def

Here, r ij = r i ’ r j , hence the minus sign in the reaction ¬eld term. Note

that the inclusion of the reaction ¬eld simply modi¬es the force function; in

the tin-foil or conducting boundary condition (µr = ∞ or κ = ∞ : f = 1)

the modi¬cation produces a shifted force, which is continuous at the cut-o¬

radius. This shifted force is well approximated for media with high dielec-

tric constant, such as water with f (µr , κ) = 0.981. Because of the smooth

force the RF modi¬cation yields acceptable dynamics, even in cases where

a reaction ¬eld is not appropriate because of anisotropy and inhomogeneity

of the medium.

The potential energy function that generates the forces of (6.50) is ob-

tained by integration:

⎡ ¤

⎢ ⎥

⎢ 2⎥

1 f (µr , κ) 2

’1 ’1

’ (rij ’ rc )⎥ .

qi ⎢

V (r) = qj (rij rc ) + qj

⎣ ¦

3

4πµ0 2rc

j>i; R ij ¤r c

i j≥i

Rij ¤rc

(i,j)∈exclst

/

(6.51)

The forces are not the exact derivatives of this potential because the trun-

cation of F ij is based on distances Rij between reporter positions of groups,

which depend not only on rij , but also on the position of other particles.

Since Rij ≈ rij , this e¬ect is small and is neglected. In addition, there are

discontinuities in the reaction-¬eld energies when dipoles cross the cut-o¬

boundary, which lead to impulsive contributions to the forces. Since these

are not incorporated into the forces, the e¬ective (i.e., integrated-force) po-

13 See Tironi et al. (1995) and H¨ nenberg and van Gunsteren (1998) for the full multipole

u

equations. Note that force ¬elds with explicit dipoles need to consider the reaction ¬eld gradient

as well!

6.3 Force ¬eld descriptions 167

tentials will slightly deviate from the real reaction-¬eld potentials. This

situation is similar to the e¬ect of truncated long-range dispersion poten-

tials as discussed above (page 159).

The reaction ¬eld addition does not account for the polarization e¬ects in

the medium beyond the cut-o¬ due to charges rather than dipoles. Naively

one may say that the Born reaction potential of a charge (see (13.72) on

page 345) is a constant potential with zero gradient, which does not lead to

extra forces on the charges and will therefore not in¬‚uence the equilibrium

distributions. Therefore the Born correction can be applied afterwards to

the results of a simulation. However, this reasoning disregards the discon-

tinuities that occur in the potential when two charges enter or leave the

interaction range. When the impulsive derivative of such potentials are not

incorporated into the force, the e¬ective potential “ which is the integral of

the forces “ does not equal the Born correction. For charges these e¬ects are

more severe than for dipoles. Let us consider the simple case of two charges.

The Born-corrected, rc -truncated, interaction energy for a set of charges is:

⎡ ¤

⎢ g(µ, κ) ⎥

qj

1 1

qi ⎢ ⎥,

’

U= qj (6.52)

⎣ rc ¦

4πµ0 rij 2

i j>i j

rij ¤rc rij ¤rc

1

def

g(µ, κ) = 1 ’ . (6.53)

µr (1 + κrc )

Applying this to two charges q1 and q2 at a distance r, we ¬nd for the

potential energy:

112 2 g(µ, κ)

r > rc : V (r) = V∞ = ’ (q1 + q2 ) ,

4πµ0 2 rc

1 q1 q2 1 g(µ, κ)

r ¤ rc : V (r) = ’ (q1 + q2 )2

4πµ0 r 2 rc

1

q1 q2 (r’1 ’ g(µ, κ)rc ).

’1

= V∞ + (6.54)

4πµ0

The e¬ective potential obtained from the integrated truncated force would

simply yield a shifted potential:

1

q1 q2 (r’1 ’ rc ).

’1

Ve¬ (r) = (6.55)

4πµ0

It is interesting that in conducting-boundary conditions (g(µ, κ = 1) the

potential equals a shifted potential, plus an overall correction V∞ equal to

the sum of the Born energies of the isolated charges.

While reaction ¬elds can be included in the force ¬eld with varying degrees

168 Molecular dynamics

of sophistication,14 they are never satisfactory in inhomogeneous systems

and systems with long-range correlation. The latter include in practice all

systems with explicit ions. The polarization in the medium is not simply

additive, as is assumed when reaction ¬elds are included per particle. For

example, the ¬eld in the medium between a positive and a negative charge

at a distance larger than rc is strong and induces a strong polarization, while

the ¬eld between two positive charges cancels and produces no polarization.

The total polarization energy is then (in absolute value) larger, resp. smaller

than predicted by the Born-corrected force ¬eld of (6.52), which is indi¬erent

for the sign of the charge. One is faced with a choice between inaccurate

results and excessive computational e¬ort.

There are three ways out of this dilemma, neither using a reaction ¬eld

correction. The ¬rst is the application of continuum corrections, treated

below. The second is the use of the fast multipole method (FMM), which

relies on a hierarchical breakdown of the charges in clusters and the evalua-

tion of multipole interactions between clusters. The method is implemented

in a few software packages, but is rather complex and not as popular as

the lattice summation methods. See the discussion in Section 13.9 on page

362. The third, most recommendable method, is the employment of e¬cient

lattice summation methods. These are, of course, applicable to periodic sys-

tems, but even non-periodic clusters can be cast into a periodic form. There

are several approaches, of which the accurate and e¬cient smooth-particle

mesh-Ewald (SPME) method of Essmann et al. (1995) has gained wide

popularity. These methods are discussed at length in Chapter 13, Section

13.10 on page 362, to which the reader is referred.

Continuum correction methods are due to Wood (1995).15 Consider a

“model world” in which the Hamiltonian is given by the force ¬eld used,

with truncated, shifted or otherwise modi¬ed long-range Coulomb inter-

action, and possibly periodic boundary conditions. Compare with a “real

world” with the full long-range interactions, possibly of in¬nite extension

without periodic boundary conditions. Now assume that the di¬erence in