Tutorial on the smooth particle-mesh Ewald algo-
rithm
These notes explain the smooth particle-mesh Ewald algorithm (SPME)
in detail. The accompanying library libpme6 is a fully-featured imple-
mentation of SPME that can be used for experimentation and reference.
This content can also be found in PDF format here.
Introduction
Motivation and overview of the SPME algorithm
Consider two atoms, labeled 1 and 2, contained in a parallelepiped
with sides determined by the (positively oriented) vectors `
1
, `
2
, `
3
∈
R
3
(see Figure ). Then the matrix L = [`
1
, `
2
, `
3
] ∈ R
3×3
is such
Figure 1: Simulation box spanned by
three vectors.
that V = det L > 0. Now suppose that atom i is at a position r
i
=
L x
i
where x
i
∈ [0, 1)
3
for i = 1, 2 and that r
ij
= kr
i
− r
j
k (see
Figure ). Between these two atoms there is a force that depends on
Figure 2: Two atoms at a distance r
ij
.
their distance and whose corresponding potential energy, U
ij
=
U
ij
(r
ij
), is shown in Figure . If instead of just two atoms we had a
total of N > 2 atoms (see Figure ), then the total potential energy
would be equal to
H =
1
2
N
∑
i=1
N
∑
j=1
i6=j
U
ij
(r
ij
). (1)
By solving Newton’s second law
1
we would be able to simulate the
1
Force equals mass times acceleration.
tutorial on spme 2
Energy
Distance (r)
U(r) = r
−12
− r
−6
U
0
(r) = r
−12
U
1
(r) = −r
−6
Figure 3: Decomposition of the Lennard
Jones 12-6 potential into repulsive
(powers of −12) and dispersive (powers
of −6) contributions.
Figure 4: Argon gas (N = 10
3
atoms).
molecular system of N atoms that interact according to the (pairwise)
forces given by F = −∇U. In practice, whatever value of N that can
be tractable with current computers is still far away from the more
realistic figure of N
A
= 6.022 × 10
23
(Avogadro’s number) which
would capture the real physics of the system.
One trick to somehow get closer to N
A
is to use a tractable value
of N and impose periodic boundary conditions (PBCs). This means
that we now have an infinite lattice of simulation boxes with identical
copies of the particles in each box. By computing the interactions
between the atoms in the original box and their copies (only the
copies in the neighboring simulation boxes are typically taken into
account —this is called the minimum image convention—), we can
make our simulations more realistic.
tutorial on spme 3
Figure 5: The protein DHFR solvated
in water (left) and two views (center,
right) without solvent. This system has
N = 23536 atoms, so N
2
≈ 10
8
and
N log
2
N ≈ 10
5
.
The total potential energy (1) then becomes
H =
1
2
N
∑
i=1
N
∑
j=1
i6=j
U
ij
(kr
i
− r
j
k) +
1
2
∑
n∈Z
3
n6=0
N
∑
i=1
N
∑
j=1
U
ij
(kr
i
− r
j
+ Lnk).
The naive approach to PBCs requires O(N
2
) calculations and is
unfeasible, hence the interest of finding efficient ways of carrying
out this task. Ewald summation
2
takes advantage of ideas from
2
P. P. Ewald. Die Berechnung optischer
und elektrostatischer Gitterpotentiale.
Annalen der Physik, 369(3):253–287, 1921
Harmonic analysis to reduce the complexity to O(N
3/2
) operations,
which is still impractical if we want to solve the equations of mo-
tion of a large system. The Particle-Mesh Ewald (PME) method
3
3
T. Darden, D. York, and L. Pedersen.
Particle mesh Ewald: An N log(N)
method for Ewald sums in large sys-
tems. The Journal of Chemical Physics,
98(12):10089, 1993
uses the Fast Fourier Transform to bring the complexity down to
O(N log N), and this algorithm has been a staple of Molecular Simu-
lation since its appearance. The Smooth Particle-Mesh Ewald (SPME)
improves on PME by giving a sufficiently smooth energy function
whose derivative can be obtained analytically, which results in more
realistic simulations due to improved energy conservation.
Some of the people involved in the development of the SPME
algorithm appear in Figure .
Smooth Particle-Mesh Ewald
We derive the formulas for Ewald summation following
4
and then
4
D. E. Williams. Accelerated conver-
gence treatment of R
−n
lattice sums. In
Uri Shmueli, editor, International Tables
for Crystallography, volume B, pages
385–397. Kluwer Academic Publishers,
2006
we introduce the smooth Particle-Mesh Ewald method from
5
.
5
U. Essmann, L. Perera, M. L.
Berkowitz, T. Darden, H. Lee, and
L. G. Pedersen. A smooth particle mesh
Ewald method. The Journal of Chemical
Physics, 103(19):8577, 1995
Decomposition of the potential energy function
Ewald summation can be applied whenever the pairwise potential
energy U = U(r) is proportional to r
−n
for some n ∈ N (a concrete
instance of U appears in Figure ). In that case, we can write U as a
sum of the form
U(r) = U
0
(r) + U
1
(r), (2)
tutorial on spme 4
(a) Paul Peter Ewald
(1888-1985).
(b) Tom Darden and
Lee Pedersen.
(c) Mike Crowley.
Figure 6: Some of the people behind
SPME. Pictures reproduced from
http://journals.iucr.org/a/issues/
2012/01/00/wx0006/ and http://www.
psc.edu/science/Darden/.
where the short-range term, U
0
, decays fast as r goes to infinity (so
that it can be truncated) and the long-range term, U
1
, is very spread
out and smooth, so that its Fourier transform is concentrated (and
can also be truncated).
There are many ways of splitting the potential energy but
6
estab-
6
B. R. A. Nijboer and F. W. De Wette.
On the calculation of lattice sums.
Physica, 23(1-5):309–321, 1957
lishes the optimality for a class of functions involving the (incom-
plete) Gamma function. If we define
ϕ
n
(r) = 1 −
Γ(n/2, α
2
π r
2
)
Γ(n/2)
,
where α
2
is a parameter to be discussed later, then we clearly have
1
r
n
=
1
r
n
(
ϕ
n
(r) + (1 − ϕ
n
(r))
)
.
We are going to concern ourselves with the case n = 6 (see Figure ),
and to define ψ( r) = r
−6
ϕ
6
(r).
The free parameter α is determined by numerically solving the
equation
r
−6
(
1 − ϕ
6
(r
cutoff
; α)
)
= ε
with prescribed values of r
cutoff
and ε (the function 1 − ϕ
n
must be
decreasing for r > r
cutoff
). This makes the short-range part of the
computation involving effectively vanish for r > r
cutoff
.
tutorial on spme 5
0
0.05
0.1
1 2 3 4 5
Energy
r
r
−6
ϕ
6
(r)
r
−6
(1 − ϕ
6
(r))
r
−6
Figure 7: Decomposition of the function
r
−6
into short-range and long-range
terms.
Recall that in the case of N > 2 particles, we have
H =
1
2
N
∑
i=1
N
∑
j=1
i6=j
U
ij
(kr
i
− r
j
k) +
1
2
∑
n∈Z
3
n6=0
N
∑
i=1
N
∑
j=1
U
ij
(kr
i
− r
j
+ Lnk).
where r
ij
is the distance between atoms i and j. We can split the
pairwise potential U
ij
above as in (2).
Ewald summation
The summation of the long-range terms (with b
i
> 0),
U
1
=
1
2
N
∑
i=1
N
∑
j=1
i6=j
b
i
b
j
ψ(kr
i
− r
j
k) +
1
2
∑
n∈Z
3
n6=0
N
∑
i=1
N
∑
j=1
b
i
b
j
ψ(kr
i
− r
j
+ Lnk)
(3)
can be rearranged as
U
1
=
1
2
∑
n∈Z
3
N
∑
i=1
N
∑
j=1
b
i
b
j
ψ(kr
i
− r
j
+ Lnk) −
1
2
N
∑
i=1
b
2
i
ψ(0), (4)
where ψ(0) = lim
r→0
+
ψ(r).
Lemma 1 (Poisson’s summation formula). If D(r) =
∑
n∈Z
3
δ(r − Ln),
then
b
D(h) =
1
V
∑
d∈Z
3
δ(h − L
−T
d).
Lemma 2. Consider the distributions ρ(r) =
∑
N
i=1
b
i
δ(r − r
i
), η(r) =
ρ(−r), and D(r) =
∑
n∈Z
3
δ(r − Ln). The following identity holds:
∑
n∈Z
3
N
∑
i=1
N
∑
j=1
b
i
b
j
ψ(kr
i
− r
j
+ Lnk) =
Z
ψ(r)
(
ρ ? η ? D
)
(r) dr,
where r = krk.
tutorial on spme 6
−4 −2 0 2 4
-4 -2 0 2 4
D
b
D
Figure 8: Example of Poisson’s summa-
tion formula.
Proof. First,
ρ ? η (r) =
Z
ρ(s) η(r − s) ds
=
Z
ρ(s) ρ(s − r) ds
=
N
∑
i=1
N
∑
j=1
b
i
b
j
Z
δ(s − r
i
) δ(r − s − r
j
) ds
=
N
∑
i=1
N
∑
j=1
b
i
b
j
δ(r − (r
i
− r
j
)).
By the associative property of the convolution, we have
ρ ? η ? D =
(
ρ ? η
)
? D,
so
ρ ? η ? D (r) =
Z
ρ ? η (s) D(s − r) ds