Ewald summation for systems with slab geometry
In-Chul Yeh and Max L. Berkowitza)
Department of Chemistry, CB# 3290, University of North Carolina, Chapel Hill, North Carolina 27599
͑Received 25 March 1999; accepted 21 May 1999͒
We propose a modiﬁcation in the three-dimensional Ewald summation technique for calculations of
long-range Coulombic forces for systems with a slab geometry that are periodic in two dimensions
and have a ﬁnite length in the third dimension. The proposed method adds a correction term to the
standard Ewald summation formula. To test the current method, molecular dynamics simulations on
water between Pt͑111͒ walls have been carried out. For a more direct test, the calculation of the pair
forces between two point charges has been also performed. An excellent agreement with the results
from simulations using the rigorous two dimensional Ewald summation technique were obtained.
We observed that a signiﬁcant reduction in computing time can be achieved when the proposed
modiﬁcation is used. © 1999 American Institute of Physics. ͓S0021-9606͑99͒70331-4͔
I. INTRODUCTION
An accurate treatment of long-range Coulombic interactions
in computer simulations of charged particles is of great
importance. For a typical three-dimensionally periodic system,
the Ewald summation technique1,2
is the most widely
used and accepted method for this purpose. Extensive optimization
techniques such as the smooth particle mesh Ewald
͑PME͒ method3,4
have been developed to perform fast and
reliable simulations of large systems using Ewald summa-
tion.
For a slab geometry which occurs frequently in surface
and interfacial systems, the conventional Ewald summation
technique cannot be used directly since there is no periodicity
in the one of three dimensions ͑let us say along the z
direction͒. For this type of two-dimensionally periodic ͑2DP͒
systems which have a ﬁnite length along the third dimension,
various methods have been proposed for the treatment of
long-range forces.5–17
Recently, some of the methods have
been compared to each other in terms of computational speed
and accuracy.16,17
A two-dimensional Ewald summation
͑EW2D͒ technique ﬁrst introduced by Parry5
and later independently
rederived by Heyes, Barber, and Clarke6
and by de
Leeuw and Perram7
is found to be the most accurate. But the
direct use of EW2D formula is known to be computationally
very expensive.8,16,17
A simple solution to this problem is to
use a precalculated table.8
The EW2D method with a precalculated
table has been successfully applied in our recent
simulation studies of water next to metal surfaces.18–20
Another
widely used approach is to apply the conventional
three-dimensional Ewald summation ͑EW3D͒ technique to a
simulation cell elongated in z direction so that a sufﬁciently
large empty space between periodic replicas in z direction is
created.21
The inclusion of empty space into the unit cell is
done to avoid an artiﬁcial inﬂuence from the periodic images
in z direction. This method was applied to simulations of
various interfacial systems.22–24
Shelley and Patey21
studied the effect of boundary conditions
for water conﬁned between planar walls. They compared
various boundary conditions such as minimum image,
spherical and cylindrical cutoffs with the Ewald summation
method. The Ewald summation was found to be the safest
choice in the calculation. Systems considered in their study
were symmetric and, therefore, the possible effect of asymmetry
and net polarization of the system on the Ewald summation
was not considered. Spohr8
compared the results
from the simulations that used EW3D method with those
from the EW2D method. He concluded that results for
EW3D converge to those for EW2D when the length of the
simulation cell in z direction (Lz) used in simulations with
EW3D was large. At the same time he noticed that even
when Lz was ﬁve times larger than the length of the simulation
cell in x or y directions ͑Lx or Ly , respectively͒, the
convergence was not satisfactory.
Recently, there has been a renewed interest in the treatment
of long-range forces in polar systems or systems carrying
a net charge.25–29
A better understanding of conditions
under which Ewald summation is performed followed from
this work. Still one issue that may be important even for
some systems that are charged neutral, but have a total dipole
moment, remains not clear. The issue is the use of the surface
term in the EW3D formula, the term that depends on the
shape of the system, its total dipole moment and the dielectric
constant of the surrounding medium.1,28–31
The surface
term is due to the fact that the electrostatic energy of the
ionic crystal is composed of two parts. First part is shapeindependent
but depends on the structure of the crystal lattice
concerned and the distribution of ions within a unit cell.
Second part depends on the shape of the piece of the crystal
and the dipole moment of a unit cell. The expression for the
shape-independent part of energy has the same mathematical
form as in the regular EW3D formula. Smith32
derived an
expression for the shape-dependent part suitable for the slab
geometry, which is of particular interest here.
In this study we propose that an EW3D method which
includes the shape-dependent correction term introduced by
Smith32
can be used for the calculation of the long-range
electrostatic forces in simulations of 2DP systems. To distin-a͒
Electronic mail: maxb@net.chem.unc.edu
JOURNAL OF CHEMICAL PHYSICS VOLUME 111, NUMBER 7 15 AUGUST 1999
31550021-9606/99/111(7)/3155/8/$15.00 © 1999 American Institute of Physics
Copyright ©2001. All Rights Reserved.
guish it from the regular EW3D method, it will be abbreviated
as EW3DC ͑the three-dimensional Ewald summation
with the correction term͒. This term is similar to the surface
term in the standard EW3D formula but has a different proportionality
factor and involves only the z component of the
total dipole moment. As far as we know, the only attempt to
include this correction term was done in a simulation on
water between two ideal classical walls performed by Hautman
et al.10
No detailed comparison with the rigorous
EW2D has been made in that work.
To test the EW3DC method we performed simulations
on water next to the charged Pt͑111͒ walls. Recently, we
used this system to calculate the dielectric constant of water
at high electric ﬁelds.19
The long-ranged forces were treated
using the EW2D method in these simulations. Below we
compare the results from the simulations performed with the
EW3DC method to those obtained with the EW2D method.
For more detailed comparison, test calculations involving
only two point charges also have been carried out. By adding
the extra correction term in the EW3D formula, both accuracy
and speed of the calculation are shown to be greatly
improved.
II. ELECTROSTATIC INTERACTIONS
Consider a system of N point charges qi at positions ri
which satisfy the charge neutrality condition ⌺i
N
qiϭ0 in a
rectangular simulation cell with lengths in x, y, and z directions
of Lx , Ly , and Lz , respectively.
The electrostatic energy can be written as
Uϭ
1
2 ͚n
Ј͚iϭ1
N
͚jϭ1
N
qiqj
͉rijϩn͉
, ͑1͒
where the sum over n is the sum over all lattice points, n
ϭ(nxLx ,nyLy ,nzLz) with nx ,ny ,nz integers. The primed
sum indicates the omission of the iϭj term when nϭ0. The
factor 1/(4␲⑀0) is omitted for simplicity ͑⑀0 is the vacuum
permittivity͒. This sum is conditionally convergent, which
means that the result depends on the order or the summation
geometry in which we add up the terms. The outcome is also
dependent on the dielectric constant of the surrounding medium
(⑀s). The electrostatic energy calculated by the Ewald
sum is given by1,2,32
Uϭ
1
2 ͚iϭ1
N
͚jϭ1
N
͚͉n͉ϭ0
ϱ
Јqiqj
erfc͑␣͉rijϩn͉͒
͉rijϩn͉
ϩ
1
2␲V ͚iϭ1
N
͚jϭ1
N
͚k 0
qiqjͩ4␲2
k2 ͪexpͩϪ
k2
4␣ͪ
ϫcos͑k•rij͒Ϫ
␣
ͱ␲
͚iϭ1
N
qi
2
ϩJ͑M,P͒. ͑2͒
J(M,P) is a shape-dependent term and depends on the
summation geometry ͑P͒.32
M is the total dipole moment of
the unit simulation cell and is given by ⌺iϭ1
N
qiri . V is the
volume of the unit simulation cell given by LxϫLyϫLz and
k is a reciprocal lattice vector given by
((2␲nxЈ)/Lx ,(2␲nyЈ)/Ly ,(2␲nzЈ)/Lz) with nxЈ ,nyЈ ,nzЈ integers.
In the calculations ␣ and number of n and k vectors are
adjustable parameters and are typically chosen for the optimum
computational efﬁciency.
When the spherical geometry is used for summation (P
ϭS), J(M,S) is given by the following relation:2,25,32
J͑M,S͒ϭ
2␲
͑2⑀sϩ1͒V
͉M͉2
. ͑3͒
If the surrounding medium has an inﬁnite dielectric constant
͑⑀sϭϱ, metal͒, J(M,S) vanishes. This boundary condition
is commonly called conducting ͑‘‘tinfoil’’͒ boundary
condition. If ⑀s has a ﬁnite value, this term cannot be simply
ignored. In particular, if there is no surrounding medium
͑⑀sϭ1, vacuum͒, this term becomes
J͑M,S͒ϭ
2␲
3V
͉M͉2
. ͑4͒
The contribution to the force from this term is
FiϭϪٌiJ͑M,S͒ϭϪ
4␲qi
3V
M. ͑5͒
Clearly, this term is not necessarily zero especially when
there is a net polarization in the simulation cell (M 0). It
can be ignored if the system is isotropic, Mϭ0. Ignoring this
term is equivalent to adopting the tinfoil boundary condition.
Implications of adopting tinfoil boundary condition in place
of vacuum boundary condition for polar systems have been
discussed by Roberts and Schnitker.30,31
The total Coulomb energy for a system periodic in two
dimensions and ﬁnite in the third dimension according to the
EW2D method is
Uϭ
1
2 ͚iϭ1
N
͚jϭ1
N
͚͉m͉ϭ0
ϱ
Јqiqj
erfc͑␣͉rijϩm͉͒
͉rijϩm͉
ϩ
␲
2A ͚iϭ1
N
͚jϭ1
N
͚h 0
qiqj
cos͑h•rij͒
h
ϫͭexp͑hzij͒erfcͩ␣zijϩ
h
2␣ͪϩexp͑Ϫhzij͒erfc
ͩϪ␣zijϩ
h
2␣ͪͮϪ
␲
A ͚iϭ1
N
͚jϭ1
N
qiqjͭzij erf͑␣zij͒
ϩ
1
␣ͱ␲
exp͑Ϫ␣2
zij
2
͒ͮϪͩ ␣
ͱ␲
͚ͪiϭ1
N
qi
2
, ͑6͒
where m is a lattice vector for the 2DP system and is given
by (mxLx ,myLy,0) with mx ,my integers. h is a reciprocal
lattice vector for the 2DP system and is given by
((2␲mxЈ)/Lx ,(2␲myЈ)/Ly,0) with mxЈ ,myЈ integers. A is the
area of the unit cell in x and y directions given by LxϫLy .
The primed sum again indicates the omission of the iϭj
term when mϭ0.
If two particles with charges qi and qj are sufﬁciently
apart in z direction and are periodically repeated in x and y
directions, the electric ﬁeld acting on the particle i due to
3156 J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 I.-C. Yeh and M. L. Berkowitz
Copyright ©2001. All Rights Reserved.
particle j can be considered as that due to the uniformly
charged sheet with a surface charge density ␴ϭqj /(LxLy)
and is given by8,33
Ejϭ
␴
2⑀0
ϭ
qj
2⑀0LxLy
͑7͒
The x, y, and z components of the force are given by the
following expression:8
FxϭFyϭ0, ͑8͒
Fzϭ
qiqj
2⑀0LxLy
. ͑9͒
By using the EW2D formula, Spohr8
found that the
above ‘‘parallel plate capacitor’’ approximation can be used
with a sufﬁcient accuracy ͑relative error less than 0.001͒ if
the z separation distance between particles i and j (͉zij͉) is
larger than Lx or Ly . The important feature of this approximation
is that no ͉zij͉ distance dependency appears in the
formulas. Therefore, if there is an empty space in the EW3D
simulation cell along the z direction with a length of at least
Lx or Ly , contributions to the forces from the image cells in
z direction should be negligible due to the neutrality of the
system (⌺j
N
Ejϭ(1/2⑀0LxLy)⌺j
N
qjϭ0). Nevertheless, the
comparison between the results from simulations using
EW3D and EW2D methods shows that the EW3D method
produces unsatisfactory results even when the length of the
empty space in z direction is signiﬁcantly larger than Lx or
Ly .8
This apparent contradiction suggests that the Ewald
summation using tinfoil boundary condition and spherical
summation geometry which is commonly used in the Ewald
summation may not be suitable for the calculation in the slab
geometry.
Smith32
showed that if a geometry of a rectangular plate
(PϭR) or disk which are inﬁnitely thin in z direction are
used as a summation geometry in the Ewald summation, then
the shape-dependent term J(M,R) is given by
J͑M,R͒ϭ
2␲
V
Mz
2
, ͑10͒
where Mz is the z component of the total dipole moment of
the simulation cell. This is equivalent to adopting a planewise
sum method34
͑summing x and y directions ﬁrst then
progressing in z direction͒. The contribution to the force
from this term is
Fx,iϭFy,iϭ0, ͑11͒
Fz,iϭϪ
‫ץ‬J͑M,R͒
‫ץ‬zi
ϭϪ
4␲qi
V
Mz . ͑12͒
Note that Eqs. ͑5͒ and ͑12͒ differ by a factor of 3 and
that only the z component of the dipole moment appears in
Eq. ͑12͒. The planewise sum combined with a sufﬁcient
empty space in z direction enables us effectively utilize the
parallel plate capacitor approximation to eliminate contributions
from image cells in z direction which are unwanted in
2DP systems. Therefore, we can use the regular EW3D
method but with the correction term given by the Eq. ͑12͒ to
calculate long-range Coulombic forces for systems which are
periodic in two dimensions and are ﬁnite in the third. This
constitutes the EW3DC method which as we show in what
follows is indeed capable of reproducing the EW2D results.
III. METHODOLOGY
For simulations of water between Pt͑111͒ walls, we use
the same interaction potentials that we used in our previous
studies,18,19,35
namely, the extended simple point charge
͑SPC/E͒ model of water36
and water–Pt͑111͒ surface
potential.37
A rectangular prism is chosen for the unit cell of the
simulations. Lx and Ly of 22.175 and 19.204 Å, respectively,
have been used to satisfy the periodicity requirement of
Pt͑111͒ surface. Simulations are performed for 512 water
molecules in the unit simulation cell. Water molecules are
conﬁned in z direction by Pt͑111͒ walls and the density of
water at the center of the simulation cell is 1.0 g/cc. This is
achieved by changing the value of a parameter zm in water–
Pt͑111͒ potential.37
External electric ﬁelds (E0) of 0–4 V/Å
are applied along the z direction.
For the calculation of long-range Coulombic interactions
we have used EW2D, EW3D, and EW3DC methods. For
EW3D and EW3DC calculations, three-dimensional periodic
boundary conditions were applied. Lz was chosen to be 90 Å
which is approximately four times larger than Lx . Ewald
parameter ␣ of 0.3 ÅϪ1
was used. As customary in the
EW3D method, only minimum image (͉n͉ϭ0) was considered
with the real space cut-off radius of 9.5 Å. The reciprocal
space cut-off radius of 1.75 ÅϪ1
has been used.
Since the reciprocal space term in the EW2D is computationally
very expensive, it is advantageous to have less
terms for the reciprocal space part and more terms for the
real space part. Therefore, for the EW2D calculation, the
parameter ␣ was chosen to be 0.104 ÅϪ1
and the real space
terms which satisfy the condition mx
2
ϩmy
2
р4 and the reciprocal
space terms which satisfy mxЈ2
ϩmyЈ2
р9 have been
considered. The direct use of the EW2D formula is so time
consuming that a precalculated table of potential energy,
forces and second derivatives on a three-dimensional grid
with the size of 0.2ϫ0.2ϫ0.2 Å3
has been constructed as
suggested by Spohr.8
The EW2D calculations have been performed
by interpolation of the table. For smaller distances
(rijϽ3.3 Å), the exact EW2D formula has been used. For
͉zij͉ϾLy , the parallel plate capacitor approximation given in
Eqs. ͑8͒ and ͑9͒ has been utilized.
The Verlet algorithm has been used to propagate the
trajectories with a time step of 2.5 fs and SHAKE routine has
been used to preserve the rigidity of the water molecules.1
The coordinates were saved at every 10 time steps for further
analysis. The temperature of the run was kept at 300 K by
coupling the system to the heatbath using the algorithm of
Berendsen et al.38
200 ps run after at least 50 ps equilibration
has been used for the analysis.
In our previous work19
where we determined the dielectric
constant of water as a function of the electric ﬁeld we
performed simulations on water lamina embedded between
two surfaces. To ﬁnd the value of the dielectric constant we
needed to calculate the total electric ﬁeld inside the lamina.
3157J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 Ewald summation for a slab
Copyright ©2001. All Rights Reserved.
This total electric ﬁeld ͓E(z)͔ was obtained using the following
relation:
E͑z͒ϭE0ϩEp͑z͒, ͑13͒
where Ep(z) is the electric ﬁeld due to orientational polarization
of water sample and is given by the following
expression:39,40
Ep͑z͒ϭ
͐ϪLz/2
z
␳q͑zЈ͒dzЈ
⑀0
, ͑14͒
where ␳q(z) is a charge density distribution.
The relative permitivity ͑dielectric constant ⑀͒ of water
in the external electric ﬁeld E0 and the total electric ﬁeld E
can be obtained from the following equation:
⑀ϭ
E0
E
. ͑15͒
To ﬁnd the permittivity we calculated the total ﬁeld in
the vicinity of the center of the slab. In previous simulation
studies on water–Pt interface it was found that the inﬂuence
of the walls on the orientational distributions of water is
almost absent for water about 10 Å away from the closest
approach of the water molecule to the Pt wall.18,35
Therefore,
it is expected that water–wall interaction has no inﬂuence on
the estimate of E when using the current method. To estimate
the total electric ﬁeld E around the center of the sample cell,
10 Å interval averages have been taken across the sample
cell in the z direction. Then to eliminate any local ﬂuctuation
of the electric ﬁeld around the center, further 5 Å interval
averages were taken. The resulting average electric ﬁeld at
zϭ0 ͑the center of the simulation cell͒ is taken to be the total
electric ﬁeld E that enters the equation determining the dielectric
constant.
The electrostatic potential ͑␾͒ has been calculated by the
following equation:
FIG. 1. Comparison of simulation results for water between Pt͑111͒ walls at
zero external electric ﬁeld obtained by using different methods to treat the
long-range forces. Top: Oxygen density distributions. Middle: Charge density
distributions. Bottom: Electrostatic potential.
FIG. 2. Comparison of simulation results for water between Pt͑111͒ walls at
E0 of 1.0 V/Å calculated by different Ewald methods. Top: Oxygen density
distributions. Middle: Charge density distributions. Bottom: Electric ﬁeld Ep
calculated by Eq. ͑14͒.
FIG. 3. Comparison of simulation results for water between Pt͑111͒ walls at
E0 of 2.0 V/Å calculated by different Ewald methods. Top: Oxygen density
distributions. Middle: Charge density distributions. Bottom: Electric ﬁeld Ep
calculated by Eq. ͑14͒.
3158 J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 I.-C. Yeh and M. L. Berkowitz
Copyright ©2001. All Rights Reserved.
␾͑z͒ϭϪ
1
⑀0
͵ϪLz/2
z
␳q͑zЈ͒͑zϪzЈ͒dzЈ, ͑16͒
where ␳q is the charge density.
IV. RESULTS AND DISCUSSION
A. Water–Pt„111… interface
Oxygen density proﬁles, charge density distributions,
and electrostatic potential proﬁles for water between Pt͑111͒
walls at zero external electric ﬁeld calculated with different
methods to treat the long-ranged Coulombic forces ͑EW2D,
EW3D, EW3DC, and the spherical cutoff͒ are compared in
Fig. 1. For the spherical cutoff method, all the interactions
with the distance larger than the cutoff radius of 9.5 Å are
ignored. Oxygen density proﬁles and charge density distributions
are almost identical regardless of methodology. But the
value of the electrostatic potential at the center of the lamina
calculated with spherical cutoff is signiﬁcantly larger than
when other methods are used. This is consistent with the
results from previous simulation studies.8,24
This result
shows that the spherical cutoff should be avoided even for a
symmetric system. Therefore the spherical cutoff results are
not considered any longer in the following comparisons. The
electrostatic potential proﬁles calculated with EW2D is in
good agreement with those obtained by EW3D and EW3DC.
The correction term has very little inﬂuence on the calculation
results since the system is symmetric in z direction and
there is no net dipole moment in z direction (Mzϭ0). From
Fig. 1 we also notice that for the value of Lzϭ90 Å an empty
spacing with a length in z direction of at least 56 Å exists in
the system, which is about 2.5 Lx .
Figure 2 shows oxygen density proﬁles, charge density
distributions, electric ﬁeld proﬁles for water between Pt͑111͒
walls at E0ϭ1 V/Å when we use EW2D, EW3D, and
EW3DC methods. For the EW3D calculations, two different
values of Lz ͑90 and 180 Å͒ have been used to see the effect
of the value of Lz . Oxygen density proﬁles show very little
difference. The charge density distribution in the water
lamina obtained from simulations with the EW3D method
are signiﬁcantly different from those obtained by EW2D and
EW3DC methods. A good agreement between the results by
EW2D and EW3DC methods has been found. Due to the
difference in the charge density distributions obtained from
simulations that use EW2D and EW3D methods, the corresponding
Ep(z) distributions are also different. From the ﬁgure
we notice that around bulk region (Ϫ5 ÅϽzϽ5 Å), the
ﬁeld EЈ, which is the sum of two components EЈϭEp
ϩE0 , has a negative value when EW3D method is used.
While in simulations of the lamina using EW2D method EЈ
would be the total ﬁeld, since it is negative in simulations
with EW3D it indicates that overcompensation of the external
ﬁeld E0 occurred in this case. This type of overcompensation
of the electric ﬁeld in thermodynamically stable system
is highly unlikely and it is a consequence of an improper
treatment of electrostatic forces using the EW3D method
without considering the slab geometry. Imposing tinfoil
boundary condition and thus neglecting the correction term
seem to introduce an extra external electric ﬁeld contribution
in addition to the initially assigned value of E0 . As a matter
of fact, using Eq. ͑12͒ the external ﬁeld E0Ј that molecules of
water experience in simulations which employ the EW3D
method can be estimated by the following relation:
E0ЈϭE0ϩ
4␲͗Mz͘
V
, ͑17͒
where ͗Mz͘ is the average of the z-component of the dipole
moment of the sample. The second term in this equation is
the contribution due to images of the lamina from the central
cell. Figure 3 shows results similar to the one from Fig. 2,
but for E0 of 2 V/Å. Again, the result from EW3D without
the correction term deviates signiﬁcantly from that of
EW2D, while the result from EW3DC agrees well with that
from EW2D.
Simulations at E0 of 3 and 4 V/Å have been also carried
out. In our recent simulations on water–Pt͑111͒ interface, a
phase transition to proton ordered cubic ice was observed
when E0 was somewhere between 3 and 4 V/Å with the
EW2D method18,19
and between 2 and 3 V/Å with the EW3D
method.35
In this study, a phase transition to a proton ordered
cubic ice has been observed in a region between 3 and 4 V/Å
with EW3DC. This again conﬁrms that water feels stronger
‘‘effective’’ external electric ﬁeld when EW3D is used, resulting
in a lower critical E0 for a phase transition. Table I
presents a summary of the results for the calculations done
on water lamina in the external electric ﬁeld when EW3D is
used. This table shows that the data are consistent with the
previously obtained results. Thus, for example, in calculations
with EW3D when the external ﬁeld E0 is 3 V/Å the
effective external ﬁeld E0Ј is actually ϳ4.15 V/Å and the
total ﬁeld E is 1.21 V/Å. As our previous simulations
showed,19
for water to undergo the transition, the total ﬁeld
is supposed to be somewhat smaller than 1 V/Å. This is why
TABLE I. The ‘‘effective’’ external electric ﬁeld E0Ј and total electric ﬁeld E for water between Pt͑111͒ walls
calculated by the EW3D method. Units of ͗Ep͘ and (4␲͗Mz͘)/V are V/Å. Errors in parentheses are estimated
from standard deviations of the results of four consecutive 50 ps runs.
E0 ͑V/Å͒ Lz ͑Å͒ ͗Ep͘ (4␲͗Mz͘)/V E0ЈϭE0ϩ(4␲͗Mz͘)/V EϭE0Јϩ͗Ep͘
1 180 Ϫ1.2060͑Ϯ0.0006͒ 0.2257͑Ϯ0.0004͒ 1.2257͑Ϯ0.0004͒ 0.020͑Ϯ0.001͒
1 90 Ϫ1.5535͑Ϯ0.0054͒ 0.5860͑Ϯ0.0009͒ 1.5860͑Ϯ0.0009͒ 0.032͑Ϯ0.006͒
2 90 Ϫ2.7025͑Ϯ0.0032͒ 1.0511͑Ϯ0.0005͒ 3.0511͑Ϯ0.0005͒ 0.349͑Ϯ0.004͒
3 90 Ϫ2.9429͑Ϯ0.0008͒a
1.1535͑Ϯ0.0005͒ 4.1535͑Ϯ0.0005͒ 1.211͑Ϯ0.001͒
a
The average is taken over the 10.6 Å interval to take into account the periodic structure in z direction
introduced by the phase transition.
3159J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 Ewald summation for a slab
Copyright ©2001. All Rights Reserved.
the phase transition was observed in a ﬁeld with E0 below 3
V/Å when the EW3D method was used.
We also calculated ⑀ by Eq. ͑15͒ and average
cos ␪␮,E0
(͗cos ␪␮,E0
͘) where ␪␮,E0
is the angle between a water
dipole ␮ and the direction of the external electric ﬁeld
E0 . Values for ⑀ and ͗cos ␪␮,E0
͘ obtained by using EW2D
and EW3DC have been compared in Table II. An excellent
agreement has been observed.
B. Forces between two particles
Calculation of the pair forces is the most direct and accurate
measure to compare different methods for the treatment
of long-range forces.8,16,17
Spohr8
compared the pair
forces calculated by EW2D and EW3D. He attributed the
deviations of the EW3D results from the EW2D results to
the coupling between the periodic replicas of the interface.
The same geometry of the simulation cell used in
Spohr’s work has been used here ͑LxϭLyϭ18 Å and Lz
ϭ3ϫLx or 5ϫLx͒. ␣ of 0.45 ÅϪ1
and the reciprocal space
cut off radius of 5.0 ÅϪ1
have been used. The parameters
have been chosen for better accuracy because we want to
eliminate any unwanted error due to numerical artifacts
when we compare the results for the EW2D and EW3D. The
choice of parameters could be relaxed for better computational
efﬁciency. For the EW2D calculation, the parameter ␣
was chosen to be 0.111 ÅϪ1
and the real space terms which
satisfy the condition mx
2
ϩmy
2
р4 and the reciprocal space
terms which satisfy mxЈ2
ϩmyЈ2
р9 have been considered.
We calculated the pair forces between two oppositely
charged unit point charges at ͑0,0,0͒ and (0,0,z) using
EW2D, EW3D, and EW3DC methods. z changes from 0 to
about Lz/2 in our tests. This type of arrangement of charges
is also used in Spohr’s work8
and creates a dipole moment in
z direction. The results are shown in Fig. 4. The deviations of
the EW3D results without the correction term from the
EW2D results reported in Spohr’s work8
have been reproduced.
At the same time the results from the EW3DC are
almost identical to the EW2D results. It is to be noted that
there are empty spaces with lengths of at least 1.5ϫLx and
2.5ϫLx in z direction for Lz of 3ϫLx and 5ϫLx , respectively,
for the results shown in Fig. 4. This clearly shows that
to use EW3D in 2DP systems, the empty space need not be
huge but the correction term should be included.
Figure 5 shows the pair forces experienced by two oppositely
charged unit point charges at ͑0,0,0͒ and ͑4.5 Å, 4.5
Å, z͒. For EW3D and EW3DC calculations, Lz of 5ϫLx ͑90
Å͒ has been used. z again changes from 0 to Lz/2. This type
of arrangement of charges creates a dipole moment not only
in z direction but also in x and y directions. The z components
of forces are shown and an excellent agreement between
EW2D and EW3DC results is evident. Note that even
for the cases when Mx and My are not zero, the correction
term ͓Eq. ͑12͔͒ which involves only Mz is sufﬁcient.
Figure 6 considers a rather extreme case when there is a
FIG. 4. Comparison of the z component of the force acting on the unit point
charge at (0,0,z) by another oppositely charged unit point charge ﬁxed at
͑0,0,0͒ in two-dimensionally periodic systems calculated by different Ewald
methods. The logarithmic scale is used for better comparison with the result
in Spohr’s work ͑Ref. 8͒.
FIG. 5. Comparison of the z component of the force acting on the unit point
charge at ͑4.5 Å, 4.5 Å, z͒ by another oppositely charged unit point charge
ﬁxed at ͑0,0,0͒ in two-dimensionally periodic systems calculated by different
Ewald methods.
TABLE II. The dielectric constant ⑀ and ͗cos ␪␮,E0
͘ for water between
Pt͑111͒ walls. ␪␮,E0
is the angle between a water dipole ␮ and the direction
of the external electric ﬁeld E0. Errors in parentheses are estimated from
standard deviations of the results of four consecutive 50 ps runs.
E0 ͑V/Å͒ Method ⑀ ͗cos ␪␮,E0
͘
1 EW2D 59.9͑Ϯ11͒ 0.330͑Ϯ0.001͒
1 EW3DC 67.5͑Ϯ7.1͒ 0.329͑Ϯ0.001͒
2 EW2D 38.1͑Ϯ1.3͒ 0.664͑Ϯ0.002͒
2 EW3DC 37.1͑Ϯ1.7͒ 0.664͑Ϯ0.002͒
3 EW2D 9.2͑Ϯ0.1͒ 0.905͑Ϯ0.001͒
3 EW3DC 9.3͑Ϯ0.1͒ 0.907͑Ϯ0.002͒
FIG. 6. Comparison of the z component of the force acting on the unit point
charge at (0,0,z) by another oppositely charged unit point charge ﬁxed at
͑0,0,Ϫ44 Å͒ in two-dimensionally periodic systems calculated by different
Ewald methods.
3160 J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 I.-C. Yeh and M. L. Berkowitz
Copyright ©2001. All Rights Reserved.
very little empty space between periodic replicas in z direction.
Forces between two oppositely charged unit point
charges at ͑0,0,Ϫ44 Å͒ and (0,0,z) are calculated. z is
changed from Ϫ43 to 44 Å. Lz used in EW3D and EW3DC
is 90 Å, which means that one of the charges is ﬁxed near the
end of the simulation cell in z direction at ͑0,0,Ϫ44 Å͒. Figure
6 shows that the result from EW3D without the correction
term deviates from the EW2D result in almost the entire
range of z values. The result from EW3DC agrees well with
EW2D result for z values up to ϳ28 Å where interaction
with the adjacent charges in periodic replicas in z direction
are beginning to dominate the interaction force. At z
ϭ28 Å, an empty space with a length of Lx in z direction is
present in the box. This shows that to get good numerical
results, there should be an empty space at least with a length
of Lx or Ly even when EW3DC is used.
C. Speed and accuracy
Computational speeds of the direct use of the EW2D
formula, the EW2D method with a precalculated table and
the EW3DC calculations have been checked for water–
Pt͑111͒ system. Only interactions between water molecules
have been considered for more direct comparisons. Calculations
have been carried out on a Silicon Graphics Origin
2000 workstation. The geometry of the simulation cell and
the Ewald parameters are the same as the ones used in this
work for water–Pt͑111͒ systems. Since the actual computing
time may strongly depend on the number of particles, the
choice of Ewald parameters, the machine used, and the degree
of the optimization, the numbers given here should be
interpreted with caution.
Computational time per time step has been estimated for
each Ewald method and is summarized in Table III. Since
the value of ␣ is small, a reciprocal space term with h 0
͓the second term in Eq. ͑6͔͒ could be neglected with very
little error as suggested by the previous studies.6,16
This optimization
of the EW2D method results in a signiﬁcant reduction
of the computing time in the direct use of EW2D
formula as shown in Table III, but does not reduce the computing
time signiﬁcantly in the EW2D calculation with the
table. The inclusion of the correction term into the EW3D
formula ͑EW3DC͒ hardly introduces any difference in the
computing time compared with EW3D. The computing time
for the EW3DC is ϳ10 times faster than the time for EW2D
with the table.
We also compared the accuracy of the calculations. The
components of forces acting on randomly selected atoms calculated
by EW2D, EW2D with the table, and EW3DC from
a given conﬁguration of 512 water molecules have been
compared. We found that EW3DC shows better accuracy
than the EW2D with the precalculated table compared with
the direct EW2D calculation. We think that an accumulation
of errors introduced by the interpolation of the table and the
parallel plate capacitor approximation results in a loss of
accuracy when the EW2D method with the precalculated
table is used.
V. CONCLUSION
In this study, Ewald summation techniques for systems
periodic in two dimensions with a ﬁnite length in the third
dimension have been tested. In terms of speed and accuracy,
the usual three-dimensional Ewald method with a correction
term for the slab geometry ͑EW3DC͒ seems to be the best
choice. The inclusion of the correction term does not introduce
any signiﬁcant computational difﬁculty and can be easily
incorporated in the standard EW3D program. Various optimization
techniques available for the EW3D technique
could be readily applicable to the EW3DC.3,4
Recently, there have been considerable interests in the
application of the ﬁrst-principles simulations using the Car–
Parrinello scheme41
to the interfacial systems such as water–
metal interfaces.42
This very promising approach may play
an important role in the future simulations of various interfacial
systems. The correct treatment of long-range forces in
this approach is required to come to any reliable conclusions
from the ﬁrst-principles simulations, especially on charged
systems. The EW3D with the correction term can be conveniently
used for the calculation of long-range forces in the
ﬁrst-principles simulations.
In this study, we clearly demonstrated how simulations
of interfacial systems using the EW3D method without the
correction term can lead to erroneous results. Numerous
simulation studies on interfacial systems in the literature
have been carried out by using the EW3D method without
considering the correction term. We do not know how serious
the errors are for each studied system and which properties
are most likely to be affected by the negligence of the
correction term. More careful studies need to be done to
resolve this matter and we plan to carry out further simulation
studies in this direction. The safest and economical
choice for the calculation of the long-range forces in most of
simulations of interfacial systems seems to be the EW3D
with the correction term.
ACKNOWLEDGMENT
We are grateful to the Ofﬁce of Naval Research for the
support.
1
M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids ͑Oxford
University, New York, 1987͒.
2
S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. London,
Ser. A 373, 27 ͑1980͒; 388, 177 ͑1983͒.
3
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G.
Pedersen, J. Chem. Phys. 103, 8577 ͑1995͒.
4
A. Y. Toukmaji and J. A. Board Jr., Comput. Phys. Commun. 95, 73
͑1996͒.
5
D. E. Parry, Surf. Sci. 49, 433 ͑1975͒; 54, 195 ͑1976͒.
6
D. M. Heyes, M. Barber, and J. H. R. Clarke, J. Chem. Soc., Faraday
Trans. II 73, 1485 ͑1977͒.
TABLE III. Comparison of computing times per time step calculated by
EW2D and EW3DC methods.
Method Time ͑s͒
EW2D 59.9
EW2D ͑optimized͒ 18.1
EW2D ͑table͒ 7.52
EW3DC 0.71
3161J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 Ewald summation for a slab
Copyright ©2001. All Rights Reserved.
7
S. W. de Leeuw and J. W. Perram, Mol. Phys. 37, 1313 ͑1979͒.
8
E. Spohr, J. Chem. Phys. 107, 6342 ͑1997͒.
9
F. E. Harris, Int. J. Quantum Chem. 68, 385 ͑1998͒.
10
J. Hautman, J. W. Halley, and Y.-J. Rhee, J. Chem. Phys. 91, 467 ͑1989͒.
11
Y.-J. Rhee, J. W. Halley, J. Hautman, and A. Rahman, Phys. Rev. B 40,
36 ͑1989͒.
12
E. Wasserman, J. R. Rustad, A. R. Felmy, B. P. Hay, and J. W. Halley,
Surf. Sci. 385, 217 ͑1997͒.
13
J. Hautman and M. L. Klein, Mol. Phys. 75, 379 ͑1992͒.
14
G. Aloisi, M. L. Foresti, R. Guidelli, and P. Barnes, J. Chem. Phys. 91,
5592 ͑1989͒.
15
J. Lekner, Physica A 176, 485 ͑1991͒.
16
S. Y. Liem and J. H. R. Clarke, Mol. Phys. 92, 19 ͑1997͒.
17
A. H. Widmann and D. B. Adolf, Comput. Phys. Commun. 107, 167
͑1997͒.
18
M. L. Berkowitz, I.-C. Yeh, and E. Spohr, to appear in Interfacial Electrochemistry:
Theory, Experiment, and Applications, edited by A. Wieckowski
͑Marcel Dekker, New York͒ in press.
19
I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 110, 7935 ͑1999͒.
20
I.-C. Yeh and M. L. Berkowitz, Chem. Phys. Lett. 301, 81 ͑1999͒.
21
J. C. Shelley and G. N. Patey, Mol. Phys. 88, 385 ͑1996͒.
22
J. Alejandre, D. J. Tildesly, and G. A. Chapela, J. Chem. Phys. 102, 4574
͑1995͒.
23
J. C. Shelley, G. N. Patey, D. R. Be´rard, and G. M. Torrie, J. Chem. Phys.
107, 2122 ͑1997͒.
24
S. E. Feller, R. W. Pastor, A. Rojnuckarin, S. Bogusz, and B. R. Brooks,
J. Phys. Chem. 100, 17011 ͑1996͒.
25
G. Hummer, N. Grønbech-Jensen, and M. Neumann, J. Chem. Phys. 109,
2791 ͑1998͒.
26
G. Hummer, L. R. Pratt, and A. E. Garcı´a, J. Phys. Chem. A 102, 7885
͑1998͒.
27
P. H. Hu¨nenberger and J. A. McCammon, J. Chem. Phys. 110, 1856
͑1999͒.
28
S. Boresch and O. Steinhauser, Ber. Bunsenges. Phys. Chem. 101, 1019
͑1997͒.
29
S. Bogusz, T. E. Cheatham III, and B. R. Brooks, J. Chem. Phys. 108,
7070 ͑1998͒.
30
J. E. Roberts and J. Schnitker, J. Chem. Phys. 101, 5024 ͑1994͒.
31
J. E. Roberts and J. Schnitker, J. Phys. Chem. 99, 1322 ͑1995͒.
32
E. R. Smith, Proc. R. Soc. London, Ser. A 375, 475 ͑1981͒.
33
R. P. Feynman, R. B. Leighton, and M. Sands, The Feynman Lectures on
Physics ͑Addison-Wesley, Reading, 1964͒, Vol. 2.
34
J. A. Hernando, Phys. Rev. A 44, 1228 ͑1991͒.
35
I.-C. Yeh and M. L. Berkowitz, J. Electroanal. Chem. 450, 313 ͑1998͒.
36
H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91,
6269 ͑1987͒.
37
K. Raghavan, K. Foster, K. Motakabbir, and M. Berkowitz, J. Chem.
Phys. 94, 2110 ͑1991͒.
38
H. J. C. Berendsen, J. P. M. Postma, W. van Gunsteren, A. DiNola, and J.
R. Haak, J. Chem. Phys. 81, 3684 ͑1984͒.
39
M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Phys. Chem. 91, 4873
͑1987͒.
40
J. N. Glosli and M. R. Philpott, Electrochim. Acta 41, 2145 ͑1996͒.
41
R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 ͑1985͒.
42
J. W. Halley, A. Mazzolo, Y. Zhou, and D. Price, J. Electroanal. Chem.
450, 273 ͑1998͒.
3162 J. Chem. Phys., Vol. 111, No. 7, 15 August 1999 I.-C. Yeh and M. L. Berkowitz
Copyright ©2001. All Rights Reserved.