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 `

, `

∈

(see Figure ). Then the matrix L = [`

, `

] ∈ 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

L x

where x

∈ [0, 1)

for i = 1, 2 and that r

= kr

− r

k (see

Figure ). Between these two atoms there is a force that depends on

Figure 2: Two atoms at a distance r

their distance and whose corresponding potential energy, U

), 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 =

∑

i=1

∑

j=1

i6=j

). (1)

By solving Newton’s second law

we would be able to simulate the

Force equals mass times acceleration.

tutorial on spme 2

Energy

Distance (r)

U(r) = r

−12

− r

−6

(r) = r

−12

(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

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 ﬁgure of N

= 6.022 × 10

(Avogadro’s number) which

would capture the real physics of the system.

One trick to somehow get closer to N

is to use a tractable value

of N and impose periodic boundary conditions (PBCs). This means

that we now have an inﬁnite 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

≈ 10

and

N log

N ≈ 10

The total potential energy (1) then becomes

H =

∑

i=1

∑

j=1

i6=j

(kr

− r

k) +

∑

n∈Z

n6=0

∑

i=1

∑

j=1

(kr

− r

+ Lnk).

The naive approach to PBCs requires O(N

) calculations and is

unfeasible, hence the interest of ﬁnding efﬁcient ways of carrying

out this task. Ewald summation

takes advantage of ideas from

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

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 sufﬁciently 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

and then

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

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

(r) + U

(r), (2)

tutorial on spme 4

(a) Paul Peter Ewald

(1888-1985).

(b) Tom Darden and

Lee Pedersen.

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

, decays fast as r goes to inﬁnity (so

that it can be truncated) and the long-range term, U

, 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

estab-

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 deﬁne

(r) = 1 −

Γ(n/2, α

π r

)

Γ(n/2)

where α

is a parameter to be discussed later, then we clearly have

(

(r) + (1 − ϕ

(r))

)

We are going to concern ourselves with the case n = 6 (see Figure ),

and to deﬁne ψ( r) = r

−6

(r).

The free parameter α is determined by numerically solving the

equation

−6

(

1 − ϕ

cutoff

; α)

)

= ε

with prescribed values of r

cutoff

and ε (the function 1 − ϕ

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.05

0.1

1 2 3 4 5

Energy

−6

(r)

−6

(1 − ϕ

(r))

−6

Figure 7: Decomposition of the function

−6

into short-range and long-range

terms.

Recall that in the case of N > 2 particles, we have

H =

∑

i=1

∑

j=1

i6=j

(kr

− r

k) +

∑

n∈Z

n6=0

∑

i=1

∑

j=1

(kr

− r

+ Lnk).

where r

is the distance between atoms i and j. We can split the

pairwise potential U

above as in (2).

Ewald summation

The summation of the long-range terms (with b

> 0),

∑

i=1

∑

j=1

i6=j

ψ(kr

− r

k) +

∑

n∈Z

n6=0

∑

i=1

∑

j=1

ψ(kr

− r

+ Lnk)

(3)

can be rearranged as

∑

n∈Z

∑

i=1

∑

j=1

ψ(kr

− r

+ Lnk) −

∑

i=1

ψ(0), (4)

where ψ(0) = lim

r→0

ψ(r).

Lemma 1 (Poisson’s summation formula). If D(r) =

∑

n∈Z

δ(r − Ln),

then

D(h) =

∑

d∈Z

δ(h − L

−T

d).

Lemma 2. Consider the distributions ρ(r) =

∑

i=1

δ(r − r

), η(r) =

ρ(−r), and D(r) =

∑

n∈Z

δ(r − Ln). The following identity holds:

∑

n∈Z

∑

i=1

∑

j=1

ψ(kr

− r

+ Lnk) =

ψ(r)

(

ρ ? η ? D

)

(r) dr,

where r = krk.

tutorial on spme 6

−4 −2 0 2 4

-4 -2 0 2 4

Figure 8: Example of Poisson’s summa-

tion formula.

Proof. First,

ρ ? η (r) =

ρ(s) η(r − s) ds

ρ(s) ρ(s − r) ds

∑

i=1

∑

j=1

δ(s − r

) δ(r − s − r

) ds

∑

i=1

∑

j=1

δ(r − (r

− r

)).

By the associative property of the convolution, we have

ρ ? η ? D =

(

ρ ? η

)

? D,

ρ ? η ? D (r) =

ρ ? η (s) D(s − r) ds

∑

i=1

∑

j=1

δ(s − (r

− r

)) D(s − r) ds

∑

i=1

∑

j=1

D(s − (r

− r

))

∑

n∈Z

∑

i=1

∑

j=1

δ(s − (r

− r

+ Ln)).

Finally,

ψ(r)

(

ρ ? η ? D

)

(r) dr

ψ(r)

∑

n∈Z

∑

i=1

∑

j=1

δ(s − (r

− r

+ Ln)) dr

∑

n∈Z

∑

i=1

∑

j=1

ψ(kr

− r

+ Lnk).

tutorial on spme 7

At this point we can go back to (4) and rewrite it as

ψ(r)

(

ρ ? η ? D

)

(r) dr −

∑

i=1

ψ(0), (5)

Furthermore, using Plancherel’s Theorem

and the fact that ρ ? η ? D

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

ψ(r)

(

ρ ? η ? D

)

(r) dr =

ψ(h)

(

ρ ? η ? D

)

(h) dh, (6)

where h = khk. Consequently, equation (5) becomes

ψ(h)

(

ρ ? η ? D

)

(h) dh −

∑

i=1

ψ(0), (7)

By the properties of the Fourier transform of a convolution

and

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) =

ρ(h)

ρ(−h)

∑

d∈Z

δ(h − L

−T



ρ(h)



∑

d∈Z

δ(h − L

−T

d).

From the above and (6), we have

(

ρ ? η ? D

)

ψ(h)



ρ(h)



∑

d∈Z

δ(h − L

−T

∑

d∈Z

ψ(kL

−T

dk)



ρ(L

−T



∑

ψ(m)



ρ(m)



∑

m6=0

ψ(m)



ρ(m)



lim

m→0

ψ(m)



ρ(0)



∑

m6=0

ψ(m)



ρ(m)



lim

m→0

ψ(m)

∑

i=1

where m = L

−T

d and m = kmk. So, to sum up, we have arrived at

the identity:

∑

i=1

∑

j=1

i6=j

ψ(kr

− r

k) +

∑

n∈Z

n6=0

∑

i=1

∑

j=1

ψ(kr

− r

+ Lnk)

∑

m6=0

ψ(m)



ρ(m)



ψ(0)

∑

i=1

−

ψ(0)

∑

i=1

tutorial on spme 8

Fast Fourier Transform and Particle-Mesh method

Ewald summation using optimal parameters requires O(N

3/2

)

operations

but it can be modiﬁed so that it involves only about

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):27–56,

October 1980

O(N log N) operations by using the Fast Fourier Transform

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

ψ(m)



ρ(m)



. (8)

By the deﬁnition of the (inverse) Fourier Transform, we see that

ρ(m) =

∑

j=1

2πi m

where r

= L x

for some x

∈ [0, 1)

. Moreover,

ρ(m) =

∑

j=1

2πi m

∑

j=1

2πi d

−1

L x

∑

j=1

2πi d

We are going to approximate the complex exponential terms above

by means of Cardinal B-Splines

. An order p B-spline function is

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

p−1

, their

−1 0 1 2 3 4 5

Order 1

Order 2

Order 3

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

, K

∈ N

and let

K =





For each r = Lx, we set u = Kx and we have

2πi m

= e

2πi d

= e

2πi d

−1

≈

∑

`∈Z

W(u − `) e

2πi d

−1

∑

n∈Z

W(u − (k + Kn)) e

2πi d

−1

(k+Kn)

∑

n∈Z

W(u − (k + Kn)) e

2πi d

−1

2πi d

| {z }

∑

n∈Z

W(u − (k + Kn)) e

2πi d

−1

where we have used the Euclidean algorithm to write

` = k + Kn with k

∈

{

0, 1, . . . , K

− 1

}

and n ∈ Z

We also deﬁne the forward (unnormalized) DFT of the 3D array

= DFT[Q

]

−1

∑

−1

∑

−1

∑

−2πi d

−1

And, likewise, the backward DFT is

IDFT[

]

−1

∑

−1

∑

−1

∑

2πi d

−1

Using the previous results, we write

ρ(m) =

∑

j=1

2πi d

≈

∑

j=1

∑

n∈Z

W(u − (k + Kn))

| {z }

2πi d

−1

∑

2πi d

−1

Remark. It is important to notice that since the splines have compact

support, the inﬁnite series above is actually a ﬁnite sum.

Remark. The following diagram depicts the relationships between the

different objects involved in the SPME method:

≈

−−−−→ Q

IFT



IDFT

≈

−−−−→

tutorial on spme 10

Going back to (8), we carry out the following approximation:

∑

m6=0

ψ(m)



ρ(m)



∑

m6=0

ψ(m)

ρ(−m)

ρ(m) ≈

∑

d6=0

ψ(kL

−T

dk)

−d

Now deﬁne







ϕ(kL

hk), if h 6= 0,

0, if h = 0.

with

(d) =







, if 0 ≤ d

− d

, if

≤ d

< K

(note that −

≤ h

− 1). Then,

−1

∑

=−

−1

∑

=−

−1

∑

=−

ψ(kL

−T

dk)

−d

−1

∑

−1

∑

−1

∑

−d

Applying the following property of the Fourier transform (see

T. H. Wolff. Lectures on harmonic

analysis, volume 29 of University Lecture

Series. American Mathematical Society,

Providence, RI, 2003

f dg =

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)

ρ =

η ρ =



C ? ρ



where we have used that

f g =

f ?

g and that, for f = f (x),

f (x) =

f (−x) (again, see

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:

→

−d

= DFT(Q)

→ A

· DFT(Q)

→ IDFT(A · DFT(Q))

= (IDFT(A) ? Q)

→ (IDFT(A) ? Q)

· Q

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.

For absolutely continuous measures µ, ν in R

. We have

T. H. Wolff. Lectures on harmonic

analysis, volume 29 of University Lecture

Series. American Mathematical Society,

Providence, RI, 2003

µ dν =

ν 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):27–56, 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):253–287,

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):309–321, 1957.

[Sch73] I. J. Schoenberg. Cardinal Spline Interpolation. Society for

Industrial and Applied Mathematics, 1973.

[Wil06] D. E. Williams. Accelerated convergence treatment of

−n

lattice sums. In Uri Shmueli, editor, International

Tables for Crystallography, volume B, pages 385–397. Kluwer

Academic Publishers, 2006.

[Wol03] T. H. Wolff. Lectures on harmonic analysis, volume 29 of

University Lecture Series. American Mathematical Society,

Providence, RI, 2003.