Real single ion solvation free energies with quantum mechanical simulation

Abstract

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 pK_a values,¹ and activity/osmotic coefficients² can be determined from ion–ion interaction free energies, whereas surface tensions,³ surface forces,⁴ colloidal/protein stability⁵ and surface potentials⁶ 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⁻¹. 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 qϕ. 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–24Fig. 1 illustrates the different approaches to calculating these quantities.

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.²⁷ 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 phase²⁸ 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 calculations¹⁵ and the unphysically large attraction of large anions to the air–water interface observed for polarizable water models.²⁹ 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 performed^38,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.

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.⁴⁰ 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:here we refer to this quantity as the ‘real’ solvation free energy. Ref. 10 provides a detailed derivation and description of it. U_XS is the ion–water interaction energy and is defined^41,42 as U_XS = U_X,Ns − U_Ns where U_X,Ns is the total energy of the solute and solvent system including the electronic energy of the ion and U_Ns 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⁻¹ 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:

U_XS = U_cav + U_XS − U_cav (2)

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

U_XS = U_cav + U_PC + U_XS − U_PC − U_cav (3)

where U_PC = U_PC,Ns − U_Ns 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:

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 (−E^vac_X) 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.

where p₀ (R_cav) is the probability of finding a cavity of size R_cav in bulk water. We have provided a calculation of this term in ref. 44.

The evaluation of the point charge term 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:

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; is a correction associated with the free energy of placing the neutralized hydrogen nucleus in water (discussed in the ESI†); and so 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 is the difference in energy for a point charge in water versus a real ion in water.¹⁵ 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:

Or its inverse:

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

and where 〈δ[U]²〉 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.

U_BM = A exp^−br (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:

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.¹⁵ It is necessary to write it in the inverse form.

where x₀ (R_cav) 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:

Here we have replaced 〈…〉₀ 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 U^exact_Ns with U^DFT_Ns 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.⁴⁴ 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 halides³² and divalent cations,³⁷ but around alkali cations non-trivial discrepancies between simulation and experimental X-ray scattering and spectroscopy results have been observed.⁴⁷

The total solvation free energy can therefore be written as:

where

To estimate the exact ion binding energy (U^exact_XS) 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⁻¹

Method	Li^+	F^−	LiF
This work	−498 ± 3	−507 ± 3	−1005 ± 4
This work	−501 ± 4	−475 ± 3	−976 ± 5
Experiment^9	−520.1	−454.1	−974.2

---

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 value⁴⁰ (≈−0.5 V compared with ≈ 4 V). This means that many calculations of single ion solvation free energies using classical-MD^25,50–53 are not comparable with experiments as they rely on an inherently arbitrary choice for the zero of the electrostatic potential.¹⁰

Table 2 Calculated values for the different types of solvation free energies in kJ mol⁻¹

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)¹⁰ 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

Fig. 4 Contributions to the ‘real’ solvation free energies for the fluoride and lithium ions in kJ mol⁻¹.

Table 3 Contributions to the ‘real’ solvation free energy for different ions in kJ mol⁻¹

Contribution	Li^+	F^−
	5.3 ± 0.2	13.6 ± 0.2
qϕ D	46.3	−46.3
qϕ C	−28.0	42.5
	−538.6 ± 3	−585.9 ± 2
	25.7 ± 1.4	77.3 ± 1.9
	−9.0 ± 1.4	−7.9 ± 0.7
	−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 and terms are quite small and can be reasonably estimated by . In contrast, fluoride has a larger charging energy than the substantially smaller lithium ion, which is then cancelled by a much larger 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.⁵⁴ 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.⁵⁴

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.⁴⁴ 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-MD¹⁰ 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⁻¹) 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.⁴⁷ 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⁻¹. 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 work^a	‘Real’	DFT-MD	−1075 ± 3^a
This work^a	Intrinsic	DFT-MD	−1122 ± 3^a
This work^a	Intrinsic-2^b	DFT-MD	−1108 ± 3^a
This work^a	Bulk	DFT-MD	−1086 ± 8^a
Hünenberger and Reif^9	‘Real’	Lit. Av.	−1095.0
Hünenberger and Reif^9	Intrinsic	Lit. Av.	−1108.0
Tissandier et al.^55–57	—^c	Cluster exp. (CPA)	−1112.5
Marcus^58	Bulk	TATB	−1064.0
Zhan and Dixon^17	—^c	Cluster theory^d	−1105.8
Asthagiri et al.^13	—^c	Cluster theory (QCT)^e	−1065.2
Pollard and Beck^59	‘Real’	Mix	−1105.4
Pollard and Beck^59	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⁻¹ and −14 kJ mol⁻¹ 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.⁴⁰ The difference between the ‘real’ and intrinsic values reported in this study is therefore much larger than the difference recommended by Hünenberger and Reif⁹ who use ϕ_D = 0.13 V.

This 20 kJ mol⁻¹ disagreement for the ‘real’ solvation free energies is particularly concerning as Table 5.15 of Hünenberger and Reif⁹ shows that the estimates of this quantity in the literature shows relatively small variation with a root-mean-square deviation of ≈15 kJ mol⁻¹ and a standard error of ≈3 kJ mol⁻¹. 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 quantity⁴⁰ 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 model⁶¹ 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.⁶² 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.¹⁰ 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⁻¹ and brings the theory into good agreement with the intrinsic solvation free energy reported by Hünenberger and Reif.⁹ The cluster pair approximation (CPA) is one of the more widely accepted approaches for determining single ion solvation free energies.⁵⁵ It is desirable to know whether the CPA estimate is reliable and what type of solvation free energy it is estimating. Hünenberger and Reif⁹ 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 model⁶⁴ indicate that this effect may be on the order of 4 kJ mol⁻¹. 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 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 functionals⁴⁰ 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 analysis^27,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.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 (http:www.cp2k.org) with the QuickStep module for the DFT calculations.⁶⁹ Shorter range double zeta basis sets optimized for the condensed phase⁷⁰ were used in conjunction with GTH pseudopotentials⁷¹ 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.⁷² The revised Perdew, Burke, and Ernzerhof (revPBE)^73,74 functional with the D3 dispersion correction due to Grimme^75,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:

For lithium R_cav 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⁻¹ for each step, which implies the fluctuations are Gaussian to a reasonable approximation.

ORCA⁷⁸ 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⁻¹ of the 24 water molecule cluster. For the MP2 calculation aug-cc-pVDZ basis sets^79,80 were used for the hydrogen, oxygen and fluoride atoms. The cc-pCVDZ⁸¹ 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⁻¹). Values for every term in the calculation that contributes to the solvation free energies are given in the ESI.†

Acknowledgements

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 MS³ (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.

References

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

---

Footnote

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

---