k = 2π(m1 a— + m2 b— + m3 c— ), (13.164)

with m1 , m2 , m3 ∈ Z, which enumerate the Fourier terms, and a— , b— , c— the

reciprocal lattice vectors. See Section 12.9 on page 331 for a description of

the reciprocal lattice and corresponding Fourier transforms. For the eval-

uation of the scalar product k · r it is generally easier to use the relative

coordinates (ξ, ·, ζ) of r:

k · r = 2π(m1 ξ + m2 · + m3 ζ), (13.165)

with ξ = r · a— , etc. Now construct the Fourier transform (often called

structure factors) of the ensemble of point charges:

qj e’ik·r j , j ∈ unit cell.

def

Qk = (13.166)

j

Since we have broadened the point charges with the spread function, the

charge density ρl (r) (13.163) is the convolution of the point charge distribu-

tion and the spread function. The Fourier transform Pk of ρl (r) therefore

l

is the product of the charge structure factors and the Fourier transform Wk

of the spread function:

1

ρl (r)e’ik·r dr

def

l

Pk = (13.167)

VV

= Qk Wk , (13.168)

where

1

w(r + Tn)e’ik·(r+Tn) dr

def

Wk = (13.169)

V V

n

1

w(r)e’ik·r dr

=

V all space

∞ π

1

dθ 2πr2 w(r) sin θ e’ikr cos θ

= dr

V 0 0

∞

1 sin kr

= 4πrw(r) dr (13.170)

V k

0

Here V means integration over one unit cell, and we have used the fact

that the spread function is isotropic. We see that Wk depends only on the

13.10 Potentials and ¬elds in periodic systems of charges 367

absolute value k of k. The validity of (13.168) can easily be checked by

evaluating (13.167) using (13.163).

The Poisson equation (13.163) in reciprocal space reads

’k 2 µ0 ¦l = ’Pk ,

l

(13.171)

k

and thus

Qk Wk

¦l = ; k = 0. (13.172)

k

µ0 k 2

Note that k = 0 must not be allowed in this equation and ¦0 is therefore not

de¬ned; indeed for electroneutral systems Q0 = 0. Electroneutrality there-

fore is required; if the system is charged the potential does not converge and

electroneutrality must be enforced by adding a homogeneous background

charge of opposite sign. This enforces Q0 = 0. The real-space potential

φl (r) follows up to a constant by Fourier transformation:

Qk Wk ik·r

φl (r) = ¦l eik·r = e . (13.173)

k

µ0 k 2

k=0 k=0

The total energy can be expressed in terms of a sum over wave vectors

1 1

k ’2 Qk Q’k Wk ’ Uself .

l l

UC = qi φ(r i ) = (13.174)

2 2µ0

k=0

i

The self-energy contained in the long-range energy is a constant given by

the j = i, n = 0 part of (13.162):

1

qi lim [r’1 ’ •s (r)]

l 2

Uself = (13.175)

4πµ0 r’0

i

(for •s see (13.159)). Similarly, the interaction energy between excluded

pairs ij for which the long-range energy must be corrected is

qi qj ’1

[rij ’ •s (rij )].

l

Uexcl = (13.176)

4πµ0

The long-range force on particle i can be evaluated from the gradient of

the potential:

Qk Wk

F l = qi E l (r i ) = ’qi (∇φl (r))r i = ’qi ikeik·r i . (13.177)

i 2

µ0 k

k=0

The sum is real since for every k-term there is a complex conjugate ’k-

term. The self-energy term does not produce a force; the ij exclusion term

produces a long-range force on particle i (and the opposite force on j):

qi qj ’2 r ij

[rij ’ f s (rij )]

Fli,excl = (13.178)

4πµ0 rij

368 Electromagnetism

(for f s see (13.161)) which should be subtracted from the long-range force

evaluation.

13.10.3 Gaussian spread function

In the special case a Gaussian distribution is chosen for the spread function,

the solution is expressed in sums of analytical functions and the classical

Ewald summation is obtained (Ewald, 1921). The advantage of a Gaussian

function is that its Fourier transform is also a Gaussian function, and both

the real-space and reciprocal-space functions taper o¬ quickly and can be

restricted to a limited range. The Gaussian function contains an inverse

width parameter β (the variance of the Gaussian distribution is 1/2β 2 ); if

β is small, the spread function is wide and the real-space summation has

many terms. The reciprocal functions, on the other hand, then decrease

rapidly with increasing k. For large β the inverse is true. Therefore β can

be tuned for the best compromise, minimizing the computational e¬ort for

a given error margin.

Noting that

x

2

e’u du,

def 2

erf (x) = √ (13.179)

π 0

∞

2

e’u du,