374 Electromagnetism

restrict the spatial solutions to lattice points. Interpolation is then needed to

obtain the energies and forces acting on charges. They scale as N log N and

form the long-range methods of choice, e.g., as implemented in the particle“

mesh“Ewald (PME) method of Darden et al. (1993) who use a Gaussian

spread function and a Lagrange interpolation, or “ preferably “ the more

accurate smooth particle“mesh“Ewald (SPME) method of Essmann et al.

(1995), who use a B-spline interpolation.18 The advantage of using B-spline

interpolation is that the resulting potential function is twice continuously

di¬erentiable if the order of the spline is at least four; smooth forces can

therefore be immediately obtained from the di¬erentiated potential. With

Lagrange interpolation the interpolated potential is only piecewise di¬eren-

tiable and cannot be used to derive the forces. The SPME works as follows,

given a system of N charges qi at positions ri within a unit cell of base

vectors a.b, c:

• Choose three integers K1 , K2 , K3 that subdivide the unit cell into small,

reasonably isotropic, grid cells. Choose an appropriate Ewald parame-

ter β (see (13.181)), a cuto¬-radius rc for the short-range interaction,

which should not exceed half the length of the smallest base vector, and

a cut-o¬ radius in reciprocal space. Choose the order p of the B-spline

interpolation, which should be at least 4 (cubic spline). A cuto¬ of 4

times the standard deviation of the Gaussian spread function implies that

√

βrc = 2 2. A grid size a/K1 , b/K2 , c/K3 of about 0.3/β and a recipro-

cal cut-o¬ of 1.5β would be reasonable for a start. The optimal values

depend on system size and density; they should be adjusted for optimal

computational e¬ciency, given a desired overall accuracy.

• Compute the structure factors using the exponential spline interpolation

as explained in Chapter 19, using (19.87) and (19.88) on page 554. Note

that for odd order p the value m = ±K/2 must be excluded

For actual programs the reader is referred to the original authors. A full

description of the SPME algorithm can be found in Griebel et al. (2003).

These methods use a splitting between short-range and long-range poten-

tials and solve the long-range potential on a grid; they are really variants

of the PPPM particle“particle particle“mesh) method developed earlier by

Hockney and Eastwood (1988). In the PPPM method the charges are dis-

tributed over grid points; the Poisson equation is solved and the potentials

and ¬elds are interpolated, using optimized local functions that minimize the

total error. Deserno and Holm (1998a) have compared the accuracies and

18 For details on B-splines see Chapter 19, Section 19.7, on page 548.

13.10 Potentials and ¬elds in periodic systems of charges 375

e¬ciencies of various methods and evaluated the error in a PPPM scheme

in a second paper (1998b).

It seems that SPME methods have not yet been applied to other charge

spread functions than Gaussian ones, although short-range force functions

that go exactly and smoothly to zero at the cut-o¬ radius would have the

advantage of avoiding cut-o¬ noise in the short-range forces. A suitable

candidate would be the cubic function, discussed on page 369.

13.10.7 Potentials and ¬elds in periodic systems of charges and

dipoles

Certain force ¬elds describe charge distributions not only with a set of

charges, but also with a set of dipoles, or even higher multipoles. The dipoles

may be permanent ones, designed to describe the charge distribution with

a smaller number of sites. The may also be induced dipoles proportional to

the local ¬eld, as in certain types of polarizable force ¬elds. In the latter

case the induced dipoles are determined in an iterative procedure until con-

sistency, or they may be considered as variables that minimize a free energy

functional. In all cases it is necessary to determine potentials and ¬elds,

and from those energies and forces, from a given distribution of charges and

dipoles.

The methods described above, splitting interactions into short- and long-

range parts, can be extended to include dipolar sources as well. One must

be aware that such extensions considerably complicate the computation of

energies and forces, as the dipolar terms involve higher derivatives than are

required for charges. It could well be advantageous to avoid dipoles “ and

certainly higher multipoles “ if the problem at hand allows formulation in

terms of charges alone. Here we shall review the methods in order to give

the reader a ¬‚avor of the additional complexity, referring the reader to the

literature for details.

Ewald summations including multipolar sources have been worked out by

Smith (1982); based on this work, Toukmaji et al. (2000) extended PME

methods to dipolar sources. The e¬ect of adding dipole moments μi to the

sources qi is that qi is replaced by qi +μi ·∇i , which has consequences for the

structure factors as well as for the short- and long-range force and energy

terms. Consider the energy U12 between two charges q1 at r 1 and q2 at r 2 :

1 1

U12 = q1 ¦(r 1 ) = q1 q2 , (13.200)

4πµ0 r12

376 Electromagnetism

where r12 = |r 1 ’ r 2 |.19 When dipoles are present this modi¬es to (cf

(13.123))

U12 = q1 ¦(r 1 ) ’ μ1 · E(r 1 ) = (q1 + μ1 · ∇1 )¦(r 1 ), (13.201)

with the potential given by (see (13.136))

q2 1 1

’ μ2 · ∇1 = (q2 + μ2 · ∇2 ) .

4πµ0 ¦(r 1 ) = (13.202)

r12 r12 r12

The interaction thus changes from q1 q2 /r12 to

1

U12 = (q1 + μ1 · ∇1 )(q2 + μ2 · ∇2 ) . (13.203)

r12

We note that this represents the total electrostatic interaction; for the en-

ergy of polarizable systems one should add the energy it costs to create the

induced dipole moments, which is a quadratic form in μ (like i μ2 /2±i )

i

for the case of linear polarizability.

This replacement works through all equations; for example, the structure

factor Qk (13.166) now becomes

(qj + μj · ∇j )e’ik·r j

Qk =

j

(qj ’ 2πiμj · k)e’ik·r j , j ∈ unit cell,

= (13.204)

j

with consequences for the spline interpolation procedure. The reader is

referred to Toukmaji et al. (2000) for details.

Exercises

13.1 Consider an electron as a (classical) sphere with homogeneous charge

distribution. What would its radius be when the total ¬eld energy

equals the relativistic rest energy mc2 ?

13.2 What is the amplitude in vacuum of E and B in a circular laser

beam of 50 mW monochromatic radiation (» = 632.8 nm), when the

beam has a diameter of 2 mm?

13.3 Consider an in¬nite planar surface with vacuum on one side and a

dielectric medium with relative dielectric constant µr on the other

side, with a charge q situated on the vacuum side at a distance

d from the plane. Show that the following potential satis¬es the

boundary conditions (13.56): on the vacuum side the direct Coulomb

19 See the beginning of Section 13.8 on page 353 for a discussion of charge interactions and where

a factor 2 should or should not appear.

Exercises 377

potential of the charge plus the vacuum potential of an image charge

qim = ’q(µr ’ 1)/(µr + 1) at the mirror position; on the medium side

the direct Coulomb potential of the charge, divided by a factor µe¬ .

Express µe¬ in µr .

13.4 Compare the exact solvation free energy of two charges q1 and q2 ,

both at a distance d from the planar surface that separates vacuum

and medium (µr ) as in Exercise 13.3, and separated laterally by a

distance r12 , with the generalized Born expression (13.103) using

Still™s expression (13.105).

13.5 Verhoeven and Dymanus (1970) have measured the quadrupole mo-

ment of D2 O. They report the values:

θa = 2.72(2), θb = ’0.32(3), θc = ’2.40,

in 10’26 esu.cm2 , for the diagonal components in a coordinate system