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
nZ
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):253287, 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
385397. 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):309321, 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
nZ
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
nZ
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
nZ
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
r0
+
ψ(r).
Lemma 1 (Poisson’s summation formula). If D(r) =
nZ
3
δ(r Ln),
then
b
D(h) =
1
V
dZ
3
δ(h L
T
d).
Lemma 2. Consider the distributions ρ(r) =
N
i=1
b
i
δ(r r
i
), η(r) =
ρ(r), and D(r) =
nZ
3
δ(r Ln). The following identity holds:
nZ
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
=
N
i=1
N
j=1
b
i
b
j
Z
δ(s (r
i
r
j
)) D(s r) ds
=
N
i=1
N
j=1
b
i
b
j
D(s (r
i
r
j
))
=
nZ
3
N
i=1
N
j=1
b
i
b
j
δ(s (r
i
r
j
+ Ln)).
Finally,
Z
ψ(r)
(
ρ ? η ? D
)
(r) dr
=
Z
ψ(r)
nZ
3
N
i=1
N
j=1
b
i
b
j
δ(s (r
i
r
j
+ Ln)) dr
=
nZ
3
N
i=1
N
j=1
b
i
b
j
ψ(kr
i
r
j
+ Lnk).
tutorial on spme 7
At this point we can go back to (4) and rewrite it as
U
1
=
1
2
Z
ψ(r)
(
ρ ? η ? D
)
(r) dr
1
2
N
i=1
b
2
i
ψ(0), (5)
Furthermore, using Plancherel’s Theorem
7
and the fact that ρ ? η ? D
7
T. H. Wolff. Lectures on harmonic
analysis, volume 29 of University Lecture
Series. American Mathematical Society,
Providence, RI, 2003
is real-valued, we can write
Z
ψ(r)
(
ρ ? η ? D
)
(r) dr =
Z
b
ψ(h)
\
(
ρ ? η ? D
)
(h) dh, (6)
where h = khk. Consequently, equation (5) becomes
U
1
=
1
2
Z
b
ψ(h)
\
(
ρ ? η ? D
)
(h) dh
1
2
N
i=1
b
2
i
ψ(0), (7)
By the properties of the Fourier transform of a convolution
8
and
8
T. H. Wolff. Lectures on harmonic
analysis, volume 29 of University Lecture
Series. American Mathematical Society,
Providence, RI, 2003
Poisson’s summation formula 1,
\
(
ρ ? η ? D
)
(h) =
b
ρ(h)
b
ρ(h)
1
V
dZ
3
δ(h L
T
d)
=
b
ρ(h)
2
1
V
dZ
3
δ(h L
T
d).
From the above and (6), we have
Z
b
ψ
\
(
ρ ? η ? D
)
=
Z
b
ψ(h)
b
ρ(h)
2
1
V
dZ
3
δ(h L
T
d)
!
dh
=
1
V
dZ
3
b
ψ(kL
T
dk)
b
ρ(L
T
d)
2
=
1
V
m
b
ψ(m)
b
ρ(m)
2
=
1
V
m6=0
b
ψ(m)
b
ρ(m)
2
+
1
V
lim
m0
b
ψ(m)
b
ρ(0)
2
=
1
V
m6=0
b
ψ(m)
b
ρ(m)
2
+
1
V
lim
m0
b
ψ(m)
N
i=1
b
i
!
2
,
where m = L
T
d and m = kmk. So, to sum up, we have arrived at
the identity:
U
1
=
1
2
N
i=1
N
j=1
i6=j
b
i
b
j
ψ(kr
i
r
j
k) +
1
2
nZ
3
n6=0
N
i=1
N
j=1
b
i
b
j
ψ(kr
i
r
j
+ Lnk)
=
1
2V
m6=0
b
ψ(m)
b
ρ(m)
2
+
1
2V
b
ψ(0)
N
i=1
b
i
!
2
1
2
ψ(0)
N
i=1
b
2
i
,
tutorial on spme 8
Fast Fourier Transform and Particle-Mesh method
Ewald summation using optimal parameters requires O(N
3/2
)
operations
9
but it can be modified so that it involves only about
9
S. W. de Leeuw, J. W. Perram, and
E. R. Smith. Simulation of Electrostatic
Systems in Periodic Boundary Con-
ditions. I. Lattice Sums and Dielectric
Constants. Proceedings of the Royal
Society A: Mathematical, Physical and
Engineering Sciences, 373(1752):2756,
October 1980
O(N log N) operations by using the Fast Fourier Transform
10
.
10
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; and R. W.
Hockney and J. J. W. Eastwood. Com-
puter Simulation Using Particles. Adam
Hilger, 1988
We now concern ourselves with summing (and later taking deriva-
tives of) the expression
m6=0
b
ψ(m)
b
ρ(m)
2
. (8)
By the definition of the (inverse) Fourier Transform, we see that
b
ρ(m) =
N
j=1
b
j
e
2πi m
T
r
j
,
where r
j
= L x
j
for some x
j
[0, 1)
3
. Moreover,
b
ρ(m) =
N
j=1
b
j
e
2πi m
T
r
j
=
N
j=1
b
j
e
2πi d
T
L
1
L x
j
=
N
j=1
b
j
e
2πi d
T
x
j
.
We are going to approximate the complex exponential terms above
by means of Cardinal B-Splines
11
. An order p B-spline function is
11
I. J. Schoenberg. Cardinal Spline
Interpolation. Society for Industrial and
Applied Mathematics, 1973
the p-fold convolution of the uniform density in the unit interval (see
Figure ). By construction, each order p B-spline is of class C
p1
, their
0
1
1 0 1 2 3 4 5
0
1
1 0 1 2 3 4 5
0
1
1 0 1 2 3 4 5
0
1
1 0 1 2 3 4 5
y
x
Order 1
y
x
Order 2
y
x
Order 3
y
x
Order 4
Figure 9: Basic B-splines
integral equals 1, and they have compact support.
Remark. Observe that these splines play the role of approximations to
the identity when used in a meshing of the unit cell. It is in this sense that
they play a role in approximating the Dirac delta functions in ρ.
tutorial on spme 9
Now, to discretize the simulation box we consider K
1
, K
2
, K
3
N
and let
K =
K
1
K
2
K
3
.
For each r = Lx, we set u = Kx and we have
e
2πi m
T
r
= e
2πi d
T
x
= e
2πi d
T
K
1
u
`Z
3
W(u `) e
2πi d
T
K
1
`
=
k
nZ
3
W(u (k + Kn)) e
2πi d
T
K
1
(k+Kn)
=
k
nZ
3
W(u (k + Kn)) e
2πi d
T
K
1
k
e
2πi d
T
n
| {z }
=1
=
k
nZ
3
W(u (k + Kn)) e
2πi d
T
K
1
k
,
where we have used the Euclidean algorithm to write
` = k + Kn with k
i
{
0, 1, . . . , K
i
1
}
and n Z
3
.
We also define the forward (unnormalized) DFT of the 3D array
Q
k
as
b
Q
d
= DFT[Q
k
]
d
=
K
1
1
k
1
=0
K
2
1
k
2
=0
K
3
1
k
3
=0
Q
k
e
2πi d
T
K
1
k
.
And, likewise, the backward DFT is
IDFT[
b
Q
d
]
k
=
K
1
1
d
1
=0
K
2
1
d
2
=0
K
3
1
d
3
=0
b
Q
d
e
2πi d
T
K
1
k
.
Using the previous results, we write
b
ρ(m) =
N
j=1
b
j
e
2πi d
T
x
j
k
N
j=1
b
j
nZ
3
W(u (k + Kn))
| {z }
=Q
k
e
2πi d
T
K
1
k
=
k
Q
k
e
2πi d
T
K
1
k
=
b
Q
d
.
Remark. It is important to notice that since the splines have compact
support, the infinite series above is actually a finite sum.
Remark. The following diagram depicts the relationships between the
different objects involved in the SPME method:
ρ
Q
IFT
y
y
IDFT
b
ρ
b
Q
tutorial on spme 10
Going back to (8), we carry out the following approximation:
m6=0
b
ψ(m)
b
ρ(m)
2
=
m6=0
b
ψ(m)
b
ρ(m)
b
ρ(m)
d6=0
b
ψ(kL
T
dk)
b
Q
d
b
Q
d
Now define
A
d
=
b
ϕ(kL
T
hk), if h 6= 0,
0, if h = 0.
with
h
i
(d) =
d
i
, if 0 d
i
<
1
2
K
i
,
K
i
d
i
, if
1
2
K
i
d
i
< K
i
.
(note that
K
i
2
h
i
<
K
i
2
1). Then,
K
1
2
1
d
1
=
K
1
2
K
2
2
1
d
2
=
K
2
2
K
3
2
1
d
3
=
K
3
2
b
ψ(kL
T
dk)
b
Q
d
b
Q
d
=
K
1
1
d
1
=0
K
2
1
d
2
=0
K
3
1
d
3
=0
A
d
b
Q
d
b
Q
d
.
Applying the following property of the Fourier transform (see
12
),
12
T. H. Wolff. Lectures on harmonic
analysis, volume 29 of University Lecture
Series. American Mathematical Society,
Providence, RI, 2003
Z
b
f dg =
Z
b
g d f ,
we write (using the notation of the previous section and keeping in
mind that the result works in the distributional sense too)
Z
b
ψ
b
η
b
ρ =
Z
c
b
ψ
b
η ρ =
Z
b
C ? ρ
ρ
where we have used that
c
f g =
b
f ?
b
g and that, for f = f (x),
b
b
f (x) =
f (x) (again, see
13
).
13
T. H. Wolff. Lectures on harmonic
analysis, volume 29 of University Lecture
Series. American Mathematical Society,
Providence, RI, 2003
Therefore, the algorithm boils down to:
Q
k
b
Q
d
= DFT(Q)
d
A
d
· DFT(Q)
d
IDFT(A · DFT(Q))
k
= (IDFT(A) ? Q)
k
(IDFT(A) ? Q)
k
· Q
k
,
Assorted results in Harmonic analysis
Theorem 3 (The convolution is commutative).
Theorem 4 (The Fourier transform of a convolution is the product of
the Fourier transforms).
Theorem 5.
14
For absolutely continuous measures µ, ν in R
n
. We have
14
T. H. Wolff. Lectures on harmonic
analysis, volume 29 of University Lecture
Series. American Mathematical Society,
Providence, RI, 2003
Z
ˆ
µ dν =
Z
ˆ
ν dµ
Theorem 6 (Plancherel’s theorem).
tutorial on spme 11
References
[dLPS80] S. W. de Leeuw, J. W. Perram, and E. R. Smith. Simulation
of Electrostatic Systems in Periodic Boundary Conditions.
I. Lattice Sums and Dielectric Constants. Proceedings of
the Royal Society A: Mathematical, Physical and Engineering
Sciences, 373(1752):2756, October 1980.
[DYP93] T. Darden, D. York, and L. Pedersen. Particle mesh Ewald:
An N log(N) method for Ewald sums in large systems. The
Journal of Chemical Physics, 98(12):10089, 1993.
[EPB
+
95] 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.
[Ewa21] P. P. Ewald. Die Berechnung optischer und elektrostatis-
cher Gitterpotentiale. Annalen der Physik, 369(3):253287,
1921.
[HE88] R. W. Hockney and J. J. W. Eastwood. Computer Simulation
Using Particles. Adam Hilger, 1988.
[ND57] B. R. A. Nijboer and F. W. De Wette. On the calculation of
lattice sums. Physica, 23(1-5):309321, 1957.
[Sch73] I. J. Schoenberg. Cardinal Spline Interpolation. Society for
Industrial and Applied Mathematics, 1973.
[Wil06] D. E. Williams. Accelerated convergence treatment of
R
n
lattice sums. In Uri Shmueli, editor, International
Tables for Crystallography, volume B, pages 385397. Kluwer
Academic Publishers, 2006.
[Wol03] T. H. Wolff. Lectures on harmonic analysis, volume 29 of
University Lecture Series. American Mathematical Society,
Providence, RI, 2003.