664 Energies and forces
The second technical problem associated with the Ewald sum is the high com-
putational overhead associated with the calculation of the struct∑ ure factorS(g) =
iqiexp(ig·ri) for large systems. AsLincreases, the number ofg-vectors satisfying
the criterion|g| ≤gmaxbecomes quite large, and the particle sum inS(g), which
also increases, must be carried out at eachg-vector. Thesmooth particle-mesh Ewald
(SPME) method, an important advance introduced by Essmannet al.(1995), substan-
tially reduces the cost of the Ewald sum and improves the scaling with system size.
The SPME approach can be viewed as a kind of “smearing” of the charges, which are
located at arbitrary spatial pointsri, over a finite set of points on a regular rectan-
gular lattice. In practice, this “smearing” is realized mathematically by employing an
interpolation formula to express the exponential exp(ig·ri) in terms of the exponential
factors that depend on the coordinates of the grid points ratherthan the arbitrary
particle coordinateri. Essmannet al.achieved this interpolation using a set of spline
functions known as thecardinal B-splinefunctions. The functions, denotedMn(x), are
defined as follows:M 2 (x) = 1−|x− 1 |for 0≤x≤2 andM 2 (x) = 0 forx <0 and
x >2. Forn >2,Mn(x) are defined via the recursion relations
Mn(x) =x
n− 1Mn− 1 (x) +n−x
n− 1Mn− 1 (x−1). (B.23)The cardinal B-spline functions have several important properties. First, they have
compact support, meaning thatMn(x) is zero outside the interval 0≤x≤n. Second,
Mn(x) isn−2 times continuously differentiable. Third, the derivative dMn/dxcan
be obtained fromMn− 1 (x) using another recurrence relation
d
dxMn(x) =Mn− 1 (x)−Mn− 1 (x−1). (B.24)Two other useful properties are:
Mn(x) =Mn(n−x)∑∞
j=−∞Mn(x−j) = 1. (B.25)A plot of the first few cardinal B-spline functions is shown in Fig. B.2.
We now introduce the scaled particle coordinatessi=V−^1 /^3 ri(which reduce to
si=ri/Lfor a cubic box), so that exp(ig·ri) = exp(2πin·si). Let the lattice be
cubic withNlpoints along each direction and define new coordinatesui =Nlsi,
such thatui,α∈[0,Nl],α=x,y,z. Then, exp(2πin·si) = exp(2πin·ui/Nl). We
can then approximate the exponential exp(2πinαui,α/Nl) using a cardinalB-spline
interpolation formula as
e^2 πinαui,α/Nl≈bn(nα)∑∞
k=−∞Mn(ui,α−k)e^2 πinαk/Nl, (B.26)where