PHYSICAL REVIEW E 84, 016707 (2011)

Coulomb interactions in charged fluids

Graziano Vernizzi,1,* Guillermo Iván Guerrero-Garcı́a,2 and Monica Olvera de la Cruz2,3,4

1
Department of Physics and Astronomy, Siena College, Loudonville, New York 12211, USA
2
Department of Materials Science, Northwestern University, Evanston, Illinois 60208, USA
3
Department of Chemical Engineering, Northwestern University, Evanston, Illinois 60208, USA
4
Department of Chemistry, Northwestern University, Evanston, Illinois 60208, USA
(Received 10 August 2010; revised manuscript received 22 May 2011; published 25 July 2011)
The use of Ewald summation schemes for calculating long-range Coulomb interactions, originally applied
to ionic crystalline solids, is a very common practice in molecular simulations of charged fluids at present.
Such a choice imposes an artificial periodicity which is generally absent in the liquid state. In this paper we
propose a simple analytical O(N 2 ) method which is based on Gauss’s law for computing exactly the Coulomb
interaction between charged particles in a simulation box, when it is averaged over all possible orientations of a
surrounding infinite lattice. This method mitigates the periodicity typical of crystalline systems and it is suitable
for numerical studies of ionic liquids, charged molecular fluids, and colloidal systems with Monte Carlo and
molecular dynamics simulations.

DOI: 10.1103/PhysRevE.84.016707 PACS number(s): 05.10.−a, 82.70.Dd, 41.20.Cv, 61.20.Gy

I. INTRODUCTION means that a system of N charges qi (i = 1, . . . ,N) each at

position bi is endowed with an infinite set of identical copies
A challenging problem that still plagues modern numer-
(images) that stem from identifying the simulation box as the
ical simulations is how to include long-range electrostatic
primitive cell of an infinite regular lattice. In mathematical
interactions in a tractable way, a question that in its most
terms, let a1 , a2 , and a3 be three vectors that define the three-
general formulation lies undefeated [1]. In infinite periodic
dimensional primitive cell, then any other cell of the periodic
crystalline systems this problem has been overcome in the last
lattice is identified by the lattice vector  = n1 a1 + n2 a2 +
century by using an approach originally proposed by Ewald
n3 a3 , ni ∈ Z. The set {ri } of the positions of all charges is
[2], which allows the calculation of electrostatic energies per
ri = bi + ,∀, i = 1,2, . . . N. The total electrostatic energy
unit cell and of the respective Madelung constants [3,4]. This
of such a periodic system is infinite. It is therefore meaningful
motivated the use of such methodology in ionic liquids, and
to consider the electrostatic energy per unit cell:
nowadays Ewald-like schemes are frequently used to take into
account long-range interactions not only in charged liquids, 1   qi qj
U= , ri = bi , rj = bj + , (1)
but also in the modeling of complex biological molecules 2 ,i,j |ri − rj |
in solution that interact electrostatically [5]. The importance
of such Coulomb systems has prompted the improvement of where the N charges interact purely via Coulomb forces within
the original performance of O(N 2 ) operations, advancing to an homogeneous isotropic dielectric environment. The sum
O(N log N ) in [6], or even to O(N ) [7,8], with sophisticated over i is for the charges in the main cell ( = 0, ri = bi )
computational schemes (for an extensive list, see [9]). Despite while the sum over j and  is for all charges in the lattice.
their significant success, it must be noted that Ewald-like The prime symbol indicates that pairs such that |ri − rj | = 0
schemes impose an artificial ordering typical of crystalline are excluded from the sum.
systems, which is not expected in the liquid state. Thus As we mentioned earlier, Eq. (1) is exact for systems with
a desirable advance in the simulation of charged liquids long-range order, such as ionic crystalline solids. However,
would be the ability of computing electrostatic long-range for ionic fluids it is just an approximation [5]. The advan-
interactions correctly, while including the natural isotropy and tage of imposing an artificial periodicity is that Eq. (1) is
homogeneity of fluids at the same time. This is the main issue tractable analytically (besides legitimate questions about its
we address in the present work. mathematical meaning, being a conditionally convergent sum
[10]). A most praised equivalent formula was introduced by
II. THEORETICAL APPROACH Ewald [2], who recast the summation into a form where the
charges interact with each other via an effective interaction
In numerical simulations, the root of the problem stated that includes all the periodic image contributions. Namely,
in the Introduction is that forcing long-range electrostatics long-range terms are expressed in the reciprocal lattice Q
inside a computational box leaves the splintery question of (defined by Q ·  = 2π m, m ∈ Z), and the total energy
what to do with the interactions at the box boundaries, where reads
the long-range nature of electrostatics does not end obviously.
One of the most popular choices, even in simulations of fluids, U = UQ + U + U0 + Ud , (2)

is to adopt periodic boundary conditions. Simply stated, it where the Fourier part is UQ = 2π Q=0 q̂(Q)q̂(−Q)
exp(−Q

2
/4α 2 )/V Q2 , the real space part is U =
,i,j qi qj erfc(Dα)/2D, the constant contribution is
 √
*
gvernizzi@siena.edu U0 = −α i qi2 / π , and the dipolar contribution is

1539-3755/2011/84(1)/016707(7) 016707-1 © 2011 American Physical Society

GRAZIANO VERNIZZI et al. PHYSICAL REVIEW E 84, 016707 (2011)


Ud = 2π ( i qi bi )2 /3V [11]. The cell volume is V = a1 ·
(a2 × a3 ), D = |bi − bj − |, and α is a parameter that in
practical applications
 can be fixed by requiring ∂U/∂α = 0
[12] (if the Q and  are not truncated, then U is always
independent from α). A derivation of Eq. (2) is presented in
the Appendix. A crucial requirement for the validity of Eq. (2)
is that the system must be globally
 electroneutral. In the rest
of this paper we too assume i qi = 0.
In order to tame the unwanted artifacts in simulations of
ionic fluids due to the existence of preferential crystalline
directions and planes of symmetry, a number of authors
suggested performing a spherical average of U over all
possible orientations of the lattice, or considered analogous FIG. 1. (Color online) The positions of the periodic images of the
isotropic radial potentials [13–15]. Such a strategy leads to main cell (square with thick borders) depend on the orientation of the
new effective potentials between charges, and it has been lattice. The arrangement of charges in the rotated periodic lattice (b)
differs from the one in the original lattice (a). Only the charges in the
generalized to other types of interactions such as dipolar,
main cell are not affected by the rotation of the lattice.
multipolar, power law, Lennard-Jones, and exponential po-
tentials [15]. Numerical studies that explored the validity of
The average is over all orthogonal 3 × 3 matrices representing
such computational schemes seem to endorse the effectiveness
proper rotations (det O = 1), that is, the group SO(3). The
of such an approximation (see, e.g., [16,17]; the list is not
volume of the group SO(3) is 8π 2 . We remind that since
exhaustive). Despite the optimistic outlook, the recipe of taking
each rotation matrix can be parametrized by three parameters,
a spherical average of a periodic structure might seem certainly
the integral in Eq. (4) is a three-dimensional integral. By
a circuitous approach at first, since a heterogeneous system
substituting Eq. (1) into Eq. (4), and under the assumption
cannot be both infinitely periodic and spherically symmetric
that the integral and the summation commute, we obtain
at the same time. This relates to the common assumption of
 π
previous approaches, that the simulation box is itself somewhat 1 
homogeneous. U  = qi dω2(1 − cos ω)
2 ,i,j 0
The main focus of this paper is to (re)analyze such a 
problem from a slightly more fundamental standpoint based qj 1
× d 2 , (5)
on the application of Gauss’s law. We obtain a new formula for 8π |bi − bj − O(ω,)|
the spherically averaged periodic potential as we show below.
We begin by assuming that since a liquid is characterized by where we introduced the convenient axis-angle representation
a finite correlation length among its constituents, all long- of a generic SO(3) matrix, i.e., the angle ω around the rotation
range ordering is inconspicuous. Therefore if one imposes axis n̂ ∈ ,  being the unit sphere. The factor 2(1 − cos ω)
a periodicity on the primitive cell (which is larger than the is the density of rotation matrices in the axis-angle parameter
fluid correlation length), the actual position of the periodic space. The key remark now is that each  integral in the
images should be of little relevance in the computation of summation can be interpreted as the potential at bi generated
physical quantities, and consequently they should have little or by a spherical charge distribution centered at bj , with radius
no influence on the structure, dynamics, and thermodynamics |O| = || and with charge density qj /(4π ||2 ) [18]. Such
of the charges within the simulation box. Ergo the orientation a simple observation allows us to evaluate all integrals exactly.
of the lattice must be irrelevant for a liquid system with a In fact, according to Gauss’s law the electrostatic potential
typical correlation length smaller than the simulation box size, V generated by a density of charges distributed uniformly
and different, rotated lattices should lead to the same physics over a sphere of radius R, with charge density σ (hence with
(see Fig. 1). For the sake of preciseness, we define a rotated total charge Q = 4π R 2 σ ), and vanishing at infinity is V (r) =
lattice by Q/R for r  R and V (r) = Q/r for r > R, where r is the
distance from the center of the sphere. We use the compact
notation V (r) = Q + ( Qr − Q )θ (r − R) where θ (x) is the unit

3 R R
(O) = O = ni ai (O) , ai (O) = Oai , (3) step function: θ (x) = 0 for x  0 and θ (x) = 1 for x > 0.
i=1
In the degenerate case with R = 0 the potential is simply
V (r) = Q/r. We thus obtain
where O is a 3 × 3 rotation matrix. The new positions are 
1  qi qj 1  1
ri (O) = bi + (O), meaning that all charges in the main cell U  = + qi qj
2 i,j |bi − bj | 2 =0,i,j ||
are unaffected by the rotation and only the periodic images
  
are actually rotated (see Fig. 1). Let U (O) from Eq. (1) be 1 1
the electrostatic energy of the rotated lattice. We then define + − θ (|bi − bj | − ||) , (6)
|bi − bj | ||
the spherical average over all possible orientations of the
lattice as where we isolated the contribution from  = 0. The first
 term in the second summation symbol is identically
 zero
1 because of the electroneutrality condition q = 0.
U  = dO U (O) . (4) i i
8π 2 SO(3) Moreover, the distance |bi − bj | cannot be larger than

016707-2
COULOMB INTERACTIONS IN CHARGED FLUIDS PHYSICAL REVIEW E 84, 016707 (2011)

3 C
2

1
B
δ
0.2 0.2 0.4
1
A
2

3
FIG. 2. (Color online) (a) Three different kinds of nearest
neighbors of a cubic cell. (b) Regions of the main cell where the 4
corrections to the Coulomb potential are relevant.
FIG. 3. (Color online) The red dots are numerical rotational
averages of Ewald Eq. (2). The thick blue line is obtained by Eq. (7).
the size of the main cell, which in the conventional√case The dashed curve is the energy of a pure Coulomb interaction (without
of a cubic simulation box of size L is |bi − bj |  3L. any corrective term).
In this case we also have || = L n21 + n22 + n23 =
√ √ √ √ √ √
0,L, 2L, 3L, 4L, 5L, 6L, 8L, . . . . However, only III. NUMERICAL TESTS
two terms survive the limits √ imposed by the θ function, that We promptly test our result with a simple example. Let
is, when || = L and || = 2L. us consider a cell of size L = 1 with N = 2 charges q1 = 1
They correspond to two shells of periodic cells surrounding and q2 = −1, at positions b1 = {−1/2, − 1/2, − 1/2} and
the simulation box (see Fig. 2), and they are precisely the 6 first b2 = {δ,δ,δ}, respectively. The parameter δ ∈ [−1/2,1/2]
nearest neighbors with a face in common with the main cell can be used to test different regions inside the main cell.
(type I), and the 12 nearest neighbors with an edge in common In particular, it enters the region where the first correction
(type II). We note explicitly that the 8 nearest neighbors with [yellow sphere in Fig. 2(b), i.e., r = |b1 − b2 | = 1]
activates √
a corner in common with the main cell (type III) do not at δ = (2 3 − 3)/6 ≈ 0.08 and the region where the second
contribute to the average energy. The final form of Eq. (6), corrections
which is also our main result, therefore reads √ activates too√ [blue spherical domain in Fig. 2(b),
i.e., r = 2] at δ = (2 6 − 3)/6 ≈ 0.32.
In Fig. 3 the numerical rotational averages of the Ewald
1  qi qj sum Eq. (2) at different values of δ are represented by red
U  = + qi qj
2 i,j |bi − bj | disks. They compare well with the plot of Eq. (7) (continuous
   thick curve). The plot is characterized by three regions: the first
1 1 region where the two particles and the effective interaction is
× 6 − θ (|bi − bj | − L) + 12
|bi − bj | L a pure Coulomb potential (up to point A), the second region
  
1 1 √ where the two particles are at a distance such that the effective
× −√ θ (|bi − bj | − 2L) . (7) potential contains a contribution from the first θ function (from
|bi − bj | 2L A to B), and the third region where both θ function corrections
contributes to the total energy (from B to C). For comparison,
The first term in the sum is the standard Coulomb energy we also plot the curve corresponding to a pure Coulomb
of the main cell, while the subsequent corrective terms are interaction for all values of δ (dashed curve). From the plot
due to averaging over all lattice orientations. From Eq. (7) it it is evident that our method is different from the so-called
is straightforward to extract the effective force between two minimum image convention [19]. In particular, the growth
charges at distance r, and we find that √ it has innoffensive of the potential at large δ is quite a natural phenomenon in
finite discontinuities at r = L and r = 2L. The θ functions ionic systems where the dipolar term in the Ewald summation
activates only when one of the charges is sufficiently close method is not neglected [20]. The deviations between our result
to the box boundaries. Figure 2(b) illustrates the regions that and the numerical integration of the Ewald formula are simply
are relevant for the corrections to the Coulomb potential. A due to numerical errors. While such errors can be reduced by
necessary condition for the first corrective term [proportional increasing the numerical accuracy of the integration, they are,
to θ (r − L)] to be nonzero is that a charge lies in the however, never completely avoidable since the infinite sums in
region outside a sphere of radius L/2 (the yellow sphere in the Ewald formula must be truncated somewhere necessarily,
Fig. 2 indicated with a thick arrow). Analogously a necessary while our result Eq. (7) is finite and exact, and does not require
condition
√ for the second corrective factor [proportional to any truncation.
θ (r − 2L)] to be nonzero is that a charge lies in the region In order to further illustrate the applicability of Eq. (7) in
near√ the corners of the main cell, at a distance larger than molecular simulations, we present here the result of a series of
L/ 2 from the center of the main cell (outside the blue regions Monte Carlo (MC) simulations. We consider the N V T ensem-
indicated with a dashed arrow). ble of a 1:1 restricted primitive model electrolyte for several

016707-3
GRAZIANO VERNIZZI et al. PHYSICAL REVIEW E 84, 016707 (2011)

TABLE I. Average electrostatic energy per particle and heat capacity for a 1:1 electrolyte, at
different molar concentrations.

This work This work Ewald sums Ewald sums

ρ (M) U ∗ /N kB T Cv /N kB U ∗ /N kB T Cv /N kB

2 −0.6277 ± 0.0028 0.149 ± 0.020 −0.6556 ± 0.0045 0.139 ± 0.020

1 −0.5268 ± 0.0029 0.147 ± 0.020 −0.5493 ± 0.0030 0.141 ± 0.020
0.5 −0.4353 ± 0.0023 0.134 ± 0.018 −0.4544 ± 0.0022 0.138 ± 0.017
0.1 −0.2582 ± 0.0021 0.111 ± 0.015 −0.2703 ± 0.0029 0.104 ± 0.015
0.05 −0.1990 ± 0.0020 0.083 ± 0.011 −0.2089 ± 0.0019 0.087 ± 0.006
0.01 −0.1010 ± 0.0017 0.049 ± 0.007 −0.1067 ± 0.0016 0.051 ± 0.008
0.005 −0.0736 ± 0.0013 0.037 ± 0.005 −0.0785 ± 0.0006 0.039 ± 0.005
0.001 −0.0339 ± 0.0010 0.018 ± 0.003 −0.0372 ± 0.0010 0.017 ± 0.003
0.0005 −0.0244 ± 0.0007 0.012 ± 0.002 −0.0267 ± 0.0007 0.012 ± 0.002
0.0001 −0.0112 ± 0.0003 0.005 ± 0.001 −0.0127 ± 0.0006 0.005 ± 0.001
0.00005 −0.0079 ± 0.0004 0.004 ± 0.001 −0.0091 ± 0.0004 0.003 ± 0.001
0.00001 −0.0037 ± 0.0003 0.001 ± 0.001 −0.0044 ± 0.0002 0.001 ± 0.001

molar concentrations, (i) in the bulk, and (ii) around a charged approach are of the same order of the standard Ewald
macroparticle. All simulations were performed by using a summation method, which is a limit case in absence of angular
standard Metropolis scheme [19,21], in a cubic simulation rotations. Morover, the RDFs in our approximation are similar
box of length L with periodic boundary conditions. Due to the to the ones obtained by the standard Ewald sum method, as
spherical average, our system is not translational invariant and one can see in Fig. 4 for a 2-M monovalent salt.
therefore in the calculation of the radial distribution functions We consider also the case of an electrical double
(RDFs) we consider only particles that are inside a spherical layer around a charged macroparticle. We use a macropar-
shell S of radius L/2, centered in the cubic simulation box. ticle’s diameter D = 20 Å at different valences zM =
Since any periodic charged particle which is located outside 4,8,12,16,20,24,28,32, for two electrolyte’s concentrations:
the spherical shell S is smeared on a homogeneous shell, it 0.01 and 2 M. An important quantity characteristic of this kind
cannot be used in the calculation of the RDFs. Thus for the of Coulombic system is the mean electrostatic potential, which
sake of simplicity we place one ion (or one macroparticle) at is defined by [22,23]
the center of the simulation box and use a cutoff L/2 to obtain   
the corresponding RDFs. The parameters for the 1:1 electrolyte e ∞ 
2
t2
simulations are the diameter of equal-sized ions a = 4.25 Å, ψ(r) = zi ρi gi (t) t − dt. (10)
0 r r
the dielectric constant = 78.5, the temperature T = 298 K, i=1

the valency z+ = |z− | = 1, and the total number of monovalent For both electrolyte concentrations, we calculated the mean
ions N ≈ 512. The length L is fixed by the bulk electrolyte electrostatic potential at the macroparticle’s surface, ψ0 ,
concentration. A Monte Carlo cycle consists of N attempts to and at the Helmholtz plane, ψH P (which is defined as the
move an arbitrary ion. The thermalization process takes about closest approach distance between ions and the macroparticle).
5 × 104 MC cycles; all canonical averages are computed over
a 2 × 105 MC sweep at full thermal equilibrium.
The electrostatic energy per particle, U ∗ /NkB T , and the
2
heat capacity Cv are defined via
U 
U∗ = , (8) 1.5
4π 0
gij(r)

2
U ∗2 − U ∗ 1
Cv = , (9)
kB T 2
where U ∗ corresponds to the total energy per cell in our 0.5
approach, (. . .) is the canonical ensemble average, and kB is
the Boltzmann constant. Such quantities are shown in Table I
for several electrolyte concentrations. For comparison, we list 0
0 1 2 3 4 5 6 7
our results next to values obtained by the standard Ewald r/a
summation technique with conducting boundary conditions,
and without taking spherical averages over all orientations of FIG. 4. (Color online) Radial distribution functions correspond-
the infinite periodic lattice. From Table I, it is evident that ing to a 1:1, 2-M bulk electrolyte in the restricted primitive model.
the thermodynamic quantities calculated using our proposed Circles correspond to our approach and squares to Ewald sums.

016707-4
COULOMB INTERACTIONS IN CHARGED FLUIDS PHYSICAL REVIEW E 84, 016707 (2011)

TABLE II. Mean electrostatic potential at the surface (ψ0 ) and 10

at the Helmholtz plane (ψH P ) around a charged macroparticle of
9
valence zM , in presence of a 1:1 electrolyte at a concentration
0.01 M. 8
7
This work Ewald sums This work Ewald sums
ψ0 (mV) ψ0 (mV) ψH P (mV) ψH P (mV) 6

gMj(r)
zM (±2) (±2) (±2) (±2) 5

4 54.90 54.91 42.14 42.15 4

8 104.69 104.78 79.21 79.3 3
12 145.21 145.67 107.11 107.56 2
16 176.60 177.42 126.04 125.82
20 202.59 203.08 139.65 140.18 1
24 225.70 225.22 150.56 150.03 0
28 243.59 242.98 156.39 155.03 0 1 2 3 4 5 6 7
32 262.47 261.52 163.17 162.54 r/a

FIG. 5. (Color online) Radial distribution functions correspond-

ing to a 1:1, 2-M electrolyte in the restricted primitive model around a
charged macroparticle. The valence and diameter of the macroparticle
are D = 20 Å and zM = 32, respectively. Circles correspond to our
These results are summarized in Tables II and III, which approach and squares to Ewald sums.
exhibit similar values between the two approaches. The RDFs
corresponding to a macroparticle of valence zM = 32 are
IV. CONCLUSIONS
displayed in Fig. 5 for a 2-M monovalent salt, showing how
the two different approximations have a similar behavior. In The results in the last section show that the effective
addition, the values of ψ0 and ψH P as a function of the total potential we present in this paper is suitable to application
number N of particles inside the simulation box for the same in simulations of ionic fluids. A major limitation of our result
macroparticle are plotted in Fig. 6(a). Here, it is seen that our is that it is valid only for three-dimensional systems. One
approach and the standard Ewald sums method display very could easily derive an equivalent expression of Eq. (7) in two
similar values of the mean electrostatic potential and the same
asymptotic behavior when the number of charged particles 150
increases. Although the proposed method is O(N 2 ) (as the
standard Ewald summation method) its simplicity foresees a (a)
better performance in terms of computational speed. This is 100
ψ[mV]

observed in Fig. 6(b), where the number of MC cycles per

second performed as a function of the total number of charged
particles are collated. For the number of particles reported, 50
it is found that our approach can be from two to five times
faster than the standard Ewald summation under analogous 0
conditions.
MC cycles/second

100 (b)
10

TABLE III. Mean electrostatic potential at the surface (ψ0 ) and 1

at the Helmholtz plane (ψH P ) around a charged macroparticle of
0.1
valence zM , in presence of a 1:1 electrolyte at a concentration 2 M.
0 1000 2000 3000 4000 5000 6000
This work Ewald sums This work Ewald sums N
ψ0 (mV) ψ0 (mV) ψH P (mV) ψH P (mV)
FIG. 6. (Color online) (a) Mean electrostatic potential at the
zM (±1) (±1) (±1) (±1)
macroparticle’s surface (circles) and at the Helmholtz plane (squares)
4 17.09 17.07 4.42 4.40 around a macroparticle as a function of the total number N of charged
8 33.87 34.02 8.54 8.60 particles inside the simulation box. The macroparticle of valence and
12 50.50 51.04 12.51 13.05 diameter D = 20 Å and zM = 32, respectively, is immersed in a 1:1,
16 67.27 67.58 16.63 16.94 2-M electrolyte in the restricted primitive model. The empty and filled
20 84.09 84.37 20.80 21.08 symbols correspond to our proposed approach and Ewald summation
24 100.66 101.01 24.76 25.10 method, respectively. (b) Number of Monte Carlo (MC) cycles per
28 117.38 117.74 28.86 29.22 second as a function of the total number N of charged particles inside
32 133.86 134.34 32.74 33.23 the simulation box. The diamonds and the triangles correspond to our
proposed approach and Ewald summation method, respectively.

016707-5
GRAZIANO VERNIZZI et al. PHYSICAL REVIEW E 84, 016707 (2011)

dimensions, but only for a logarithmic potential. In general, where erfc(x) is the complementary error function:
 ∞
our arguments can be repeated for long-range potentials in 2
dt e−t ,
2

d dimensions of the type 1/r d−2 . We also emphasize that erfc(x) = √ (A4)
π x
the dipolar term Ud in the Ewald summation formula Eq. (2)
we obtain
must be included and cannot be dropped in order to obtain a  qi qj  α α  2
dt e−D t − √
2 2
perfect match with Eq. (7). This is interesting because several U = √ q
works [19,21] discuss how the dipolar term can be dropped ,i,j
π 0 π i i
according to what kind of dielectric boundary conditions 
 qi qj
one imposes, while our formula has been derived without + erfc(αD). (A5)
making any particular choice about the boundary conditions. ,i,j
2D
Furthermore, we notice that in contrast with previous similar
methods [14], ours does not rely on the homogeneity of In Eq. (A5) we added and subtracted the term with D = 0
the system within the main cell, which makes it suitable to in the first sum. That allows the application of the Poisson
study anisotropic systems frequently occurring in molecular summation formula:
 1  ix·Q ˆ
biology and fluid state physics. As a final remark, Fig. 2(b) f (x + ) = e f (Q), (A6)
suggests that the effective energy in Eq. (7) is amenable 
V Q
to local computational algorithms, meaning that only the
where the reciprocal lattice is defined by Q ·  = 2π m, m ∈
particles within the simulation box are sufficient for computing
Z, V is the volume of the primitive cell, and fˆ(k) is the Fourier
the total electrostatic energy of the rotational average of an
transform of f (x): 
infinite periodic lattice, despite the long-range nature of the
interactions. Computational-intensive resources for evaluat- fˆ(k) ≡ d 3 x e−ik·x f (x). (A7)
ing special functions, or reciprocal lattices, are no longer R3
necessary. Although our O(N 2 ) method is outperformed We have
by, for instance, a O(N log N ) optimized Ewald summation  qi qj  α
dt e−D t
2 2
scheme, its simplicity should be particularly appealing for √
,i,j
π 0
implementations with parallel computing at an intermediate
  α
number of charges [17]. 1  i(bi −bj )·Q −iQ·x
dt e−x t
2 2
= √ qi qj e 3
d xe
V π Q,i,j R3 0
ACKNOWLEDGMENTS 
π  α
1
qi qj ei(bi −bj )·Q dt 3 e−Q /4t
2 2
G.V. and M.O.C. acknowledge the support of NSF Grant =
V Q,i,j 0 t
No. DMR-0907781, and G.I.G.-G. acknowledges the support
of NSF Grant No. DMR-05020513. 2π  e−Q /4α
2 2

= q̂(Q)q̂(−Q) , (A8)
V Q Q2
APPENDIX: ELEMENTARY DERIVATION OF THE EWALD
SUMMATION FORMULA where 
The energy per cell of a infinite periodic ionic lattice is q̂(Q) ≡ qj e−ibj ·Q . (A9)
j
defined as
1   qi qj In the first line of Eq. (A8) we used the Poisson summation
U= , D ≡ |bi − bj − |. (A1) formula, then in the second line we evaluated the Gaussian
2 ,i,j D
integral over dx3 , and finally we used the definition Eq. (A9).
It is convenient to introduce the following integral representa- The final expression for U is
2π  e−Q /4α α  2
2 2
tion of the Coulomb potential:
U = q̂(Q)q̂(−Q) −√ q
 ∞ V Q Q 2 π i i
1 2
dt e−D t .
2 2
=√ (A2)   qi qj
D π 0 + erfc(Dα). (A10)
2D
The large-D behavior of Eq. (A2) is controlled by the integra- ,i,j
tion over small-t values, and the small-D behavior is controlled The term Q = 0 in the first sum requires special care. In fact,
by the integration over large-t values. Such an observation is the asymptotic expansion at small Q of the argument in the
relevant for the application of the Ewald method [2], which sum is
basically consists of decomposing the small-distance behavior  2
e−Q /4α
2 2
i qi
of the integration domain from the large-distance behavior, q̂(Q)q̂(−Q) 2
∼
and afterwards by tuning the movable boundary separating Q Q2
   2
m,n qm qn i(bm − bn ) · Q
the two domains to achieve equal convergence rates for the
i qi
two sums. Namely, by decomposing the integration region + −
∞ α ∞ Q 4α 2
0 = 0 + α , and by using the integral
√   
 ∞ 1 Q 2
dt e −D 2 t 2
=
π
erfc(αD), (A3) + qm qn (bm − bn ) · + O(Q). (A11)
2D 2 m,n Q
α

016707-6
COULOMB INTERACTIONS IN CHARGED FLUIDS PHYSICAL REVIEW E 84, 016707 (2011)

2 Eq. (A10) as
The divergence of order O(1/Q ) is removed by the elec-
troneutrality condition i qi = 0, which also expunges the
2π  e−Q /4α α  2
2 2

divergence of order O(1/Q) and the first term of the O(1) U = q̂(Q)q̂(−Q) − √ q
constant term. The remaining term is anisotropic. We therefore V Q=0 Q2 π i i
take its spherical average (for a more rigorous mathematical  2
derivation of the same result, see [11])   qi qj 2π 
+ erfc(Dα) + qm bm , (A13)
 π  2π
1 1 ,i,j
2D 3V m
d cos θ dφ qm qn (bm − bn )2 cos2 θ
4π 0 0 2 m,n which is the final expression for the Ewald sum we used in this
1 1 paper. Equation (A13) is independent from α. However, for
= qm qn (bm − bn )2 = − qm qn bm · bn numerical estimates the infinite sums must be truncated, and
6 m,n
3 m,n such a truncation yields a function U (α), which is α dependent.
 2 However, there is an interval of values for α where U is roughly
1  constant. In our numerical utilizations of the Ewald method, we
=− qm bm , (A12)
3 m implemented the feature for α to be automatically selected in
the region where ∂U (α)/∂α = 0 [12]. Once α has been fixed,
which shows that the total dipole moment of the cell the Ewald method provides an efficient way to estimate the
contributes to the Ewald sum. We can then write total electrostatic energy per cell of the periodic ionic lattice.

[1] M. Karttunen, J. Rottler, I. Vattulainen, and C. Sagui, Comp. [14] E. Yakub and C. Ronchi, J. Chem. Phys. 119, 11556 (2003);
Model. Membr. Bil. 60, 49 (2008). J. Low Temp. Phys. 139, 633 (2005); E. Yakub, J. Phys. A 39,
[2] P. P. Ewald, Ann. Phys. (Leipzig) 64, 253 (1921). 4643 (2006).
[3] C. Kittel, Introduction to Solid State Physics, 8th ed. (Wiley, [15] X. Wu and B. R. Brooks, J. Chem. Phys. 122, 044107 (2005);
New York, 2004). 129, 154115 (2008).
[4] E. Madelung, Phys. Zs. 19, 524 (1918). [16] E. Yakub, C. Ronchi, and D. Staicu, J. Chem. Phys. 127, 094508
[5] C. Sagui and T. A. Darden, Annu. Rev. Biophys. Biomol. Struct. (2007); J. Nucl. Mater. 389, 119 (2009).
28, 155 (1999). [17] P. K. Jha, R. Sknepnek, G. I. Guerrero-Garcı́a, and M. Olvera de
[6] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and la Cruz, J. Chem. Theor. Comput. 6, 3058 (2010).
L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). [18] D. Jackson, Classical Electrodynamics, 3rd ed. (John Wiley &
[7] C. Sagui and T. Darden, J. Chem. Phys. 114, 6578 (2001). Sons, New York, 1998).
[8] A. C. Maggs and V. Rossetto, Phys. Rev. Lett. 88, 196402 [19] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids
(2002). (Clarendon, Oxford, 1989).
[9] M. M. Reif, V. Kräutler, M. A. Kastenholz, X. Daura, and P. H. [20] C. S. Hoskins and E. R. Smith, Mol. Phys. 41, 243 (1980).
Hünenberger, J. Phys. Chem. B 113, 3112 (2009). [21] D. Frenkel and B. Smit, Understanding Molecular Simulation
[10] E. V. Kholopov, Phys. Usp. 47, 965 (2004). (Academic, San Diego, 2002).
[11] S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. [22] E. González-Tovar and M. Lozada Cassou, J. Phys. Chem. 93,
London, Ser. A 373, 27 (1980). 3761 (1989).
[12] G. Hummer, Chem. Phys. Lett. 235, 297 (1995). [23] G. I. Guerrero-Garcı́a, E. González-Tovar, M. Lozada Cassou,
[13] D. J. Adams and G. S. Dubey, J. Comput. Phys. 72, 156 and F. de J. Guevara-Rodrı́guez, J. Chem. Phys. 123, 034703
(1987). (2005).

016707-7