Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Real single ion solvation free energies with quantum mechanical simulation

Timothy T. Duignan*, Marcel D. Baer, Gregory K. Schenter and Christopher J. Mundy
Physical Science Division, Pacific Northwest National Laboratory, P.O. Box 999, Richland, Washington 99352, USA. E-mail:; Tel: +1 509 3756940

Received 12th May 2017 , Accepted 26th May 2017

First published on 4th July 2017

Single ion solvation free energies are one of the most important properties of electrolyte solutions and yet there is ongoing debate about what these values are. Only the values for neutral ion pairs are known. Here, we use DFT interaction potentials with molecular dynamics simulation (DFT-MD) combined with a modified version of the quasi-chemical theory (QCT) to calculate these energies for the lithium and fluoride ions. A method to correct for the error in the DFT functional is developed and very good agreement with the experimental value for the lithium fluoride pair is obtained. Moreover, this method partitions the energies into physically intuitive terms such as surface potential, cavity and charging energies which are amenable to descriptions with reduced models. Our research suggests that lithium's solvation free energy is dominated by the free energetics of a charged hard sphere, whereas fluoride exhibits significant quantum mechanical behavior that cannot be simply described with a reduced model.

1 Introduction

A grand challenge in theory, simulation and modeling is to accurately predict the interaction free energies between ions and other species in water. These free energies determine the density distributions of ions in equilibrium, which in turn determine a huge range of important properties of electrolyte solution. For example, absolute pKa values,1 and activity/osmotic coefficients2 can be determined from ion–ion interaction free energies, whereas surface tensions,3 surface forces,4 colloidal/protein stability5 and surface potentials6 are directly related to ion–surface interaction free energies.

These free energies are determined by a subtle balance of contributions, but one of the most important is the change in the ion–water interaction energy. For example, as an ion approaches another ion or an interface there is a significant energy cost associated with removing water from the ion's hydration layers. Theoretical models therefore need to be carefully tested to ensure that they are correctly reproducing these ion–water interactions. Ionic solvation free energies, the free energy change associated with transferring an ion from vacuum to water, are the most direct experimental measurement of ion–solvent interactions. This is why molecular dynamics with classical interaction potentials (classical-MD) and continuum solvent models are often parameterized or tested by comparison with measured solvation free energies.7,8 These free energies are also important in their own right as they determine the partitioning of ions between different phases.

Due to the electro-neutrality requirement, single ion solvation free energies are one of the only examples of a solvation free energy that is not directly experimentally accessible. A number of ‘extra thermodynamic assumptions’ have been hypothesized in order to provide a convenient estimate of the single ion solvation free energy. (See ref. 9 or ref. 10 for a discussion of some of these approaches.) Unfortunately, none of these have proven sufficiently compelling for the community to agree on, necessitating the use of theoretical methods to resolve this question. Theory has proven inadequate at this task so far, with estimates varying by more than 50 kJ mol−1. Because of the importance of these energies to physical chemistry, their conclusive determination would be a significant achievement in its own right. Additionally, the ability to reliably and accurately compute free energies of molecules in solution is a central problem of physical chemistry. A methodology capable of doing so would have a broad range of very exciting potential applications.

Another challenge is to partition the ionic solvation free energy into separate, physically meaningful terms, such as cavity formation and electrostatic interaction energies. Coarse-grained models which reproduce these separate contributions would not suffer from problems associated with error cancellation. This partitioning is also useful as it will enable us to identify which terms show a linear response and so can potentially be treated with reduced models. The quasi-chemical theory (QCT)11–15 is useful for this purpose. Ref. 15 in particular applies QCT to perform this partitioning with the AMOEBA water model. One particularly important contribution is associated with moving the ion across the surface potential at the air–water interface. This term is of the form . There are several different definitions of ϕ, which correspond to different definitions of the single ion solvation free energy. The expressions for these quantities are provided in the ESI. A full discussion is beyond the scope of this article but is provided in ref. 10 and the references therein.

One important approach to calculating solvation free energies of ions is the cluster continuum method. This approach combines quantum chemistry calculations on small ion–water clusters in the minimum energy geometry with a continuum solvent model.13,16–21 This approach relies on several approximations, namely anharmonicity is neglected; the contribution of the surface potential is ambiguous; and the effect of the surrounding solvent on the water structure is neglected. The validity of these approximations has been discussed extensively elsewhere.10,20–24 Fig. 1 illustrates the different approaches to calculating these quantities.

image file: c7sc02138k-f1.tif
Fig. 1 Schematic depicting the two different approaches to calculating single ion solvation free energies with quantum mechanics. The cluster continuum model is the most widely used, but it relies on several approximations and has no bulk air–water interface and so it is unclear what the surface potential contribution is. We will show how to use DFT-MD to calculate these energies including the contribution of the surface potential at the distant air–water interface.

Attempts at using classical-MD to address this problem have been made.25,26 There are significant challenges with this approach however. For example, it has recently been shown that AMOEBA relies on substantial cancellation of errors to reproduce ion–water dimer binding energies.27 This undermines the notion that these parameters and functional forms are transferable to the condensed phase. In addition, properties such as ionic polarizability are known to vary significantly from the gas phase to the condensed phase28 compounding this issue. As a result, problems have arisen such as the over-polarization of the chloride anion by a factor of 2 with AMOEBA compared with ab initio calculations15 and the unphysically large attraction of large anions to the air–water interface observed for polarizable water models.29 It remains to be seen whether a new generation of polarizable models can overcome these problems.27,30,31

A number of recent studies have determined that density functional theory interaction potentials combined with molecular dynamics simulation (DFT-MD) can provide an accurate description of the water structure around simple ions.32–37 Given the accuracy achieved in determining the local water structure around an ion, it is surprising that very few attempts to calculate solvation free energies with DFT-MD have been performed38,39 particularly given the importance of free energies in determining a range of experimentally relevant properties. It therefore remains to be seen whether this accurate structural description translates into accurate solvation and interaction free energies. In recent years significant advances in the protocols necessary to apply DFT-MD have brought us much closer to resolving this question. Herein, we establish the simulation protocols necessary to calculate single ion solvation free energies and apply these to the lithium and fluoride ions. To this end we use a modified version of QCT that goes beyond the harmonic approximation and includes the important fluctuations beyond the first hydration shell in determining accurate solvation free energies. We proceed by first calculating the solvation free energy of creating a cavity in revPBE-D3 water. We then calculate the free energy of turning a charge on in that cavity using thermodynamic integration. The free energy of replacing this charged hard sphere with a full quantum mechanically treated ion is then estimated using a free energy perturbation method, the free energy of relaxing the hard sphere repulsion and corrections associated with the use of periodic boundary conditions and a small system size are also included. Finally we estimate a correction associated with the error in the revPBE-D3 functional. These contributions are depicted in Fig. 2. This allows us to arrive at real single ion solvation free energies that compare well with experiment. This methodology has the added advantage that it partitions the solvation free energies up into physically intuitive terms that can be mapped onto reduced theories for solvation. Our results suggest that lithium's solvation free energy is dominated by the free energetics of a charged hard sphere whereas fluoride exhibits behavior that requires a quantum mechanical description.

image file: c7sc02138k-f2.tif
Fig. 2 Schematic of the partitioning of the single ion solvation free energies used here. The contributions are the cavity formation, point charge, quantum mechanical and hard sphere relaxation terms.

Our research highlights the importance of using DFT-MD to provide estimates for both the dipolar surface potential due to the presence of a distant air–water interface and the Bethe potential.40 These surface potentials are essential for comparing our predictions with other published theoretical and experimental values in the literature. Ref. 40 and ref. 10 provide a comprehensive discussion of these surface potential.

2 Theory

The goal is to calculate the excess chemical potential of an ion (X) in water at infinite dilution:
image file: c7sc02138k-t1.tif(1)
here we refer to this quantity as the ‘real’ solvation free energy. Ref. 10 provides a detailed derivation and description of it. UXS is the ion–water interaction energy and is defined41,42 as UXS = UX,NsUNs where UX,Ns is the total energy of the solute and solvent system including the electronic energy of the ion and UNs gives the total energy of a given water structure with only the water molecules present. The asterisk indicates that the ‘point to point’ or ‘Ben-Naim’ standard state convention is used. These values differ by −7.95 kJ mol−1 from the 1 atm to 1 M standard state often used.

Following the application of QCT in ref. 15 it is useful to partition the interaction energy of an ion in solution into a hard sphere repulsion, which creates a cavity for the ion to occupy, which is then relaxed after the ion is solvated:

UXS = Ucav + UXSUcav (2)

Ucav is a hard sphere repulsion term, which pushes only on the oxygen atoms. We use this as it allows for a simple determination of the cavity formation energy.

However, instead of placing the real ion in the cavity in a one step process as is done in ref. 15, we break the process up into smaller steps. This is because in contrast to ref. 15 the placement of the ion with a DFT-MD appears to be characterized by non-Gaussian fluctuations. This implies that the free energy cannot reliably be estimated using only equilibrium simulations with the ion present and not present. Instead we must break the process down into smaller steps that can be shown to be approximately Gaussian. Breaking the process up into smaller steps has the added advantage that we can identify the contributions that exhibit linear response behavior as was done in previous studies.10,43 For these reasons we add an additional term to the partitioning which amounts to placing a point charge in the center of the hard cavity that is gradually turned on and then swapped out for the real ion. Because this charging can be performed incrementally, the steps can be made small enough so that the assumption of Gaussian fluctuations is accurate.

UXS = Ucav + UPC + UXSUPCUcav (3)
where UPC = UPC,NsUNs is the energy change on inserting a point charge into a water structure. This partitioning is depicted in Fig. 2.

We can then write the free energy of solvation as:

image file: c7sc02138k-t2.tif(4)

image file: c7sc02138k-t3.tif gives the quantum mechanical contributions to the solvation free energy, i.e., the chemical, dispersion and exchange contributions. It accounts for the difference between the real quantum mechanically treated ion and a charged hard sphere. The electronic vacuum energy (−EvacX) is included in this QM term.

We can estimate the cavity formation energy directly from simulation for cavities up to 3–4 Å by observing the probability of cavity formation at equilibrium.

image file: c7sc02138k-t4.tif(5)
where p0 (Rcav) is the probability of finding a cavity of size Rcav in bulk water. We have provided a calculation of this term in ref. 44.

The evaluation of the point charge term image file: c7sc02138k-t5.tif was carried out in ref. 10 where an extensive discussion of the complexities associated with the correct treatment of the surface potential terms was provided. For the purposes of this paper we calculate the Ewald solvation free energies and then make the appropriate corrections to determine estimates for the intrinsic, bulk and ‘real’ solvation free energies. The definitions of these quantities are provided in the ESI and details on how to calculate them are also provided in ref. 10.

The point charge term can be broken into three separate contributions:

image file: c7sc02138k-t6.tif(6)
where ϕD is the potential created by the dipolar orientation of water molecules at a distant air–water interface; ϕC is the potential created by the orientation of the water molecules surrounding the neutral cavity; image file: c7sc02138k-t7.tif is a correction associated with the free energy of placing the neutralized hydrogen nucleus in water (discussed in the ESI); and so image file: c7sc02138k-t8.tif is the free energy associated with the response of the water to the charging of the ion. To model the point charge a hydrogen nucleus with no basis functions and with a scaled charge is used.

This quantum mechanical term image file: c7sc02138k-t9.tif is the difference in energy for a point charge in water versus a real ion in water.15 The complex electrostatic corrections will mostly cancel as they only depend on the charge and we can therefore simply take the difference in total energy when the point charge is replaced by the real ion. This is given by:

image file: c7sc02138k-t10.tif(7)

Or its inverse:

image file: c7sc02138k-t11.tif(8)

We can expand the averages out with a cumulant expansion up to second order by assuming Gaussian fluctuations and performing the integral analytically.

image file: c7sc02138k-t12.tif(9)
image file: c7sc02138k-t13.tif(10)
where 〈δ[U]2〉 simply indicates the standard deviation squared. We can use both of these expressions and take the average to get a best estimate of this term.

There is one complication, which is that the Bethe potential of the cell (trace of the quadrupole moment) is not precisely the same with the real ion present versus the point charge present. It is therefore necessary to include a small correction associated with the change in the Bethe potential given by qΔϕB when calculating the ‘real’ solvation free energies (see ref. 10 and the ESI for details). We include this correction in the charging energy term.

This method seems to work well for the lithium cation without modification. This suggests that a charged hard sphere is a good model for a lithium ion. For fluoride however, this is not the case. The anion has a large diffuse electron cloud that pushes weakly on the water molecules over a larger range so a hard sphere repulsion is a very poor model for it and so this step is a non-linear/non-Gaussian process. In order to make the charged hard sphere similar to the real ion we use a Born–Mayer type repulsion that acts on the oxygen atoms.

UBM = A[thin space (1/6-em)]expbr (11)

Here r is the ion to oxygen distance and A and b are ion specific parameters. The final solvation free energy should not depend on the choice of these parameters. This process can be performed by rewriting the QM term as:

image file: c7sc02138k-t14.tif(12)

We can then break the first term up into smaller increments, gradually turning on the Born–Mayer repulsion potential so that the Gaussian approximation is accurate. The direct and inverse estimates for all of the contributions are given in the ESI, showing that the Gaussian approximation is reasonable.

The last term in eqn (4) is just the energy of relaxing the hard sphere repulsion. If the hard sphere wall is put just inside the peak in the ion–oxygen radial distribution function, then this term is quite small.15 It is necessary to write it in the inverse form.

image file: c7sc02138k-t15.tif(13)
where x0 (Rcav) is the probability of there being no oxygen atoms inside the hard sphere radius around the ion when the sampling is performed with the real ion–water interactions.

Finally we account for any errors associated with the DFT functional we are using. We can do this by writing the free energy at the exact level as:

image file: c7sc02138k-t16.tif(14)

Here we have replaced 〈…〉0 with 〈…〉UNs to indicate that the sampling is performed with the solvent–solvent interactions turned on. Currently computational limitations mean that the simulation must be performed with DFT level interactions, which means we must replace UexactNs with UDFTNs in the sampling. There is substantial evidence that, although it benefits from cancellation of errors,45,46 revPBE-D3 does a good job describing the structure of pure water, which indicates that this is a reasonable assumption.44 The second term then becomes the solvation free energy determined with DFT plus the DFT vacuum energy. To estimate the first term we take advantage of the same approximation and use structures extracted from the DFT simulation. Note however that the exact expression uses DFT sampling for the ion–water interaction energy. Hence, we do not need to assume that the ion–water interaction energy is described perfectly by the DFT level of theory as any error in this energy will be corrected for assuming it has Gaussian fluctuations. This is an important point as the ability of revPBE-D3 to reproduce bulk water structure has been well tested.44–46 Its ability to reproduce ion–water interactions however, is much less certain. There is strong evidence that GGA functionals accurately describes the water structure around halides32 and divalent cations,37 but around alkali cations non-trivial discrepancies between simulation and experimental X-ray scattering and spectroscopy results have been observed.47

The total solvation free energy can therefore be written as:

image file: c7sc02138k-t17.tif(15)
image file: c7sc02138k-t18.tif(16)

To estimate the exact ion binding energy (UexactXS) we use the MP2 level of theory (the details are discussed below). As we currently lack the capability to calculate the binding energy with MP2 for the full 96 water molecule system in periodic boundary conditions we extract approximately forty ion–water clusters from the simulation and compute the difference in ion–water binding energy for both methods with non-periodic boundary conditions. This term converges reasonably well as the cluster size increases indicating that distant water molecules only interact electrostatically.

3 Results and discussion

3.1 Total single ion solvation free energies

Table 1 and Fig. 3 show that the final theoretical solvation free energy for the lithium fluoride pair is calculated with chemical accuracy and agrees with the experimental value within the statistical uncertainty in the calculation.
Table 1 Values for the ‘real’ solvation free energies. The experimental values are taken from ref. 9. The division of the experimental free energy of the lithium fluoride pair into separate contributions is uncertain due to the difficulty of determining this split experimentally. All energies are given in units of kJ mol−1
Method Li+ F LiF
This work image file: c7sc02138k-t25.tif −498 ± 3 −507 ± 3 −1005 ± 4
This work image file: c7sc02138k-t26.tif −501 ± 4 −475 ± 3 −976 ± 5
Experiment9 −520.1 −454.1 −974.2

image file: c7sc02138k-f3.tif
Fig. 3 Values for the ‘real’ solvation free energies. The spread in the experimental estimates is indicated with the double sided arrow. The statistical uncertainty in the theoretical calculation is much smaller than the spread of experimental values, highlighting why theory is useful for resolving this problem.

These simulations were performed under bulk periodic boundary conditions using Ewald summation. Under these conditions the zero of the electrostatic potential is set so that the average potential over the cell is zero. Thus, the raw solvation free energies are not referenced to the potential in vacuum and they neglect the surface potential created by the real air–water interface.10,40 We refer to these values as the Ewald solvation free energies and Table 2 shows that these values when computed with quantum mechanics are implausible. It is well established that Ewald based free energies are not an experimentally measurable quantity due to the fact that they include a contribution from the large Bethe potential of water which has been extensively discussed.9,10,40,48 Ewald solvation free energies must be carefully corrected to account for role of the Bethe potential as well as finite cell size effects.10,49 These corrections are used to determine the ‘real’, intrinsic, and bulk solvation free energies, as defined in ref. 10 and the ESI. These values are much more in line with experimental estimates than the Ewald values. Unfortunately, these corrections are rarely made in the context of classical-MD. This is likely due to the fact that the Bethe potential calculated with classical-MD is normally much smaller than the quantum mechanical value40 (≈−0.5 V compared with ≈ 4 V). This means that many calculations of single ion solvation free energies using classical-MD25,50–53 are not comparable with experiments as they rely on an inherently arbitrary choice for the zero of the electrostatic potential.10

Table 2 Calculated values for the different types of solvation free energies image file: c7sc02138k-t27.tif in kJ mol−1
Ion ‘Real’ Intrinsic Bulk Ewald
Li+ −501.4 −547.7 −519.7 −873.7
F −474.9 −428.6 −471.0 −91.2

As stated above, the methods employed herein afford a detailed partitioning of the solvation free energy, allowing us to connect with reduced models of solvation. Fig. 4 and Table 3 give the contributions to the single ion solvation free energy for lithium and fluoride. We can see that the free energy is dominated by the charging energy, as is to be expected from a simple Born model. Furthermore, we have added the contributions from the surface dipole potential (ϕD) and the multipolar cavity potential (ϕC)10 that have been discussed in detail in previous publications.10,40 ϕD and ϕC have been demonstrated to exhibit a large dependence on the form of molecular interaction and the corresponding local solvation structure around the ion. Moreover, these electrostatic potentials play a necessary role in defining the important contributions to the solvation free energy. For the case of DFT-MD, ϕC and ϕD largely cancel for both ions resulting in a small net potential.10,40

image file: c7sc02138k-f4.tif
Fig. 4 Contributions to the ‘real’ solvation free energies for the fluoride and lithium ions in kJ mol−1.
Table 3 Contributions to the ‘real’ solvation free energy image file: c7sc02138k-t28.tif for different ions in kJ mol−1
Contribution Li+ F
image file: c7sc02138k-t29.tif 5.3 ± 0.2 13.6 ± 0.2
D 46.3 −46.3
C −28.0 42.5
image file: c7sc02138k-t30.tif −538.6 ± 3 −585.9 ± 2
image file: c7sc02138k-t31.tif 25.7 ± 1.4 77.3 ± 1.9
image file: c7sc02138k-t32.tif −9.0 ± 1.4 −7.9 ± 0.7
image file: c7sc02138k-t33.tif −3.1 ± 1.5 31.8 ± 0.7

An important finding of our research that can be gleaned from examining Table 3 strongly suggests that lithium resembles a simple charged hard sphere, i.e., lithium's image file: c7sc02138k-t19.tif and image file: c7sc02138k-t20.tif terms are quite small and can be reasonably estimated by image file: c7sc02138k-t21.tif. In contrast, fluoride has a larger charging energy than the substantially smaller lithium ion, which is then cancelled by a much larger image file: c7sc02138k-t22.tif term. This is not unexpected, fluoride is known to have a significantly larger exchange energy than similarly sized cations due to the diffuse nature of the wave function which overlaps substantially with the water molecules.54 This cancellation effect would be even larger if the QM term was divided into dispersion and exchange terms as these substantially cancel each other as well.54

Although the exact partitioning of the ion solvation free energy used here has not been applied in the case of classical-MD, it is possible to arrive at some general conclusions based on previous studies. First, the cavity energy is fairly similar with both classical-MD and DFT-MD.44 The relaxation energy could be similar, assuming the classical-MD properly reproduces the structure of water around the ion. The charge hydration asymmetry is much larger with the DFT-MD10 however and in order to compensate for this classical-MD will necessarily underestimate the charge asymmetry in the quantum mechanical term. In particular the large exchange repulsion for the fluoride ion is almost certainly not properly captured by the simple fitted Lennard-Jones potential often used with these models.

The correction associated with the DFT functional used is non-trivial but within the expected accuracy of dispersion corrected GGAs. The small size of the correction term for lithium is slightly misleading as it implies that the revPBE-D3 functional is very accurate for the lithium ion. This is not strictly the case. If we look at the correction for small lithium water clusters of the size of eight water molecules the correction is much larger, (≈−17 kJ mol−1) but for the larger clusters this correction is much smaller indicating that there is a significant cancellation of errors between the ion–water first shell and second shell interactions. This suggests first solvation shell water molecules are too weakly bound whereas the more distant ones bind too strongly. This has been observed in other contexts, namely DFT functionals cannot precisely reproduce the X-ray determined peak position in the sodium–oxygen and potassium–oxygen radial distribution functions.47 Practically speaking, we observe that the cluster correction for lithium converges slower as a function of cluster size than for fluoride. This indicates that ion-specific interactions with the second hydration layer can be important and difficult to reproduce. This highlights a potential problem with classical-MD simulations that are often fitted to reproduce only small ion–water cluster energies and only consider first hydration layer water molecules.25,30,31 The importance of second hydration layer effects has already been established on the basis of cluster-continuum calculations.18,20 These results also highlight the importance of accurately modeling the full condensed phase environment and its fluctuations in order to obtain good estimates of solvation free energies.

To aid comparison with other studies we can use the experimentally well accepted difference in the solvation free energy between the lithium and hydrogen ions, given in ref. 9, in order to arrive at an estimate of the proton solvation free energy. The free energy of the proton is often used as a standard to compare different approaches and models of solvation free energies. Table 4 provides values for the proton solvation free energy using the definitions of the ‘real’, intrinsic and bulk solvation free energies provided in ref. 10 and the ESI.

Table 4 Estimates of the proton solvation free energy (μ*(H+)) in kJ mol−1. A few relevant examples from the literature are also provided for comparison. See Tables 5.15 and 5.19 of ref. 9 for a more complete list. Note that the ‘point to point’ or ‘Ben-Naim’ standard state convention is used
Source Type Method μ*(H+)
a Error is estimated from statistical error in simulation.b Estimated using the center of nuclear charge as the molecular center.c It is unclear how the cluster based values map onto the definitions provided here.d Ref. 20 provides a discussion of cluster-continuum theory methodology generally.e Ref. 60 provides a discussion of cluster-continuum QCT calculations.
This worka ‘Real’ DFT-MD −1075 ± 3a
This worka Intrinsic DFT-MD −1122 ± 3a
This worka Intrinsic-2b DFT-MD −1108 ± 3a
This worka Bulk DFT-MD −1086 ± 8a
Hünenberger and Reif9 ‘Real’ Lit. Av. −1095.0
Hünenberger and Reif9 Intrinsic Lit. Av. −1108.0
Tissandier et al.55–57 c Cluster exp. (CPA) −1112.5
Marcus58 Bulk TATB −1064.0
Zhan and Dixon17 c Cluster theoryd −1105.8
Asthagiri et al.13 c Cluster theory (QCT)e −1065.2
Pollard and Beck59 ‘Real’ Mix −1105.4
Pollard and Beck59 Bulk Mix −1066.8

This table shows our estimate of the ‘real’ and intrinsic solvation free energy differ from the average experimental estimates by 20 kJ mol−1 and −14 kJ mol−1 respectively.

The correction from the intrinsic to the ‘real’ solvation free energy is determined by the dipolar surface potential, ϕD. In this study, we use the dipolar surface potential of 0.48 V calculated with DFT-MD.40 The difference between the ‘real’ and intrinsic values reported in this study is therefore much larger than the difference recommended by Hünenberger and Reif9 who use ϕD = 0.13 V.

This 20 kJ mol−1 disagreement for the ‘real’ solvation free energies is particularly concerning as Table 5.15 of Hünenberger and Reif9 shows that the estimates of this quantity in the literature shows relatively small variation with a root-mean-square deviation of ≈15 kJ mol−1 and a standard error of ≈3 kJ mol−1. The reason for this disagreement could be associated with the calculation of the air–water surface potential. In order to correctly estimate ‘real’ solvation free energies it is necessary to know the dipolar surface potential of the air–water interface. The revPBE-D3 and BLYP-D2 functionals give a value of 0.48 V for this quantity40 if the oxygen atom is chosen to be the center of the water molecule. This value does not show any significant basis set or system size dependency, but there could be quantitative errors in this quantity associated with the use of generalized gradient approximation functionals. The MB-Pol water model61 has a dipolar surface potential of ≈0.3 V, using the oxygen atom as the molecular centre. This model has been shown to agree well with SFG measurements of the air–water interface.62 These measurements are very sensitive to the orientational structure of water molecules at the interface and so this provides some indication that the real dipolar surface potential is ≈0.2 V smaller than the revPBE-D3 value, which would explain this discrepancy.

The intrinsic solvation free energies do not depend on the properties of the air–water interface and so the discrepancy of these values with ref. 9 cannot be explained by an incorrect surface potential. The intrinsic solvation free energies do however depend on the inherently arbitrary choice of the oxygen atom as the center of the water molecule.10 An alternative choice for the origin of the water molecule will result in a different value for the intrinsic solvation free energy and so any comparison of this quantity with experiment must be treated with caution. For example, we can make a potentially more reasonable choice for the center of the water molecule, namely the center of nuclear charge. This definition increases the Bethe potential by 0.14 V and reduces ϕD by 0.14 V. This alters the intrinsic solvation free energy (Intrinsic-2) by 14 kJmol−1 and brings the theory into good agreement with the intrinsic solvation free energy reported by Hünenberger and Reif.9 The cluster pair approximation (CPA) is one of the more widely accepted approaches for determining single ion solvation free energies.55 It is desirable to know whether the CPA estimate is reliable and what type of solvation free energy it is estimating. Hünenberger and Reif9 argue that their intrinsic solvation free energy is consistent with the CPA. The agreement of out Intrinsic-2 value with this value therefore suggests that the solvation free energies determined with the CPA are equivalent to the solvation free energies assuming that the water molecules at a distant air–water interface are isotropically oriented about the center of nuclear charge, namely ϕD = 0.

Additional evidence for this interpretation of the CPA can be inferred from ref. 63. This paper calculated the free energies of forming small ion water clusters with the SPC/E water model and then combined these energies with the CPA to estimate single ion solvation free energies. These calculations showed that the surface potential that was consistent with the SPC/E based CPA estimate of the solvation free energy was 0.16 V less than the actual surface potential of SPC/E water. This difference is similar to the dipolar surface potential of the SPC/E water, which is +0.15 V if the center of nuclear charge is used to determine the water molecules origin.

The estimate for the proton solvation free energy provided here does not depend on any fitting to experiment or adjustable parameters (other then the single parameter used in the development of the revPBE functional). More importantly, it does not rely on the harmonic approximation. The methods used in this study include complex electron correlation effects and do not require any unjustified assumptions about the structure of water around ions such as hydration numbers, as these are self-consistently determined by DFT-MD. The most significant approximation is the use of revPBE-D3 for the contribution of water–water interactions to the structure of water around the ions.

3.2 Uncertainty and future work

These error bars should be considered approximate as they are mainly determined from the differences between the direct and inverse forms of the PDT formula. Relying on direct thermodynamic integration combined with longer trajectories would be required for a precise determination of the error bars. The inverse and direct estimates agree quite closely for all of the terms and these error bars can be combined, assuming they are independent, to determine the error in the total solvation free energies. The resulting uncertainty indicates that this estimate is close to ‘chemical accuracy’.

This does not account for the unquantifiable error that arises from assuming that revPBE-D3 structures are accurate. The neglect of quantum nuclear effects is an issue that should also be addressed in future. Path integral simulations with a classical water model64 indicate that this effect may be on the order of 4 kJ mol−1. This is similar in size to the uncertainty in the estimates reported here. Surprisingly this correction is positive for lithium and negative for fluoride, resulting in a much smaller correction for the salt value.

Improving the estimate of the image file: c7sc02138k-t23.tif correction will be important to test the values determined here. In particular, the CCSD(T) level of theory should be combined with larger basis sets and larger cluster sizes. Additionally, as discussed above the calculation of the surface potential of the air–water interface relies on GGA functionals40 and it is not possible to easily determine the error associated with using DFT for this term, as it depends on the water structure. An improved level of theory could therefore lead to a different value, which would change the cation–anion split reported here. This would not change the experimental agreement of the salt values however as this term makes compensating contributions for the cation and the anion. Performing sampling at the RI-MP2 or RI-RPA sampling level is therefore an important goal.65,66

In future the method should be generalized to other ions such as water's self ions, potassium, sodium, cesium, iodide, divalent ions, tetra-phenyl arsonium and tetra-phenyl borate. An energy decomposition analysis27,67,68 should be used to partition the quantum mechanical energy into dispersion, exchange, induction etc. The MP2 correction for ion–ion and ion–surface PMFs should be estimated. There are also complexities associated with treating an electrolyte solution in the limit of infinite dilution, which are discussed in ref. 10. Finally, coarse grained models should be fitted to reproduce the contributions so that more complex systems can be modeled cheaply.

4 Conclusion

We have used DFT-MD to calculate the ‘real’ solvation free energy of the lithium and fluoride ions including a correction that accounts for the error in the ion–water interaction calculated with DFT. The resulting salt values show excellent experimental agreement and the intrinsic single ion solvation free energies agree well with experimental values based on the CPA, provided that the center of nuclear charge is chosen to be the molecular origin for the water molecule. This calculation moves beyond older approaches that rely on the harmonic approximation and only explicitly consider interactions with the first solvation layer.

This work has important implications for simple models of electrolyte solutions that we believe should be parametrized to reproduce the values calculated herein. We have shown that the lithium ion is reasonably well approximated as a charged hard sphere because the corrections associated with the quantum mechanical nature of the ion are relatively small. In contrast, the fluoride anion has a large quantum mechanical correction that compensates for the large charging contribution. By using a simple, well defined correction to the DFT-MD single ion solvation free energies based in MP2, our research suggests the exquisite sensitivity to ion-specific interactions with water molecules in the second hydration layer that are not properly described with gradient corrected functionals such as revPBE-D3.

5 Calculation details

The system contained 96 water molecules with the ion located at the center of a 14.33 Å3 supercell. Details for the charging free energy are given in ref. 10. Born–Oppenheimer NVT simulations (at 300 K) were performed under PBC using the CP2K simulation suite ( with the QuickStep module for the DFT calculations.69 Shorter range double zeta basis sets optimized for the condensed phase70 were used in conjunction with GTH pseudopotentials71 and a 400 Ry cutoff for the auxiliary plane wave basis and a 0.5 fs time step. A Nosé–Hoover thermostat was attached to every degree of freedom to ensure equilibration.72 The revised Perdew, Burke, and Ernzerhof (revPBE)73,74 functional with the D3 dispersion correction due to Grimme75,76 was used. The energies were accumulated for ≈12 ps after ≈2 ps of equilibration for each simulation. The details of the Bethe potential calculation are provided in ref. 10.

For the hard sphere repulsion we use a potential of the form:

image file: c7sc02138k-t24.tif(17)

For lithium Rcav was set to 2 Å and for fluoride it was set to 2.6 Å. The cavity formation energy was calculated by rasterizing the cell for large revPBE-D3 slab calculations (details given in ref. 44). The relaxation energies were determined from the cumulative radial distribution functions from the ion to the closest oxygen atom (see ESI for details). For the parameters for the Born–Mayer potential we use the value for b determined in ref. 77 and we choose several different A values so that the inverse and direct estimates agree to within 2 kJ mol−1 for each step, which implies the fluctuations are Gaussian to a reasonable approximation.

ORCA78 was used to calculate the cluster correction to the revPBE-D3 functional at the MP2 level of theory. For the revPBE-D3 calculation CP2K was used with the periodicity none option and a larger cell size to remove any box size dependence. Otherwise the same parameters, basis sets etc. as the simulation were used. Clusters of 24 water molecules were used in the cluster correction calculation with 3 frames per picosecond. These calculations showed good convergence with the 32 water molecule cluster being within 1 kJ mol−1 of the 24 water molecule cluster. For the MP2 calculation aug-cc-pVDZ basis sets79,80 were used for the hydrogen, oxygen and fluoride atoms. The cc-pCVDZ81 basis set was used for the lithium ion. The 1s orbitals on the fluoride and oxygen atoms were frozen for the MP2 level calculations. There is some error associated with the basis set size and the MP2 level of theory that was used for the cluster calculations. To correct for this, the binding energies of ions to the nearest four water molecules were calculated with the CCSD(T) level of theory and with MP2 using quadruple zeta basis sets and the counter poise correction. The average differences compared to the MP2 double zeta level calculations were used to estimate corrections for these two issues. These corrections were relatively small (≈6 kJ mol−1). Values for every term in the calculation that contributes to the solvation free energies are given in the ESI.


We would like to thank David Dixon, Thomas Beck, Shawn Kathmann, Liem Dang, Philippe Hünenberger, Sotiris Xantheas, Richard Remsing and John Weeks for helpful discussions. Computing resources were generously allocated by PNNL's Institutional Computing program. This research also used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. TTD, GKS and CJM were supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences. MDB was supported by MS3 (Materials Synthesis and Simulation Across Scales) Initiative, a Laboratory Directed Research and Development Program at Pacific Northwest National Laboratory (PNNL). PNNL is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy.


  1. K. S. Alongi and G. C. Shields,Theoretical Calculations of Acid Dissociation Constants: A Review Article, Elsevier Masson SAS, 2010, vol. 6, ch. 8, pp. 113–138 Search PubMed.
  2. T. T. Duignan, M. D. Baer and C. J. Mundy, Curr. Opin. Colloid Interface Sci., 2016, 23, 58–65 CrossRef CAS.
  3. T. T. Duignan, D. F. Parsons and B. W. Ninham, J. Phys. Chem. B, 2014, 118, 8700–8710 CrossRef CAS PubMed.
  4. D. F. Parsons, M. Boström, T. J. Maceina, A. Salis and B. W. Ninham, Langmuir, 2010, 26, 3323–3328 CrossRef CAS PubMed.
  5. Y. Zhang and P. S. Cremer, Curr. Opin. Chem. Biol., 2006, 10, 658–663 CrossRef CAS PubMed.
  6. R. I. Slavchov, J. K. Novev, T. V. Peshkova and N. A. Grozev, J. Colloid Interface Sci., 2013, 403, 113–126 CrossRef CAS PubMed.
  7. D. Horinek, S. I. Mamatkulov and R. R. Netz, J. Chem. Phys., 2009, 130, 124507 CrossRef PubMed.
  8. T. T. Duignan, D. F. Parsons and B. W. Ninham, J. Phys. Chem. B, 2013, 117, 9421–9429 CrossRef CAS PubMed.
  9. P. Hünenberger and M. Reif, Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities, The Royal Society of Chemistry, 2011 Search PubMed.
  10. T. T. Duignan, M. D. Baer, G. K. Schenter and C. J. Mundy, submitted for publication, 2017, arXiv:1702.05203.
  11. G. Hummer, L. R. Pratt and A. E. García, J. Phys. Chem. A, 1998, 102, 7885–7895 CrossRef CAS.
  12. S. B. Rempe, L. R. Pratt, G. Hummer, J. D. Kress, R. L. Martin and A. Redondo, J. Am. Chem. Soc., 2000, 112, 966–967 CrossRef.
  13. D. Asthagiri, L. R. Pratt and H. S. Ashbaugh, J. Chem. Phys., 2003, 119, 2702–2708 CrossRef CAS.
  14. T. L. Beck, M. E. Paulaitis and L. R. Pratt, The Potential Distribution Theorem and Models of Molecular Solutions, Cambridge University Press, 2006 Search PubMed.
  15. D. M. Rogers and T. L. Beck, J. Chem. Phys., 2010, 132, 014505 CrossRef PubMed.
  16. G. J. Tawa, I. A. Topol, S. K. Burt, R. A. Caldwell and A. A. Rashin, J. Chem. Phys., 1998, 109, 4852–4863 CrossRef CAS.
  17. C.-G. Zhan and D. A. Dixon, J. Phys. Chem. A, 2001, 105, 11534–11540 CrossRef CAS.
  18. C.-G. Zhan and D. A. Dixon, J. Phys. Chem. A, 2004, 108, 2020–2029 CrossRef CAS.
  19. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093 CrossRef CAS PubMed.
  20. V. S. Bryantsev, M. S. Diallo and W. A. Goddard, J. Phys. Chem. B, 2008, 112, 9709–9719 CrossRef CAS PubMed.
  21. D. Sabo, D. Jiao, S. Varma, L. R. Pratt and S. B. Rempe, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2013, 109, 266–278 RSC.
  22. S. Kathmann, G. Schenter and B. Garrett, J. Phys. Chem. C, 2007, 111, 4977–4983 CAS.
  23. D. M. Rogers and S. B. Rempe, J. Phys. Chem. B, 2011, 115, 9116–9129 CrossRef CAS PubMed.
  24. S. Merchant, P. D. Dixit, K. R. Dean and D. Asthagiri, J. Chem. Phys., 2011, 135, 054505 CrossRef PubMed.
  25. A. Grossfield, P. Ren and J. W. Ponder, J. Am. Chem. Soc., 2003, 125, 15671–15682 CrossRef CAS PubMed.
  26. H. Yu, T. W. Whitfield, E. Harder, G. Lamoureux, I. Vorobyov, V. M. Anisimov, A. D. Mackerell and B. Roux, J. Chem. Theory Comput., 2010, 6, 774–786 CrossRef CAS PubMed.
  27. Y. Mao, O. Demerdash, M. Head-Gordon and T. Head-Gordon, J. Chem. Theory Comput., 2016, 12, 5422–5437 CrossRef CAS PubMed.
  28. A. Serr and R. R. Netz, Int. J. Quantum Chem., 2006, 106, 2960–2974 CrossRef CAS.
  29. M. D. Baer and C. J. Mundy, J. Phys. Chem. Lett., 2011, 2, 1088–1093 CrossRef CAS.
  30. D. J. Arismendi-Arrieta, M. Riera, P. Bajaj, R. Prosmiti and F. Paesani, J. Phys. Chem. B, 2016, 120, 1822–1832 CrossRef CAS PubMed.
  31. M. Riera, A. W. Götz and F. Paesani, Phys. Chem. Chem. Phys., 2016, 18, 30334–30343 RSC.
  32. J. L. Fulton, G. K. Schenter, M. D. Baer, C. J. Mundy, L. X. Dang and M. Balasubramanian, J. Phys. Chem. B, 2010, 114, 12926–12937 CrossRef CAS PubMed.
  33. M. D. Baer, V.-T. Pham, J. L. Fulton, G. K. Schenter, M. Balasubramanian and C. J. Mundy, J. Phys. Chem. Lett., 2011, 2, 2650–2654 CrossRef CAS.
  34. J. L. Fulton, E. J. Bylaska, S. Bogatko, M. Balasubramanian, E. Cauët, G. K. Schenter and J. H. Weare, J. Phys. Chem. Lett., 2012, 3, 2588–2593 CrossRef CAS PubMed.
  35. S. Bogatko, E. Cauët, E. Bylaska, G. Schenter, J. Fulton and J. Weare, Chem.–Eur. J., 2013, 19, 3047–3060 CrossRef CAS PubMed.
  36. M. D. Baer, I.-F. W. Kuo, D. J. Tobias and C. J. Mundy, J. Phys. Chem. B, 2014, 118, 8364–8372 CrossRef CAS PubMed.
  37. M. D. Baer and C. J. Mundy, J. Phys. Chem. B, 2016, 120, 1885–1893 CrossRef CAS PubMed.
  38. K. Leung, S. B. Rempe and O. A. von Lilienfeld, J. Chem. Phys., 2009, 130, 204507 CrossRef PubMed.
  39. V. Weber and D. Asthagiri, J. Chem. Phys., 2010, 133, 141101 CrossRef PubMed.
  40. R. C. Remsing, M. D. Baer, G. K. Schenter, C. J. Mundy and J. D. Weeks, J. Phys. Chem. Lett., 2014, 5, 2767–2774 CrossRef CAS PubMed.
  41. A. Ben-Naim, J. Phys. Chem., 1978, 82, 792–803 CrossRef CAS.
  42. D. Ben-Amotz, F. O. Raineri and G. Stell, J. Phys. Chem. B, 2005, 109, 6866–6878 CrossRef CAS PubMed.
  43. R. C. Remsing and J. D. Weeks, J. Phys. Chem. B, 2016, 120, 6238–6249 CrossRef CAS PubMed.
  44. M. Galib, T. T. Duignan, Y. Misteli, M. D. Baer, G. K. Schenter, J. Hutter and C. J. Mundy, J. Chem. Phys., 2017 DOI:10.1063/1.4986284 , arXiv:1703.06940.
  45. L. Ruiz Pestana, N. Mardirossian, M. Head-Gordon and T. Head-Gordon, Chem. Sci., 2017, 8, 3554–3565 RSC.
  46. O. Marsalek and T. E. Markland, J. Phys. Chem. Lett., 2017, 8, 1545–1551 CrossRef CAS PubMed.
  47. M. Galib, M. D. Baer, L. B. Skinner, C. J. Mundy, T. Huthwelker, G. K. Schenter, N. Govind and J. L. Fulton, J. Chem. Phys., 2017, 146, 084504 CrossRef CAS PubMed.
  48. S. M. Kathmann, I.-F. W. Kuo, C. J. Mundy and G. K. Schenter, J. Phys. Chem. B, 2011, 115, 4369–4677 CrossRef CAS PubMed.
  49. M. A. Kastenholz and P. H. Hünenberger, J. Chem. Phys., 2006, 124, 224501 CrossRef PubMed.
  50. S. Rajamani, T. Ghosh and S. Garde, J. Chem. Phys., 2004, 120, 4457–4466 CrossRef CAS PubMed.
  51. G. Lamoureux and B. Roux, J. Phys. Chem. B, 2006, 110, 3308–3322 CrossRef CAS PubMed.
  52. J. P. Bardhan, P. Jungwirth and L. Makowski, J. Chem. Phys., 2012, 137, 124101 CrossRef PubMed.
  53. F. Sedlmeier and R. R. Netz, J. Chem. Phys., 2013, 138, 115101 CrossRef PubMed.
  54. T. P. Pollard and T. L. Beck, Curr. Opin. Colloid Interface Sci., 2016, 23, 110–118 CrossRef CAS.
  55. M. D. Tissandier, K. A. Cowen, W. Y. Feng, E. Gundlach, M. H. Cohen, A. D. Earhart, J. V. Coe and T. R. Tuttle, J. Phys. Chem. A, 1998, 102, 7787–7794 CrossRef CAS.
  56. D. M. Camaioni and C. A. Schwerdtfeger, J. Phys. Chem. A, 2005, 109, 10795–10797 CrossRef CAS PubMed.
  57. C. P. Kelly, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2006, 110, 16066–16081 CrossRef CAS PubMed.
  58. Y. Marcus, Ion Solvation, Wiley, 1985 Search PubMed.
  59. T. P. Pollard and T. L. Beck, J. Chem. Phys., 2014, 141, 18C512 CrossRef PubMed.
  60. D. Asthagiri, P. D. Dixit, S. Merchant, M. E. Paulaitis, L. R. Pratt, S. B. Rempe and S. Varma, Chem. Phys. Lett., 2010, 485, 1–7 CrossRef CAS PubMed.
  61. G. R. Medders, V. Babin and F. Paesani, J. Chem. Theory Comput., 2014, 10, 2906–2910 CrossRef CAS PubMed.
  62. G. R. Medders and F. Paesani, J. Am. Chem. Soc., 2016, 138, 3912–3919 CrossRef CAS PubMed.
  63. L. Vlcek, A. A. Chialvo and J. M. Simonson, J. Phys. Chem. A, 2013, 117, 11328–11338 CrossRef CAS PubMed.
  64. D. M. Wilkins, D. E. Manolopoulos and L. X. Dang, J. Chem. Phys., 2015, 142, 064509 CrossRef PubMed.
  65. M. Del Ben, J. Hutter and J. VandeVondele, J. Chem. Phys., 2015, 143, 054506 CrossRef PubMed.
  66. M. Del Ben, O. Schütt, T. Wentz, P. Messmer, J. Hutter and J. Vandevondele, Comput. Phys. Commun., 2015, 187, 120–129 CrossRef CAS.
  67. K. U. Lao, R. Schäffe, G. Jansen and J. M. Herbert, J. Chem. Theory Comput., 2015, 11, 2473–2486 CrossRef CAS PubMed.
  68. P. Horn, Y. Mao and M. Head-Gordon, Phys. Chem. Chem. Phys., 2016, 18, 23067–23079 RSC.
  69. J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128 CrossRef CAS.
  70. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105 CrossRef PubMed.
  71. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710 CrossRef CAS.
  72. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643 CrossRef.
  73. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
  74. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890 CrossRef CAS.
  75. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473 CrossRef CAS PubMed.
  76. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104 CrossRef PubMed.
  77. T. T. Duignan, D. F. Parsons and B. W. Ninham, J. Phys. Chem. B, 2013, 117, 9412–9420 CrossRef CAS PubMed.
  78. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78 CrossRef CAS.
  79. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023 CrossRef.
  80. R. A. Kendall, T. H. Dunning Jr and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806 CrossRef CAS.
  81. D. E. Woon and T. H. Dunning Jr, J. Chem. Phys., 1995, 103, 4576–4585 CrossRef.


Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sc02138k

This journal is © The Royal Society of Chemistry 2017