A Continuum Model of Solvation Energies Including Electrostatic

Article pubs.acs.org/JPCB

A Continuum Model of Solvation Energies Including Electrostatic, Dispersion, and Cavity Contributions Timothy T. Duignan, Drew F. Parsons,* and Barry W. Ninham Department of Applied Mathematics, Research School of Physical Sciences and Engineering, Australian National University, Canberra ACT 0200, Australia ABSTRACT: Physically accurate continuum solvent models that can calculate solvation energies are crucial to explain and predict the behavior of solute particles in water. Here, we present such a model applied to small spherical ions and neutral atoms. It improves upon a basic Born electrostatic model by including a standard cavity energy and adding a dispersion component, consistent with the Born electrostatic energy and using the same cavity size parameter. We show that the well-known, puzzling differences between the solvation energies of ions of the same size is attributable to the neglected dispersion contribution. This depends on dynamic polarizability as well as size. Generally, a large cancellation exists between the cavity and dispersion contributions. This explains the surprising success of the Born model. The model accurately reproduces the solvation energies of the alkali halide ions, as well as the silver(I) and copper(I) ions with an error of 12 kJ mol−1 (±3%). The solvation energy of the noble gases is also reproduced with an error of 2.6 kJ mol−1 (±30%). No arbitrary fitting parameters are needed to achieve this. This model significantly improves our understanding of ionic solvation and forms a solid basis for the investigation of other ion-specific effects using a continuum solvent model.

I. INTRODUCTION The interpretation of a suite of standard measurements and our understanding of electrolyte solutions depends inextricably on classical continuum solvent theories. It is sufficient to just mention the words Born or Debye−Hückel theory to conjure up the foundations that inform our intuition. Whether it be the setting and prescribing of pH via buffers, the measurements of interfacial tensions, osmotic pressures, ion binding, membrane and zeta potentials, solubility and viscosity, colloidal stability, or protein precipitation, Born and Debye are not too far removed.1 That electrolyte solutions occupy center stage in nature and industry is a trite observation. Less obvious is the fact that 200 years since Berthelot’s then counterintuitive observations of soda lime rocks in the receding Nile river,1 and 150 years since Hofmeister’s work,2 we still remain bemused and stymied by specific ion effects. The only certain thing is that any progress has to begin with a quantitative theory of the solvation free energies of ions.3 These properties are important in their own right, e.g, in ion exchange, purification, and physiological applications. But they are also important in laying the foundations on which to build better models of electrolyte solutions. Background. Solvation energies measure the change in free energy of a solute in a vacuum compared to that of the solute immersed in water. A large and venerable literature3,4 is dedicated to their calculation for ions (and molecules), for different sizes and shapes, using a variety of methods. The free energy change arises from three sources: The first two are the electrostatic and dispersion interaction energy of a solute molecule with water. The third is the change in the free energy of the water due to the rearrangement and displacement of the © 2013 American Chemical Society

water molecules around the solute particle. This includes the cavity formation energy. Significance. Determination of the free energy changes for many processes first requires an understanding of solvation free energies, in the simplest case of bulk water. This is because the change in the ion−water and water−water interactions makes up a significant contribution to the total free energy change in many processes. For instance, knowing the free energy change on moving an ion close to a surface or to another particle allows the prediction of ionic distribution profiles which determine many important macroscopic properties. Some examples are an ion moving closer to air or a mineral surface, colloidal particle, protein, or another ion, which allows the prediction of surface tensions and forces, colloidal stability, protein precipitation, and activity and osmotic coefficients. The calculation of these energies is an important and long-term goal, but is hampered by an incomplete understanding of the ion−water interaction. As a result, ad hoc corrections need to be made with parameters that are adjusted to fit theory to experiment. “Hydration forces” are an example of this. They are not accurately modeled, and must remain an ad hoc correction as long as the theory omits major contributing forces.5 Modeling Approaches. A recent book3 gives an overview of the situation. In this section we provide some context to the model presented in this paper by discussing some specific models antecedent to it. This is followed by a general discussion of the current approaches found in the literature. Received: April 11, 2013 Revised: June 20, 2013 Published: July 9, 2013 9421

dx.doi.org/10.1021/jp403596c | J. Phys. Chem. B 2013, 117, 9421−9429

The Journal of Physical Chemistry B

Article

Born Model. The first and most important is the Born model.6 This model gives the electrostatic contribution to the solvation energy due to polarization of the solvent from the ion’s charge. The ion is modeled as a point charge placed in a cavity inside a dielectric medium. The input parameters it requires are the static dielectric constant, the ionic valency, and the cavity size.7 It includes both an enthalpic and entropic contribution, as the dielectric constant and possibly the cavity size are temperature dependent. Its extension to include ionic correlations gives Debye−Hückel theory.8 The Born equation is G Born = −

q2 ⎛ 1⎞ ⎜1 − ⎟ εr ⎠ 8πεoR cav ⎝

(MSA), which treats the solvent on a molecular level. This model can also be used to derive a value of Radd = 0.54 Å.13 There is some debate regarding the best experimental solvation energy values to compare against. This will be discussed below. The experimental values used here are those of Tissandier et al.15,16 Cavity Size. The cavity radius represents the distance at which the polarization distribution of the water effectively begins. Various physically based arguments for why this Radd parameter is different for cations and anions have been put forward.3,12,14,17 In brief, the argument is that the distribution of water dipoles should be considered to be further away from cations than anions due to the orientation of the water molecules in the first solvation layer. This difference in solvation energy of similarly sized anions and cations is a problem even in more complex continuum models. For instance Marcus’s solvation energy model18 included an arbitrary correction which shifted the solvation energy of cations and anions a constant amount in opposite directions to improve agreement with experiment. The need for this correction was again attributed to the difference in orientation of the water molecules in the first layer. This fact is often cited as evidence of the importance of water structure on electrolyte properties. It is argued that the specific size and orientation of the water molecules play an important role. Hence it is concluded that explicit solvent models are necessary to understand the properties of electrolyte solutions.3,17,19 An objection to this argument can be raised based on several observations. The first is clear if we look at the sodium and fluoride ions. These ions have almost identical solvation energies.16 However, if we look at the distance from the center of these ions to the center of the oxygen atoms in the nearest water molecules we see that the water molecules are substantially (0.33 A) further away from the fluoride ion. So, clearly the difference in orientation is not allowing the water molecules to approach any closer. The second is that the electrostatic potential created around a water molecule, given in ref 20, is in fact roughly symmetric around a water molecule. This is true in the sense that at a distance of 2.75 Å from the central oxygen atom in a water molecule the maximum positive and negative potentials are equal in magnitude (±100 kJmol−1). This implies that the electrostatic interaction energy of an ion of fixed size with a water molecule should be the same irrespective of whether it is positively or negatively charged. This means Radd should be the same for cations and anions. An additional problem with attributing different Radd parameters on the basis of different water orientation is that it does not explain the variation in this parameter that occurs with cations. Consider the solvation energy of a monovalent silver ion and of a sodium ion. Both have essentially the same size,21 but silver has a solvation energy 65 kJmol−1 more negative then sodium. The water molecules will have the same orientation in both cases. So the use of a different Radd cannot be justified with this argument. This then raises the question of how to account for the clear difference in the solvation energies of two ions of the same size. A common answer is to argue that “chemical” interactions occur between the ions and water. In the context of a many body system like water, it is not at all clear what the distinction between chemical and physical interaction means. It would be unfortunate if it were true. Presumably, it would mean that ab initio explicit water simulations will be necessary in order to

(1)

Here Rcav is normally interpreted as the size of the cavity in water that the ion creates. Although this model is much criticized and rather simplistic, there is some justification for believing it is relatively accurate in capturing the essence of electrostatic interactions, at least for monovalent ions.9−11 Given its simplicity and neglect of other contributions, this model is surprisingly successful at reproducing the experimental solvation energies. Indeed, several calculations of solvation energies based solely on this model have been put forward. They adduce various justifications for choices of the cavity radius. These are typically related to the crystal radius of the ion, referred to as RI, by taking Rcav = RI + Radd. Here, Radd is some additional cavity increment which may be different depending on whether it is for a cation or anion, with some physical rationalization.7,12,13 Or alternatively the distance to the peak (RS) in the ion-oxygen radial distribution function (RDF) is used instead of the crystal radius: Rcav = RS − Radj.14 The success of this schema can be seen in Figure 1 which plots the solvation energy as a function of radius, where two

Figure 1. Born model solvation energy of the alkali halide ions compared to experiment. Plotted as a function of the crystal radius of the ion RI. The ionic cavity radius is Rcav = RI + Radd.

fitted values of Radd have been chosen. With Radd = 0.35 Å the Born model roughly follows the anionic experimental solvation energies. With Radd = 0.55 Å it follows the cations almost exactly. The need for different values of Radd for cations and anions can be clearly seen from the smaller anions, which have a significantly larger solvation energy then a cation of the same size would. It is also possible to derive the Born equation from the apparently more sophisticated mean spherical approximation 9422



Article

solvent separation. This significantly complicates the situation when ions are concerned. Cations have a missing electron. This significantly reduces their polarizabilities relative to their size. Anions have an extra electron correspondingly increasing their polarizabilities. It is true that the electron cloud will have a corresponding change in size too, but this effect is not significant enough to counter the change in polarizability. An additional complication is that ions produce an electric field that means water molecules are pulled much closer to ions than they are to neutral solutes of the same size. The dispersion interaction increases very quickly at small separations and so this results in a significant increase in the dispersion energy. These complications illustrate the need for an independent and accurate calculation of the dispersion contribution to the solvation energy. To test that the cavity and dispersion contributions are being modeled accurately, it is important that any solvation model used to account for the solvation energies of ions can also give a sensible interpretation of the solvation energies of neutral molecules as well. Dispersion Model. We have recently developed a method which calculates the dispersion contribution to the solvation energy of a molecule at the same level of theory as the Born model. Significantly, this model includes higher order multipole moment contributions up to the octupole moment.27 The water is treated as a polarizable continuum with a spherical cavity occupied by a solute particle. An experimentally determined dielectric susceptibility function for water is used outside the cavity, and the properties of the solute (size and multipole polarizability) are determined independently, from a combination of experiment and ab initio calculations.27 As this model uses the same assumptions as the Born model, it is consistent to use the same cavity size for both Born and dispersion calculations. Model Classification. Modern models can be classified by making two distinctions, following Zuo et al.28 The first is the obvious distinction between explicit and continuum solvent models, for example, molecular dynamics simulation29 compared with dielectric continuum models.30 The second classification is between physically based and fitted models. This distinction relates to the determination of the parameters of the model. Physically based models can determine the parameters in two ways. The first way is from fundamental methods such as a quantum chemical calculations of the solute charge distribution.30 A second way is from independent experimental measurements such as the dielectric susceptibility and crystal ionic size data.18 In contrast, fitted models determine the parameters by direct comparison with solvation energies; for example, the Lennard-Jones parameters of classical MD simulation29 and the surface tension parameters in continuum solvent models.30 Clearly, the physically based models are more desirable than fitted models. If they agree with experiment this is much stronger evidence that the physics has been correctly modeled. They can thereby provide a foundation for modeling more complex phenomena. However, as mentioned above, some authors argue that for ions in solution accurate physically based continuum water models are impossible. This would mean that modeling must resort to either fitted models or computationally costly explicit water models, which can be difficult to interpret. Goal. This paper presents a continuum solvent model of solvation energies which is based on sound physical arguments. The parameters input into the model have good independent

model even the simplest properties of ions in water. We therefore seek simpler models and explanations of these solvation energies without attributing them to unquantifiable “chemical” interactions. Additional Contributions: Dispersion and Cavity Energies. A problem with using the Born model to calculate solvation energies is that it leaves out other known contributions. These contributions can be roughly categorized as cavity formation and dispersion energy contributions. Both have been shown to be significant.22,23 They need to be modeled separately and included along with the Born contribution. The cavity contribution has an opposite sign to the dispersion contribution, and there is reason to believe they will be of comparable magnitude too,24 and so will largely cancel. In fact this cancellation is very likely the reason for the surprising successes of the Born model. This means that both of these contributions must be taken into account to improve experimental agreement. Both contributions increase with solute size. This is consistent with the fact that the Born model is less accurate for the larger anions compared with the smaller cations. With the above as background, we now attempt to remedy the deficiencies of continuum solvent theory. We present a model which shows that the failures of the Born model can be fixed without having to resort to arguments regarding specific orientation of water molecules. This is achieved by including an accurate calculation of the dispersion and cavity contributions. Qualitatively the argument is that the larger solvation energy of fluoride and silver compared with potassium and sodium, respectively, is explained by the fact that fluoride and silver are much more polarizable than potassium and sodium. Hence they will have larger dispersion contributions to the solvation energy. Incorporating an accurate description of the dispersion interaction of ions with the surrounding water validates this argument and removes the need to vary the Radd parameters of the ions. This is a much more satisfactory explanation for the discrepancy. It is conceptually simple, and in principle it does not necessarily require time-consuming and highly nontrivial explicit water simulation. However, these methods may be useful in determining the size parameters where X-ray data is not available. Further, it does not depend on arbitrarily adjusted parameters to fix the difference in solvation energy between ions of similar size. Conceptually, the model we develop is unavoidable. A model of electrolyte solutions that ignores half of the forces operating, those due to quantum mechanics, in favor of electrostatics is clearly deficient. Simple Models. One approach in the literature to account for the cavity and dispersion contributions is to assume that they are functions of the size of the solute. For instance, some models18,25 use the noble gas solvation energies. They argue that with a consistent definition of ion and gas particle size the contribution to the ionic solvation energy from these two contributions can be determined by interpolation of the solvation energy of the noble gases as a function of size. Other models include the dispersion and cavity energies in a surface tension parameter.26 This approach is sensible for some contributions; for instance, we will argue below that a surface tension parameter is a sensible way to incorporate the energy cost of cavity formation. However, for the dispersion energy this approach is not reasonable. The dispersion energy is a function of the polarizability, the size of the solute particle, and the solute− 9423



Article 10 ⎛ ⎞ 15k T (bR S)k GQ = −⎜⎜1 − ∑ exp−bRS ⎟⎟ B5 k! ⎝ ⎠ 2R cav k=0 ∞ ε(iξ ) − 1 ∑ ⎛⎝⎜1 − 1 δn0⎞⎠⎟α2(iξn) 3 n 2 ε(iξ ) + 1

physical rationalization and are not arbitrarily adjusted in order to reproduce the experimental values. Our model incorporates the newly developed calculation of the dispersion contribution to the solvation energy, which includes quadrupole and octupole moment contributions.27 Our goal is to accurately reproduce and explain the solvation energy of ions and atoms, and hence provide physical insight, and a stepping stone for models of more complex ion-specific effects like activities and membrane potentials, while also conserving computational time and conceptual simplicity. A continuum model is used for several reasons. The first is that it is easy to implement, as it has the obvious advantage of simplicity and potentially low computational demand. The second is that accurately calculating electrostatic and dispersion interactions is feasible. Specifically, a Lifshitz style model can be used to take into account all of the complex many body interactions using experimentally derived dielectric data. The electrostatic interaction is treated at the same level of theory for consistency, reducing the number of additional parameters. Although this is not an explicit solvent model, the methods presented and discussed here still have implications for the modeling of dispersion interactions in explicit solvent approaches. They may be useful in a hybrid approach or may provide a qualitative indication of what to expect from explicit water models.

n=0

n=0

n

14 ⎛ ⎞ 28k T (bR S)k GO = −⎜⎜1 − ∑ exp−bRS ⎟⎟ B7 k! ⎝ ⎠ 3R cav k=0 ∞ ε(iξ ) − 1 ∑ ⎛⎝⎜1 − 1 δn0⎞⎠⎟α3(iξn) 4 n 2 ε(iξ ) + 1 n=0

3

n

(3)

(4)

Cavity Contribution. Calculating the free energy cost of creating a cavity in the solvent for a solute to occupy is a much studied, but not well resolved, problem.23,36 Evidence from molecular dynamics simulations shows that for small hard sphere solutes, of the size considered in this paper, the effective surface tension reduces to zero roughly linearly with the radius. This means that for these small solutes the energy of cavity formation scales with the volume of the cavity.37,38 The hard sphere solutes used in these calculations are very similar to the noble gas solutes because the dispersion interaction does not significantly alter the structure of the surrounding water. They might then be expected to be a good model for the cavity formation energy of noble gas solutes. Because these models use a simple hard sphere solute to model the cavity, the volume of the cavity is defined as the region where the centers of the water molecules are excluded. This means that the most appropriate choice of radius for this cavity is RS, as this gives the average distance to the center of the closest water molecules. We therefore model the cavity formation energy of the noble gas elements with the following expression:

II. THEORY Electrostatic Contribution. We include the Born model in this calculation to account for the electrostatic contribution to the solvation energy.6 The expression is eq 1. Many modifications to this model have been proposed31 to incorporate dielectric saturation,32 electrostriction,33 or other such higher order effects.34 These methods often involve a spatially dependent dielectric function. However, there is no clear consensus in the literature regarding the best method for incorporating these effects or how important they are. For this reason, as well as to keep the model simple, and to minimize the number of fitted parameters, we use the unmodified Born model. Dispersion Contribution. Here, we apply a model of the dispersion contribution to the solvation energy which we have recently developed.27 The model uses the frequency dependent multipole polarizabilities of solute particles (αn(iξn)) in order to calculate the multipole dispersion interaction of a solute particle with the surrounding water up to the octupole moment. The water is modeled as a dielectric function (ε(iξn)) with a spherical cavity of size Rcav occupied by the solute. TangToennies damping functions35 are included which reduce the size of the dispersion energy in order to account for solute− solvent wave function overlap. These functions require the Born-Mayer (b) parameters as input, which are a measure of the spatial extent of the molecular wave functions, as well as the ion−solvent separation parameter, RS. This calculation is based on the same assumptions as the Born model, allowing for their combination. The equations for dipole, quadrupole, and octupole, respectively, are 6 ⎛ ⎞ 6k T (bRS)k GD = −⎜⎜1 − ∑ exp−bRS ⎟⎟ B3 k! ⎝ ⎠ R cav k=0 ∞ ε(iξ ) − 1 ∑ ⎛⎝⎜1 − 1 δn0⎞⎠⎟α(iξn) n 2 2ε(iξn) + 1

2

G NG − Cav = σNGR S3

(5)

Here σNG is a parameter determined from explicit water simulation.38 We use the same value as ref 39, σNG = 0.73 kJmol−1 Å−3. However the situation for the ions is quite different. The reason the effective surface tension of a hard sphere solute cavity reduces to zero is that the water molecules around very small sized cavities can form hydrogen bonded networks.37 This means that the cavities do not break the structure of the water, and hence there is a smaller energy cost to their formation. However, for small ions this is not the case. The ion’s electric field reorients the water molecules, breaking the hydrogen bonded structure and costing energy. This is similar to the effect that very large solutes have on water, where the bulk surface tension parameter is accurate.37 This means that it is more appropriate to use the bulk surface tension in determining the cavity formation energy of ions. We therefore use the expression

G Ion − Cav = σIon 4πR S2

(6)

For σIon, we use the surface tension of bulk water σion = 0.434 kJ mol−1 Å−2.3 Input Data. We compare the solvation energies calculated with this model against the experimental ionic solvation energies taken from ref 16 and noble gas solvation energies taken from ref 40. Cavity Size. The cavity size used for electrostatic and dispersion contributions is Rcav = RS − 0.84 Å. RS is the distance

(2) 9424



Article

Proton Solvation Energy. The values for the experimental solvation energies of individual ions used here are the so-called ‘benchmark experimental values’49 of Tissandier et al.15,16 They are significantly shifted from the values of Marcus50 which are those most often used in the specific ion effects literature. This difference is due to the difficult nature of splitting the measured salt pair solvation energies into two contributions from the individual ions. Fixing the value of the solvation energy for one ion (normally the proton) allows all the other ionic solvation energies to be determined from the pair data. There are two important implications of these newer experimental values. One is that the gap between anions and cations of the same size is now significantly reduced. This makes it more plausible that the discrepancy is due to dispersion interactions. It also means that justifications for different choices of the Born radius need to be reconsidered as they have been chosen to agree with experimental data which may not be correct. Another important implication is that these new experimental results seem to show a qualitatively different trend from the Born model predictions as seen in Figure 1. There is no choice of Radd that gives adequate agreement with the experimental halide ion solvation energies. This indicates that there are additional contributions to the solvation energy which the Born model neglects. Surface Potential. A complicating factor in comparing theoretical and experimental solvation energies is the contribution from the surface potential of water. This is the potential difference over a water interface due to a charge imbalance. There is ongoing debate in the literature as to whether Tissandier’s values include some contribution from the surface potential.51 If there is some contribution of this type to Tissandier’s values, then this would need to be corrected for before the values could be compared with a continuum type model such as this one, which neglects any such contribution. In this paper we compare with the uncorrected estimates of Tissandier’s values given in ref 16. We agree with the argument of Tissandier, Hunenberger, and Reif52 that these values do not include any surface potential contribution. An accurate value can be obtained from clusters containing only one water molecule. The notion that there will be a significant contribution from the surface potential of bulk water at this small a cluster size seems unlikely. In addition, any surface

between the centers of the solute and the water molecules in the first layer, taken as the position of the peak in the soluteoxygen radial distribution function. The determination of these values using ab initio and experimental methods is given in ref 27. The values are given in Table 1. The values Radj = 0.84 Å is Table 1. Solute Crystal Size of Coordination Number 6 (RI), Solute−Oxygen Distance (RS), Born-Mayer Parameter (b), and Cavity Size (Rcav) atom/ion +

Li Na+ K+ Rb+ Cs+ Cu+ Ag+ F− Cl− Br− I− He Ne Ar Kr

RI(Å) a

0.76 1.02a 1.38a 1.52a 1.67a 0.77a 1.15a 1.33a 1.81a 1.96a 2.20a 1.08b 1.21b 1.64b 1.78b

RS(Å) 2.06 ± 2.35 ± 2.79 ± 2.92c 3.11 ± 2.25h 2.41 ± 2.68 ± 3.20 ± 3.35 ± 3.64 ± 3.08i 3.21 ± 3.55 ± 3.73 ±

c

0.02 0.02c 0.03c 0.04c 0.01c 0.06c 0.01c 0.03c 0.01c 0.01d 0.01d 0.03d

b(Å−1)

Rcav(Å)

5.49e 5.00e 4.33e 4.06e 3.77e 5.47e 4.76e 4.42e 3.47e 3.18e 2.71e 4.89f 4.64f 3.84f 3.52f

1.22g 1.51g 1.95g 2.08g 2.27g 1.41g 1.57g 1.84g 2.36g 2.51g 2.80g 2.24g 2.37g 2.71g 2.89g

a

ref 41. bref 42. cAverage of values from refs 43−45. dref 46, 47. Interpolation of b versus. RI for noble gas particles. fref 48. gRS − 0.84 Å. hab initio calculation. iExtrapolated by linear relationship between RI and RS. e

derived from the expression Rcav = RI + 0.54 Å from MSA calculations of ref 13. Combined with the approximate empirical relationship, RS = RI + 1.38 Å.21 It is important to distinguish between the Rcav parameter which is used in the Born and dispersion energy calculations and the parameter used to determine the cavity formation energy which is RS. Rcav represents the region where the water is not significantly polarized. RS gives the average distance between the centers of the solute and the water molecules in the first solvation shell. It is therefore a more appropriate measure of the region from which the water molecules are excluded.

Table 2. Contributions to the Solvation Energies of Some Small Ions and Noble Gas Atomsa

a

atom/ion

Gexp

Gtheory

Li+ Na+ K+ Rb+ Cs+ Cu+ Ag+ F− Cl− Br− I− He Ne Ar Kr

−537.2 −431.8 −359.8 −337.2 −314.2 −591.2 −496.6 −436.8 −311.7 −285.8 −250.6 11.5 11.2 8.4 6.9

−546.1 −445.6 −355.4 −344.9 −324.3 −571.3 −524.9 −434.1 −325.1 −305.4 −258.2 14.8 12.2 4.5 4.5

electrostatic −562.1 −454.2 −351.7 −329.7 −302.1 −486.4 −436.8 −372.7 −290.6 −273.2 −244.9 0 0 0 0

(105%) (105%) (98%) (98%) (96%) (82%) (88%) (85%) (93%) (96%) (98%)

cavity 23.1 30.1 42.4 46.5 52.7 27.6 31.6 39.1 55.8 61.1 72.2 21.3 24.1 32.7 37.9

(−4%) (−7%) (−12%) (−14%) (−17%) (−5%) (−6%) (−9%) (−18%) (−21%) (−29%) (185%) (216%) (390%) (547%)

dispersion −7.1 −21.5 −46.2 −61.6 −74.9 −112.5 −119.8 −100.5 −90.3 −93.3 −85.4 −6.5 −12.0 −28.2 −33.4

(1%) (5%) (13%) (18%) (24%) (19%) (24%) (23%) (29%) (33%) (34%) (−57%) (−107%) (−336%) (−482%)

The rounded percentage contributions to the total are given in brackets. The free energies are given in units of kJ mol−1. 9425



Article

also agrees well with the noble gas solvation energies where the error is 2.6 kJ mol−1 (30%). The percentage error is quite large due to the significant cancellation of the cavity and dispersion contributions. This cancellation can be seen in Figure 3 and Figure 4 which give the contributions to the solvation energy.

potential contribution should change as the cluster size increases resulting in Tissandier’s estimate of the proton solvation energy varying with cluster size, but the estimate is independent of cluster size up to six water molecules.16 This implies that there is no surface potential contribution included in this estimate. Hunenberger and Reif53 also give a comprehensive review of the literature regarding the proper choice of the intrinsic solvation energy of the proton. They conclude, based on multiple lines of evidence, that the best choice is a value basically the same as the one used here. If this contribution were important it would not be clear at all how to correct the values to account for it. The correct surface potential of bulk water is still debated.54 In addition to this the surface potential contribution to the solvation energy will have a contribution from both the bulk air−water interface and the ion−water interface which are likely to significantly cancel.55 This means that knowing the bulk air−water interface potential would be insufficient to go forward. Finally, even if this contribution had been accurately calculated, there is little reason to believe that it will be the same in the very small ion− water clusters used in Tissandier’s experimental estimates.

Figure 4. Calculated contributions to the total ionic solvation energy.

III. RESULTS AND DISCUSSION The solvation energy values calculated with this method are given in Table 2. Figure 2 shows that the model gives good

In addition, Figure 3 and Figure 2 both show that the model reproduces all of the qualitative features of the solvation energies. Including the non-monotonic jumps due to Cu+ (RI = 0.77) and Ag+ (RI = 1.15). The large cancellation of dispersion and cavity formation energies has also been accurately captured. The fact that only one value for the Radj parameter is needed to fit the solvation energy of both cations and anions is significant. This indicates that differences in water orientation around the two types of ions do not play a significant role. The explanation of the unusually large solvation energies of the copper and silver ions with no adjustments to the model is an additional success. Unfortunately, as with any model, there is a degree of subjectivity involved. This means it is possible that prior knowledge of the experimental data has had some influence on the choice of parameters. However, we emphasize that this model shows quantitative agreement with solvation energies for eleven ions and four noble gas atoms, without the use of any parameters fitted to reproduce the data. The parameters all have good independent physical rationalizations. This addresses a common criticism of continuum solvent models, that they are simply equations with parameters fitted to the data, with little physical basis or predictive power. The independent origin of the parameters used in this model renders such objections impotent. Uncertainty. The minimum value for the error in the experimental values is 2 kJ mol−116 This is significantly smaller than the error associated with the theoretical values. One important source of error in the theoretical values is the poor convergence of the multipole expansion of the dispersion solvation energy for the halide ions. For fluoride even the octupole dispersion energy is quite significant. This means that it may be necessary to take into account even higher order multipole moments, which is difficult, or use an entirely different approach. With a simple linear extrapolation we can infer that the next higher order contribution for fluoride may be roughly −8 kJ mol−1. This is 2% of the solvation energy and is quite significant. In the case of the pairwise dispersion interaction between neutral atoms, one approach is to cut off the energy expansion at the octupole moment even if the energy has not fully

Figure 2. Calculated solvation energies compared with experiment for some ions.

agreement with experimental solvation energies. The mean absolute error is 12 kJ mol−1 (3%). Figure 3 shows the model

Figure 3. Calculated contributions to the solvation energy of the noble gases. Including comparison to the experiment values. 9426



Article

have much larger polarizabilities than similarly sized alkali cations, and hence they have a much larger dispersion contribution to the solvation energy. This result implies that we will profit from further attention focused on improving the description of dispersion energies. This task is already well under way. The continuum approach is clearly most suited to this challenge due to its ability to take into account complex many body interactions through experimentally determined dielectric properties. One possibility is to combine continuum model calculations of the dispersion interaction with explicit solvent calculations of the electrostatic interaction. This would be particularly useful in the case of ab initio explicit solvent simulations using DFT, which give notoriously inaccurate dispersion energies.57 Explicit Water. The success of this model has some hopeful implications. The first is that it may be possible to build successful models of ion-specific effects using continuum solvent models that have predictive and explanatory power. Here a relatively simple continuum model reproduces cation, anion, and neutral atom solvation energies. This means that many contributions which are often cited as being the cause of ion specific effects, such as specific chemical interactions with water, or specific orientational structures, cage-like or otherwise, that the water molecules form, can possibly be neglected. Indeed, the fact that the solvation energies show such a simple regular relationship with ion size already implies this fact. The consistency and completeness of the model is supported by the fact that the cation and anion discrepancy disappears once the proper proton solvation energy is chosen and dispersion and cavity contributions are taken into account. Explaining the discrepancy between cation and anion solvation is an important step in understanding ion-specific effects. This discrepancy is often used as a justification of the need for explicit water models,3,17,19 but our result implies that it is due to the dispersion interaction, as the solvation energy has been modeled in exactly the same way for both types of ions, and the dispersion energy is the only contribution which does not solely scale with the size of the solute. This suggests the continuum model should be successful in other applications. This statement is supported by successful continuum solvent models which have been developed in recent years. A model developed by Levin et al.39,58 can reproduce electrolyte surface tensions accurately. A dielectric model developed by Lund et al.19 qualitatively explains the law of matching water affinity. To a certain degree, however, the model developed here is still consistent with the claim that explicit water structure and orientation are important. These effects may play a role in the RS parameter and hence the size of the cavity. This is particularly the case for the silver(I) and copper(I) ions where there is a 0.1 Å deviation from RS = RI + 1.38 Å. The determination of this RS parameter may also depend on explicit water simulation if the ab initio calculations or experimental estimates are not sufficiently accurate. However, if the model could be generalized to calculate ion−ion and ion−surface interactions in water then the model would still have a significant computational advantage, as it is these calculations that can require significant computational time at an explicit solvent level. Further Applications. This solvation energy calculation can serve as a starting point for improving models of more complex phenomena. For instance, the size of the cavity, the magnitude of the cavity formation energy, and the ion−solvent

converged. This is because this gives a good approximation of the total energy. Beyond the octupole moment contribution there are various terms of opposite signs which will partially cancel.56 This argument justifies the approximation in this case as well. Another important source of error in the theoretical calculations is the experimental value of the ionic cavity size in water. The calculations are quite sensitive to the choice of these values. They were determined from the average solute− water distances derived from experimental results27 given in the literature.43−45 The literature values generally vary over a reasonably large range. Although the average value can be determined accurately, some systematic error in the measurements due to the finite concentration of the solutions or other effects is possible. The experimental radii were supported by ab initio calculations which largely agreed with the experimental values,27 but more accurate estimates of these values may explain some of the error in the model. This is particularly true for the noble gas atoms where the radii were based on only a few molecular dynamics simulations. One issue with the model is that the Tang-Toennies functions and its parameters (b and Rs) are not rigorously derived for use in this case. It is reasonable to think that equally plausible arguments could be made for other choices of the key parameters, and indeed this has been done for the cavity size parameter.12,14 However, the final justification for the use of these functions and parameters is the good agreement with experiment. However, this does mean that there is a small degree of implicit empirical parametrization in the model. There may also be additional contributions which this model neglects such as a repulsive wave function overlap energy or a chemical ion−water interaction. These contributions are very hard to accurately calculate and the success of this model implies they are not important. Implications. There are two important implications of this model. The first and most important is the importance of the dispersion interaction. The second is that the structure and orientation of water molecules can largely be ignored in calculating the solvation energy. Dispersion Interaction. Figure 4 clearly shows that for the anions as well as the silver and copper cations the dispersion contribution to the solvation energy is very significant. One third of the iodide ion solvation energy comes from the dispersion energy. In addition, and we believe crucially, the dispersion interaction is not simply a quantitative correction, but plays a role in a qualitative description of the solvation energies. The use of a different cavity size definition for cations and anions in previous models can be attributed to inaccurate modeling due to the neglect of the dispersion interaction. The larger solvation energy of an anion, compared with a cation of the same size, is due to the anion’s larger dispersion energy. In addition the dispersion energy explains the observation that the solvation energy of the silver(I) and copper(I) ions is significantly larger than expected based on the Born model and their size. The most common explanation of the discrepancy between the solvation energy of similarly sized cations and anions is based on the different orientation of the water molecules.3 This explanation cannot explain the unusually large solvation energies of copper(I) and silver(I). They are monovalent cations; the orientation of water molecules around them should not be very different than around a similarly sized alkali cation. This model provides a simple and consistent explanation, due to the fact that they 9427



Article

explicit solvent model. This is shown by the fact that no traditional MD simulation can reproduce individual ionic solvation energies and entropies simultaneously.29 The qualitative explanation for the difficulty can be understood from the fact that around noble gas solutes the water molecules are significantly restricted in the orientations they can have while still maintaining their (presumed) strong hydrogen bonded network. The larger the gas solute, the bigger this effect, and hence the lower the entropy. This effect is often attributed to the formation of clathrates which are the hydrogen bonded cages of water around a solute. With an ionic solute, however, the water molecules are oriented toward the charge. They do not form a low entropy network, but are simply individually strongly oriented in one position and hence have a low entropy. As the ion size increases they are less strongly oriented due to the weaker field and hence have a less negative entropy. This is opposite to the behavior of the noble gas atoms. This effect can be seen from the temperature derivative of the surface tension as a function of cavity size,64 which is negative for large cavities and positive for small ones indicating an order breaking effect of large cavities and an order making effect for small cavities. In addition, the entropies are a relatively minor contributor to the ion−water interaction. For instance, 97% of the iodide solvation free energy comes from the enthalpy,15 so the entropic effect is very small. Despite these problems, progress has been made generalizing this model to calculate ionic solvation entropies. Partial Molar Volumes. An additional important task which is currently being studied is calculating the partial molar volumes from this model. This task suffers similar problems to the entropy calculation. In order to calculate partial molar volumes the pressure dependence of the cavity surface tension parameter should be determined, which is difficult. Ref 3 points out that no models have considered this in calculating the partial molar volumes. In addition to this, the pressure dependence of the surface tension of water is very difficult to calculate,65 but attempts to do so3,65 have resulted in an unphysically large positive contribution to the partial molar volumes. The small negative contribution from the pressure dependence of the electrostatic and dispersion contribution is not enough to counter this. It is also worth pointing out that the partial molar volumes can be explained qualitatively by a simple model which calculates the volume of the ion’s cavity. This model does, however, slightly overestimate the volume change due to the neglect of electrostriction. This implies that including a correction from the electrostatic and dispersion terms may provide quantitative agreement. Again, progress has been made on this calculation using this model. However, complications due to dissolved atmospheric gas at low salt concentrations may inevitably keep us busy.

dispersion interaction are all important in determining the free energy change as an ion approaches a surface. These free energies determine the distribution profile of ions at the surface and hence the surface tensions and surface forces of electrolyte systems, which are central problems in physical chemistry. Ions at the Air−Water Interface. The propensity of large anions to adsorb at the air−water surface observed in simulation and experiment has been attributed to the large polarizability of these ions.59 However, an alternative case has been made that polarizability may not be so important,60 and instead, it is the change in the cavity contribution to the solvation energy which explains the observed surface affinity.58 The model of Levin et al.39,58 can reproduce electrolyte surface tensions while neglecting the significant dispersion repulsion of ions from a surface,61 which cancels out a significant part of the cavity energy. The approach used here may be useful in shedding light on this puzzle. The model of Levin et al. uses the hard sphere solute model to calculate the cavity formation energy, which we have argued significantly underestimates the actual ionic cavity formation energy. In addition, it uses the Born cavity size. This cavity size represents the region where the polarization of the water molecules is insignificant. This is different from the region where the centers of the water molecules cannot penetrate which is the definition used in a hard sphere solute model. The resulting underestimation of the cavity size likely explains why the model can neglect the repulsive dispersion contribution and still achieve reasonable results. Other Properties. The model can also, in principle, be applied to ion−ion interactions relevant to the determination of activity and osmotic coefficients. Likewise, it applies to ion− mineral or ion−protein interactions, determining ion distribution profiles and hence explaining a wide range of experimental phenomena. For instance, Lund et al.19 have shown that the law of matching water affinity can largely be explained by the change in ionic solvation energy as two ions approach. Hence, an accurate continuum model of the solvation energy is essential to improve our understanding of ionic osmotic and activity coefficients. In addition, it should be possible to generalize the model to nonaqueous solvents, as a further test and to explore the role of water.62 Entropy. An important task is to calculate the temperature dependence and hence entropy of solvation. There are some difficulties. One is that the cavity radius will vary with temperature, which would make a significant contribution to the entropy. The approach given in ref 13 does provide a method of calculating this contribution. An alternative is to look at how the measured RS parameter varies with temperature, which may correlate with the variation in Rcav with temperature. An additional problem is that the cavity contribution is likely to have a significant entropic component and it is not necessarily true that this is the same for microscopic and macroscopic surfaces. Providing more accurate calculation of this contribution and its temperature dependence is a key goal of future research. It is possible that an explicit solvent treatment is necessary to calculate entropies, as they depend on the number of ways the solvent particles can be arranged. In other words, tracking the positions and orientations of the water molecules may be necessary. Promising work has been carried out on this front.63 Nevertheless, this is a far from simple problem, even with an

IV. CONCLUSION We have presented a model which reproduces the solvation free energies of eleven monovalent spherical ions and four noble gas atoms with good accuracy using a continuum solvent model, with no fitted parameters. This model shows that the discrepancy between anion and cation solvation energy which was previously attributed to solvent orientation effects can instead be explained by the larger dispersion interaction of the much more polarizable halide ions. This same explanation also predicts the unusually large solvation energy of the copper(I) and silver(I) ions, which cannot be explained by differing water 9428



Article

(28) Zuo, C.-S.; Wiest, O.; Wu, Y.-D. J. Phys. Chem. A 2009, 113, 12028−12034. (29) Horinek, D.; Mamatkulov, S. I.; Netz, R. R. J. Chem. Phys. 2009, 130, 124507. (30) Liu, J.; Kelly, C. P.; Goren, A. C.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G.; Zhan, C.-G. J. Chem. Theory Comput. 2010, 6, 1109− 1117. (31) Hunenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; The Royal Society of Chemistry, 2011; p 491. (32) Buckingham, A. D. Discuss. Faraday Soc. 1957, 24, 151−157. (33) Wood, R. H. J. Phys. Chem. 1981, 85, 3944−3949. (34) Basilevsky, M. V.; Parsons, D. F. J. Chem. Phys. 1998, 108, 9107. (35) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726. (36) Huang, D. M.; Geissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704−6709. (37) Chandler, D. Nature 2005, 437, 640−7. (38) Rajamani, S.; Truskett, T. M.; Garde, S. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 9475−80. (39) Levin, Y.; dos Santos, A. P.; Diehl, A. Phys. Rev. Lett. 2009, 103, 257802. (40) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016−2027. (41) Shannon, R. D. Acta Crystallogr. 1976, A32, 751. (42) Zhang, Y. O.; Xu, Z. Am. Mineral. 1995, 80, 670−675. (43) Ohtaki, H.; Radnai, T. Chem. Rev. 1993, 93, 1157−1204. (44) Marcus, Y. Chem. Rev. 2009, 109, 1346−70. (45) Fulton, J. L.; Pfund, D. M.; Wallen, S. L.; Newville, M.; Stern, E. A.; Ma, Y. J. Chem. Phys. 1996, 105, 2161. (46) Irudayam, S. J.; Henchman, R. H. J. Phys.: Condens. Matter 2010, 22, 284108. (47) Filipponi, A.; Bowron, D. T.; Lobban, C.; Finney, J. L. Phys. Rev. Lett. 1997, 79, 1293−1296. (48) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 2003, 118, 4976. (49) Camaioni, D. M.; Schwerdtfeger, C. A. J. Phys. Chem. A 2005, 109, 10795−10797. (50) Marcus, Y. Ion Properties; CRC Press: Jerusalem, 1997. (51) Arslanargin, A.; Beck, T. L. J. Chem. Phys. 2012, 136, 104503− 12. (52) Hunenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; The Royal Society of Chemistry, 2011; p 451. (53) Hunenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; The Royal Society of Chemistry, 2011; p 469. (54) Kathmann, S. M.; Kuo, I.-F. W.; Mundy, C. J.; Schenter, G. K. J. Phys. Chem. B 2011, 115, 4369−77. (55) Harder, E.; Roux, B. J. Chem. Phys. 2008, 129, 234706. (56) Sheng, X. W.; Li, P.; Tang, K. T. J. Chem. Phys. 2009, 130, 174310. (57) Grimme, S.; Antony, J.; Schwabe, T.; Mück-Lichtenfeld, C. Org. Biomol. Chem. 2007, 5, 741−58. (58) dos Santos, A. P.; Diehl, A.; Levin, Y. Langmuir 2010, 26, 10778−10783. (59) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2002, 106, 6361− 6373. (60) Kunz, W. Curr. Opin. Colloid Interface Sci. 2010, 15, 34−39. (61) Boström, M.; Williams, D. R. M.; Ninham, B. W. Langmuir 2001, 17, 4475−4478. (62) Peruzzi, N.; Ninham, B. W.; Lo Nostro, P.; Baglioni, P. J. Phys. Chem. B 2012, 116, 14398−405. (63) Beck, T. L. J. Phys. Chem. B 2011, 115, 9776−81. (64) Ashbaugh, H. S.; Asthagiri, D.; Pratt, L. R.; Rempe, S. B. Biophys. Chem. 2003, 105, 323−338. (65) Vavruch, I. J. Colloid Interface Sci. 1995, 169, 249−250.

orientations. This result underscores the important role dispersion interactions play in ion specific effects, showing that specific water structure and chemical interactions can be neglected. We expect that the method of calculating dispersion interactions applied in this paper will be useful in explaining other properties of electrolyte solutions, such as surface tensions and activity and osmotic coefficients.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the NCI National Facility at the ANU. We would like to thank Stefan Scheel, Stefan Buhmann, Angus Grey Weale, and Pierandrea Lo Nostro for very useful discussions.

■

REFERENCES

(1) Ninham, B. W.; Lo Nostro, P. Molecular Forces and Self Assembly; Cambridge University Press: Cambridge, 2010. (2) Kunz, W.; Henle, J.; Ninham, B. W. Curr. Opin. Colloid Interface Sci. 2004, 9, 19−37. (3) Hunenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; The Royal Society of Chemistry, 2011. (4) Kollman, P. Chem. Rev. 1993, 93, 2395−2417. (5) Ninham, B. W.; Duignan, T. T.; Parsons, D. F. Curr. Opin. Colloid Interface Sci. 2011, 16, 612−617. (6) Born, V. M. Z. Phys. 1919, 1, 45−48. (7) Rashin, A. A.; Honig, B. J. Phys. Chem. 1985, 89, 5588−5593. (8) Robinson, R. A.; Stokes, R. H. Electrolyte Solutions; Butterworth & Co.: Devon, 1959. (9) Jayaram, B.; Fine, R.; Sharp, K.; Honig, B. J. Phys. Chem. 1989, 93, 4320−4327. (10) Sandberg, L.; Edholm, O. The J. Chem. Phys. 2002, 116, 2936. (11) Israelachvili, J. N. Intermolecular and Surface Forces, 3rd ed.; Academic Press: San Diego, 2011, p 98. (12) Babu, C. S.; Lim, C. J. Phys. Chem. A 2001, 105, 5030−5036. (13) Fawcett, W. R. J. Phys. Chem. B 1999, 103, 11181−11185. (14) Schmid, R.; Miaha, A. M.; Sapunovb, V. N. Phys. Chem. Chem. Phys. 2000, 2, 97−102. (15) Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R. J. Phys. Chem. A 1998, 102, 7787−7794. (16) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2006, 110, 16066−81. (17) Netz, R. R.; Horinek, D. Annu. Rev. Phys. Chem. 2012, 63, 401− 418. (18) Marcus, Y. J. Chem. Soc., Faraday Trans. 1991, 87, 2995. (19) Lund, M.; Jagoda-Cwiklik, B.; Woodward, C. E.; Vácha, R.; Jungwirth, P. J. Phys. Chem. Lett. 2010, 1, 300−303. (20) Niu, S.; Tan, M.-L.; Ichiye, T. J. Chem. Phys. 2011, 134, 134501. (21) Marcus, Y. Chem. Rev. 1988, 88, 1475−1498. (22) Boström, M.; Ninham, B. W. J. Phys. Chem. B 2004, 108, 12593−12595. (23) Pierotti, R. A. Chem. Rev. 1976, 76, 717−726. (24) Weijo, V.; Mennucci, B.; Frediani, L. J. Chem. Theory Comput. 2010, 6, 3358−3364. (25) Abraham, H.; Liszi, J. J. Chem. Soc., Faraday Trans. 1: Phys. Chem. Condens. Phases 1978, 74, 1604−1614. (26) Cramer, C. J.; Truhlar, D. G. Chem. Rev. 1999, 99, 2161−2200. (27) Duignan, T. T.; Parsons, D. F.; Ninham, B. W. J. Phys. Chem. B 2013, DOI: 10.1021/jp403595x. 9429


A Continuum Model of Solvation Energies Including Electrostatic

Recommend Documents