static continuum theory, since the di¬erence concerns long-range e¬ects on

a scale much coarser than atomic detail, and at such distances from real

14 See for advanced reaction ¬elds, e.g., Hummer et al. (1996), Bergdorf et al. (2003) and Christen

et al. (2005).

15 These methods are based on earlier ideas of Neumann (1983), who applied continuum methods

to interpret simulation results on dielectric behavior. Several authors have made use of con-

tinuum corrections, and most applications have been reviewed by Bergdorf et al. (2003), who

considered the e¬ects of truncation, reaction ¬eld functions and periodic boundary conditions

on ionic hydration and on the interaction between two ions in a dielectric medium.

6.3 Force ¬eld descriptions 169

charges that the dielectric response can be assumed to be linear and ¬eld-

independent. Now correct the simulation results of the model world with

the di¬erence obtained by continuum theory. Separate corrections can be

obtained for charge“dipole and dipole“dipole cuto¬s, and for the e¬ect of

periodic boundary conditions. The principle is as simple as that, but the

implementation can be quite cumbersome and not applicable to all possible

cases. The principle of the method is as follows.

Consider a system of sources, taken for simplicity to be a set of charges

qi at positions r i , in a dielectric medium with linear (i.e., ¬eld-independent)

local dielectric constant µ = µr µ0 and without electrostriction (i.e., ¬eld-

dependent density). The local polarization P (r) (dipole density) is given

by a product of the (local) electric susceptibility χ and the electric ¬eld E:16

χ = µr ’ 1,

P = µ0 χE, (6.56)

The ¬eld is determined by the sum of the direct ¬eld of the sources and the

dipolar ¬elds of the polarizations elsewhere:

G(r ’ r i ) + T (r ’ r )P (r ) dr .

E(r) = (6.57)

i

Here G(r) is the ¬eld produced at r by a unit charge at the origin, and

T (r) is the tensor which “ multiplied by the dipole vector “ yields the ¬eld

at r due to a dipole at the origin. The vector function G equals minus the

gradient of the potential ¦(r) by a unit charge at the origin. For G and

T one can ¬ll in the actual truncated or modi¬ed functions as used in the

simulation, including periodic images in the case of periodic conditions. For

example, for a truncated force, these in¬‚uence functions are:

1r

for r ¤ rc ,

G(r) =

4πµ0 r3

= 0 for r > rc , (6.58)

1

(3x± xβ ’ r2 δ±β )r’5 for r ¤ rc ,

=

T ±β

4πµ0

= 0 for r > rc . (6.59)

The ¬rst term in (6.57) is the vacuum ¬eld of the set of charges. The

integral equation (6.56) with (6.57) can then be numerically solved for P

and E and from there energies can be found. The total electrostatic energy

is 1 µE 2 dr, which diverges for point charges and must be corrected for

2

the vacuum self-energy (see (13.42) on page 340); the polarization energy is

16 See Chapter 13 for the basic electrostatic equations.

170 Molecular dynamics

Table 6.3 Continuum corrections to the simulated hydration free energy of

a sodium ion in 512 water molecules in a periodic box, simulated with

i’w w’w for

Coulomb forces truncated at rc for ion-water interactions and rc

water“water interactions. Column 3 gives the simulation results

(Straatsma and Berendsen, 1988), column 4 gives the Born correction

resulting from the ion-water cuto¬, column 5 gives the Born-corrected

results, column 6 gives the continuum correction resulting from water-water

cuto¬ (Wood, 1995) and column 7 gives the corrected results. Column 8

gives the corrections resulting from periodicity

1 2 3 4 5 6 7 8

i’w

i’w w’w

rc rc ”Gsim CBorn 3+4 Cw’w 5+6 CPBC

nm nm kJ/mol kJ/mol kJ/mol kJ/mol kJ/mol kJ/mol

’424 ’76.2 ’500 ’479 ’0.1

0.9 0.9 +21.0

’444 ’65.3 ’509 ’477 ’0.2

1.05 0.9 +32.5

’461 ’57.2 ’518 ’476 ’0.2

1.2 0.9 +42.5

’404 ’76.2 ’480 ’477 ’0.1

0.9 1.2 +03.6

’429 ’57.2 ’486 ’472 ’1.3

1.2 1.2 +13.9

given by the interaction of each charge with the induced dipoles:

1

Upol = ’ P (r) · G(r ’ r i ) dr.

qi (6.60)

2

i

The factor 1 is implicit in the basic equations (see Chapter 13, Section 13.6

2

on page 339), but can also be derived from changing the charges from zero

to the full charge and integrating the work it costs to increment the charges

against the existing polarization ¬eld. While the solution of the integral

equation (which can be transformed into a matrix equation, to be solved

by a suitable iterative method, see Bergdorf et al., 2003) is needed for the

“model world,” the reference calculation for the “real world” can be cast

into a Poisson equation and solved by more elegant methods. Finally, the

di¬erence in solvation energy should be used to correct “model-world” sim-

ulations. The continuum approximations are rather poor in the immediate

neighborhood of an ion, where speci¬c solvent structures and dielectric sat-

uration may occur; such regions must be covered correctly by simulation

with atomic detail. These regions, however, do not contribute signi¬cantly

to the di¬erence between “model world” and “real world” because they fall

well within the cuto¬ range of the model potentials.

Wood (1995) applied continuum corrections to a series of MD simulations

by Straatsma and Berendsen (1988) on the free energy of hydration of an

6.3 Force ¬eld descriptions 171

ion, with various values of water“water and water“ion cuto¬ radii. Because

these early simulations were limited in size (one ion and 512 water molecules

in a cubic periodic box, with cuto¬s up to 1.2 nm), they serve well for

demonstration purposes because the corrections are substantial. Table 6.3

shows the results of the free energy gained by charging a neon atom to

become a Na+ ion. Di¬erent values for the ion“water and the water“water

cuto¬ were used (columns 1 and 2). When the simulation results (column

3) are corrected for the Born energy (column 4), using the ion“water cuto¬,

the results (column 5) appear to be inconsistent and depend strongly on the

water“water cuto¬. The reason for this is that the ion polarizes the water

with the result that there is a substantial water“water interaction beyond the

water“water cuto¬ radius. When the continuum correction is applied to the

water“water cuto¬ as well (column 6), the result (column 7) is now consistent

for all simulation conditions, i.e., within the accuracy of the simulations

of about 3 kJ/mol. The correction for the periodic condition (column 8)

appears to be rather small and to increase with the cuto¬ radius. This is

consistent with the more detailed calculations of Bergdorf et al. (2003), who

observe that the corrections for periodicity are quite small when applied to

truncated potentials, but much larger when applied to full lattice sums. This

can be interpreted to mean that the use of full lattice sums enhances the

periodicity artifacts. On the other hand, it prevents the “ generally worse “

artifacts due to truncation.

6.3.6 Polarizable force ¬elds

The dominant and most relevant omission in the usual force ¬elds is the

incorporation of electronic polarizability. The electron distribution around

a given nuclear con¬guration depends on the presence of external electric

¬elds; in ¬rst order a dipole moment μ is induced proportional to the (ho-

mogeneous) electric ¬eld E:

μ = ±E, (6.61)