Inferring properties of macroscopic solutions from molecular simulations is complicated by the limited size of systems that can be feasibly examined with a computer. When long-ranged electrostatic interactions are involved, the resulting finite size effects can be substantial and may attenuate very slowly with increasing system size, as shown by previous work on dilute ions in bulk aqueous solution. Here we examine corrections for such effects, with an emphasis on solvation near interfaces. Our central assumption follows the perspective of Hünenberger and McCammon [J. Chem. Phys. 110, 1856 (1999)]: Long-wavelength solvent response underlying finite size effects should be well described by reduced models like dielectric continuum theory, whose size dependence can be calculated straightforwardly. Applied to an ion in a periodic slab of liquid coexisting with vapor, this approach yields a finite size correction for solvation free energies that differs in important ways from results previously derived for bulk solution. For a model polar solvent, we show that this new correction quantitatively accounts for the variation of solvation free energy with volume and aspect ratio of the simulation cell. Correcting periodic slab results for an aqueous system requires an additional accounting for the solvent’s intrinsic charge asymmetry, which shifts electric potentials in a size-dependent manner. The accuracy of these finite size corrections establishes a simple method for a posteriori extrapolation to the thermodynamic limit and also underscores the realism of dielectric continuum theory down to the nanometer scale.
Molecular simulations of ions near the boundaries of liquid solutions1–31 have challenged and reshaped traditional microscopic perspectives on solvation thermodynamics.32 In the case of aqueous interfaces, they have inspired a new generation of experiments10,15,28,29,33–43 and theories4,8–10,19,24,30,44 to identify basic driving forces and their implications across a broad range of fields such as electrochemistry,45 aerosols,46,47 and protein surfaces.48,49 The great strength of such simulations—resolving atomistically detailed structure on ultrafast time scales—is counterbalanced, however, by the necessary use of imperfect potential energy models and by the necessary restriction to small system sizes. This paper addresses the latter problem, aiming to extract from microscopically finite simulations the behavior of macroscopically large solutions and their interfaces.
Efforts to correct for finite size effects in molecular simulation are nearly as old as the simulation methods themselves. Most notably, the use of periodic boundary conditions (PBC) mitigates the unrealistically high surface-to-volume ratios, and degrees of interfacial curvature of small systems. Even for systems that interact with relatively short-ranged interactions, however, careful treatment of interactions between distant particles is often required to obtain well-converged thermodynamic properties.50–64 In the case of ion solvation, the very slow decay of Coulomb interactions with distance demands an even more careful treatment. Indeed, a simulation cell bearing a net charge cannot be periodically replicated without incurring a super-extensive energetic cost. Whether implicitly or explicitly, neutralizing countercharge must be introduced to ensure a well-defined thermodynamic limit when an ion is subject to periodic boundary conditions. Here we will exclusively consider systems that are periodic in three dimensions, whose unit cells contain a single charged solute and a spatially uniform neutralizing background charge density, in addition to a collection of electroneutral solvent molecules. Due to the long range of electrostatic forces, periodic replicas of these charged constituents can generate substantial fields and potentials, which in turn induce solvent response that decays only gradually with increasing system size. These are the contributions we wish to assess and correct.
While far from macroscopic, computationally accessible systems are large enough to capture much of the complicated molecular-scale response that is missing from the classic linear response theory for solvent polarization, i.e., dielectric continuum theory (DCT). Our approach asserts that DCT accurately describes solvent response at all larger scales. As such, finite size effects in DCT should serve as useful estimates for those of more complicated molecular models. This perspective was employed by Hünenberger and McCammon65 to investigate the nature and magnitude of periodicity-induced artifacts for bulk solutions, building on previous studies by Hummer et al.66,67 The development described in Sec. III casts the resulting finite size corrections in a different form that facilitates analysis of interfacial systems. In Sec. III B, we confirm that our approach is equivalent to that of Ref. 65 for a charged solute in bulk solution. We also verify Hummer’s demonstration that simulations, performed using methods described in Sec. II, closely follow the predicted approach to the macroscopic limit.
Extended interfaces between coexisting phases can also be simulated using PBC. The widely used “slab geometry” includes regions of each phase (e.g., liquid and vapor) that span the periodic simulation cell in two Cartesian directions (x and y), separated by a planar interface perpendicular to the third direction (z). We show in Sec. II that solvation free energies computed from simulations with this geometry exhibit a strong size dependence that differs significantly from the bulk case. Results are presented both for the SPC/E model of water and also for a simple polar solvent comprising diatomic molecules with a pair of opposite charges. In Sec. III C, we derive a thermodynamic finite size correction for ions solvated by an ideal dielectric solvent in the slab geometry. This correction quantitatively reconciles the size dependence observed for the simple polar solvent.
For the case of liquid water coexisting with vapor in the slab geometry, the finite size correction obtained from DCT accounts only partly for the size dependence observed in simulations. We attribute the remaining discrepancy to the intrinsic asymmetry of water’s molecular charge distribution. In the SPC/E model, positive and negative charges are not equivalently distributed within each water molecule, in contrast to the diatomic molecules in our simple polar solvent. This asymmetry creates thermodynamic biases in the solvation of cations versus anions and also underlies the net orientation of molecules at the air-water interface. In Sec. III D, we show that charge asymmetry in the periodic slab geometry generates a size-dependent electric potential that can be easily extrapolated to the macroscopic limit. The corresponding finite size correction, added to the dielectric response correction, yields ion solvation free energies from aqueous slab simulations that are very nearly independent of unit cell size and aspect ratio. A similar charge asymmetry correction is applied to results of bulk periodic simulations in Sec. III E, yielding good agreement with the limiting behavior of slab simulations.
The size dependence of ion solvation free energies near simulated air-water interfaces thus closely follows predictions of an idealized continuum model, once effects of charge asymmetry under PBC have been reckoned. This success of DCT extends to simulation unit cells with dimensions as small as 1 nm, indicating a surprising robustness of linear response theory at nearly molecular length scales. In Sec. IV, we discuss implications of this result for microscopic mechanisms and theories of ion solvation near interfaces. In Sec. V, we conclude.
II. SIMULATION METHODS
The general system of interest, shown schematically in Fig. 1, comprises a set of periodically replicated liquid slabs with width and dielectric constant ϵ, separated by regions of vacuum in the direction normal to the interface, which we take to be along the z-direction. The length of the primary unit cell in the z-direction is Lz. The liquid/vapor interface lies in the xy-plane and has a surface area A = LxLy, where Lx and Ly are the lengths of the simulation cell in x and y, respectively. The total volume of the simulation cell is = LxLyLz. A periodically replicated solute is situated in the liquid slab at a distance d from one of the liquid/vapor interfaces. A bulk solution lacking liquid/vapor interfaces is obtained by setting = Lz.
We have numerically examined finite size effects on ion solvation through molecular simulations of two model solvents. The simpler of these solvents, diatomic molecules comprising oppositely charged particles separated by a short distance, closely resembles the Stockmayer fluid. The other, a simple point charge model of water, adds the complication of asymmetry under charge inversion, whose implications for finite size effects are analyzed in Sec. III D.
Our model solute is the same in both sets of simulations, namely, a Lennard-Jones particle with a point charge at its center. Parameters for the Lennard-Jones interaction, u(r) = 4ε[(σ/r)12 − (σ/r)6] as a function of distance r, are identical to those of the SPC/E model for water, i.e., ε = 0.1553 kcal/mol and σ = 3.166 Å. This choice of σ corresponds roughly to the size of a Ca2+ ion.66 (The leading order correction considered in this work is, however, independent of solute size.) For consistency with standard notation, we use the symbol σ here to denote the Lennard-Jones diameter. The same symbol will later denote a different quantity, namely, the surface charge that results from dielectric response. The distinction should be clear from context.
In water studies, we employed the SPC/E model for interactions involving solvent molecules.68 This choice, made for simplicity, affects the exact values of computed solvation free energies, but our conclusions about finite size effects should apply equally well to more elaborate molecular models. Bulk liquid water simulations included a single solute and a number N of water molecules ranging from 64 to 756 such that the total number density was 0.033 33 Å−3. Simulations of air-water coexistence included one solute and N = 128 water molecules in a simulation cell with lateral dimensions Lx = Ly = 12.429 Å and varying separation between periodic liquid slabs. (As the sole exception, in Fig. 6, we show results for Lx = Ly = 31.1 Å as well.) The resulting aspect ratios Lz/Lx ranged from 4 to 30. To suppress interfacial instabilities at very high aspect ratios, we introduced an external confining potential, specifically walls that effectively constrain molecules to a region slightly larger than the neat liquid slab occupies. These walls interact with the oxygen atom of each water molecule through a potential usolv-wall(z), which inside the walls (|z| < |zwall|) has a Weeks-Chandler-Andersen form:69 for |δz| < 21/6σsolv-wall, where δz = z − zwall, and zwall = ±15 Å. Outside the walls (|z| > |zwall|), usolv-wall(z) is infinite. Otherwise (|δz| > 21/6σsolv-wall), usolv-wall(z) vanishes. The coordinate −Lz/2 < z < Lz/2 is the molecule’s position within the primary simulation cell. Away from the interfaces, the density of the liquid varied approximately between 0.032 00 Å−3 and 0.033 33 Å−3, depending on the aspect ratio and the charge of the solute.
The simpler polar fluid we consider consists of “dumbbell” molecules, each a pair of Lennard-Jones spheres (ε = 0.1553 kcal/mol and σ = 3.166 Å) separated by a rigid bond of length lD = 0.25 Å. Opposite charges of magnitude qD = 0.9e reside at the centers of these spheres. We consider a liquid state density of 0.025 Å−3. Bulk simulations included a number N of dumbbells ranging from 48 to 384. Slab simulations included N = 96 molecules. As in aqueous simulations, lateral dimensions were fixed at Lx = Ly = 12.429 Å, while the aspect ratio Lz/Lx varied from 4 to 30. Walls were again used to confine solvent molecules to the range −15.0 Å < z < 15.0 Å. However, in this case, the confining potential interacts with both Lennard-Jones spheres of each dumbbell molecule.
Some of the finite size corrections detailed below involve the static dielectric constant ϵ. In the case of SPC/E water, the value of ϵ has been established by previous work and is sufficiently large that required factors of (ϵ − 1)/ϵ are nearly unity. The dumbbell fluid is considerably less polarizable so that a more precise estimate of ϵ is needed. From a 5 ns bulk simulation of 512 dumbbell molecules (and no solute), we estimated the dielectric constant to be ϵ ≈ 7.1, obtained from the fluctuation-dissipation theorem,70
where M is the total dipole moment of the simulation cell and (kB is Boltzmann’s constant). We append the subscript “E = 0” to emphasize that this formula holds with the use of tin foil boundary conditions (see, e.g., Refs. 70–72). Similar dumbbell models have been examined in previous studies,73 which reported similar values of ϵ for comparable molecular dipole moments.
All simulations were performed with the LAMMPS simulation package.74 Dynamics were propagated at constant volume V and temperature T = 298 K using Langevin dynamics as implemented in LAMMPS,75,76 with a time step of 1 fs. The SHAKE algorithm was used to constrain all bonds and angles.77,78 Long range electrostatic interactions were evaluated using Ewald summation with tin foil boundary conditions. (The solvent’s net dipole would experience different forces under a different choice of boundary conditions.70,79–84 The solvent response of interest, however, does not directly involve this dipole, so we expect the choice of boundary conditions to have minor significance.) Ewald sums were calculated using the particle-particle particle-mesh solver,85 with parameters chosen such that the RMS error in the forces were a factor 105 smaller than the force between two unit charges separated by a distance of 1.0 Å.86
The solvation free energies we consider are precisely defined by
where ϕsolv is the electric potential generated by the solvent, evaluated at the center of a solute; P(L)(ϕsolv; q) is its probability distribution in the presence of a solute charge q; and denotes an average over P(L)(ϕsolv; q). The superscript (L) specifies the dimensions L = (Lx, Ly, Lz) of the simulated unit cell. In the limit L → ∞ of large simulation cell size, approaches the reversible work of charging a solute at infinite dilution. Evaluating the integral in Eq. (2), which is a statement of the potential distribution theorem,87 requires knowledge of P(L)(ϕsolv; 0) in its extreme wings, which we obtained from umbrella sampling. Specifically, we introduced a series of solute charges q/e = −1.0, −0.9, …, 0, …, +0.9, +1.0, which bias typical solvent potential fluctuations to a corresponding series of ranges centered on , where
Samples from the probability distributions were then reweighted and combined according to the MBAR algorithm.88 The resulting statistics of extreme solvent fluctuations required to compute could as well be obtained from alternative methods such as thermodynamic integration, reversible work calculation, or irreversible transformations via Jarzynski’s identity.66,89,90 The biasing method we used is especially convenient for molecular dynamics simulation, as it can be achieved simply by modulating the charge of a single particle.
The value of ϕsolv(RN) for a given configuration RN of solvent molecules was evaluated by inserting a solute test charge, and then subtracting contributions due to the test charge itself. In detail, we compute the total potential energy u(RN; qtest) of a system including solvent, a periodic collection of test charges qtest, and a charge neutralizing background. From this, we first subtract the energy u(RN; 0) of interactions among solvent molecules and of Lennard-Jones interactions with the solute. We then subtract the energy of interactions among periodic test charges and the neutralizing background, qtestϕwig/2, where ϕwig is the canonical Wigner potential discussed in Sec. III B and Appendix B. Dividing the remainder by qtest yields the solvent potential, .
Charging free energies computed with these methods are presented in Figs. 2 and 3, for bulk and interfacial systems, respectively. A strong dependence on periodicity L is evident for both model solvents considered. For cubic simulation cells, bulk free energies vary by 10s of kBT as cell dimensions grow from roughly 1 to 2 nm. Bulk results for anisotropic cells differ even more dramatically, by 100s of kBT as the aspect ratio Lz/Lx grows from 2 to 12, and appear not to converge as Lz → ∞ with fixed Lx. Similarly large finite size effects are exhibited by interfacial simulations. In this case, however, charging free energies do appear to converge in the limit of infinite aspect ratio, pointing to qualitative differences in the system size dependence of periodic bulk and slab geometries. Such nontrivial convergence behavior also warns that extrapolating results naively from ever-larger systems may not achieve physically meaningfully asymptotic properties.
III. CORRECTING FOR FINITE SIZE EFFECTS ON IONIC SOLVATION FREE ENERGIES
This work is by no means the first to address the role of finite size effects in ionic solvation free energies obtained from molecular simulation,65–67,91,92 and indeed, the central assumption—that DCT adequately describes long-wavelength dielectric response—follows the perspective of Ref. 65. The majority of previous work that addresses finite size effects has, however, focused on a scenario that lacks macroscopic dielectric boundaries, such as the liquid/vapor interface. Here we present a general formalism that is also applicable in such cases.
The quantity that we seek to calculate is the free energy to introduce a charge q to the center of a neutral cavity embedded in a solvent with dielectric constant ϵ, at infinite dilution,
F(∞)(q) is the free energy of a macroscopic system in the presence of a single solute with charge q. Molecular simulations cannot directly access , but can instead provide corresponding free energies under PBC,
As in Sec. II, the superscript (L) indicates spatial periodicity L = (Lx, Ly, Lz) in three Cartesian directions. Our goal is to connect the charging free energies for periodic and macroscopic systems, i.e., to calculate a correction
that can be applied to simulation results as an extrapolation to infinite dilution.
We consider two distinct contributions to ΔF(L; q). The first accounts for size dependence of a solvent’s structural response to solute charging, which we estimate using dielectric continuum theory. The second considers biases on ion solvation that exist even in the unperturbed solvent, specifically electric potentials experienced by a neutral solute, which are also size-dependent.
A. Size-dependent dielectric response
DCT is defined by a simple linear response relationship between the electric field at a point r and the average induced polarization density m(r) at that location. Together with Poisson’s equation, this relation connects spatial variations in m(r) to the density ρ(r) of free charge,
(Throughout this formulation, we work in a unit system in which 4πϵ0 = 1, where ϵ0 is the permittivity of free space.) In our case, the free charge represents periodically arranged solutes and their compensating uniform background of charge density ,
Here, ρ0(r) is the charge distribution of a single solute, localized at position r0 inside the dielectric and with net charge q = ∫drρ0(r). The solute’s periodic replicas are separated by lattice vectors b = (nxLx, nyLy, nzLz), with nx, ny, and nz taking on all integer values. If dielectric boundaries are present, we take them to have the same spatial periodicity as the solutes, as well as a clear dependence on L. For the periodic slab geometry, these boundaries are a pair of infinite parallel planes with fixed separation w, repeated in the z direction with periodicity Lz (see Fig. 1).
To calculate the charging free energy for this theory, we exploit a generic property of linear response, namely, that the change in free energy due to an external field is half the average interaction energy between the system and field. For a dielectric material interacting with periodic solutes and their neutralizing background,
where ϕsolv(r) is the electric potential generated by the polarization field m(r). Because the solvent is electrically neutral, its interaction with the uniform background charge does not contribute to Eq. (9). Determining a finite size correction ΔFDCT(L; q) from dielectric continuum theory thus amounts to calculating the average solvent potential at the center of a solute.
We obtain the solvent potential by integrating over the volume Ω occupied by the solvent,
This volume is periodic by construction but it need not be connected. If Ω comprises disconnected regions, then a periodic series of surfaces ∂Ω bounds these regions. The unit normal vector on these boundaries we denote by . Integrating by parts and using the divergence theorem gives
where R is a point on the set of boundaries ∂Ω. The normal component of solvent polarization at a boundary, , serves as an effective surface charge density due to polarization. Using the basic linear response relation in Eq. (7), we finally obtain
The correction ΔFDCT(L; q) compares response at finite L with the macroscopic result. In the limit L → ∞, all solute replicas described by ρ(r) are irrelevant except for the central charge at r0, and the background charge density vanishes. As a result,
where Δσ(R) = σ(R; L) − σ(R; ∞) is a difference between periodic and macroscopic systems. Motivated by the slab geometry of interest, we have assumed that there is a clear correspondence between boundaries for finite and infinite L. The L → ∞ slab possesses just two planar surfaces (neglecting the solute’s excluded volume93), bounding the solvent from above and below; these two surfaces are also present at finite L. On all other boundaries included in ∂Ω, there is no macroscopic contribution, i.e., σ(R; ∞) = 0. A more complicated set of periodic boundaries might not permit this simplification, requiring instead an explicit subtraction of surface integrals evaluated on different sets of boundaries.
For any choice of boundaries, evaluating Eq. (13) requires knowledge of the function σ(R; L). We have not solved a general dielectric boundary value problem, but have instead rephrased it in a convenient way. For the slab geometry, σ(R; L) can be worked out exactly in terms of a series of image charges, as detailed in Appendix A. More importantly, Δσ(R) is amenable to a greatly simplifying approximation, which we will show to be very accurate in the case of periodic slabs. For a bulk system lacking boundaries, Eq. (13) does constitute a full solution. In Sec. III B, we show that this result is consistent with previous work on ion solvation in bulk solvents.
B. Application to ion solvation in the absence of interfaces
Understanding the behavior of ΔF(L; q) for “bulk” periodic systems—those lacking any macroscopic dielectric boundary—has been the focus of much research. In an early study on the subject, Hummer et al.66 argued that estimates of charging free energy should include interactions of the ion with its periodic replicas and background charge. For simulations of liquid water using Ewald summation and a cubic cell of side length L, they showed that incorporating these “self-interactions” can yield free energies that are essentially independent of system size. The finite size correction is thus well approximated in this case by
where ϕwig/q ≈ −2.837 297/L (see, e.g., Ref. 94, and note that atomic units are implied in the numerical value). The Wigner potential ϕwig/q is defined as the electric potential at the site of a unit point charge due to all of its periodic replicas and a homogeneous background charge that acts to neutralize the primitive cell. (ϕwig/q is commonly referred to as ξEW in the literature.) Based on a dimensional analysis, Figueirido et al.91 proposed the following correction for finite dielectric constant ϵ:
which has the correct limiting behavior for ϵ → 1. Later studies extended these corrections to also account for the finite size of the ion.65,67 We do not consider such higher order corrections in this article.
In the absence of interfaces, our DCT expression in Eq. (13) simplifies to
For a cubic primitive cell (Lx = Ly = Lz = L), this result is equivalent to Figueirido’s, as expected from the DCT analysis of Ref. 65. Figure 2 exemplifies the remarkable effectiveness of this correction, echoing the conclusions of Hummer et al. Here we have added ΔFDCT(L; q) to the values of obtained from bulk molecular simulations of water (d) and the simple polar fluid (b). The DCT correction reconciles solvation free energies from simulations with significantly different periodicity, reducing discrepancies of ∼25kBT down to ∼1kBT. Residual finite size effects beyond ΔFDCT(L; q) are noticeable only for the smallest simulations of the simple polar liquid. We have confirmed numerically that these differences can be largely removed with higher-order DCT corrections that account for the solute’s excluded volume.65,67
For an anisotropic cell, Eq. (16) is a very straightforward generalization of the standard correction. To evaluate it numerically, we have used Ewald summation as outlined in Appendix B. In Fig. 2, we include simulation results for cuboidal simulation cells with Lz/Lx = 2 and 12, demonstrating that the DCT correction remains accurate even as the cell’s aspect ratio becomes large. The interesting aspect of this comparison is that ΔFDCT(L; q) diverges as Lz/Lx → ∞ with LxLy fixed. The accuracy of our correction then suggests that diverges in the opposite sense with growing aspect ratio. Indeed, simulation results with 756 molecules and Lz/Lx = 12 extend far beyond the plotted energy scale , hinting at such a divergence.
C. Application to periodic slabs
Applying the general result in Eq. (13) to systems with boundaries is considerably more difficult, as it requires calculation of the polarization surface charge σ(R; L) [as well as integration of ρ(r)/|r − r0| over disconnected volumes]. The spatial variation of σ(R; L) encodes important driving forces of solvation. A point ion near a dielectric boundary experiences an electric field that grows sharply in magnitude as it approaches the surface. For a semi-infinite dielectric with a planar boundary, this field is equivalent in DCT to the force generated by an image charge opposite to the boundary, accompanied by strong spatial variation in the polarization surface charge. A poor description of this position dependence would yield a poor approximation to . This feature of σ(R; L), however, is predominantly local in space and therefore only weakly sensitive to PBC. Moreover, such local aspects of solvation are likely not well described by DCT, which lacks microscopic complexities that can dominate on small scales. The difference quantity Δσ(R; L), we argue, should have much weaker spatial dependence and should be well described by DCT.
For the specific case of periodic slabs, a solution for σ(R; L) can be obtained by summing contributions from an infinite series of effective image charges, as described in Appendix A. This solution exhibits the rapid spatial variation described above, particularly as the solute ion nears the dielectric interface. But as anticipated, the image charges that contribute most strongly to this position dependence are identical in value and placement when L → ∞. The finite size effects of interest to our work originate instead in the long-wavelength features of polarization surface charge. We thus propose a greatly simplifying approximation, namely, that is constant along the planar boundaries of each dielectric slab. (This argument relies on the ion residing within the solvent’s liquid phase. As discussed in Appendix A, the approximation is inappropriate when r0 lies in the vapor phase instead.)
An appropriate value for can be determined by integrating over a single periodic replica ∂ω of the slab’s boundaries,95
where ω denotes a single periodic replica of the slab (see Fig. 1). In using the divergence theorem, we have exploited the fact that symmetry and continuity of m(r) within the slab requires that at the lateral interfaces between periodic replicas. The use of Eqs. (7) and (8) then yields
where is the slab’s thickness and A = LxLy is the area of a single periodic replica of the upper (or lower) dielectric boundary. This value of is precisely what is required for ΔFDCT(L; q) to be finite.
The finite size correction that results from this uniform surface charge approximation, together with Eq. (18),
has a simple physical interpretation. The quantity in square brackets in Eq. (19) corresponds to an electric potential generated by (i) periodically arranged unit point charges (excepting the central replica); (ii) a partially neutralizing charge density that is uniform within the dielectric slabs and vanishes outside; and (iii) uniformly charged plates at the dielectric boundaries, which achieve overall neutralization. It can be regarded as a generalization of the Wigner potential, in which the uniform compensating charge outside the dielectric has been moved to the dielectric boundaries. Interestingly, in the limit Lz → ∞ with A and fixed, all of the compensating charge moves to these boundaries. This scenario—a single dielectric slab, with ions periodically replicated in the x and y directions, bounded above and below by uniformly charged neutralizing plates—is markedly similar to the neutralization scheme proposed by Ballenegger et al.96 for molecular dynamics simulations of charged slab systems that are aperiodic in the z direction.
The modification of the Wigner potential described above can be calculated exactly in closed form, as shown in Appendix B. The DCT finite size correction can then be written simply as
In the limit that Lx, Ly, and Lz all become infinite, with still fixed, , an important property of the full solution for ΔFDCT that is preserved by the uniform surface charge approximation. Note that does not depend upon the solute’s position within the slab.
Applying this correction to simulation results for periodic slabs of the simple polar fluid, we can account almost completely for the substantial variation of charging free energy with system size. Figure 3(b) shows results for Lx = Ly ≈ 12.4 Å, slab width ≈ 27.5 Å, and various aspect ratios Lz/Lx, with the solute located equidistant from the liquid slab’s two planar boundaries. Corrected free energies for these systems vary by less than 1.05 kBT.
Aqueous simulations exhibit a more complicated size dependence, which is only partially captured by . Charging free energies corrected according to Eq. (20) are shown in Fig. 3(d) for Lx = Ly ≈ 12.4 Å, slab width ≈ 29.2 Å, and various aspect ratios Lz/Lx, with the solute placed equidistant from the two interfaces. The variation of with Lz is greatly reduced by , but a systematic size dependence remains, leaving differences as large as 13.4 kBT unexplained. A breakdown of the uniform surface charge approximation cannot alone explain this failure since the remnant size dependence is different for cations and anions—a charge asymmetry that cannot be described by DCT. The source of this charge-asymmetric finite size effect is considered in Sec. III D.
For the sake of simplicity, Fig. 3 presents results only for the case d = /2, i.e., with the ion situated midway between the liquid slab’s boundaries. In Sec. III D, we show that the corrections derived above are equally accurate for other values of d, provided the solute remains in the liquid slab.
D. Finite size effects due to charge asymmetry
The simple polar fluid we have considered is charge symmetric, in the sense that the positive and negative charge is equivalently distributed within each dumbbell molecule. By contrast, the SPC/E model of water is charge asymmetric, with positive charges situated farther from the center of volume exclusion. Such asymmetry can lead to a physically meaningful discrimination between cations and anions. In simulations with PBC, it can also generate unphysical differences. For example, a solvent comprising volume-excluding molecules whose internal charge distribution is spherically symmetric cannot respond to the charging of a solute, yet a bulk simulation of such particles under PBC can yield a nonzero charging free energy that is different for cations and anions.97,98 This artifact reflects a sensitivity of ϕsolv(r0) to the existence of interfaces even when they are arbitrarily far away.
Bulk simulations with PBC lack interfaces regardless of periodicity so that artifacts in persist even in the limit L → ∞. The slab geometry is an intermediate case, with explicit interfaces but also artifacts and finite size effects of charge asymmetry that vanish in the L → ∞ limit. To demonstrate this fact, we examine the average electric potential ϕneut experienced by a neutral solute in a simulation of periodic slabs. The dependence of this “neutral cavity potential” on system size provides a way to extrapolate computed charging free energies to the macroscopic limit.
The neutral cavity potential can be written exactly as
where denotes a single unit cell of the periodic slab simulation. is the solvent’s average charge density at position r, with a neutral solute at position r0. Due to volume exclusion, ⟨ρsolv(r)⟩ vanishes inside the solute (and its periodic replicas). If the solute resides within the slab, this constraint requires that ⟨ρsolv(r)⟩ depends on x and y as well as z, a significant feature for contributions to ϕneut from the primary simulation cell (b = 0). Our interest lies, however, in contributions from more distant unit cells, which dictate the system size dependence of ϕneut. These distant contributions are much less sensitive to the solvent’s electrostatic inhomogeneity in x and y. We thus expect that
is a good approximation, where const is independent of L, and the quantity has been spatially averaged over x and y. Figure 4(b) shows this laterally averaged charge density, obtained from simulations of SPC/E water with a neutral solute. Results for two different solute locations (one at the middle of the slab and the other at the liquid/vapor interface) are nearly indistinguishable. Provided the solute’s cross-sectional area πσ2/4 is not a substantial fraction of LxLy, we expect the laterally averaged charge density to be insensitive to the solute’s presence at all.
Removal of the dependence on the lateral coordinates allows for great simplification. As shown in Appendix B, exploiting the symmetries of yields
The difference quantity
can now be written compactly as
where we have eliminated terms that vanish due either to electroneutrality of the solvent or to the zero net dipole moment of a pure liquid slab. Numerical results for Δϕneut(L) are shown in Fig. 4(a), obtained directly from aqueous simulations and also from Eq. (25). They indicate that neglecting the lateral dependence of is indeed a very good approximation.
Equation (25) describes an intrinsic bias on ion solvation in a system of periodic slabs, one that vanishes in the limit Lz → ∞. We take this bias, in the absence of any solute charge, as a baseline for dielectric response, and therefore propose a full finite size correction for the slab geometry,
Equations (20) and (25) provide accurate, inexpensive, and easily implemented approximations for ΔFDCT(L; q) and Δϕneut(L), respectively. Charging free energies corrected according to Eq. (26) are shown in Fig. 5(a) for an ion at the center of a periodic slab of liquid water. The remnant size dependence of evident in Fig. 3(d) is captured almost exactly by the neutral cavity potential. Equation (26) accounts equally well for finite size effects observed for an ion placed at the air/water interface, as shown in Fig. 5(b).
Several features of the neutral cavity potential are noteworthy. First, within the approximation leading to Eq. (25), Δϕneut(L) is independent of the solute’s location, even for points well outside the slab. This invariance suggests that the bias we are correcting is dominated by very distant contributions, in effect a consequence of omitting each slab’s boundary at x, y = ∞ under PBC. Indeed, Eq. (25) resembles the conventional surface potential , where zliq and zvap denote locations on either side of a neat liquid/vapor interface (zvap < zliq). These potentials are simply related within a molecular multipole expansion up to quadrupole order, assuming that the solvent’s average polarization field is nonzero only at the interface and that its quadrupole field is uniform within the liquid phase. The neutral cavity potential then evaluates simply to Δϕneut(L) ≈ ϕsurf/Lz. From this relationship, our numerical results for Δϕneut suggest a surface potential of roughly −0.54 V for the SPC/E model, comparable to published values.12,98,100–102 Like the conventional surface potential, Δϕneut(L) may therefore have little to do with the molecular physics of solvation, as underscored by the possibility that Δϕneut(L) ≠ 0 even for spherically symmetric solvent molecules. In computing the solvent potential from a lattice sum method, we effectively adopt the so-called “P-summation scheme” (see, e.g., Refs. 73, 89, and 103) such that the surface potential contains unphysical contributions from the internal charge distribution of the solvent molecules.
For bulk liquid simulations, such artifacts of PBC are not removed by extrapolation to the macroscopic limit. For the slab geometry, Eq. (22) indicates that these artifacts are completely removed as Lz → ∞. Specifically, if the neutral solute resides in the vapor phase,104 r0 = rvap, then Eq. (22) is exact, with a vanishing constant term. Integrating Eq. (23) then gives . In other words, solvent potentials are properly referenced to the vapor phase when the spacing between slabs is infinite. Based on this fact, one convenient way to evaluate Δϕneut(L) in slab simulations is simply to compute the average electric potential at a point within the vapor phase (and far from the interface). Also shown in Fig. 5(a) are the results from a slab simulation employing the popular Yeh-Berkowitz method99 for simulating slab systems, in conjunction with the neutralization scheme outlined in Ref. 96. (In this case Δϕneut = 0, and the generalized Wigner potential has been evaluated numerically in the limit Lz → ∞.) Good agreement is observed between these results and those obtained with the a posteriori correction derived in this work.
E. Reconciling bulk and slab results
For the periodic slab geometry, with ions at the slabs’ centers, the limit L → ∞ we have considered corresponds to an individual ion in the middle of an individual slab with infinite lateral extent. This system of course still retains a finite dimension, namely, the slab’s thickness . Only as grows should we expect to approach the charging free energy for a macroscopic solution at infinite dilution. For the simple polar fluid, this remaining finite size effect is evident in Fig. 3(b), where deviations are visible between for periodic slab and bulk systems. A correction for finite can be estimated from DCT, as described in Appendix A, yielding an estimate of the macroscopic charging free energy away from any interface
where is the result of periodic slab simulations with the corresponding finite periodicity correction ΔFDCT + qΔϕneut applied. As shown in Fig. 6(a), this adjustment brings slab and bulk simulation results into good agreement for the simple polar solvent.
The same finite thickness correction should apply to the aqueous slab simulations as well. It is not sufficient, however, to reconcile differences between slab and bulk results in this case. Indeed, it is clear from Fig. 5(a) that these deviations are charge asymmetric, and thus cannot be accounted by DCT. We argue that the remaining discrepancy is due not to finite slab thickness but instead to artifacts of bulk periodic simulations. As discussed in Sec. III D, such bulk periodic systems lack interfaces even in the limit L → ∞, which can generate artifacts correctible only through interfacial considerations. (See Ref. 105 for a recent overview of this issue and of strategies to correct for it.) Specifically, a neutral solute at the center of a macroscopic spherical liquid droplet experiences an average electric potential,
where rliq denotes a point in the liquid far from both the solute and the liquid/vapor interface. The first term in Eq. (29) should be well described by a bulk periodic simulation as L → ∞, while the latter contribution is absent. Following the reasoning of Sec. III D, we account for the interfaces missing in a bulk periodic simulation by adding qϕsurf to the computed charging free energy. Our estimate of from bulk periodic simulations is thus
where is the result of periodic bulk simulations with the corresponding finite periodicity correction ΔFDCT applied. Numerically integrating for all of our neutral cavity simulations (using the “P-summation” convention as described above), we obtain values of ϕsurf ranging from −0.53 V to −0.58 V, comparable to previously reported surface potential values12,98,100–102 for point charge models, as well as our estimate based on Δϕneut. Figure 6(b) shows the macroscopic estimates for SPC/E water obtained from Eq. (30) with ϕsurf = −0.58 V, a value that yields particularly close alignment of bulk and slab results. This agreement suggests that the charging free energies in Fig. 6(b) are robust assessments of ion solvation in the physically realistic scenario of a macroscopic solution with boundaries. The full macroscopic solvation free energy would require adding the thermodynamic cost to introduce a neutral solute, which should not face significant finite size effects.
Many simulation studies have examined small ions near aqueous interfaces using a slab geometry. For system sizes typical of these studies, we have demonstrated finite size effects that can shift solvation free energies by 10s of kBT. While substantial in this sense, these effects are unlikely to change many qualitative conclusions from previous work. In particular, the finite size corrections of Eqs. (20) and (25) do not depend on the ion’s location within the liquid. They are therefore insensitive to a solute’s approach to the interface. Previously computed density and free energy profiles should therefore be unaffected for solute positions on the liquid side of the interface; in this case, our corrections simply amount to resetting the zero of energy. Our finite size corrections do not apply well when an ion resides outside the liquid, both because the uniform surface charge approximation of Eq. (19) fails and because non-dielectric response such as interface deformation is pronounced. Empirically, we find in this case that aligning simulation results from different system sizes requires a correction that does depend on the ion’s position relative to the interface. Capturing that dependence is a challenge that likely requires an approach more nuanced than DCT.
Our results for two model molecular liquids demonstrate that DCT can predict the system size dependence of interfacial ion solvation thermodynamics with striking accuracy. This success is remarkable in light of the length and energy scales involved and the complexities of interfacial response. The smallest systems we have investigated are periodic on a scale of just a few molecular diameters, well below the scale on which the liquid appears continuous. The scale of charging free energies in Fig. 6 underscores that the fluctuations governing solvation of a fully charged ion are extraordinarily rare (with relative probability <e−100) in its absence. In the case of water, the overall response to charging a solute is markedly nonlinear. Both in bulk water and in aqueous slabs, Fchg(q) deviates strongly from the parabolic form characteristic of linear response. [Most notably, the curvature d2Fchg/dq2 differs significantly at positive and negative q so that anions are far more favorably solvated than cations.66,106–109 A contribution linear in q, which favors cation solvation, is much weaker for appreciably charged solutes (|q| ≳ 0.5).73] The softness of liquid/vapor interfaces and the altered hydrogen bond network structure at the liquid’s boundary contribute additional nonlinearities that shape ions’ association with the surface.10,18,31,110,111 And yet a continuum linear response theory describes nm-scale thermodynamic finite size effects to an accuracy of ≲kBT. From the non-parabolic shape of Fchg(q), we know that this success of DCT cannot hold on all length scales. Our results indicate that charge-asymmetric nonlinear response in such water models is confined to the very near field of one or two solvation shells.
Our formulation of finite periodicity corrections for inhomogeneous geometries is motivated by a specific interest in how liquid water accommodates solutes near interfaces, and it builds directly on previous work from the molecular physics community. But similar finite size effects have been investigated in other areas, e.g., condensed matter physics (see, e.g., Refs. 112–116). The results and methods described here may thus be useful in a range of contexts involving polarization response, such as charged defect formation in semiconductors, and obtaining reliable reference potentials in computational electrochemistry.117,118
In this article, we have advanced a physical perspective and a mathematical framework for assessing the impact of spatial periodicity on ion solvation free energies computed from molecular simulation. Our general result Eq. (13) is a finite size correction from dielectric continuum theory, suitable for application to systems with boundaries where the local polarizability changes sharply. It requires as input long-wavelength features of induced surface polarization charge. For the case of a charged solute in a dielectric slab under periodic boundary conditions, taking this boundary charge to be uniform is an excellent approximation, which yields Eq. (20) as an easily computable correction.
For solvent molecules that are charge-asymmetric (such as water), an additional consideration outside the scope of DCT is needed to describe finite size effects in the slab geometry. Related to the surface potential, the resulting correction in Eq. (25) removes artifacts associated with the absence of distant interfaces under PBC. As Lz → ∞, this absence becomes inconsequential. By contrast, for periodic bulk simulations, the absence of interfaces is problematic even in the limit L → ∞. Using the surface potential to effectively restore those interfaces, we find that fully corrected bulk and slab results approach a consistent macroscopic limit.
We regard these successes as evidence that the electrostatic response of aqueous environments is well described by DCT down to length scales little larger than a molecular diameter. Theoretical work by others on the association of ions with the air/water interface has asserted that the realism of DCT extends even further in that context,6,8,9,19,27,119 in essence to arbitrarily small scales. The methods we have described facilitate a careful scrutiny of this hypothesis, which we pursue in ongoing work.
We would like to thank Dayton G. Thorpe and Layne B. Frechette for many fruitful discussions. The work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, through the Chemical Sciences Division (CSD) of Lawrence Berkeley National Laboratory (LBNL), under Contract No. DE-AC02-05CH11231.
APPENDIX A: METHOD OF IMAGES APPLIED TO SLAB SYSTEMS
Calculating the dielectric response to a point charge near a planar boundary between two semi-infinite polarizable media is a textbook exercise.120 Focusing on solvation near liquid-vapor interfaces, we take one of the regions to be vacuum (ϵ = 1) and the other region Ω to have a large dielectric constant ϵ. For a point charge q that is located in the high-dielectric medium at a point r0, lying a distance d from the interface, Poisson’s equation can be easily solved using the method of images. The average electric field within Ω is constructed with an effective external point charge, whose field is divergence-free inside Ω and enforces dielectric boundary conditions at the interface ∂Ω. The electric potential ϕ(r) is then specified by the requirement that ϕ vanishes at infinite distance from the real charge. For a point r within Ω,
is the potential generated by an image charge aq at a point that mirrors the location of the real charge, and a = (ϵ − 1)/(ϵ + 1). The unit vector is perpendicular to the interface, pointing outwards from the dielectric.
Solvation of a point charge in a dielectric slab, bounded by two planar interfaces with vapor, can be similarly solved. In this case, satisfying boundary conditions on both interfaces requires an infinite series of image charges, as depicted in Fig. 7. Here the real charge lies a distance d1 from its closest interface (the top interface in Fig. 7), and a distance d2 = − d1 from the other interface (the lower interface in Fig. 7), where is the slab thickness. Above the slab reside image charges: (i) aq, a3q, a5q, … at positions , and (ii) a2q, a4q, a6q, … at positions , where is an outward normal for the upper interface. Below the slab reside image charges: (i) aq, a3q, a5q, … at positions , and (ii) a2q, a4q, a6q, … at positions . Since the electric potential must again vanish as r → ∞, ϕ(r) is still given by Eq. (A1), now with
Finite size corrections in Sec. III C involve the potential ϕ(r) evaluated at the solute’s position. For this case of a single ion in a single slab, ϕ(r0) can be written exactly in terms of the Hurwitz-Lerch transcendent ,
The potential ϕself due to the point ion itself is singular, but it is independent of the boundaries’ geometry and therefore cancels in all quantities of interest in the main text. In the special case that the ion is equidistant from the interfaces (d1 = d2 = /2),
In Fig. 8, we plot the non-singular part of this potential, ϕimage(r0), for various slab widths.
The case of central interest in this paper—the periodic slab simulation geometry depicted in Fig. 1—is more complicated in two respects. First, image charges that determine the electric field do not follow a simple pattern. Second, ϕ(r) need not vanish as r → ∞ so that additional information is needed to fully specify a boundary value problem for the electric potential. As a result, the collection of image charges used to solve determine ϕ(r) only up to an additive constant. Our approach circumvents this issue [see Eq. (12)] by calculating the potential directly from the polarization surface charge density σ(R), which is in turn completely specified by the electric field,
A full solution in this scheme would first determine image charges that enforce boundary conditions for , calculate the resulting surface charge as a function of position R on each boundary, and then integrate as per Eq. (12) to obtain the electric potential. This procedure would be tedious but feasible. For the physical systems we are considering, which feature an ion inside a high-dielectric solvent that coexists with vapor, an accurate approximation can be computed much more easily, as explained below.
For a periodic collection of point charges in dielectric slabs, the method of images gives a series of effective charges exemplified by Fig. 9. For a given real charge, one or more of the nearest images are identical to the non-periodic case of Fig. 7 (assuming Lz ≥ 3/2). The charges and positions of more distant images are complicated by the system’s periodic structure, although simple recursion relations can be derived when Lz is an integer multiple of . The infinite collection of real charges also necessitates introducing a neutralizing background charge, which generates a series of slabs of image charge density (not shown). The background charge and its images contribute importantly to and thus also to σ(R). Those background contributions, however, are spatially uniform in x and y, i.e., they are independent of R on each boundary. We exploit this fact by splitting into a zero-wavevector component and a component δσ(R) that varies with R and integrates to zero. The average surface charge density can be calculated simply and without reference to image charges, as described in Sec. III C. The remainder δσ(R) receives no contribution from the background charge or its images.
Results presented in Sec. III C were obtained neglecting the lateral variation of polarization surface charge, in effect setting δσ to zero. The justifications for this uniform surface charge approximation are several-fold. First, the most significant contributions to Fchg through δσ come from the nearest image charges, which do not depend on Lz and thus cancel in the calculation of ΔFDCT. Contributions to σ(R) from more distant images are smaller in magnitude and, due to the greater distance, significantly dominated by their impact on the average charge density, for which already accounts. Finally, contributions to the surface polarization charge from nearby pairs of real and/or image charges that are equidistant from a boundary are typically of order 1/ϵ for large ϵ. For example, the physical charge q is largely offset by an image charge aq opposite to the boundary, with the residual (1 − a)q = 2q/(ϵ + 1) ∼ 1/ϵ. By contrast, the average charge density is of order unity.
The uniform surface charge approximation is less defensible when the solute ion resides instead in the low-dieletric vapor phase. Such a solute can experience significant solvation forces from multiple periodic replicas of the liquid slab; contributions from only one of these replicas will cancel in the calculation of ΔFDCT. Contributions from nearby charges also do not offset to order 1/ϵ in this case. For example, the image charge nearest the interface contributes to δσ with magnitude (ϵ − 1)/(ϵ + 1) ∼ 1 and does not have an offsetting counterpart opposite the boundary. Empirically, we find that finite size scaling predicted from this approximation is indeed not closely followed by simulations. The discrepancy likely has multiple sources. In addition to poorly approximating ΔFDCT, we expect that linear dielectric response is a much cruder model in this case. As small ions leave the liquid phase, they can deform the interface substantially,110,111,121 a response that lies distinctly outside the scope of DCT.
APPENDIX B: GENERALIZING THE WIGNER POTENTIAL
Consider a periodic charge distribution c(r) that is (i) electroneutral, drc(r) = 0, (ii) non-dipolar, drc(r) r = 0, and (iii) spatially uniform in two directions x and y, c(r) = c(z). The net potential ϕ(r0) generated by this distribution at observation point r0 is
Exploiting the lateral uniformity of c(r), this potential may be written as a Fourier series in just one dimension,
where , kz = 2πn/Lz, and the sum runs over all nonzero integers n. This sum can be rewritten as
in terms of a function J(z) that can be expressed in closed form,
The second equality in Eq. (B4) can be shown by twice differentiating J(z), noting that , and solving the ordinary differential equation that results. Doing so is aided by recognizing that d2|z|/dz2 = 2δ(z), that J(0) = πLz/3, and that J(z) is an even function of z. The assumed symmetries of c(z) allow further simplification,
Because the neutralizing charge distributions encountered in Sec. III C are uniform in x and y, their contributions to ΔFDCT can be simplified using Eq. (B5). In particular, the generalized Wigner potential that arises for periodic dielectric slabs differs from the conventional Wigner potential ϕwig only in the form of the neutralizing background. Defining
we can regard Δϕ* as the potential due to periodic, uniformly charged plates with charge density , together with a compensating background that is uniform outside the dielectric slabs and zero inside. This charge distribution has precisely the form assumed for c(r) above, with
where hout(z) = 1 outside the dielectric slabs (i.e., for −Lz/2 < z < − /2 and for /2 < z < Lz/2) and hout(z) = 0 inside (i.e., for −/2 < z < /2). For an observation point inside one of the dielectric slabs (−/2 < z0 < z/2), substituting Eq. (B7) into Eq. (B5) yields, after simple integration,
Note that this result is independent of the observation point z0, provided it lies within the dielectric. Note also that Δϕ* vanishes when Lz = so that the standard Wigner potential is recovered for a bulk periodic system.
The conventional Wigner potential, arising from periodic point charges and a homogeneous background charge ,
is most conveniently evaluated using Ewald summation,
For cubic simulation cells, the value of ϕwig/q ≈ −2.837 297/L is known to high precision.94 The finite size corrections described in this paper generally require ϕwig for anisotropic cells, which we computed from Eq. (B10) (with a large value of κLx that justifies neglecting the first sum entirely). For highly anisotropic cells, the magnitude of ϕwig becomes very large (see Fig. 10), appearing to diverge linearly with growing aspect ratio, . This divergence precisely matches that of Δϕ* in Eq. (B8) so that the generalized Wigner potential remains finite even for very large aspect ratios.
Acknowledging that the dielectric does not penetrate the solute is essential for local polarization response, which would otherwise be singular. This local response, however, should be consistent across different system sizes; a comparison among them eliminates the singularity and justifies treating the solute as a point charge within the dielectric domain Ω.
We take the net surface charge on the slab’s two boundaries to be equal. For a single aperiodic slab, this symmetry holds regardless of the ion’s position in the slab, as can be shown by summing image charges above and below each interface. For the periodic case, we assume that symmetrically placed replicas do not upset this balance (Figs. 7 and 9).
We have in mind that the solute is sufficiently far from the interface that capillary fluctuations are not disrupted.