Interfacial ion solvation: Obtaining the thermodynamic limit from molecular simulations

Interfacial ion solvation: Obtaining the thermodynamic limit from molecular simulations

Stephen J. Cox Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, United States.    Phillip L. Geissler Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, United States. Department of Chemistry, University of California, Berkeley, CA 94720, United States. geissler@berkeley.edu
July 4, 2019
Abstract

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.

I Introduction

Molecular simulations of ions near the boundaries of liquid solutions Jungwirth and Tobias (2001, 2002, 2006); Ben-Amotz (2016); Baer and Mundy (2011); Tobias et al. (2013); Netz and Horinek (2012); Levin (2009); Levin, dos Santos, and Diehl (2009); Otten et al. (2012); Archontis and Leontidis (2006); Arslanargin and Beck (2012); Baer et al. (2014); Beck (2011); McCaffrey et al. (2017); Caleman et al. (2011); Herce et al. (2005); Kumar, Knight, and Voth (2013); Levin and dos Santos (2014); Ou et al. (2013); Ou and Patel (2013); Stern et al. (2013); Tse et al. (2015); Vaikuntanathan, Shaffer, and Geissler (2013); Whitmer et al. (2014); Yagasaki, Saito, and Ohmine (2010); Baer et al. (2012); Mucha et al. (2005); Petersen et al. (2005); Archontis, Leontidis, and Andreou (2005); Noah-Vanhoucke and Geissler (2009) have challenged and reshaped traditional microscopic perspectives on solvation thermodynamics Onsager and Samaras (1934). In the case of aqueous interfaces they have inspired a new generation of experiments Otten et al. (2012); Piatkowski et al. (2014); Petersen and Saykally (2004, 2006); Verreault, Hua, and Allen (2012); Mishra et al. (2012); Liu et al. (2004); Mucha et al. (2005); Raymond and Richmond (2004); Tarbuck, Ota, and Richmond (2006); Viswanath and Motschmann (2007); Shultz et al. (2000); Petersen et al. (2005); McCaffrey et al. (2017); Weber et al. (2004) and theories Archontis, Leontidis, and Andreou (2005); Vrbka et al. (2004); Levin (2009); Levin, dos Santos, and Diehl (2009); Levin and dos Santos (2014); Otten et al. (2012); Vaikuntanathan, Shaffer, and Geissler (2013); Ben-Amotz (2016) to identify basic driving forces and their implications across a broad range of fields, such as electrochemistry Schmickler and Santos (2010), aerosols Finlayson-Pitts (2003); Knipping et al. (2000) and protein surfaces Baldwin (1996); Jungwirth and Winter (2008). 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 curvatures of small systems. For particles that interact solely through short-ranged interactions, using PBC and calculating interactions between particles’ nearest periodic replicas can yield well converged thermodynamic properties with only tens of particles Frenkel and Smit (2002); Allen and Tildesley (1987). In the case of ion solvation, the very slow decay of Coulomb interactions with distance demands 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 McCammon Hünenberger and McCammon (1999) to investigate the nature and magnitude of periodicity-induced artifacts for bulk solutions, building on previous studies by Hummer et al.Hummer, Pratt, and García (1996, 1997) 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.2 we confirm that our approach is equivalent to that of Ref. Hünenberger and McCammon, 1999 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 ( and ), separated by a planar interface perpendicular to the third direction (). 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.3 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.4 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.5, 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

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

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, as a function of distance , are identical to those of the SPC/E model for water, i.e.,  kcal/mol and  Å. 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 moleculesBerendsen, Grigera, and Straatsma (1987). 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 of water molecules ranging from 64 to 756, such that the total number density was  Å. Simulations of air-water coexistence included one solute and water molecules in a simulation cell with lateral dimensions  Å and varying separation between periodic liquid slabs. (As the sole exception, in Fig. 6 we show results for  Å as well.) The resulting aspect ratios 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 , which inside the walls () has a WCA form: for , where , and  Å. Outside the walls (), is infinite. Otherwise (), vanishes. The coordinate is the molecule’s position within the primary simulation cell. Away from the interfaces, the density of the liquid varied approximately between 0.03200 Å and 0.03333 Å, 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 ( kcal/mol and  Å) separated by a rigid bond of length  Å. Opposite charges of magnitude reside at the centers of these spheres. We consider a liquid state density of  Å. Bulk simulations included a number of dumbbells ranging from 48 to 512. Slab simulations included molecules. As in aqueous simulations, lateral dimensions were fixed at  Å, while the aspect ratio varied from 4 to 30. Walls were again used to confine solvent molecules to the range . 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 . For SPC/E water the value of has been established by previous work, and is sufficiently large that required factors of 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 , obtained from the fluctuation-dissipation theorem: Neumann (1983)

(1)

where is the total dipole moment of the simulation cell, and ( is Boltzmann’s constant). We append the subscript ‘’ to emphasize that this formula holds with the use of tin foil boundary conditions (see e.g. Refs. Zhang, Hutter, and Sprik, 2016; De Leeuw, Perram, and Smith, 1986; Neumann, 1983).

All simulations were performed with the LAMMPS simulation package Plimpton (1995). Dynamics were propagated at constant volume and temperature  K using Langevin dynamics as implemented in LAMMPS Schneider and Stoll (1978); Dünweg and Paul (1991), with a time step of 1 fs. The SHAKE algorithm was used to constrain all bonds and angles Ryckaert, Ciccotti, and Berendsen (1977); Andersen (1983). 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 conditionsde Leeuw, Perram, and Smith (1980a, b); Smith (1981); Redlack and Grindlay (1975); Neumann (1983); Kantorovich and Tupitsyn (1999); Ballenegger (2014). 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 Hockney and Eastwood (1988), with parameters chosen such that the RMS error in the forces were a factor smaller than the force between two unit charges separated by a distance of 1.0 ÅKolafa and Perram (1992).

The solvation free energies we consider are precisely defined by

(2)

where is the electric potential generated by the solvent, evaluated at the center of a solute; is its probability distribution in the presence of a solute charge ; and denotes an average over . The superscript specifies the dimensions of the simulated unit cell. In the limit of large simulation cell size, approaches the reversible work of charging a solute at infinite dilution. Evaluating the integral in Eq. 2 requires knowledge of in its extreme wings, which we obtained from umbrella sampling. Specifically, we introduced a series of solute charges , biasing typical solvent potential fluctuations to a corresponding series of ranges centered on , where

(3)

Samples from the probability distributions were then reweighted and combined according to the MBAR algorithm Shirts and Chodera (2008).

The value of for a given configuration 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 of a system including solvent, a periodic collection of test charges , and a charge neutralizing background. From this we first subtract the energy 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, , where is the canonical Wigner potential discussed in Sec. III.2 and Appendix B. Dividing the remainder by yields the solvent potential, .

Charging free energies computed with these methods are presented in Figs. 1 and 3, for bulk and interfacial systems, respectively. A strong dependence on periodicity is evident for both model solvents considered. For cubic simulation cells bulk free energies vary by 10s of as cell dimensions grow from roughly 1 to 2 nm. Bulk results for anisotropic cells differ even more dramatically, by 100s of as the aspect ratio grows from 2 to 12, and appear not to converge as with fixed . 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 Hummer, Pratt, and García (1996); Figueirido, Del Buono, and Levy (1997); Hummer, Pratt, and García (1997); Hünenberger and McCammon (1999); Hummer, Pratt, and García (1998), and indeed, the central assumption—that DCT adequately describes long-wavelength dielectric response—follows the perspective of Ref. Hünenberger and McCammon, 1999. The majority of previous work 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 to the center of a neutral cavity embedded in a solvent with dielectric constant , at infinite dilution:

(4)

is the free energy of macroscopic system in the presence of a single solute with charge . Molecular simulations cannot directly access , but can instead provide corresponding free energies under PBC:

(5)

As in Sec. II the superscript indicates spatial periodicity in three Cartesian directions. Our goal is to connect the charging free energies for periodic and macroscopic systems, i.e., to calculate a correction

(6)

that can be applied to simulation results as an extrapolation to infinite dilution.

We consider two distinct contributions to . 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.

iii.1 Size-dependent dielectric response

DCT is defined by a simple linear response relationship between the electric field at a point and the average induced polarization density at that location. Together with Poisson’s equation, this relation connects spatial variations in to the density of free charge,

(7)

(Unless otherwise specified, atomic units are implied throughout this development.) In our case the free charge represents periodically arranged solutes and their compensating uniform background of charge density ,

(8)

Here, is the charge distribution of a single solute, localized at position inside the dielectric and with net charge . The solute’s periodic replicas are separated by lattice vectors , with , and 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 . For the periodic slab geometry these boundaries are a pair of infinite parallel planes with fixed separation , repeated in the direction with periodicity (see Fig. 2).

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 system and field. For a dielectric material interacting with periodic solutes and their neutralizing background,

(9)

where is the electric potential generated by the polarization field . Because the solvent is electricially neutral, its interaction with the uniform background charge does not contribute to Eq. 9. Determining a finite size correction 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:

(10)

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:

(11)

where 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:

(12)

The correction compares response at finite with the macroscopic result. In the limit , all solute replicas described by are irrelevant except for the central charge at , and the background charge density vanishes. As a result,

(13)

where 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 . The slab possesses just two planar surfaces (neglecting the solute’s excluded volume111Acknowledging 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 .), bounding the solvent from above and below; these two surfaces are also present at finite . On all other boundaries included in , there is no macroscopic contribution, i.e., . 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 . We have not solved a general dielectric boundary value problem, but have instead rephrased it in a convenient way. For the slab geometry, can be worked out exactly in terms of a series of image charges, as detailed in Appendix A. More importantly, 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 the next section we show that this result is consistent with previous work on ion solvation in bulk solvents.

iii.2 Application to ion solvation in the absence of interfaces

Understanding the behavior of 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. Hummer, Pratt, and García (1996) 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 , 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

(14)

where (see e.g. Ref. Nijboer and Ruijgrok, 1988, and note that atomic units are implied in the numerical value). The Wigner potential 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. ( is commonly referred to as in the literature.) Based on a dimensional analysis, Figueirido et al. Figueirido, Del Buono, and Levy (1997) proposed the following correction for finite dielectric constant :

(15)

which has the correct limiting behavior for . Later works extended these corrections to also account for the finite size of the ion Hummer, Pratt, and García (1997); Hünenberger and McCammon (1999). We do not consider such higher order corrections in this article.

In the absence of interfaces, our DCT expression in Eq. 13 simplifies to

(16)

For a cubic primitive cell (), this result is equivalent to Figueirido’s, as expected from the DCT analysis of Ref. Hünenberger and McCammon, 1999. Fig. 1 exemplifies the remarkable effectiveness of this correction, echoing the conclusions of Hummer et al. Here we have added 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 down to . Residual finite size effects beyond 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 Hummer, Pratt, and García (1997); Hünenberger and McCammon (1999).

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. 1 we include simulation results for cuboidal simulation cells with and , demonstrating that the DCT correction remains accurate even as the cell’s aspect ratio becomes large. The interesting aspect of this comparison is that diverges as with 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 extend far beyond the plotted energy scale (), hinting at such a divergence.

Figure 1: System size dependence of solute charging free energy , as a function of solute charge , for periodic bulk simulations. (a) Results for the simple polar fluid (with 48, 192 and 384 molecules in cubic simulation cells) exhibit strong finite size effects. These variations with cell size are largely reconciled in (b) by the dielectric correction from Eq. 15. Higher-order corrections due to the finite size of the solute Hummer, Pratt, and García (1997); Hünenberger and McCammon (1999) have been neglected. (c) Results for SPC/E water also vary substantially with cell size, and also with its aspect ratio. The filled circles were obtained with 64, 256 and 512 molecules in cubic simulation cells, while the empty squares were obtained with 128 and 756 molecules in cuboidal cells ( and , respectively). Note that the results for 756 molecules extend far beyond the data range shown (e). (d) Corrected free energies for SPC/E water are nearly indistinguishable on this energy scale.

iii.3 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 (as well as integration of over disconnected volumes). The spatial variation of 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 the boundary, accompanied by strong spatial variation in polarization surface charge. A poor description of this position dependence would yield a poor approximation to . This feature of , 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 , 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 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 . 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 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:222We 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).

(17)

where denotes a single periodic replica of the slab (see Fig. 2). In using the divergence theorem, we have exploited the fact that symmetry and continuity of within the slab requires that at the lateral interfaces between periodic replicas. Use of Eqs. 7 and 8 then yields

(18)

where is the slab’s thickness and is the area of a single periodic replica of the upper (or lower) dielectric boundary. This value of is precisely what is required for to be finite.

The finite size correction that results from this uniform surface charge approximation, together with Eq. 18,

(19)

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 with and fixed, all of the compensating charge moves to these boundaries. This scenario – a single dielectric slab, with ions periodically replicated in the and directions, bounded above and below by uniformly charged neutralizing plates – is markedly similar to the neutralization scheme proposed by Ballenegger et al. Ballenegger, Arnold, and Cerdà (2009) for molecular dynamics simulations of charged slab systems that are aperiodic in the 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

(20)

In the limit that , , and all become infinite, with still fixed, , an important property of the full solution for that is preserved by the uniform surface charge approximation. Note that does not depend upon the solute’s position within the slab.

Figure 2: Schematic of the slab geometry typically used in molecular simulations of interfacial systems. In the case of a liquid coexisting with its vapor, a slab of width (blue) is separated from its periodic replica normal to the interface (taken to be along ) by a region of vacuum (orange). The length of the primary unit cell in the -direction is . The liquid/vapor interface lays in the -plane, and has a surface area , where and are the lengths of the simulation cell in and , respectively. The total volume of the simulation cell is , and the volume of the slab per cell is .

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. Fig. 3 (b) shows results for  Å, slab width  Å, and various aspect ratios , 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 .

Figure 3: System size dependence of solute charging free energy for periodic slab simulations. In all cases, the solute’s position was restrained with a harmonic bias potential whose minimum coincides with the slab’s center of mass. The lateral dimensions  Å  of these simulation cells are comparable to the cubic cell length for 64 water or 48 simple polar molecules (see Fig. 1). (a) Results for the simple polar fluid (with 96 molecules and  Å) depend strongly on . (b) Applying the dielectric correction from Eq. 20 removes nearly all size dependence for this charge-symmetric solvent. For comparison, we also include the bulk results obtained with 48 molecules (empty circles). (c) Results for SPC/E water (with 128 molecules and  Å) also depend strongly on . The correction from Eq. 20 in this case removes much of the observed size dependence, but a systematic dependence on remains.

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  Å, slab width  Å, and various aspect ratios , with the solute placed equidistant from the two interfaces. The variation of with is greatly reduced by , but a systematic size dependence remains, leaving differences as large as 13.4 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 the next section.

iii.4 Finite size effects due to charge asymmetry

The simple polar fluid we have considered is charge symmetric, in the sense that 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 anionsWilson, Pohorille, and Pratt (1989); Harder and Roux (2008). This artifact reflects a sensitivity of 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 . The slab geometry is an intermediate case, with explicit interfaces but also artifacts and finite size effects of charge asymmetry that vanish in the limit. To demonstrate this fact, we examine the average electric potential 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:

(21)

where denotes a single unit cell of the periodic slab simulation. is the solvent’s average charge density at position , with a neutral solute at position . Due to volume exclusion, vanishes inside the solute (and its periodic replicas). If the solute resides within the slab, this constraint requires that depends on and as well as , a significant feature for contributions to from the primary simulation cell (). Our interest lies, however, in contributions from more distant unit cells, which dictate the system size dependence of . These distant contributions are much less sensitive to the solvent’s electrostatic inhomogeneity in and . We thus expect that

(22)

is a good approximation, where is independent of , and the quantity has been spatially averaged over and . Fig. 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, the other at the liquid/vapor interface) are nearly indistinguishable. Provided the solute’s cross-sectional area is not a substantial fraction of , we expect the laterally averaged charge density to be insensitive to the solute’s presence at all.

Figure 4: Collective effects of charge asymmetry of an SPC/E water molecule. (a) The neutral cavity potential, relative to that of an infinite simulation cell, decays with aspect ratio as predicted by Eq. 25 (dotted line). Symbols show simulation results (with 128 water molecules,  Å, and  Å) for the solute at the center of the slab () and close to the liquid/vapor interface (). Asymptotic values of were obtained by fitting to the form , where and are free parameters. (b) Average solvent charge density profiles obtained from simulations with a neutral solute. Simulation parameters are the same as in (a). Results are shown for two solute positions, in the middle of the slab and close to the liquid/vapor interface.

Removal of the dependence on the lateral coordinates allows for great simplification. As shown in Appendix B, exploiting the symmetries of yields

(23)

The difference quantity

(24)

can now be written compactly,

(25)

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

Eq. 25 describes an intrinsic bias on ion solvation in a system of periodic slabs, one that vanishes in the limit . 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:

(26)

Eqs. 20 and 25 provide accurate, inexpensive, and easily implemented approximations for and , 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. Eq. 26 accounts equally well for finite size effects observed for an ion placed at the air/water interface, as shown in Fig. 5 (b).

Figure 5: Charging free energies for solutes in periodic water slabs, corrected for size dependence of both dielectric response and effects of charge asymmetry. Results are shown for simulations with  Å and  Å, and for two solute positions: (a) equidistant from either interface (), and (b) at the air/water interface (). For comparison, the empty circles show the bulk results obtained with 64 molecules (see Fig. 1 (d)).

Several features of the neutral cavity potential are noteworthy. First, within the approximation leading to Eq. 25, 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 under PBC. Indeed, Eq. 25 resembles the conventional surface potential , where and denote locations on either side of a neat liquid/vapor interface (). 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 . From this relationship our numerical results for suggest a surface potential of roughly  V for the SPC/E model, comparable to published values Sokhan and Tildesley (1997); Arslanargin and Beck (2012); Harder and Roux (2008); Wick et al. (2007); Lee Warren and Patel (2007). Like the conventional surface potential, may therefore have little to do with the molecular physics of solvation, as underscored by the possibility that even for spherically symmetric 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 . Specifically, if the neutral solute resides in the vapor phase 333We have in mind that the solute is sufficiently far from the interface that capillary fluctuations are not disrupted., , 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 in slab simulations is simply to compute the average electric potential at a point within the vapor phase (and far from the interface).

iii.5 Reconciling bulk and slab results

For the periodic slab geometry, with ions at the slabs’ centers, the limit 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

(27)

where is the result of periodic slab simulations with the corresponding finite periodicity correction applied. As shown in Fig. 6 (a) this adjustment brings slab and bulk simulation results into good agreement for the simple polar solvent.

Figure 6: Macroscopic charging free energies, for a solute that is distant from any interface, computed from slab and from bulk simulations. (a) For the simple polar solvent, slab results (with  Å) were corrected according to Eq. 27, yielding excellent agreement with the bulk values. (b) For water, slab results (with  Å) were again corrected according to Eq. 27. Bulk results in this case were also corrected, using the surface potential as in Eq. 29. We also show SPC/E results for a system with larger lateral periodicity (see legend), demonstrating that effects of finite have been addressed as well as those of finite . The aspect ratio is for all slab results shown.

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 the previous section, such bulk periodic systems lack interfaces even in the limit , which can generate artifacts correctible only through interfacial considerations. (See Ref. Duignan et al., 2017 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

(28)

where denotes a point in the liquid far from both the solute and the liquid/vapor interface. The first term in Eq. 28 should be well described by a bulk periodic simulation as , while the latter contribution is absent. Following the reasoning of the previous section, we account for the interfaces missing in a bulk periodic simulation by adding to the computed charging free energy. Our estimate of from bulk periodic simulations is thus

(29)

where is the result of periodic bulk simulations with the corresponding finite periodicity correction applied. Fig. 6 (b) shows the resulting macroscopic estimates for SPC/E water, where has been treated as a fitting parameter. A value of  V – comparable to previously reported surface potential values Sokhan and Tildesley (1997); Arslanargin and Beck (2012); Harder and Roux (2008); Wick et al. (2007); Lee Warren and Patel (2007) for point charge models, and our estimate based on – brings slab and bulk simulation results into good agreement. This convergence 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.

Iv Discussion

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 . 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 ) 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, deviates strongly from the parabolic form characteristic of linear response. (Most notably, the curvature differs significantly at positive and negative , so that anions are far more favorably solvated than cations Hummer, Pratt, and García (1996); Hirata, Redfern, and Levy (1988); Grossfield (2005); Bardhan, Jungwirth, and Makowski (2012); Lynden-Bell and Rasaiah (1997).) 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 Benjamin (1991); Venkateshwaran, Vembanur, and Garde (2014); Otten et al. (2012); Noah-Vanhoucke and Geissler (2009); Kumar, Knight, and Voth (2013). And yet a continuum linear response theory describes nm-scale thermodynamic finite size effects to an accuracy of . From the non-parabolic shape of 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. Makov and Payne, 1995; Freysoldt, Neugebauer, and Van de Walle, 2009; Komsa, Rantala, and Pasquarello, 2012; Komsa and Pasquarello, 2013; Andreussi and Marzari, 2014). 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 Cheng et al. (2014); Le et al. (2017).

V Conclusions

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 this absence becomes inconsequential. By contrast, for periodic bulk simulations the absence of interfaces is problematic even in the limit . 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 Levin (2009); Levin, dos Santos, and Diehl (2009); Levin and dos Santos (2014); Tobias et al. (2013); Baer et al. (2012); dos Santos and Levin (2013), in essence to arbitrarily small scales. The methods we have described facilitate a careful scrutiny of this hypothesis, which we pursue in ongoing work.

Acknowledgements.
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 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 Jackson (1999). Focusing on solvation near liquid-vapor interfaces, we take one of the regions to be vacuum () and the other region to have a large dielectric constant . For a point charge that is located in the high-dielectric medium at a point , lying a distance 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 is then specified by the requirement that vanishes at infinite distance from the real charge. For a point within ,

(30)

where

(31)

is the potential generated by an image charge at a point that mirrors the location of the real charge, and . 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 from its closest interface (the top interface in Fig. 7), and a distance from the other interface (the lower interface in Fig. 7), where is the slab thickness. Above the slab reside image charges: (i) at positions , and (ii) at positions , where is an outward normal for the upper interface. Below the slab reside image charges: (i) at positions , and (ii) at positions . Since the electric potential must again vanish as , is still given by Eq. 30, now with

Finite size corrections in Sec. III.3 involve the potential evaluated at the solute’s position. For this case of a single ion in a single slab, can be written exactly in terms of the Hurwitz-Lerch transcendent :

(33)

The potential 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 ():

(34)

In Fig. 8 we plot the non-singular part of this potential, , for various slab widths.

Figure 7: The method of images applied to a point charge in a dielectric slab. The slab (light blue) has finite width and is surrounded by vacuum (orange) on either side. The charge (white circle) is located a distance from the upper interface, and a distance from the lower interface. Dark blue circles represent image charges, which collectively generate an electric field in the dielectric that satisfies required boundary conditions on both interfaces. Their charges follow a geometric progression with increasing distance from the real charge. For odd, for images above the slab and below, as shown by the solid black lines. For even, on both sides of the slab, as shown by the dotted black lines.
Figure 8: The potential due to image charges at the location of a point charge in a slab of width , as given by Eq. 33. The point charge is located a distance from the nearest interface. Both and are reported in units of Bohr (); for convenience, in Ångstrom is given in parentheses in the legend. The solid black line shows the result for the water slab used in the main article. The solid gray line shows the result for the semi-infinite slab given by Eq. 31.

The case of central interest in this paper – the periodic slab simulation geometry depicted in Fig. 2 – is more complicated in two respects. First, image charges that determine the electric field do not follow a simple pattern. Secondly, need not vanish as , 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 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 , which is in turn completely specified by the electric field:

(35)

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 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 ). The charges and positions of more distant images are complicated by the system’s periodic structure (although simple recursion relations can be derived when 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 . Those background contributions, however, are spatially uniform in and , i.e., they are independent of on each boundary. We exploit this fact by splitting into a zero-wavevector component and a component that varies with and integrates to zero. The average surface charge density can be calculated simply and without reference to image charges, as described in Sec. III.3. The remainder receives no contribution from the background charge or its images.

Results presented in Sec. III.3 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 through come from the nearest image charges, which do not depend on and thus cancel in the calculation of . Contributions to 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 for large . For example, the physical charge is largely offset by an image charge opposite the boundary, with the residual . 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 . Contributions from nearby charges also do not offset to order in this case. For example, the image charge nearest the interface contributes to with magnitude , 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 , 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 Benjamin (1991, 1993); Venkateshwaran, Vembanur, and Garde (2014), a response that lies distinctly outside the scope of DCT.

Figure 9: The method of images applied to a periodic collection of point charges in a periodic set of dielectric slabs (the periodic generalization of Fig. 7). The image charges in the vacuum regions immediately above and below the slab have the same values and positions as in the non-periodic case. Beyond these regions, boundary conditions generate a more complicated set of values and locations for the images (depicted as open circles). Because is an integer in this case, the locations are discrete and regularly spaced, and the charge values can be written recursively.

Appendix B Generalizing the Wigner potential

Consider a periodic charge distribution that is (i) electroneutral, ; (ii) non-dipolar, ; and (iii) spatially uniform in two directions and , . The net potential generated by this distribution at observation point is

(36)

Exploiting the lateral uniformity of , this potential may be written as a Fourier series in just one dimension,

(37)

where , , and the sum runs over all nonzero integers . This sum can be rewritten

(38)

in terms of a function that can be expressed in closed form,

(39)

The second equality in Eq. 39 can be shown by twice differentiating , noting that , and solving the ordinary differential equation that results. Doing so is aided by recognizing that , that , and that is an even function of . The assumed symmetries of allow further simplification:

(40)

Because the neutralizing charge distributions encountered in Sec. III.3 are uniform in and , their contributions to can be simplified using Eq. 40. In particular, the generalized Wigner potential that arises for periodic dielectric slabs differs from the conventional Wigner potential only in the form of the neutralizing background. Defining

(41)

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