Inﬂuence of cut-off truncation and artiﬁcial periodicity of electrostatic
interactions in molecular simulations of solvated ions:
A continuum electrostatics study
Michael Bergdorf, Christine Peter, and Philippe H. Hu¨nenbergera)
Laboratorium fu¨r Physikalische Chemie, ETH Zu¨rich, CH-8093 Zu¨rich, Switzerland
͑Received 13 May 2003; accepted 6 August 2003͒
A new algorithm relying on ﬁnite integration is presented that solves the equations of continuum
electrostatics for truncated ͑and possibly reaction-ﬁeld corrected͒ solute–solvent and solvent–
solvent interactions under either nonperiodic or periodic boundary conditions. After testing and
validation by comparison with existing methods, the algorithm is applied to investigate the effect of
cut-off truncation and artiﬁcial periodicity in explicit-solvent simulations of ionic solvation and
ion–ion interactions. Both cut-off truncation and artiﬁcial periodicity signiﬁcantly alter the
polarization around a spherical ion and thus, its solvation free energy. The nature and magnitude of
the two perturbations are analyzed in details, and correction terms are proposed for both effects.
Cut-off truncation is also shown to induce strong alterations in the potential of mean force for
ion–ion interaction. These observations help to rationalize artifacts previously observed in explicit–
solvent simulations, namely spurious features in the radial distribution functions close to the cut-off
distance and alterations in the relative stabilities of contact, solvent-separated and free ion
pairs. © 2003 American Institute of Physics.
͓DOI: 10.1063/1.1614202͔
I. INTRODUCTION
Computer simulation with an explicit representation of
the solvent molecules has become a standard tool for investigating
the structure, dynamics, and function of ͑bio-͒molecules
in solution.1–6
However, due to important computational
costs, the system sizes that are accessible to
such simulations remain truly microscopic (typically
Ͻ1000 nm3
). As a direct consequence, the longest-range
(Ͼ5 nm) component of intermolecular interactions, which is
generally dominated by electrostatics, cannot be computed in
an exact manner. Unfortunately, because electrostatic interactions
are of large magnitude, many simulated observables
turn out to be highly sensitive to the treatment of these interactions
and, due to their long range, to the boundary conditions
used in the simulation ͑system size and shape, ﬁnite
versus periodic system͒. For this reason, the approximate
representation of long-range electrostatic interactions in
explicit-solvent ͑bio-͒molecular simulations is probably
nowadays one of the principal bottlenecks in the accuracy of
these methods. Uncontrolled approximations can give rise to
important artifacts ͑so-called ﬁnite-size effects͒, which may
strongly impair the reliability of many current simulations.
There is thus considerable effort in the scientiﬁc community
towards the goal of improving the representation of electrostatics
in computer simulations.
The vast majority of explicit-solvent ͑bio-͒molecular
simulations are carried out under periodic boundary
conditions6,7
͑PBC͒. In this case, the solute ͑bio-͒molecule is
placed into a computational box ͑space-ﬁlling shape, e.g.,
rectangular͒, and the empty volume is ﬁlled by solvent molecules.
The system considered in the simulation consists of
the central box surrounded by an inﬁnite array of periodic
copies of itself, which has the advantage of removing any
surface distortion associated with a solvent-vacuum boundary.
There are essentially three methods to handle electrostatic
interactions in simulations under PBC: ͑i͒ Straight
truncation of the Coulomb interactions at a convenient cutoff
distance;6,7
͑ii͒ smooth truncation of the Coulomb interactions,
e.g., by means of a switching or shifting function8–13
or by including a reaction-ﬁeld correction;14–19
͑iii͒ use of
lattice-sum methods ͑Ewald,20
P3
M,21
or PME22,23
methods͒.
Cut-off truncation reduces the computational costs and the
effect of artiﬁcial periodicity in simulations. However,
straight truncation ͑ST͒ represents a very severe approximation,
leading to heating as well as important artifacts in simulated
properties of liquids,16,24,25
solvated ions,26–35
ion
pairs,9,36–44
and biomolecules.45–48
Smooth-truncation methods
may be applied to reduce the heating caused by the application
of a cut-off, but nevertheless retain ͑and sometimes
amplify͒ a number of the undesirable effects of abrupt
truncation.9,10,27,28,41,48–51
Furthermore, these methods are
generally ad hoc and lack any physical basis. An exception is
the inclusion of a Barker–Watts reaction-ﬁeld correction14–16
͑RF͒ to the cut-off truncation. This correction scheme approximately
accounts for the mean effect of the medium beyond
the cut-off sphere of each particle by assuming that this
medium behaves as a homogeneous dielectric continuum of
permittivity equal to that of the solvent. Due to the form of
the added reaction-ﬁeld term, the correction effectively acts
as a ͑physically based͒ switching function in the limit of
a͒
Author to whom correspondence should be addressed. Telephone: ϩ41 1
632 5503; Fax: ϩ41 1 632 1039.
Electronic mail: phil@igc.phys.chem.ethz.ch
JOURNAL OF CHEMICAL PHYSICS VOLUME 119, NUMBER 17 1 NOVEMBER 2003
91290021-9606/2003/119(17)/9129/16/$20.00 © 2003 American Institute of Physics
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
high solvent permittivities. Finally, lattice-sum ͑LS͒
methods20–23,52–54
rely on Fourier series to describe the longrange
component of electrostatic interactions, i.e., they assume
that these interactions are exactly periodic within the
inﬁnite system. Although LS and RF methods rely on more
or less reasonable approximations for dealing with the longrange
component of electrostatic interactions, an approximate
treatment is certainly preferable to the plain omission
of this long-range component, as done in the ST scheme.
Nevertheless, some dependence of simulated observables on
the cut-off distance or system size has also been evidenced
for the LS,19,55–59
and RF27,28,60
methods. It is therefore of
importance to carefully investigate and compare the properties
of the three most common electrostatic schemes ͑ST, RF,
and LS͒.
A general strategy to analyze ﬁnite-size effects and improve
electrostatic schemes for explicit-solvent simulations
relies on the use of continuum electrostatics.26–28,51,52,56–64
In
the continuum-electrostatics approach, the solute is treated as
a low-dielectric cavity encompassing the solute atomic point
charges, and embedded in a dielectric continuum of permittivity
equal to that of the solvent. In the classical implementation
of the method,65–68
the electrostatic potential in the
system is computed by numerically solving the Poisson ͑or
Poisson–Boltzmann, in the presence of implicit counterions͒
equation, giving access to the electrostatic solvation
free energy of the solute. Although there is no choice of
boundary conditions that adequately mimics an inﬁnitely dilute
solution in explicit-solvent simulations, this is not the
case in continuum-electrostatics calculations. There, the
boundary conditions to solve the Poisson equation are speciﬁed
in the form of the potential at the surface of the computational
box. For a reasonably large solute–wall distance,
this potential is well approximated by the solvent-screened
Coulomb potential of the solute charges.57
In this way, continuum
electrostatics can be used to estimate, for a given
solute conﬁguration, the electrostatic solvation free energy
corresponding to exact ͑nontruncated͒ Coulomb interactions
͑CB͒ under nonperiodic boundary conditions ͑NPBC͒, a
good model for the ideal situation of a solute at inﬁnite dilution.
This suggests that artifacts linked with the use of
approximate electrostatic interactions and periodic boundary
conditions in explicit-solvent simulations could be investigated
using continuum electrostatics, provided that the
method is generalized to these modiﬁed interactions and
boundary conditions. Such generalizations have recently
been developed26–28,51,52,57–64
for nearly all types of relevant
electrostatic interaction schemes ͑CB/LS, ST, or RF͒ and
boundary conditions ͑NPBC or PBC͒, as summarized in
Table I. By comparing, for a given solute conﬁguration, the
outcome of a continuum-electrostatics calculation based on
modiﬁed interactions and boundary conditions with that of
another calculation based on CB interactions under NPBC, it
is possible to estimate the perturbation ⌬⌬Gsolv of the solvation
free energy. The corresponding perturbation ⌬Edirect
in the direct electrostatic interaction energy between solute
atomic charges is straightforward to calculate. The sum
⌬⌬Gel of the two contributions represents the perturbation
of the electrostatic free energy of the system due to the use of
approximate electrostatics and boundary conditions in the
simulation. This procedure is illustrated schematically in Fig.
1 for the speciﬁc case of an explicit-solvent simulation employing
ST or RF electrostatics under PBC. The quantity
⌬⌬Gel ͑possibly evaluated for multiple solute conﬁgurations͒
gives the required information to investigate the nature
TABLE I. Generalizations of the continuum-electrostatics approach to
modiﬁed electrostatic interactions and boundary conditions. The methods
have been classiﬁed using the codes 3D ͑problem solved in three dimensions͒,
1D ͑problem reduced to a one-dimensional equation by symmetry͒,
analytical ͑analytical solution available͒, Poisson ͑based on solving the Poisson
equation͒, FFT ͑based on the use of fast Fourier transforms͒, or direct
͑based on solving ﬁeld equations analogous to Eqs. ͑1͒ and ͑6͒ in real-
space͒.
Electrostatics NPBC PBC
CB/LSa
3D-Poisson ͑Ref. 95͒b
3D-Poisson ͑Ref. 57͒
3D-FFT ͑Refs. 51, 63, and 64͒
͑spherical ionc
͒ 1D-Analytical ͑Ref. 86͒ 1D-Analytical ͑Ref. 57͒
SC,RFd
3D-Directe
3D-Directe
3D-FFT ͑Refs. 51 and 60͒
͑spherical ionc
͒ 1D-Direct ͑Refs. 26–28͒ -
a
Interactions follow from the Coulomb ͑CB͒ potential under NPBC or the
Ewald lattice-sum ͑LS͒ potential under PBC.
b
First calculation on a biomolecule using a ﬁnite-difference algorithm ͑many
alternatives have been proposed, including, e.g., ﬁnite-element, boundaryelement
and inducible-multipole algorithms͒.
c
Solutions developed for the special case of a single spherical ion.
d
Interactions follow from the Coulomb potential truncated at a given cut-off
distance, without ͑ST͒ or with ͑RF͒ the inclusion of a reaction-ﬁeld correc-
tion.
e
Developed in the present article.
FIG. 1. Schematic illustration of the procedure used to assess, based on
continuum electrostatics, artifacts linked with approximate electrostatics and
boundary conditions in explicit-solvent simulations. Ideally, an explicitsolvent
simulation aiming to describe a solute molecule ͑symbolized by a
black sphere͒ at inﬁnite dilution should be based on a quasi-macroscopic
system under NPBC together with exact CB interactions ͑top left drawing͒.
Due to computational limitations this is not feasible in practice, and one may
simulate instead a system under PBC with ST ͑or RF͒ electrostatic interactions
͑top right drawing͒. The corresponding perturbation can be evaluated
by considering the implicit-solvent analogs of the two cases. Using continuum
electrostatics ͑for a given solute conﬁguration͒, the solvation free
energies and the direct interactions between solute charges can be computed
both under CB/NPBC ͑bottom left drawing; based on a good approximation
for the electrostatic potential at the surface of the computational volume͒
and under ST/PBC or RF/PBC ͑bottom right drawing͒. The free-energy
difference ⌬⌬Gel represents the perturbation of the electrostatic free energy
induced by the use approximate electrostatics and boundary conditions, and
is a key quantity for the analysis of ﬁnite-size effects in the explicit-solvent
simulation.
9130 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
and magnitude of the corresponding artifacts in explicitsolvent
simulations.
Since they ignore the discrete nature of the solvent,
continuum-electrostatics models have some limitations, including
an important sensitivity to empirical model parameters
͑atomic charges, atomic radii, solute permittivity, exact
deﬁnition of the solute–solvent boundary͒, the neglect of
nonlinear effects ͑electrostriction, dielectric saturation͒, and
the neglect of the detailed solvent structure around the solute
͑structure of the ﬁrst solvation shells, speciﬁc hydrogen
bonds͒. Furthermore, the electrostatic contribution to the solvation
free energy should be complemented by a nonpolar
contribution, typically assumed to be proportional to the
solvent-exposed surface area. However, since the present
method relies on the comparison of two closely related
continuum-electrostatics calculations involving the same parameters
and solute conﬁguration ͑Fig. 1, bottom drawings͒,
it is likely that errors in the short-range description of solvation
cancel out to a large extent. The difference will depend
almost exclusively on long-range effects, for which continuum
electrostatics can be expected to be accurate. Thus,
electrostatic free energy differences from continuum electrostatics
should be almost quantitatively transposable to interpret
ﬁnite-size artifacts in explicit-solvent simulations.
Inspection of Table I reveals one missing entry. There is
currently no general continuum-electrostatics method to deal
with truncated electrostatic interactions ͑ST or RF͒ under
NPBC, although a method exists in the special case of a
single spherical ion26–28
͓one-dimensional ͑1D͒-Direct
method͔. The goal of the present article is to describe and
apply a general method based on ﬁeld equations and a ﬁniteintegration
algorithm ͓three-dimensional ͑3D͒-Direct
method͔. In addition to dealing with the NPBC case, this
new method is also applicable to systems under PBC. However,
it scales rather unfavorably with the system size ͑as
Ng
6
, where Ng is the number of grid points along each Cartesian
direction͒, and can only be used for small systems.
Therefore, its application is restricted here to the investigation
of the consequences of cut-off truncation and artiﬁcial
periodicity of electrostatic interactions in molecular simulations
of ionic solvation and ion–ion interaction. These systems
are very important benchmark systems for evaluating
the accuracy of electrostatic interactions in molecular simulations
because ͑i͒ they offer the simplest context to investigate
electrostatic ﬁnite-size effects, and ͑ii͒ despite the apparent
simplicity of the problem, the accurate determination
of ionic solvation free energies9,19,26–35,55,57,69–79
and ion–ion
potentials of mean force9,36–44,57,80–84
has turned out to be a
surprisingly difﬁcult problem.
In the present article, the algorithm is described in details
and the inﬂuence of various parameters controlling its
behavior is investigated. The accuracy of the algorithm is
further tested by comparing solvation free energies computed
for a single spherical ion to values estimated through the
1D-Direct26–28
͑NPBC͒ or 3D-FFT51,60
͑PBC͒ methods. Finally,
the present 3D-Direct method is applied to investigate
the effect of cut-off truncation and artiﬁcial periodicity in
computer simulations of ionic solvation and ion–ion interac-
tions.
II. THEORY
A. Continuum electrostatics
In the continuum-electrostatics approach, a solute–
solvent system is modeled ͑for a given solute conﬁguration͒
as a set of solute atomic partial charges embedded in a polarizable
medium of heterogeneous dielectric permittivity.
Application of the laws of electrostatics within such a system
leads to the following expression24,25,51,64,85
for the electric
ﬁeld E͑r͒
E͑r͒ϭV͑r͒ϩ ͵͵͵R3
d3
rЈTᠪ ͑rϪrЈ͒P͑rЈ͒, ͑1͒
where V͑r͒ is the vacuum ﬁeld ͑electric ﬁeld generated by
the solute atomic partial charges in the absence of polarizable
medium͒, P͑r͒ the polarization ͑dipole moment density͒,
and Tᠪ (r) the dipole–dipole interaction tensor characterizing
the solvent–solvent interactions within the system.
In the application of Eq. ͑1͒, it will be assumed that both
solute–solute and solute–solvent electrostatic interactions
are described by truncated Coulomb interactions with a
Barker–Watts reaction-ﬁeld correction.14–16
In the Barker–
Watts scheme, the potential generated at r by a unit charge at
the origin is given by
␺BW͑r͒ϭ
1
4␲⑀o
H͑RϪr͒ͩrϪ1
ϩ
␣r2
2R3 Ϫ
␣ϩ2
2R ͪ, ͑2͒
where H(r) is the Heaviside function ͓H(r)ϭ1 if rϾ0,
H(r)ϭ0 otherwise͔, ⑀o the permittivity of vacuum, and R is
the cut-off distance. The parameter ␣ ͑with 0р␣р1) determined
by the relative dielectric permittivity ⑀Ј of the medium
surrounding the cut-off sphere of each particle through
␣ϭ
2͑⑀ЈϪ1͒
2⑀Јϩ1
. ͑3͒
The function ␺BW in Eq. ͑2͒ accounts both for the direct
Coulombic potential generated by the charge (rϪ1
term͒ and
for the polarization by the neighboring charges of the medium
outside its cut-off sphere (r2
term͒. In the present
work, the discussion of the general form of the Barker–Watts
interaction function ͑BW͒ will essentially focus on the cases
␣ϭ0 (⑀Јϭ1), corresponding to straight truncation of the
Coulomb interactions without reaction-ﬁeld correction ͑ST͒,
and ␣ϭ1 (⑀Ј→ϱ), corresponding to truncated Coulomb interactions
with a reaction-ﬁeld correction corresponding to a
conducting medium ͑RF͒.
When Eq. ͑2͒ is applied to the solvent–solvent interactions
under nonperiodic boundary conditions ͑NPBC͒, the
dipole–dipole interaction tensor reads
Tᠪ NPBC͑r͒ϭ
1
4␲⑀o
H͑RSSϪr͒
ϫٌ ٌͩrϪ1
ϩ
␣r2
2RSS
3 Ϫ
␣ϩ2
2RSS
ͪ
ϭ
1
4␲⑀o
H͑RSSϪr͒ͩ3r rϪr2
1ᠪ
r5 ϩ
␣1ᠪ
RSS
3 ͪ, ͑4͒
9131J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
where RSS is the solvent–solvent cut-off distance and the
notation a b has been introduced for the tensor with elements
␮␯ given by a␮b␯ . Under periodic boundary conditions
͑PBC͒, and provided that RSS is smaller than half the
smallest dimension of the computational box ͑which will be
assumed from here on͒, the dipole–dipole interaction tensor
reads
Tᠪ PBC͑r͒ϭ
1
4␲⑀o
H͑RSSϪr¯͒ͩ3r¯ r¯Ϫr¯2
1ᠪ
r¯5 ϩ
␣1ᠪ
RSS
3 ͪ, ͑5͒
where r¯ is the minimum-image vector corresponding to r.
Using the approximation of linear response, the reaction
of the polarizable medium is linear in the local electric
ﬁeld,24,25,51,64,85
i.e.,
P͑r͒ϭ⑀o͓⑀͑r͒Ϫ1͔E͑r͒, ͑6͒
where ⑀(r) the relative dielectric permittivity of the medium,
which may be heterogeneous in space. Typically, one distinguishes
between solute and solvent regions, characterized by
distinct homogeneous permittivity values.
The dipole–dipole interaction tensor Tᠪ (r) deﬁned by
Eqs. ͑4͒ or ͑5͒ is singular at the origin. However, the singularity
is integrable when applying Eq. ͑1͒. More precisely,
deﬁning ⍀(r;␳) as the sphere of radius ␳ and surface ⌺(r;␳)
centered at r, one may write
͵͵͵R3
d3
rЈTᠪ ͑rϪrЈ͒P͑rЈ͒
ϭI͑r͒ϩ lim
␳→0
͵͵͵R3‫گ‬⍀(r;␳)
d3
rЈTᠪ ͑rϪrЈ͒P͑rЈ͒. ͑7͒
The ﬁrst term can be evaluated as
I͑r͒ϭ lim
␳→0
͵͵͵⍀(r;␳)
d3
rЈTᠪ ͑rϪrЈ͒P͑rЈ͒
ϭ
1
4␲⑀o
ͫlim
␳→0
͵͵͵⍀(0;␳)
d3
sٌ ٌsϪ1
ͬP͑r͒
ϭϪ
1
4␲⑀o
ͫlim
␳→0
␳Ϫ2
͵͵⌺(0;␳)
d2
␴sϪ2
s sͬP͑r͒
ϭϪ
1
2⑀o
ͫ͵0
␲
d␪ sin ␪ cos2
␪ͬP͑r͒ϭϪ
1
3⑀o
P͑r͒.
͑8͒
The second equality follows from inserting Eqs. ͑4͒ or ͑5͒,
deﬁning sϭrЈϪr, and noting that as ␳ tends towards zero:
͑i͒ The Heaviside function evaluates to one for any ﬁnite
RSS ; ͑ii͒ the contribution proportional to the unit tensor vanishes;
͑iii͒ P(rЈ) may be approximated by P͑r͒ and factorized
from the integral. The third equality follows from applying
the gradient theorem and inserting ٌsϪ1
ϭϪsϪ3
s.
The fourth equality follows from observing that, due to symmetry,
the off-diagonal elements of the tensor vanish upon
integration, and that the diagonal elements are all equal. The
ﬁfth equality follows from evaluating one of these diagonal
elements ͑integrand sϪ2
sz
2
ϭcos2
␪) in spherical coordinates.
Combining Eqs. ͑1͒, ͑6͒, ͑7͒, and ͑8͒, the equation to be
solved for the electrostatic ﬁeld E͑r͒ reads
E͑r͒ϭV͑r͒Ϫ 1
3 ͓⑀͑r͒Ϫ1͔E͑r͒
ϩ⑀o lim
␳→0
͵͵͵R3‫گ‬⍀(r;␳)
d3
rЈTᠪ ͑rϪrЈ͒
ϫ͓⑀͑rЈ͒Ϫ1͔E͑rЈ͒, ͑9͒
with Tᠪ (r) given by Eqs. ͑4͒ or ͑5͒. If this equation can be
solved for E͑r͒, the free energy of interaction between the
solute atomic point charges and the polarizable medium is
given by
⌬GϭϪ
1
2
͵͵͵R3
d3
rV͑r͒•P͑r͒. ͑10͒
In the special case of a nonpolarizable solute, this quantity
represents the solvation free energy of the solute.
B. Discretization
To transform the solving of Eq. ͑9͒ into a computationally
tractable problem, three approximations are made. First,
the inﬁnite integration domain is reduced to a ﬁnite region of
space. More precisely, two types of computational domains
are considered: ͑i͒ A spherical volume of radius S surrounded
by vacuum under NPBC, or ͑ii͒ a cubic unit cell of edge L
under PBC. In both cases, the restriction to a ﬁnite computational
domain is expected to have limited consequences
when SӷRSS under NPBC or Lӷ2RSS under PBC ͑because
the truncation of solvent–solvent interactions largely reduces
dipole–dipole correlations at large distances͒, provided that
the vacuum ﬁeld is only active over a small region within
this domain ͑which will be the case due to truncation of the
solute–solvent interactions͒. Second, the solute is assumed to
be nonpolarizable and the solvent to be represented by a
medium of homogeneous permittivity. Thus, the computational
domain comprises two subdomains characterized by
different homogeneous dielectric permittivity values: A ͑possibly
discontinous͒ solute subdomain of relative permittivity
one, and a solvent subdomain of relative permittivity ⑀s .
Third, the problem is discretized by paving the computational
domain using N grid cells, leading to piecewiseconstant
representations Vi and Ei ͑with iϭ1,...,N) of the
vacuum and electric ﬁelds.
Within these approximations, Eq. ͑9͒ becomes
Ei␮ϭVi␮Ϫ
⑀sϪ1
3
␴iEi␮
ϩ⑀o͑⑀sϪ1͚͒jϭ1
N
vj͑1Ϫ␦ij͒␴j ͚␯ϭ1
3
T␮␯͑riϪrj͒Ej␯ ,
͑11͒
where the ␮ and ␯ indexes enumerate Cartesian components,
ri and vi are the center coordinate and volume of grid cell i,
the exterior function ␴i evaluates to one if ri is within the
solvent subdomain and zero otherwise, and the matrix elements
of the dipole–dipole interaction tensor ͓Eqs. ͑4͒ or
͑5͔͒ take the form
9132 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
TNPBC,␮␯͑r͒
ϭ
1
4␲⑀o
H͑RSSϪr͒ͩ3r␮r␯Ϫ␦␮␯r2
r5 ϩ
␣␦␮␯
RSS
3 ͪ ͑12͒
or
TPBC,␮␯͑r͒
ϭ
1
4␲⑀o
H͑RSSϪr¯͒ͩ3r¯␮r¯␯Ϫ␦␮␯r¯2
r¯5 ϩ
␣␦␮␯
RSS
3 ͪ. ͑13͒
Deﬁning the 3Nϫ3N-dimensional matrix Aᠪ as
Ai␮j␯ϭͩ1ϩ
⑀sϪ1
3
␴iͪ␦ij␦␮␯
Ϫ⑀o͑⑀sϪ1͒vj͑1Ϫ␦ij͒␴jT␮␯͑riϪrj͒, ͑14͒
Eq. ͑11͒ can be rewritten in matrix notation
Aᠪ EϭV. ͑15͒
Because the size of the matrix Aᠪ is generally very large,
it cannot be stored in memory and direct methods cannot be
used to solve Eq. ͑15͒. Therefore, an under-relaxed Jacobi
method is applied here to obtain successive approximate solutions
for the discretized electric ﬁeld. Given the approximate
solution E(k)
at iteration k, the grid E(kϩ1)
at the next
iteration is computed as
E(kϩ1)
ϭE(k)
Ϫ␭Dᠪ Ϫ1
͑Aᠪ E(k)
ϪV͒, ͑16͒
where Dᠪ is the diagonal matrix deﬁned by the diagonal elements
of Aᠪ , and ␭ ͑with 0Ͻ␭р1) is a relaxation parameter.
A reasonable initial guess for E(0)
is provided by the vacuum
ﬁeld scaled by the solvent permittivity, i.e.,
E(0)
ϭ⑀s
Ϫ1
V. ͑17͒
It is easily seen that a self-consistent solution of Eq. ͑16͒
must satisfy Eq. ͑15͒. In order to assess the convergence of
the numerical solution upon iterating, the residual ͑with units
of an electric ﬁeld͒
␶(k)
ϭ͚ͩiϭ1
N
vi␴iʈ͑Aᠪ E(k)
ϪV͒iʈ2
͚iϭ1
N
vi␴i
ͪ1/2
, ͑18͒
is introduced as a measure of accuracy.
After solving Eq. ͑15͒ for the discretized electric ﬁeld,
the solvation free energy can be evaluated as ͓Eq. ͑10͒; nonpolarizable
solute͔
⌬GsolvϭϪ 1
2 ⑀o͑⑀sϪ1͚͒iϭ1
N
vi␴iVi•Ei , ͑19͒
where Eq. ͑6͒ was used. Note that the solvation free energy
solely depends on the electric ﬁeld within the solvent subdomain.
Furthermore, due to the form of Eq. ͑14͒, the ﬁeld Ei
corresponding to a point i in the solvent subdomain does not
depend on the ﬁeld Ej at any point j within the solute subdomain.
For this reason, increased computational efﬁciency
can be achieved by omitting all grid points of the solute
subdomain from the deﬁnition of the matrix Aᠪ and the determination
of the solution of Eq. ͑15͒.
C. Application to solvated spherical ions
For a system consisting of a single spherical ion of radius
RI and charge qI centered in the computational domain,
the vacuum ﬁeld corresponding to ion–solvent interactions
described by the Barker–Watts scheme ͓Eq. ͑2͔͒ is given by
V͑r͒ϭ
qI
4␲⑀o
H͑RISϪr͒ͩr
r3 Ϫ
␣r
RIS
3 ͪ, ͑20͒
where RIS is the ion–solvent cut-off distance. This equation
is valid under both NPBC and PBC, provided that RIS is
smaller than half the smallest dimension of the computational
box ͑which will be assumed from here on͒.
Because the ion is generally small compared to the size
of the computational domain, while the variations of the
electric ﬁeld within the solvent are typically largest close to
its surface, the accuracy of the results will depend crucially
on the detailed representation of the ionic surface. For this
reason, two levels of grid resolution are used. First, a coarse
grid of spacing ⌬ is generated, that covers the entire computational
domain ͑leading to grid-cell volumes viϭ⌬3
). Second,
all cells of the coarse grid with their center closer than
()/2)⌬ from the surface of the ion are further discretized,
i.e., they are replaced by a set of ﬁner grid cells of edge ␦
ϭnϪ1
⌬ where n is a positive integer ͑leading to grid-cell
volumes viϭ␦3
). Any grid-cell center of the ﬁner grid that is
closer than RI from the ion center is discarded from the
calculation ͑solute point͒. In this representation, the vacuum
potential at any grid point i is evaluated in practice as
Viϭ
qI
4␲⑀o
I͑RISϪri͒ͩri
ri
3 Ϫ
␣ri
RIS
3 ͪ, ͑21͒
where I(RISϪri) represents the fraction of the grid cell located
within the ion–solvent cut-off sphere.64
After solving Eq. ͑15͒ for the discretized electric ﬁeld,
the radial polarization p(r) around the ion ͓Eq. ͑6͔͒ can be
computed in the form of a histogram. To avoid artifacts
linked with the use of two different grid spacings, this calculation
is based on a uniform grid of NЈ points obtained by
partitioning all cells of the coarse grid into ﬁner grid cells of
edge ␦ sharing a common value of the electic ﬁeld. Under
NPBC, the radial polarization is then computed as
pNPBC͑rn͒ϭ⑀o͑⑀sϪ1͒H͑rnϪRI͒
ϫ
͚iϭ1
NЈ ␴iw͑ri ;rn ,⌬r͒ri
Ϫ1
ri•Ei
͚iϭ1
NЈ ␴iw͑ri ;rn ,⌬r͒
, ͑22͒
where
rnϭ͑nϩ 1
2͒⌬r, nϭ0,1,...,nmax , ͑23͒
⌬r being the histogram width, and
w͑r;rn ,⌬r͒ϭͭ1 if rnϪ 1
2 ⌬rрrϽrnϩ 1
2 ⌬r
0 otherwise.
͑24͒
Under PBC, the periodic copies of the central box must be
taken into account and Eq. ͑22͒ is modiﬁed to
9133J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
pPBC͑rn͒ϭ⑀o͑⑀sϪ1͒H͑rnϪRI͒
ϫ
͚l෈Z3͚iϭ1
NЈ ␴iw͑ri,l ;rn ,⌬r͒ri,l
Ϫ1
ri,l•Ei
͚l෈Z3͚iϭ1
NЈ ␴iw͑ri,l ;rn ,⌬r͒
,
͑25͒
with ri,lϭriϩLl. In practice, the sum over l is restricted to
vectors with integer components in the range ͓Ϫlmax ;lmax͔,
ensuring a correct description of the polarization up to r
ϭ(lmaxϩ1/2)L. These polarization histograms can be compared
to the ideal ͑CB/NPBC͒ Born polarization86
pBorn͑r͒ϭ
qI
4␲
⑀sϪ1
⑀s
rϪ2
. ͑26͒
Combining Eqs. ͑19͒ and ͑21͒, the ionic solvation free
energy is evaluated ͑based on the reﬁned grid of NЈ points͒
as
⌬GsolvϭϪ
qI
8␲
͑⑀sϪ1͚͒iϭ1
NЈ
vi␴iI͑RISϪri͒
ϫͩri
ri
3 Ϫ
␣ri
RIS
3 ͪ•Ei . ͑27͒
This value can be compared to the ideal ͑CB/NPBC͒ Born
electrostatic solvation free energy86
⌬Gsolv
Born
ϭϪ
qI
2
8␲⑀o
⑀sϪ1
⑀s
RI
Ϫ1
. ͑28͒
The application to two ͑or more͒ spherical ions is
straightforward and only requires the following minor adaptations:
͑i͒ The quantity ␴i is zero ͑point discarded from the
calculation͒ for all grid cells with centers located inside any
ion, and one otherwise; ͑ii͒ the vacuum ﬁeld V ͓Eq. ͑21͔͒ is
expressed as a sum of contributions arising from each ion;
͑iii͒ the solvation free energy ⌬Gsolv ͓Eq. ͑27͔͒ is expressed
as a sum of contributions arising from each ion; ͑iv͒ ﬁne
grids ͑spacing ␦Ͻ⌬) are used to describe the close neighborhood
of all ions.
III. COMPUTATIONAL DETAILS
The solution of Eq. ͑15͒, restricted to the case of one or
two spherical ions, was implemented in a C program. The
single ion or the two ions are placed on the z axis of the
coordinate system ͑single ion at zϭ0; two ions at zϭ
Ϯd/2, where d is the interionic distance͒. Taking advantage
of the symmetry of the problem, the storage of the discretized
vacuum and electric ﬁelds is only required for one
quadrant (x,yу0) of the computational domain.
After an evaluation of the convergence properties of the
algorithm, a set of computational parameters was selected
and adopted for all subsequent calculations. The corresponding
values are as follows ͑unless otherwise speciﬁed͒. The
spacings corresponding to the coarse and ﬁne grids were set
to ⌬ϭ0.1 nm and ␦ϭ0.025 nm, respectively. The relaxation
parameter ␭ was set to 0.4. The algorithm was terminated
when the residual ␶(k)
was either below 10Ϫ3
kJ molϪ1
nmϪ1
eϪ1
or reached a minimum value. All calculations under
NPBC used a sphere of radius Sϭ4.0 nm as computational
domain. To ensure that this domain is large enough to
be representative of an inﬁnite nonperiodic system, a number
of single-ion and two-ion calculations were repeated with S
ϭ5.0 nm. The observed differences in solvation free energy
were in all cases below 0.1 kJ molϪ1
in magnitude.
The method was ﬁrst applied to solvated spherical ions.
Ionic solvation free energies ͓Eq. ͑27͔͒ and radial polarization
histograms ͓Eqs. ͑22͒ or ͑25͒ with ⌬rϭ0.025 nm and
lmaxϭ2] were computed for all combinations of the following
parameters: ionic charge qIϭ1 e, ionic radii RI
ϭ0.2 nm ͑about the size of Naϩ
) or 0.4 nm ͑about the size
of ClϪ
), cut-off radii RISϭRSSϭRCϭ0.8 or 1.2 nm, solvent
permittivities ⑀sϭ2 ͑alkanelike solvent͒ or 78 ͑water͒, ␣
ϭ0 ͑ST scheme͒ or 1 ͑RF scheme͒, and NPBC or PBC ͑with
Lϭ2.6 nm). To validate the method, the results are compared
with those of calculations employing other methods
͑Table I͒, namely: ͑i͒ 1D-Direct27,28
͑NPBC; bin size 0.005
nm, range 4.0 nm͒ or ͑ii͒ 3D-FFT60
͑PBC; 180 grid points
along each Cartesian axis͒.
The method was then used to investigate the effect of
periodicity on ionic solvation free energies in systems with
truncated electrostatic interactions ͑ST or RF͒. To this end,
ionic solvation free energies were computed under PBC using
the above combination of parameters, for cubic box
edges L ranging from 1.6 nm (RCϭ0.8 nm) or 2.4 nm (RC
ϭ1.2 nm) to 8.0 nm. The effect on periodicity can be quantiﬁed
by the relative periodicity-induced perturbation of the
ionic solvation free energy ␥(L), deﬁned as
␥͑L͒ϭ
⌬⌬Gsolv͑L͒
⌬Gsolv
NPBC
with
⌬⌬Gsolv͑L͒ϭ⌬Gsolv
PBC
͑L͒Ϫ⌬Gsolv
NPBC
. ͑29͒
Finally, the effect of cut-off truncation ͑ST or RF͒ and
periodicity on the electrostatic solvation contribution
⌬Gsolv(d) to the potential of mean force for the interaction
between two ions at distance d ͑under PBC, ions aligned
along an axis of the cubic unit cell͒ was evaluated for the
special case of charges qIϭϮqJϭ1 e, radii RIϭRJ
ϭ0.4 nm, cut-off radii RISϭRSSϭRCϭ1.2 nm and for a cubic
unit-cell of edge Lϭ6 nm ͑PBC͒. For validation, the results
under PBC were compared with those of calculations
employing the 3D-FFT60
method ͑180 grid points along each
Cartesian axis͒. The corresponding overall electrostatic contribution
⌬Gel(d) to the potential of mean force was also
evaluated as
⌬Gel͑d͒ϭ⌬Gsolv͑d͒ϩqIqJ␺BW͑d˜ ͒, ͑30͒
where d˜ϭd ͑NPBC͒ or d˜ϭmin͕d;LϪd͖ ͑PBC͒, and ␺BW is
given by Eq. ͑2͒ with ␣ϭ0 ͑ST͒ or ␣ϭ1 ͑RF͒. These proﬁles
can be compared with the expected long-range behavior
of the electrostatic potential of mean force for a Coulombic
interaction between the ions, namely
9134 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
⌬Gel
lr
͑d͒ϭ⌬Gsolv
ion
͑qI ,RI͒ϩ⌬Gsolv
ion
͑qJ ,RJ͒
ϩ
qIqJ
4␲⑀o⑀s
,dϪ1
͑31͒
where ⌬Gsolv
ion
(q,R) is the solvation free energy of an isolated
ion of radius R and charge q under NPBC when applying the
speciﬁc electrostatic scheme.
IV. RESULTS
A. Convergence properties
The convergence properties of the under-relaxed Jacobi
algorithm ͓Eq. ͑16͔͒ used to solve Eq. ͑15͒ are illustrated in
Fig. 2 for a spherical ion of charge qIϭ1 e and radius RI
ϭ0.4 nm in a solvent of permittivity ⑀sϭ78, based on the
RF scheme with a single cut-off radius RCϭ1.2 nm and using
three choices of the relaxation parameter ␭. Results for
the ST scheme are qualitatively very similar ͑data not
shown͒. Within few iterations, the residual ␶(k)
decreases
from about 20 to values below 3 kJ molϪ1
nmϪ1
eϪ1
͓Fig.
2͑a͔͒. Convergence to zero residual only occurs when the
solvent permittivity ⑀s is smaller than about 10 ͑data not
shown͒. This limited convergence is probably related to the
presence of a strong discontinuity in the system permittivity
at the ion–solvent boundary in the case of high solvent permittivities.
For ⑀s values larger than about 10, the residual
reaches a minimum after a certain number of iterations ͑typically
about 15–20 for ␭ϭ0.4) and slowly rises again afterwards.
When this situation occurred, the algorithm was terminated
at the minimum value ␶ of the residual. However,
convergence of ␶(k)
towards ␶ is associated with the simultaneous
convergence of ⌬Gsolv
(k)
to a well-deﬁned value
⌬Gsolv ͓Fig. 2͑b͒ and Table II͔. Because the values of both
the minimum residual and the associated converged solvation
free energy are essentially independent of the convergence
parameter ␭, it appears that the method is nevertheless
able to produce accurate results for high ⑀s values. This is
also supported by the observation that for ⑀sϭ2, values of
⌬Gsolv
(k)
when ␶(k)
ϭ3 kJ molϪ1
nmϪ1
eϪ1
typically differ
from the corresponding converged values (␶(k)
Ͻ10Ϫ3
kJ molϪ1
nmϪ1
eϪ1
) by less than 1% ͑data not
shown͒.
Although the convergence parameter ␭ does not inﬂuence
the ﬁnal values of ␶ and ⌬Gsolv , it has a strong impact
on the convergence rate. For ⑀sϭ78, ␭ϭ0.2 leads to slow
convergence, ␭ϭ0.6 to slow convergence and oscillatory
evolution of ⌬Gsolv
(k)
, while the algorithm fails to converge
for ␭ϭ0.8 ͑data not shown͒. In practice, it was found that
␭ϭ0.4 is the optimal choice in this case, and is also adequate
for ⑀sϭ2 ͑although a somewhat larger values may slightly
accelerate convergence͒. This value was adopted for all subsequent
calculations.
The rates of convergence to the minimum ␶ of ␶(k)
appear
to be very similar under NPBC or PBC, and for the RF
and ST ͑data not shown͒ schemes. The ﬁnal values of the
residual, however, are somewhat lower under NPBC compared
to PBC ͑Table II͒.
B. Single spherical ion
The radial polarization p(r) around around a spherical
ion ͓Eqs. ͑22͒ and ͑25͔͒ of charge qIϭ1 e and radius RI
ϭ0.4 nm in a solvent of permittivity ⑀sϭ78 is displayed in
Fig. 3 for different choices of boundary conditions and treatments
of the electrostatic interactions based on a single cutoff
RCϭ1.2 nm. The polarization corresponding to the Born
model86
͓CB/NPBC; Eq. ͑26͔͒ or to the lattice-sum case ͑LS/
PBC; computed using the 3D-FFT method64
͒, and the polarization
computed from the 1D-Direct method27,28
͑ST,RF/
NPBC͒ are also displayed for comparison. For both the ST
and RF schemes under NPBC, the agreement between the
results of the 1D-Direct and of the present 3D-Direct methods
is excellent over the whole range of distances. The only
noticeable difference is the more progressive transition of
p(r) around RC in the 3D-Direct calculation for the ST
scheme, which is due to a signiﬁcantly lower resolution and
to the smoothing of the vacuum ﬁeld at the ion–solvent cutoff
distance ͓function I in Eq. ͑21͔͒. These curves, however,
differ signiﬁcantly from the Born polarization, corresponding
to the ideal situation of a spherical ion solvated by a nonperiodic
Coulombic continuum of inﬁnite extent.
For the ST case, the polarization curve is discontinuous
at the ion–solvent cut-off distance RIS . The polarization below
RIS is consistently larger than predicted by the Born
model, whereas the polarization above is smaller, although
always positive. Underpolarization of the solvent above RIS
is easily understood since the solvent beyond this distance
does not feel directly the electrostatic ﬁeld of the ion. However,
the solvent in this region reacts indirectly to the ionic
FIG. 2. Convergence properties of the under-relaxed Jacobi algorithm ͓Eq.
͑16͔͒ used to solve Eq. ͑15͒. The residual ␶(k)
͓Eq. ͑18͔͒ is displayed as a
function of the number of iterations k ͑a͒, and the solvation free energy
⌬Gsolv
(k)
at iteration k ͓Eq. ͑27͔͒ as a function of the corresponding residual
␶(k)
͑b͒. The system consists of a single spherical ion of charge qIϭ1 e and
radius RIϭ0.4 nm in a solvent of permittivity ⑀sϭ78, and is either nonperiodic
͑NPBC; spherical domain of radius Sϭ4.0 nm) or periodic ͑PBC;
cubic unit cell of edge Lϭ2.6 nm). Electrostatic interactions correspond to
the RF scheme, with cut-off radii RISϭRSSϭRCϭ1.2 nm. Three choices of
the Jacobi relaxation parameter ␭ are compared ͓for NPBC, only the curve
corresponding to ␭ϭ0.4 is displayed in ͑a͔͒.
9135J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
ﬁeld through interactions with the polarized solvent within
the cut-off sphere of the ion, leading to the observed residual
polarization. Inside the cut-off sphere of the ion, the solvent
is overpolarized because each solvent volume element only
interacts with a fraction of the highly polarized solvent
within the cut-off sphere of the ion. This partial interaction
results in a bias of the solvent polarization towards the ion.
For the RF case, the polarization curve is continuous at
RIS , although its derivative is not. The polarization is consistently
smaller than the corresponding Born curve over the
whole range of distance, both below and above RIS . The
difference between the curves is largest at distances close to
RIS , and becomes progressively smaller at either short or
long distances from the ion. It has been shown that in the
limit of high solvent permittivities, the Barker–Watts potential
represents ͑for the solvent–solvent interactions͒ the cutoff-truncated
polynomial of second order ͑terms in rϪ1
to
r2
) that leads to the best agreement between Born and effective
polarizations.27,28
A number of additional results related
to this comparison are derived in Appendix A, namely that ͑i͒
the RF/NPBC polarization converges to the Born polarization
in the limit RIS ,RSS→ϱ; ͑ii͒ in the limit of small distances
͑compared to RC), the RF/NPBC polarization converges
towards the Born polarization; ͑iii͒ in the limit of
large distances ͑compared to RC), the RF/NPBC polarization
becomes proportional to rϪ2
, just as the Born polarization.
For both the ST and RF schemes, the polarization curves
corresponding to PBC are systematically lower ͑in the range
RI to L) compared to the polarization under NPBC. The
reason for this is that under PBC, the solvent in the reference
TABLE II. Solvation free energy ⌬Gsolv of a single spherical ion ͓Eq. ͑27͔͒ computed using the 3D-Direct
method ͑present article͒. The system consists of a single ion of charge qIϭ1 e and radius RI in a solvent of
permittivity ⑀s , and is either nonperiodic ͑NPBC; spherical domain of radius Sϭ4 nm) or periodic ͑PBC; cubic
unit cell of edge Lϭ2.6 nm). Electrostatic interactions correspond to either the ST or RF schemes, with cut-off
radii RISϭRSSϭRC . For the RF scheme, the reaction-ﬁeld contribution to the solvation free energy
͓␣-dependent term in Eq. ͑27͔͒ is reported between parentheses. The converged value ␶ of the residual is also
indicated. The solvation free energies ⌬Gsolv
corr
corrected by the inclusion of a self-energy term ͓Eq. ͑B3͔͒ are also
given. For comparison, the corresponding Born solvation free energies ͓Eq. ͑28͔͒ are Ϫ342.9 (⑀ϭ78, RI
ϭ0.2 nm), Ϫ171.4 (⑀ϭ78, RIϭ0.4 nm), Ϫ173.7 (⑀ϭ2, RIϭ0.2 nm), and Ϫ86.8 (⑀ϭ2, RIϭ0.4 nm)
kJ molϪ1
.
BC Interaction ⑀s
RI
͓nm͔
RC
͓nm͔
␶
͓kJ molϪ1
nmϪ1
eϪ1
͔
⌬Gsolv
͓kJ molϪ1
͔
⌬Gsolv
ref
͓kJ molϪ1
͔
⌬Gsolv
corr
͓kJ molϪ1
͔
NPBC ST 78 0.2 0.8 1.2 Ϫ284.8 Ϫ281.7a
Ϫ370.5
NPBC ST 78 0.2 1.2 1.2 Ϫ306.3 Ϫ303.4a
Ϫ363.4
NPBC ST 78 0.4 0.8 0.8 Ϫ100.9 Ϫ100.1a
Ϫ186.6
NPBC ST 78 0.4 1.2 0.7 Ϫ129.4 Ϫ128.8a
Ϫ186.5
PBC ST 78 0.2 0.8 3.0 Ϫ284.1 Ϫ280.3b
Ϫ369.8
PBC ST 78 0.2 1.2 2.9 Ϫ301.0 Ϫ297.9b
Ϫ355.0
PBC ST 78 0.4 0.8 2.0 Ϫ100.6 Ϫ102.4b
Ϫ186.3
PBC ST 78 0.4 1.2 1.8 Ϫ125.3 Ϫ125.7b
Ϫ182.4
NPBC RF 78 0.2 0.8 1.1 Ϫ217.0 (37.0) Ϫ213.7 (37.1)a
Ϫ345.6
NPBC RF 78 0.2 1.2 1.1 Ϫ259.2 (26.2) Ϫ256.3 (26.4)a
Ϫ344.9
NPBC RF 78 0.4 0.8 0.6 Ϫ47.6 (23.7) Ϫ46.8 (23.7)a
Ϫ176.2
NPBC RF 78 0.4 1.2 0.6 Ϫ86.4 (22.0) Ϫ85.8 (22.2)a
Ϫ172.1
PBC RF 78 0.2 0.8 2.8 Ϫ214.6 (35.8) Ϫ211.0 (35.6)b
Ϫ343.2
PBC RF 78 0.2 1.2 2.9 Ϫ253.1 (22.3) Ϫ249.5 (22.5)b
Ϫ338.8
PBC RF 78 0.4 0.8 1.6 Ϫ47.5 (22.9) Ϫ48.4 (23.0)b
Ϫ176.1
PBC RF 78 0.4 1.2 1.6 Ϫ82.3 (18.7) Ϫ82.3 (18.8)b
Ϫ168.0
NPBC ST 2 0.2 0.8 0.0 Ϫ133.4 Ϫ135.2a
Ϫ176.8
NPBC ST 2 0.2 1.2 0.0 Ϫ146.7 Ϫ148.3a
Ϫ175.6
NPBC ST 2 0.4 0.8 0.0 Ϫ45.6 Ϫ46.4a
Ϫ89.0
NPBC ST 2 0.4 1.2 0.0 Ϫ60.3 Ϫ60.8a
Ϫ89.2
PBC ST 2 0.2 0.8 0.0 Ϫ133.1 Ϫ135.0b
Ϫ176.5
PBC ST 2 0.2 1.2 0.0 Ϫ145.6 Ϫ147.6b
Ϫ174.5
PBC ST 2 0.4 0.8 0.0 Ϫ45.6 Ϫ46.3b
Ϫ89.0
PBC ST 2 0.4 1.2 0.0 Ϫ59.6 Ϫ60.2b
Ϫ88.5
NPBC RF 2 0.2 0.8 0.0 Ϫ99.8 (13.7) Ϫ101.0 (13.8)a
Ϫ164.9
NPBC RF 2 0.2 1.2 0.0 Ϫ123.1 (9.8) Ϫ124.4 (9.8)a
Ϫ166.5
NPBC RF 2 0.4 0.8 0.0 Ϫ20.4 (9.0) Ϫ20.6 (9.1)a
Ϫ85.5
NPBC RF 2 0.4 1.2 0.0 Ϫ39.2 (8.3) Ϫ39.5 (8.4)a
Ϫ82.6
PBC RF 2 0.2 0.8 0.0 Ϫ99.5 (13.8) Ϫ100.9 (13.8)b
Ϫ164.6
PBC RF 2 0.2 1.2 0.0 Ϫ122.4 (9.5) Ϫ124.0 (9.5)b
Ϫ165.8
PBC RF 2 0.4 0.8 0.0 Ϫ20.4 (9.0) Ϫ20.7 (9.1)b
Ϫ85.5
PBC RF 2 0.4 1.2 0.0 Ϫ38.9 (8.0) Ϫ39.2 (8.1)b
Ϫ82.3
a
Solvation free energies ⌬Gsolv
ref
estimated from the 1D-Direct method ͑Refs. 27 and 28͒ are given for comparison
͑ST,RF/NPBC͒.
b
Solvation free energies ⌬Gsolv
ref
estimated from the 3D-FFT method ͑Ref. 60͒ are given for comparison
͑ST,RF/PBC͒.
9136 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
unit cell is perturbed by its interaction with the solvent in
adjacent unit cells ͑itself polarized by the periodic copies of
the ion͒ and, for rϾLϪRIS , with the periodic copies of the
ion themselves. The consequences of these interactions are a
depolarization of the solvent in the reference unit cell ͑compared
to NPBC͒, and the occurence of negative polarization
values for the solvent in the neighboring unit cells. At large
distances, p(r) displays an irregular oscillatory behavior
with values close to zero at the location of the nearest neighbor
ions ͑i.e., L, &L, )L, . . . ; data not shown͒. The depolarization
of the solvent within the reference unit cell is
slightly more important in the RF case compared to the ST
case, which is probably a consequence of the larger magnitude
of the residual solvent polarization above RIS observed
in the RF/NPBC case. In both cases, however, the solvent
depolarization within the reference unit cell remains relatively
small because, due to the truncation of ion–solvent
interactions at RISϽL/2, dipoles in the reference unit cell do
not interact directly with the periodic copies of the ion. Finally,
the polarization corresponding to the RF/PBC scheme
is seen to agree reasonably well with the LS/PBC curve, the
difference being expectedly largest in the neighborhood of
the cut-off distance. There is, however, an important difference
between the two schemes. When cut-off truncation is
applied, the solvation free energy only depends on the polarization
in the range RI to RISϽL/2 ͓Eqs. ͑10͒ and ͑20͔͒ and
the effect of periodicity on the ionic solvation free energy is
expected to be relatively small. If this restriction is removed
͑e.g., when nontruncated LS interactions are considered͒, the
effect of artiﬁcial periodicity on the ionic solvation free energy
becomes dramatically more important.19,55,57
Ionic solvation free energies ⌬Gsolv computed for
spherical ions with different parameter combinations, electrostatic
schemes, and boundary conditions are listed in Table
II. The values ⌬Gsolv
ref
computed using the 1D-Direct
method27,28
͑ST,RF/NPBC͒ or the 3D-FFT method60
͑ST,RF/
PBC͒ are also listed for comparison. Note that, while the
former values are certainly very accurate, the latter values
are probably subject to errors of a similar magnitude as the
present method. The agreement between the values computed
using different methods is in general very good. The
average and maximal relative differences between the
present 3D-Direct and the reference values are 1.1% and
1.9%, respectively. Not unexpectedly, these relative differences
tend to be somewhat larger for ͑i͒ the smaller ionic
radius, ͑ii͒ the larger permittivity value, ͑iii͒ periodic boundary
conditions. The following observations can be made: ͑i͒
the solvation free energies are larger in magnitude for the
smaller ion and the larger permittivity value, in qualitative
agreement with the Born model; ͑ii͒ the solvation free energies
are larger in magnitude for the larger cut-off value, a
consequence of including a larger amount of polarized solvent
within the cut-off sphere of the ion; ͑iii͒ the solvation
free energies are more negative for the ST scheme compared
to the RF scheme, a consequence of the solvent overpolarization
within the cut-off sphere of the ion for the ST scheme
͑Fig. 3͒ and of the inclusion of an additional positive term in
the solvation free energy for the RF scheme ͓␣-dependent
term in Eq. ͑27͒; reported in Table II between parentheses͔;
͑iv͒ the solvation free energies are larger in magnitude under
NPBC compared to PBC, a consequence of the periodicityinduced
solvent depolarization within the reference unit cell
͑Fig. 3͒; ͑v͒ the periodicity-induced perturbation of the solvation
free energy (NPBC→PBC) is larger for the smaller
ion, for the larger cut off, for the higher permittivity value,
and for RF compared to ST. The latter effect is a consequence
of the larger periodicity-induced solvent depolarization
within the reference unit cell for the RF scheme ͑Fig. 3͒.
The values reported in Table II are strongly cut-offdependent
and compare poorly with the corresponding Born
solvation free energies of Ϫ342.9 and Ϫ171.4 kJ molϪ1
(⑀sϭ78, qIϭ1 e, RIϭ0.2 or 0.4 nm͒ or Ϫ173.7 and
Ϫ86.8 kJ molϪ1
(⑀sϭ2, qIϭ1 e, RIϭ0.2 or 0.4 nm͒. As
discussed in Appendix B, these large discrepancies could be
reduced by the inclusion of a charge self-energy term into the
total electrostatic energy of the system. It is also suggested
that such a self-energy term should be systematically included
in the total electrostatic energy during molecular
simulations relying on effective cut-off-based electrostatic
interaction functions to ensure the obtension of meaningful
energies. In this context, a new deﬁnition ͓Eq. ͑B4͔͒ is proposed
for the electrostatic interaction energy in simulations
employing the Barker–Watts reaction-ﬁeld scheme.
The effect of artiﬁcial periodicity on the solvation free
energy of a spherical ion computed using cut-off-based ͑ST
or RF͒ electrostatics is illustrated in Fig. 4 for an ion of
charge qIϭ1 e, a solvent of permittivity ⑀sϭ78, and for
different values of the ionic radius RI and cut-off radius RC .
The relative periodicity-induced perturbation ␥(L) of the
FIG. 3. Radial polarization p(r) around a solvated spherical ion ͓Eqs. ͑22͒
or ͑25͔͒. The system consists of a single ion of charge qIϭ1 e and radius
RIϭ0.4 nm in a solvent of permittivity ⑀sϭ78, and is either nonperiodic
͑NPBC; spherical domain of radius Sϭ4.0 nm) or periodic ͑PBC; cubic unit
cell of edge Lϭ2.6 nm). Electrostatic interactions correspond to either the
ST ͑a͒ or RF ͑b͒ schemes, with cut-off radii RISϭRSSϭRCϭ1.2 nm. In
addition to the results of the 3D-Direct method ͑present article͒, the analytical
polarization corresponding to the Born model ͓CB/NPBC; Eq. ͑26͔͒ and
the lattice-sum case ͓LS/PBC; computed using the 3D-FFT method ͑Ref.
64͔͒, and the polarization computed from the 1D-Direct method ͑Refs. 27
and 28͒ for the speciﬁc interaction scheme ͑ST,RF/NPBC͒ are also presented
for comparison. The cut-off distance as well as the
͑half-͒box edge ͑PBC only͒ are indicated by arrows.
9137J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
ionic solvation free energy ͓Eq. ͑29͔͒ is displayed as a function
of L/2 in Fig. 4͑a͒. All curves converge to a limiting
value of one when L/2ӷRC , indicating that the solvation
free energy under PBC indeed converges to its NPBC value
when the computational box becomes large compared to the
cut-off radius. For example, the magnitude of ␥(L) for L/2
ϭ4 nm is smaller than 10Ϫ4
for all parameter combinations
considered. When L/2 is only moderately larger than RC ,
artiﬁcial periodicity causes a depolarization of the solvent
͑Fig. 3͒ and a decrease in the magnitude of the solvation free
energy. As a consequence, ␥(L) becomes negative. When
L/2ϭRC , the solvation free energy is reduced by 2%–9%
compared to its NPBC value for the parameter combinations
considered. In agreement with previous observations ͑Table
II͒, the relative periodicity-induced perturbation of the solvation
free energy ͑i͒ increases in magnitude with increasing
ionic radius; ͑ii͒ increases in magnitude with increasing cutoff
radius; ͑iii͒ is larger for the RF scheme compared to the
ST scheme.
The different curves in Fig. 4͑a͒ can be adequately represented
by exponential functions. This is shown in Fig. 4͑b͒,
where the quantity log10͓Ϫ(RC /RI)␥(L)͔ is displayed as a
function of L/(2RC). The resulting data can be ﬁtted by two
straight lines, corresponding to the ST and RF schemes, with
linear correlation coefﬁcients ͑over the interval L/(2RC)
ϭ0 to 2.5͒ of Ϫ0.9998 and Ϫ0.9985, respectively. Thus,
irrespective of RI and RC , the relative periodicity-induced
perturbation appears to be approximately of the form
␥͑L͒ϷϪ
RI
RC
10␮ L/͑2RC͒ ϩ␯
. ͑32͒
Based on all available data for ␧sϭ78, the constants in Eq.
͑32͒ evaluated to ␮ϭϪ2.203 and ␯ϭ1.290 for the ST
scheme, and ␮ϭϪ1.468 and ␯ϭ0.715 for the RF scheme.
This behavior can be contrasted to the case of nontruncated
electrostatic interactions. In this case, the solvation free
energy corresponding to CB/NPBC is given by the Born
expression86
͓Eq. ͑28͔͒. A corresponding analytical
expression57
has been derived for the LS/PBC case, namely
⌬Gsolv
LS/PBC
ϭϪ
qI
2
8␲⑀o
⑀sϪ1
⑀s
ͫL
RI
ϩ␰EW
ϩ
4␲
3 ͩRI
L ͪ2
Ϫ
16␲2
45 ͩRI
L ͪ5
ͬLϪ1
, ͑33͒
with ␰EWϷϪ2.837 297. Thus, it follows from Eq. ͑29͒ that:
␥͑L͒ϭ
RI
L ͫ␰EWϩ
4␲
3 ͩRI
L ͪ2
Ϫ
16␲2
45 ͩRI
L ͪ5
ͬ. ͑34͒
In this case, the evolution of ␥(L) towards zero when L
ӷRI is in LϪ1
, i.e., much slower than the exponential
distance-dependence observed for cut-off-based schemes
͓Eq. ͑32͔͒. For example, for an ion of radius RIϭ0.4 nm,
␥(L) evaluates to Ϫ0.95 for L/2ϭ0.4 nm, Ϫ0.19 for L/2
ϭ3.0 nm, and is above Ϫ0.1 for L/2Ͼ5.7 nm ͓compare
with the smaller magnitude and faster relaxation observed
for cut-off-based schemes in Fig. 4͑a͔͒. This shows that the
application of a cut-off in the computation of ionic solvation
free energies by explicit-solvent simulation dramatically reduces
the system-size dependence of the calculated solvation
free energies compared to lattice-sum methods.19,55,57
More
generally, cut-off truncation ͑with the possible inclusion of a
reaction-ﬁeld correction͒ efﬁciently reduces the impact of
ﬁnite-size effects and artiﬁcial periodicity on the energies
and forces in any molecular dynamics simulation. However,
this is at the expense of introducing other ͑potentially more
harmful͒ artifacts related with the cut-off truncation itself.
C. Interaction between two spherical ions
The electrostatic solvation free energy proﬁles
⌬Gsolv(d) for a pair of monovalent spherical ions ͑same or
opposite charges, identical radii of 0.4 nm͒ in a solvent of
permittivity ⑀sϭ78 are displayed in Fig. 5 as a function of
the interionic distance d for different choices of boundary
conditions and treatments of the electrostatic interactions
based on a single cut-off RCϭ1.2 nm. The PBC curves correspond
to ions aligned along an axis of a cubic unit cell of
edge Lϭ6 nm. The corresponding proﬁles computed using
the 3D-FFT method60
͑ST,RF/PBC͒ are also displayed for
comparison.
As expected, the curves corresponding to the NPBC case
present a minimum ͑maximum͒ at ionic contact for ions of
identical ͑opposite͒ charges, and asymptotically converge to
a common value for a given scheme. More precisely, in the
limit of large interionic distances ͑isolated ions͒, ⌬Gsolv(d)
converges towards twice the solvation free energy ⌬Gsolv
ion
of
a single ion (⌬Gsolv
ion
ϭϪ129.4 or Ϫ86.4 kJ molϪ1
for the ST
or RF schemes; see Table II͒. In the limit d→0 ͑superimFIG.
4. Periodicity-induced perturbation of the solvation free energy of a
spherical ion. The system consists of a single ion of charge qIϭ1 e and
radius RI in a cubic periodic box of edge L ﬁlled by a solvent of permittivity
⑀sϭ78. Electrostatic interactions correspond to either the ST or RF
schemes, with cut-off radii RISϭRSSϭRC . ͑a͒ Relative periodicity-induced
shift ␥(L) in the solvation free energy ͓Eq. ͑29͔͒, displayed as a function of
L/2. ͑b͒ Logarithm of minus ␥(L) ampliﬁed by RC /RI , displayed as a
function of L/2RC . The dashed lines corresponding to a least-squares ﬁt
͓over the interval L/(2RC)ϭ0–2.5] corresponding to either the RF or the
ST schemes.
9138 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
posed ions͒, ⌬Gsolv(d) converges towards four times ⌬Gsolv
ion
͑zero͒ for ions of identical ͑opposite͒ charges.
Under PBC, symmetry and periodicity constraints impose
that the proﬁle possesses a stationary point at dϭL/2
and is symmetrical with respect to this point. This induces a
difference between the NPBC and PBC curves for d in the
interval ͓0;L/2͔. Within this range, the perturbation caused
by the introduction of periodicity is attractive ͑repulsive͒ for
ions of identical ͑opposite͒ charges. However, the magnitude
of the effect is very small. For example, the differences between
the values of ⌬Gsolv(L/2) under NPBC and PBC is
only about 0.45 kJ molϪ1
in magnitude for all cases considered.
Thus, in contrast to the case of lattice-sum methods,57
artiﬁcial periodicity has very little inﬂuence on the solvation
free energy proﬁle for pairs of small monovalents ions in a
solvent of high permittivity when cut-off truncation is applied.
Finally, the agreement between the present calculations
employing the 3D-Direct method under PBC and the
3D-FFT method60
is quite good, especially for the ST
scheme. For the RF scheme, the agreement is slightly worse,
probably due to small differences in the application of
boundary smoothing at the ion surface and ion–solvent cutoff
distance.60
The corresponding proﬁles for the overall electrostatic
contribution ⌬Gel(d) to the potential of mean force ͓Eq.
͑30͔͒ are displayed in Fig. 6 as a function of the interionic
distance d. At short distances, the curves corresponding to
NPBC and PBC are nearly identical. In the ST case, the
curves present minima ͑maxima͒ at contact and at the cut-off
distance for ions of identical ͑opposite͒ charges. The presence
of an extremum at the cut-off distance is clearly an
artifact related to the use of truncated electrostatic interactions.
These proﬁles provide an explanation for a number of
observations made in explicit-solvent simulations employing
the ST scheme: ͑i͒ The tendency for ion pairs with like
charges to cluster at distances close to the cut-off
distance;9,96
͑ii͒ the tendency for ion pairs of opposite
charges to avoid distances close to the cut-off
distance;9,42,43,96
͑iii͒ the artiﬁcially increased38,40,41,57
stability
of contact ion pairs for ions of like charges;36,37,39,41,87
͑iv͒ the artiﬁcially decreased stability of contact pairs for
ions of opposite charges.42,44,80
In the RF case, the proﬁles nearly present the expected
behavior, namely repulsion ͑attraction͒ for ions of like ͑opposite͒
charges, except for a signiﬁcant artifact in the neighborhood
of the cut-off distance. For ions of like charges, a
spurious minimum occurs just below the cut-off distances,
while for ions of opposite charges, a spurious maximum occurs
just above the cut-off distance. Although these artifacts
might affect the populations of contact pairs for ions of like
or opposite charges in simulations of ionic solution,88,89
the
magnitude of these artifacts is limited compared to the ST
scheme, in agreement with previous observations.44,90
At large distances from the ion, and for both the ST and
RF schemes, the NPBC proﬁles tend to be close to the expected
Coulombic limit ͓Eq. ͑31͔͒. However, the exact agreement
is difﬁcult to assess since positive deviations occur,
which are probably related to the limited size of the computational
domain.
V. CONCLUSION
In the present study, continuum electrostatics was used
to investigate the nature and magnitude of the perturbations
induced by cut-off truncation and artiﬁcial periodicity in
FIG. 5. Electrostatic solvation free energy proﬁle ⌬Gsolv(d) for a pair of
monovalent spherical ions. The system consists of two ions of radii RI
ϭ0.4 nm bearing identical ͑a and c͒ or opposite ͑b and d͒ charges, and
separated by a distance d in a medium of permittivity ⑀sϭ78. It is either
nonperiodic ͑NPBC; spherical domain of radius Sϭ4.0 nm) or periodic
͑PBC; cubic unit cell of edge Lϭ6 nm; ions aligned with an axis of the unit
cell͒. Electrostatic interactions correspond to either the ST ͑a and b͒ or RF ͑c
and d͒ schemes, with cut-off radii RISϭRSSϭRCϭ1.2 nm. In addition to the
results of the 3D-Direct method ͑present article͒, the solvation free energies
computed from the 3D-FFT method ͑Refs. 51 and 60͒ for the speciﬁc interaction
scheme ͑ST,RF/PBC͒ and the same value of L are also presented for
comparison.
FIG. 6. Electrostatic contribution ⌬Gel(d) to the potential of mean force for
a pair of monovalent spherical ions. The system consists of two ions of radii
RIϭ0.4 nm bearing identical ͑a and c͒ or opposite ͑b and d͒ charges, and
separated by a distance d in a medium of permittivity ⑀sϭ78. It is either
nonperiodic ͑NPBC; spherical domain of radius Sϭ4.0 nm) or periodic
͑PBC; cubic unit cell of edge Lϭ6 nm; ions aligned with an axis of the unit
cell͒. Electrostatic interactions correspond to either the ST ͑a and b͒ or RF ͑c
and d͒ schemes, with cut-off radii RISϭRSSϭRCϭ1.2 nm. The ideal longrange
limit ⌬Gel
lr
(d) is also presented for comparison ͓Eq. ͑31͔͒. It is calculated
using ⌬Gsolv
ion
ϭϪ129.4 (ST) or Ϫ86.4 (RF) kJ molϪ1
͑Table II͒.
9139J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
explicit-solvent simulations of ions and ion pairs. To this
purpose, a new algorithm relying on ﬁnite integration was
developed to solve the equations of continuum electrostatics
based on truncated ͑and possibly reaction-ﬁeld corrected͒
solute–solvent and solvent–solvent electrostatic interactions,
under either nonperiodic ͑NPBC͒ or periodic ͑PBC͒ boundary
conditions. This algorithm was tested and validated by
comparison with available methods ͑Table I͒ whenever pos-
sible.
In the context of the solvation of a single spherical ion,
the main observations can be summarized as follows:
͑A͒ The application of cut-off truncation ͑under NPBC͒
signiﬁcantly affects the solvent polarization around a
spherical ion ͑compared to the ideal Born result͒. With
straight truncation ͑ST͒ of the interactions, the solvent
is overpolarized within the cut-off sphere of the ion and
underpolarized outside this sphere. When a reactionﬁeld
͑RF͒ correction is applied, the agreement with the
Born ͑NPBC͒ or lattice-sum ͑PBC͒ polarization is signiﬁcantly
improved, the deviations being largest in the
neighborhood of the cut-off distance.
͑B͒ The introduction of artiﬁcial periodicity (NPBC
→PBC) when applying cut-off-based electrostatics
leads to a depolarization of the solvent around the ion
in the reference unit cell. This effect is caused by the
indirect ͑solvent-mediated͒ perturbation of the solvent
molecules in this reference cell by the periodic copies
of the ion. The depolarization is more signiﬁcant for
the RF scheme compared to the ST scheme.
͑C͒ The application of cut-off truncation ͑under NPBC͒ decreases
the magnitude of the ionic solvation free energy
of a spherical ion ͑compared to the ideal Born
result͒. The magnitude of this effect is highly sensitive
to the electrostatic scheme ͑ST or RF͒ and to the choice
of a cut-off radius. However, as discussed in Appendix
B, the problem could be largely ͑though approximately͒
remedied in explicit-solvent simulations by the
systematic inclusion of an appropriate self-energy term
in the total electrostatic energy of the system. Alternatively,
an exact correction term to ionic solvation free
energies computed from explicit-solvent simulations
can be obtained by the application of the present
continuum-electrostatics method under NPBC or of its
one-dimensional analog.27,28
͑D͒ The introduction of artiﬁcial periodicity (NPBC
→PBC) when applying cut-off-based electrostatics
causes a further decrease in the magnitude of the ionic
solvation free energy. In contrast to lattice-sum
methods,57
where this free-energy shift is important
even for relatively large system sizes ͑proportional to
LϪ1
, L being the edge length of the cubic unit cell͒, the
effect decays rapidly with increasing system sizes ͑proportional
to RϪ1
exp(ϪcL/R), R being the cut-off distance͒
in the case of cut-off-based electrostatics. Here
also, an approximate correction term was derived ͓Eq.
͑34͔͒ that can be applied to correct ionic solvation free
energies evaluated from explicit-solvent simulations
under PBC.
The relevance of these observations can be appreciated
by recalling the ﬁve sources of error related to the computation
of ionic solvation free energies from explicit-solvent
simulations relying on cut-off-based electrostatic interactions:
͑i͒ Incorrect polarization of the solvent around the ion
due to truncated ͑and possibly modiﬁed͒ electrostatic interactions;
͑ii͒ Cut-off- and system size-dependent perturbation
of the solvent polarization due to artiﬁcial periodic boundary
conditions; ͑iii͒ artifacts at the cut-off distance arising from
the ﬁnite size of the solvent molecules, and related to the use
of either a molecular or an atomic cut-off; ͑iv͒ artiﬁcial heating
during molecular-dynamics simulations due to the possible
presence of discontinuities in the atomic forces; ͑v͒
inaccuracy of the ion–solvent and solvent–solvent interaction
functions and parameters ͑force ﬁelds͒.
Only with the understanding of these ﬁve sources of error
and the design of appropriate correction schemes will it
be possible to obtain accurate ionic solvation free energies
from explicit-solvent simulations. The discussion ͑and correction͒
of the two ﬁrst sources of error has been the focus of
the present article. The third problem has been previously
discussed by a number of groups.55,72–75,77,78,91
Due to the
ﬁnite size of solvent molecules, the solvent-generated potential
at the ion site ͑and thus the solvation free energy͒ may
vary considerably depending on whether cut-off truncation is
applied on an atomic or on a molecular basis ͑and in the
latter case, on the choice of a molecular center͒. Although
the debate is not yet completely settled, a number of arguments
suggest that molecular truncation ͑based on the center
of van der Waals interactions for a spherical molecule77
͒ represents
the appropriate method for evaluating the solventgenerated
potential at the ion site,74,77,91
while a ͑generally
sizeable͒ correction term must be applied if atomic truncation
͑or a lattice-sum method͒ is employed instead. The
fourth problem, namely the artiﬁcial heating of molecules at
distances close to the cut-off radius, may be alleviated by the
use of an effective truncated electrostatic interaction that is
continuously differentiable at the cut-off distance, together
with atomic truncation. For example, the Barker–Watts
reaction-ﬁeld interaction14–16
causes very limited heating
provided that the permittivity of the solvent considered is
high and that an ͑unusual͒ atomic-cut-off implementation is
applied.16
Finally, in regard to the ﬁfth problem, it should be
stressed that the derivation of force-ﬁeld parameters for ion–
solvent interactions based on experimental ionic solvation
free energies makes little sense before the four other problems
are solved ͑i.e., methodology-independent ionic solvation
free energies can be be obtained from explicit-solvent
simulations͒.
In the context of the potential of mean force for the
interaction between two spherical ions, the main observations
can be summarized as follows:
͑A͒ The application of cut-off truncation ͑under NPBC͒ induces
serious artifacts in the overall electrostatic contribution
to the potential of mean force for the interaction
between two spherical ions. As previously
observed in explicit-solvent simulations, these lead to
spurious features in the radial distribution functions
9140 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
close to the cut-off distance9,42,43
and to artifacts in the
relative stability of the contact, solvent-separated and
free ion pairs.36–42,44,57,80,87
These effects are reduced
͑although not compleletly eliminated͒ by the application
of a reaction-ﬁeld correction.
͑B͒ The introduction of artiﬁcial periodicity (NPBC
→PBC) when applying cut-off-based electrostatics appears
to cause very small changes in the electrostatic
contribution to the ͑minimum-image͒ potentials of
mean force for small monovalent ions in a solvent of
high permittivity. A rather weak periodicity-induced
perturbation was also reported in this case for latticesum
methods.57
However, the causes of the limited
overall effect are different. In the cut-off case, both the
solvation free energy proﬁle and the direct ion–ion interaction
are almost unaffected by periodicity. In the
lattice-sum case, both contributions are largely affected,
but the two perturbations nearly cancel each
other.
Explicit-solvent simulations of ion solvation and ion–
ion interactions are currently in progress to conﬁrm the validity
of the above considerations derived from a continuumelectrostatics
analysis, and their compatibility with the
results of simulations employing lattice-sum methods.
ACKNOWLEDGMENTS
The Swiss National Science Foundation ͑Project 2100-
061939͒ is gratefully acknowledged for ﬁnancial support, as
well as Peter Gee for his careful reading of the manuscript
and helpful suggestions.
APPENDIX A: POLARIZATION IN THE
BARKER–WATTS REACTION-FIELD SCHEME
As shown in Fig. 3, the polarization corresponding to the
RF scheme under NPBC is very close to the corresponding
Born polarization, the deviation being largest at distances
close to the cut-off radius. Here, we derive a number of
results related to this comparison, in the more general context
of the Barker–Watts ͑BW͒ interaction function ͓Eq. ͑2͔͒,
namely that ͑i͒ the BW/NPBC and CB/NPBC ͑Born͒ polarizations
become identical in the limit RIS ,RSS→ϱ, irrespective
of the value of ␣; ͑ii͒ when RISϭRSS , the BW/NPBC
polarization converges towards the Born polarization at short
distances from the ion, provided that ␣ϭ2(⑀sϪ1)/(2⑀s
ϩ1); ͑iii͒ the BW/NPBC polarization, just as the Born polarization,
is proportional to r2
at large distances from the
ion, provided that ␣ϭ(⑀sϩ2)/(⑀sϪ1). The latter results obviously
remain approximately valid for a solvent of high
permittivity (⑀sӷ1) when ␣ is set to one ͑RF͒ or close to
one.
The derivations are based on continuum-electrostatics
results presented in a previous article,27,28
and applied to the
speciﬁc case of truncated electrostatic interactions corresponding
to the Barker–Watts potential ͓Eq. ͑2͔͒. In this case
the radial polarization around a solvated spherical ion is a
solution of the integral equation
p͑r͒ϭg͑r͒ϩ
⑀sϪ1
⑀s
͵max(RI ,͉rϪRSS͉)
rϩRSS
K͑r,rЈ͒p͑rЈ͒drЈ
for rϾRI , ͑A1͒
together with p(r)ϭ0 for rрRI . The inhomogeneous term
g(r) is related to the vacuum ﬁeld generated by the ion
g͑r͒ϭ
qI
4␲
⑀sϪ1
⑀s
H͑RISϪr͒ͩrϪ2
Ϫ
␣r
RIS
3 ͪ, ͑A2͒
and the kernel K(r,rЈ) to the form of the solvent–solvent
interactions
K͑r,rЈ͒ϭϪ͑␣ϩ2͒
r4
ϩ͓RSS
2
Ϫ͑rЈ͒2
͔2
Ϫ2r2
͓RSS
2
ϩ͑rЈ͒2
͔
16r2
RSS
3 .
͑A3͒
First, it is shown that in the limit RIS ,RSS→ϱ, the polarization
deﬁned by Eq. ͑A1͒ converges to the Born polarization
͓Eq. ͑26͔͒. Due to the form of Eq. ͑A2͒,
limRIS→ϱ͓g(r)ϪpBorn(r)͔ϭ0. It is thus sufﬁcient to prove
that the integral term in Eq. ͑A1͒ vanishes in the limit of an
inﬁnite solvent–solvent cut-off radius. In this limit, the lower
bound of integration can be set to RSSϪr. For rЈ within the
interval RSSϪr to RSSϩr, K(r,rЈ) is positive and possesses
a single maximum at r˜Јϭ(RSS
2
ϩr2
)1/2
with K(r,r˜Ј)ϭ(␣
ϩ2)/(4RSS). Assuming that the polarization is positive and
ﬁnite over the whole range of distance ͑with a maximum
value p˜), upper and lower bounds can be given to the integral
in Eq. ͑A1͒
0р ͵RSSϪr
RSSϩr
K͑r,rЈ͒p͑rЈ͒drЈ
р2rp˜K͑r,r˜Ј͒ϭ
␣ϩ2
2
p˜
r
RSS
, ͑A4͒
which shows that the integral vanishes for any ﬁnite r in the
limit RSS→ϱ. Thus, the BW/NPBC polarization converges
to the Born polarization in the limit of large cut-off radii,
irrespective of the value of ␣. This result is in particular valid
for the ST27
(␣ϭ0) and RF (␣ϭ1) schemes.
As a second result, it is shown that the polarization deﬁned
by Eq. ͑A1͒ converges to the Born polarization in the
limit of short distances when RISϭRSS , provided that ␣
ϭ2(⑀sϪ1)/(2⑀sϩ1), i.e., when Eq. ͑3͒ is used with ⑀Ј
ϭ⑀s ͑adjusted boundary conditions61
͒. For rрmin͓RIS ,RSS
ϪRI͔, the Heaviside function in Eq. ͑A2͒ can be omitted and
the lower bound in Eq. ͑A1͒ replaced by RSSϪr. Thus, one
looks for a solution of
p͑r͒ϭ
⑀sϪ1
⑀s
ͫqI
4␲ ͩrϪ2
Ϫ
␣r
RIS
3 ͪ
ϩ ͵RSSϪr
RSSϩr
K͑r,rЈ͒p͑rЈ͒drЈͬ. ͑A5͒
Using the result
͵RSSϪr
RSSϩr
K͑r,rЈ͒͑rЈ͒Ϫ2
drЈϭ
␣ϩ2
3
r
RSS
3 , ͑A6͒
9141J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
one shows easily that the function satisfying Eq. ͑A5͒ when
RISϭRSS and ␣ϭ2(⑀sϪ1)/(2⑀sϩ1) is the Born polarization
͓Eq. ͑26͔͒. Therefore, with this speciﬁc value of ␣, the
BW/NPBC polarization will converge to the Born polarization
at short distances, irrespective of the cut-off values. For
a solvent of high permittivity (⑀sӷ1), this result will remain
approximately valid when ␣ϭ1 ͑RF͒.
As a third result, it is shown that the polarization deﬁned
by Eq. ͑A1͒ possesses an rϪ2
dependence in the limit of
large distances ͑just as the Born polarization͒, provided that
␣ϭ(⑀sϩ2)/(⑀sϪ1). For rуmax͓RIS ,RSSϩRI͔, Eq. ͑A2͒
implies g(r)ϭ0 and the lower bound in Eq. ͑A1͒ can be
replaced by rϪRSS . Thus, one looks for a solution of
p͑r͒ϭ
⑀sϪ1
⑀s
͵rϪRSS
rϩRSS
K͑r,rЈ͒p͑rЈ͒drЈ. ͑A7͒
Using the result
͵rϪRSS
rϩRSS
K͑r,rЈ͒͑rЈ͒Ϫ2
drЈϭ
␣ϩ2
3
rϪ2
, ͑A8͒
one shows easily that the function satisfying Eq. ͑A7͒ when
␣ϭ(⑀sϩ2)/(⑀sϪ1) is crϪ2
, where c is a constant. Therefore,
with this speciﬁc value of ␣, the BW/NPBC polarization
will possess a rϪ2
dependence in the limit of large distances.
For solvent of high permittivities (⑀sӷ1), this result
will remain approximately valid when ␣ϭ1 ͑RF͒ or when ␣
is given by Eq. ͑3͒ with ⑀Јϭ⑀s ͑adjusted boundary condi-
tions͒.
The above observations are illustrated in Fig. 7 for the
case of an ion of radius RIϭ0.4 nm immersed in a solvent of
permittivity ⑀sϭ78. In the inset, the polarization p(r) is displayed
as a function of rϪ2
for the RF scheme (␣ϭ1), using
four different values of the cut-off radius RC . The Born polarization
is linear in rϪ2
and corresponds to the diagonal of
the graph. Besides a small artifact close to the cut-off distance
͑only visible in the graph for RCϭ1.2 and 1.6 nm͒ and
a slight offset ͑only visible in the graph for RCϭ1.2), the
curves are nearly indistinguishable from each other and from
the Born polarization. In the main graph, the differences between
p(r) and pBorn(r) are displayed as a function of rϪ2
for the same values of RC . In the short distance limit ͑right
side of the graph͒, convergence towards the Born polarization
is evident although slow. However, the differences are of
rather small magnitude. For example, the relative difference
at the surface of the ion, ͓pBorn(RI)Ϫp(RI)͔/pBorn(RI),
evaluates to 0.91, 0.29, 0.12, or 0.06% for RCϭ1.2, 1.6, 2.0
or 2.4 nm. In the long distance limit ͑left side of the graph͒,
the approximate rϪ2
evolution of p(r)ϪpBorn(r) can also be
observed. The maximal error in the polarization occurs exactly
at the cut-off distance. Increasing the cut-off distance
rapidly reduces the magnitude of this error, the difference
pBorn(RC)Ϫp(RC) evaluating to Ϫ0.021, Ϫ0.009, Ϫ0.004,
Ϫ0.003 e nmϪ2
for RCϭ1.2, 1.6, 2.0, or 2.4 nm. These
results indicate that for large-enough cut-off distances, the
RF/NPBC scheme provides an essentially correct representation
of the polarization around a spherical ion, in both the
short- and long-distance ranges.
APPENDIX B: SELF-ENERGY TERM
FOR CUT-OFF-BASED INTERACTION FUNCTIONS
Here it is shown that when an effective cut-off-based
interaction function is used to handle electrostatic interactions
in an explicit-solvent simulation, a charge self-energy
term should be included into the total electrostatic energy of
the system to ensure a fast convergence of ionic solvation
free energies towards the Born result in the limit of large
cut-off distances. Generalizing this observation to the case of
more complex molecular systems, a new deﬁnition ͓Eq.
͑B4͔͒ is proposed for the electrostatic interaction energy in
simulations employing the Barker–Watts reaction-ﬁeld
scheme.
Consider an effective cut-off-based electrostatic interaction
function where the potential generated at r by a unit
charge at the origin is given by
␺͑r͒ϭ
1
4␲⑀o
H͑RϪr͓͒rϪ1
ϩ␺˜ ͑r͔͒, ͑B1͒
where R is the cut-off radius, chosen smaller than half the
smallest dimension of the computational box. It is further
assumed that ͑i͒ the interaction function vanishes at rϭR,
i.e., ␺(R)ϭ0, ͑ii͒ ␺˜ is a sum of terms of the form RϪlϪ1
rl
with lу0, and ͑iii͒ the polarization around a spherical ion
converges to the Born polarization ͓Eq. ͑26͔͒ in the limit of
an inﬁnite cut-off distance. For example, as shown in Appendix
A, the Barker–Watts interaction function ͓Eq. ͑2͔͒ satisﬁes
the three conditions irrespective of the value of ␣. In
fact, there is some hint that the form of Eq. ͑B1͒ and the
second condition automatically imply the third one.60
FIG. 7. Radial polarization p(r) around a solvated spherical ion ͓Eqs. ͑22͔͒
compared to the Born polarization ͓Eq. ͑26͔͒. The system consists of a
single ion of charge qIϭ1 e and radius RIϭ0.4 nm in a solvent of permittivity
⑀sϭ78 under NPBC ͑spherical domain of radius Sϭ10.0 nm). Electrostatic
interactions correspond to the RF scheme, with cut-off radii RIS
ϭRSSϭRCϭ1.2, 1.6, 2.0, or 2.4 nm. The polarization is computed using the
1D-Direct method ͑Refs. 27 and 28͒ for the speciﬁc interaction scheme
͑RF/NPBC͒. In the inset, p(r) is displayed as a function of rϪ2
, while the
Born polarization pBorn(r) ͑not displayed͒ is a straight line corresponding to
the diagonal of the graph. In the main graph, the difference p(r)
ϪpBorn(r) is displayed as a function of rϪ2
.
9142 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
Under these assumptions and for large-enough cut-off
distances, the polarization within the cut-off sphere of an ion
should be quite close to the Born polarization ͓see Fig. 3͑b͒
for the RF scheme͔. Using the Born polarization as an approximation
in this case, one may estimate the corresponding
ionic solvation free energy. Combining Eq. ͑10͒ with V(r)
ϭϪqIٌ␺(r) and Eq. ͑26͒, one obtains the approximate ex-
pression
⌬GsolvϷ
qI
2
8␲
⑀sϪ1
⑀s
͵RI
R
dr4␲
r
r
•ٌ␺͑r͒
ϭ
qI
2
2
⑀sϪ1
⑀s
͓␺͑R͒Ϫ␺͑RI͔͒
ϭ⌬GBornϪ
qI
2
8␲⑀o
⑀sϪ1
⑀s
␺˜ ͑RI͒, ͑B2͒
where we used ␺(R)ϭ0 and inserted Eq. ͑28͒. In the RF
case, estimates based on this equation can be compared to
the data in Table II. For example, for ⑀sϭ78, qIϭ1 e, RI
ϭ0.2 nm and RISϭRSSϭ1.2 nm, Eq. ͑B2͒ ͓using Eq. ͑2͒
with ␣ϭ1] gives an estimate of Ϫ258.0 kJ molϪ1
͑including
27.8 kJ molϪ1
for the ␣-dependent contribution͒, to be
compared with the numerical value of Ϫ259.2 kJ molϪ1
͑including
26.2 kJ molϪ1
for the ␣-dependent contribution͒ in
Table II. If ␺˜ is a sum of terms of the form RϪlϪ1
rl
with l
у0, ␺˜ (RI) can be approximated by ␺˜ (0) in the limit of large
cut-off radii and small ions. In this case, Eq. ͑B2͒ shows that
the ionic solvation free energy computed from an explicitsolvent
simulation employing a cut-off-based interaction
function will converge signiﬁcantly faster towards the Born
result upon increasing the cut-off distance if a self-energy
term of the form
⌬Gselfϭ
qI
2
8␲⑀o
⑀sϪ1
⑀s
␺˜ ͑0͒, ͑B3͒
is included in the electrostatic energy of the system. Note
that ⌬Gself converges towards zero in the limit R→ϱ. However,
because this term is generally large and converges only
as RϪ1
, its inclusion makes a signiﬁcant difference even for
relatively large cut-off radii.
Generalizing this observation to more complex molecular
systems suggests that a charge self-energy term should be
included in explicit-solvent molecular dynamics simulations
employing any effective cut-off-based electrostatic interaction
function. Intuitively, this term may be interpreted as the
reversible work required to individually charge the atoms
when they are at inﬁnite separation. This work excludes the
͑inﬁnite͒ Coulombic self-energy, but retains the contribution
arising from the non-Coulombic term associated with ␺˜ in
Eq. ͑B1͒.
In the speciﬁc case of the Barker–Watts reaction-ﬁeld
method ͓Eq. ͑2͔͒, a reasonable expression for the total electrostatic
energy UBW could be
UBWϭ͚i
͚jϾi,j excl(i)
qiqj␺BW͑r¯ij͒
ϩ
1
4␲⑀o
͚ͭi
͚jϾi,j෈excl(i)
qiqj␺˜ BW͑r¯ij͒
ϩ
1
2
␺˜ BW͑0͚͒ͫi
qi
2
Ϫ⑀s
Ϫ1
͚ͩi
qiͪ2
ͬͮ, ͑B4͒
where r¯ij is the minimum-image vector corresponding to rij ,
excl(i) denotes the exclusion list of atom i ͑the distance
between any two excluded atoms is assumed to be smaller
than R), and
␺˜ BW͑r͒ϭ
␣r2
2R3 Ϫ
␣ϩ2
2R
, ͑B5͒
with ␣ deﬁned by Eq. ͑3͒. Note that current simulation programs
͑e.g., GROMOS
92
and GROMACS
93
͒ typically restrict the
implementation of the Barker–Watts reaction-ﬁeld method to
the ﬁrst term in Eq. ͑B4͒. The second term is explicitly included
here because so-called excluded neighbors ͑usually
ﬁrst and second covalent neighbors͒ should only be excluded
from the summation of the Coulombic (rϪ1
) contribution,
but not of the reaction-ﬁeld (␺˜ BW) contribution. The form of
the third term has been chosen for consistency in the context
of small molecules ͑compared to the cut-off radius and unitcell
size͒. For a small molecule ͑or ion͒ gathered by periodicity
around its center, r¯ij can be replaced by rij in Eq. ͑B4͒
and the Heaviside function involved in ␺BW can be omitted.
In this case, the reaction-ﬁeld ͑non-Coulombic͒ contribution
contribution to UBW can be written
UBW
rf
ϭ
1
8␲⑀o
ͭ␣
R3 ͚i
͚jϾi
qiqjrij
2
ϩ
⑀sϪ1
⑀s
␺˜ BW͑0͚͒ͩi
qiͪ2
ͮ.
͑B6͒
For a neutral molecule, one has
UBW
rf,dip
ϭϪ
1
8␲⑀o
␣␮2
R3 , ͑B7͒
where ␮ is the molecular dipole moment, which matches the
Onsager expression for a dipole in a spherical cavity94
provided
that adjusted boundary conditions ͓⑀Јϭ⑀s in Eq. ͑3͔͒
are applied. For a monoatomic ion, one has UBW
rf,ion
ϭ⌬Gself ,
i.e., the self-energy term suggested by Eq. ͑B3͒. Note that the
last term in Eq. ͑B4͒ only affects the energy of the system,
but does not induce atomic forces. However, it may be essential
to include it in free-energy calculations involving alterations
of the atomic partial charges.
In the speciﬁc case of a single ion, the inclusion of such
a self-energy term should substantially reduce the error on
ionic solvation free energies computed from explicit-solvent
simulations with ﬁnite cut-off distances. This can be seen
from the corresponding corrected values ⌬Gsolv
corr
ϭ⌬Gsolv
ϩ⌬Gself reported in Table II. In the RF case, taking the same
example as above (⑀sϭ78, qIϭ1 e, RIϭ0.2 nm and RIS
ϭRSSϭ1.2 nm), ⌬Gsolv
Born
evaluates Ϫ342.9 kJ molϪ1
, to be
compared with an estimate ⌬Gsolv
corr
of Ϫ344.9 kJ molϪ1
. The
corresponding estimate for the ST case, Ϫ363.4 kJ molϪ1
, is
signiﬁcantly less accurate. This is probably due to the poorer
9143J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Cut-off and periodicity effects in simulations of ions
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
agreement between the polarization within the cut-off sphere
of an ion and the Born polarization in this case ͓see Fig.
3͑a͔͒.
1
M. Karplus and J. A. McCammon, Nat. Struct. Biol. 9, 646 ͑2002͒.
2
X. Daura, A. Gla¨ttli, P. Gee, C. Peter, and W. F. van Gunsteren, Adv.
Protein Chem. 62, 341 ͑2002͒.
3
T. E. Cheatham, III and P. A. Kollman, Annu. Rev. Phys. Chem. 51, 435
͑2000͒.
4
L. Saiz, S. Bandyopadhyay, and M. L. Klein, Biosci Rep. 22, 151 ͑2002͒.
5
S. B. Engelsen, C. Monteiro, C. Herve´ de Penhoat, and S. Pe´rez, Biophys.
Chem. 93, 103 ͑2001͒.
6
W. F. van Gunsteren and H. J. C. Berendsen, Angew. Chem., Int. Ed. Engl.
29, 992 ͑1990͒.
7
M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids ͑Oxford
University Press, New York, 1987͒.
8
B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan,
and M. Karplus, J. Comput. Chem. 4, 187 ͑1983͒.
9
C. L. Brooks III, B. M. Pettitt, and M. Karplus, J. Chem. Phys. 83, 5897
͑1985͒.
10
K. S. Kim, Chem. Phys. Lett. 156, 261 ͑1989͒.
11
R. J. Loncharich and B. R. Brooks, Proteins: Struct., Funct., Genet. 6, 32
͑1989͒.
12
P. J. Steinbach and B. R. Brooks, J. Comput. Chem. 15, 667 ͑1994͒.
13
J. Norberg and L. Nilsson, Biophys. J. 79, 1537 ͑2000͒.
14
J. A. Barker and R. O. Watts, Mol. Phys. 26, 789 ͑1973͒.
15
I. G. Tironi, R. Sperb, P. E. Smith, and W. F. van Gunsteren, J. Chem.
Phys. 102, 5451 ͑1995͒.
16
P. H. Hu¨nenberger and W. F. van Gunsteren, J. Chem. Phys. 108, 6117
͑1998͒.
17
G. Hummer and D. M. Soumpasis, Phys. Rev. E 49, 591 ͑1994͒.
18
G. Hummer and D. M. Soumpasis, J. Phys.: Condens. Matter 6, A141
͑1994͒.
19
G. Hummer, L. R. Pratt, and A. E. Garcia, J. Phys. Chem. 100, 1206
͑1996͒.
20
P. P. Ewald, Ann. Phys. ͑Leipzig͒ 64, 253 ͑1921͒.
21
R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles
͑Institute of Physics, Bristol, 1981͒.
22
T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 ͑1993͒.
23
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G.
Pedersen, J. Chem. Phys. 103, 8577 ͑1995͒.
24
M. Neumann, Mol. Phys. 50, 841 ͑1983͒.
25
M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. 52, 97 ͑1984͒.
26
R. H. Wood, J. Chem. Phys. 103, 6177 ͑1995͒.
27
N. A. Baker, P. H. Hu¨nenberger, and J. A. McCammon, J. Chem. Phys.
110, 10679 ͑1999͒.
28
N. A. Baker, P. H. Hu¨nenberger, and J. A. McCammon, J. Chem. Phys.
113, 2510 ͑1999͒.
29
C. L. Brooks III, J. Phys. Chem. 90, 6680 ͑1986͒.
30
C. L. Brooks III, J. Chem. Phys. 86, 5156 ͑1987͒.
31
J. D. Madura and B. M. Pettitt, Chem. Phys. Lett. 150, 105 ͑1988͒.
32
T. P. Straatsma and H. J. C. Berendsen, J. Chem. Phys. 89, 5876 ͑1988͒.
33
S. G. Kalko, G. Sese, and J. A. Padro, J. Chem. Phys. 104, 9578 ͑1996͒.
34
H. Resat and J. A. McCammon, J. Chem. Phys. 104, 7645 ͑1996͒.
35
H. Resat and J. A. McCammon, J. Chem. Phys. 108, 9617 ͑1998͒.
36
B. M. Pettitt and P. J. Rossky, J. Chem. Phys. 84, 5836 ͑1986͒.
37
L. X. Dang and B. M. Pettitt, J. Phys. Chem. 94, 4303 ͑1990͒.
38
E. Gua`rdia, R. Rey, and J. A. Padro´, J. Chem. Phys. 95, 2823 ͑1991͒.
39
L. X. Dang, B. M. Pettitt, and P. J. Rossky, J. Chem. Phys. 96, 4046
͑1992͒.
40
G. Hummer, D. M. Soumpasis, and M. Neumann, Mol. Phys. 81, 1155
͑1994͒.
41
J. S. Bader and D. Chandler, J. Phys. Chem. 96, 6423 ͑1992͒.
42
R. A. Friedman and M. Mezei, J. Chem. Phys. 102, 419 ͑1994͒.
43
G. S. Del Buono, F. E. Figueirido, and R. M. Levy, Chem. Phys. Lett. 263,
521 ͑1996͒.
44
X. Rozanska and C. Chipot, J. Chem. Phys. 112, 9691 ͑2000͒.
45
H. Schreiber and O. Steinhauser, Biochemistry 31, 5856 ͑1992͒.
46
H. Schreiber and O. Steinhauser, Chem. Phys. 168, 75 ͑1992͒.
47
H. Schreiber and O. Steinhauser, J. Mol. Biol. 228, 909 ͑1992͒.
48
P. E. Smith and B. M. Pettitt, J. Chem. Phys. 95, 8430 ͑1991͒.
49
G. S. Del Buono, T. S. Cohen, and P. J. Rossky, J. Mol. Liq. 60, 221
͑1994͒.
50
L. Perera, U. Essmann, and M. L. Berkowitz, J. Chem. Phys. 102, 450
͑1994͒.
51
S. Boresch and O. Steinhauser, J. Chem. Phys. 115, 10780 ͑2001͒.
52
P. H. Hu¨nenberger, in Simulation and Theory of Electrostatic Interactions
in Solution: Computational Chemistry, Biophysics, and Aqueous Solution,
edited by G. Hummer and L. R. Pratt ͑American Institute of Physics, New
York, 1999͒, pp. 17–83.
53
M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 ͑1998͒.
54
M. Deserno and C. Holm, J. Chem. Phys. 109, 7694 ͑1998͒.
55
G. Hummer, L. R. Pratt, and A. E. Garcia, J. Phys. Chem. A 102, 7885
͑1998͒.
56
M. Kastenholz and P. H. Hu¨nenberger, J. Phys. Chem B ͑to be published͒.
57
P. H. Hu¨nenberger and J. A. McCammon, J. Chem. Phys. 110, 1856
͑1999͒.
58
P. H. Hu¨nenberger and J. A. McCammon, Biophys. Chem. 78, 69 ͑1999͒.
59
W. Weber, P. H. Hu¨nenberger, and J. A. McCammon, J. Phys. Chem. B
104, 3668 ͑2000͒.
60
C. Peter, W. F. van Gunsteren, and P. H. Hu¨nenberger, J. Chem. Phys. ͑to
be published͒.
61
S. Boresch and O. Steinhauser, J. Chem. Phys. 111, 8271 ͑1999͒.
62
S. Boresch and O. Steinhauser, J. Chem. Phys. 115, 10793 ͑2001͒.
63
M. Marchi, D. Borgis, N. Levy, and P. Ballone, J. Chem. Phys. 114, 4377
͑2001͒.
64
C. Peter, W. F. van Gunsteren, and P. H. Hu¨nenberger, J. Chem. Phys. 116,
7434 ͑2002͒.
65
J. D. Madura, M. E. Davis, M. K. Gilson, R. C. Wade, B. A. Luty, and J.
A. McCammon, in Reviews in Computational Chemistry, Vol. 4, edited by
K. B. Lipkowitz and D. B. Boyd ͑VCH, New York, 1994͒, pp. 229–267.
66
B. Honig, K. Sharp, and A.-S. Yang, J. Phys. Chem. 97, 1101 ͑1993͒.
67
M. K. Gilson, Curr. Opin. Struct. Biol. 5, 216 ͑1995͒.
68
B. Honig and A. Nicholls, Science 268, 1144 ͑1995͒.
69
F. Figueirido, G. S. Del Buono, and R. M. Levy, J. Phys. Chem. B 101,
5622 ͑1997͒.
70
G. Hummer, L. R. Pratt, and A. E. Garcia, J. Chem. Phys. 107, 9275
͑1997͒.
71
S. Sakane, H. S. Ashbaugh, and R. H. Wood, J. Phys. Chem. B 102, 5673
͑1998͒.
72
G. Hummer, L. R. Pratt, A. E. Garcia, B. J. Berne, and S. W. Rick, J. Phys.
Chem. B 101, 3017 ͑1997͒.
73
H. S. Ashbaugh and R. H. Wood, J. Chem. Phys. 106, 8135 ͑1997͒.
74
J. A˚ qvist and T. Hansson, J. Phys. Chem. B 102, 3837 ͑1998͒.
75
H. S. Ashbaugh and R. H. Wood, J. Phys. Chem. B 102, 3844 ͑1998͒.
76
G. Hummer, L. R. Pratt, A. E. Garcia, S. Garde, B. J. Berne, and S. W.
Rick, J. Phys. Chem. B 102, 3841 ͑1998͒.
77
Y. N. Vorobjev and J. Hermans, J. Phys. Chem. 103, 10234 ͑1999͒.
78
T. Darden, D. Pearlman, and L. G. Pedersen, J. Chem. Phys. 109, 10921
͑1998͒.
79
Y. Y. Sham and A. Warshel, J. Chem. Phys. 109, 7940 ͑1998͒.
80
H. Resat, M. Mezei, and J. A. McCammon, J. Phys. Chem. 100, 1426
͑1996͒.
81
F. Figueirido, G. S. Del Buono, and R. M. Levy, J. Chem. Phys. 103, 6133
͑1995͒.
82
L. R. Pratt, G. Hummer, and A. E. Garcia´, Biophys. Chem. 51, 147 ͑1994͒.
83
A. Kovalenko and F. Hirata, J. Chem. Phys. 112, 10403 ͑2000͒.
84
A. Masunov and T. Lazaridis, J. Am. Chem. Soc. 125, 1722 ͑2003͒.
85
M. Neumann and O. Steinhauser, Mol. Phys. 57, 97 ͑1986͒.
86
M. Born, Z. Phys. 1, 45 ͑1920͒.
87
J.-C. Soetens, C. Millot, C. Chipot, G. Jansen, J. G. A´ ngya´n, and B.
Maigret, J. Phys. Chem. B 101, 10910 ͑1997͒.
88
L. Degre`ve and F. L. B. da Silva, J. Chem. Phys. 110, 3070 ͑1999͒.
89
D. Sardella Vieira and L. Degre`ve, J. Mol. Struct. 580, 127 ͑2002͒.
90
G. Hummer, D. M. Soumpasis, and M. Neumann, Mol. Phys. 77, 769
͑1992͒.
91
C. Satheesan Babu, P.-K. Yang, and C. Lim, J. Biol. Phys. 28, 95 ͑2002͒.
92
W. F. van Gunsteren, S. R. Billeter, A. A. Eising, P. H. Hu¨nenberger, P.
Kru¨ger, A. E. Mark, W. R. P. Scott, and I. G. Tironi, Biomolecular Simulation:
The GROMOS96 Manual and User Guide ͑Verlag der Fachvereine,
Zu¨rich, 1996͒.
93
H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comput. Phys.
Commun. 91, 43 ͑1995͒.
94
L. Onsager, J. Am. Chem. Soc. 58, 1486 ͑1936͒.
95
J. Warwicker and H. C. Watson, J. Mol. Biol. 157, 671 ͑1982͒.
96
P. Aufﬁnger and D. L. Beveridge, Chem. Phys. Lett. 234, 413 ͑1995͒.
9144 J. Chem. Phys., Vol. 119, No. 17, 1 November 2003 Bergdorf, Peter, and Hu¨nenberger
Downloaded 13 Jun 2006 to 128.200.197.134. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp