Re-examining the tetraphenyl-arsonium/tetraphenyl-borate (TATB) hypothesis for single-ion solvation free energies

For well over 100 years, researchers have sought a comprehensive theoretical understanding of electrolyte solutions, including the specific ion effects (SIEs) observed throughout chemistry and biology.^1,2 SIEs can arise in diverse and unexpected places.^3–5 A knowledge of the single-ion solvation free energy scale is necessary for understanding, quantifying, and predicting these effects.^6,7 The free energy change upon inserting an ion in a solvent is composed of both the free energy to grow an ion-sized cavity and the subsequent free energy change to turn on the ion charge. The latter free energy includes electrostatic effects from distant interfaces. Early efforts took a purely electrostatic viewpoint (neglecting the cavity formation energy).^8–10 Later theories built on these dominating electrostatic effects to incorporate many-body dispersion forces and the cavity formation term.¹¹ These modern iterations, however, remain incomplete, and they struggle to account for the strong local effects near the ion, where field magnitudes are on the order of 1 V/Å.¹² Under the stress of such strong fields, it is perhaps not surprising to find that simple models do not accurately portray the significant distortion and even exchange of electron density between the ion and nearby solvent molecules.^1,13–20

In order to explain the myriad SIEs, it can be argued that initial helpful goals might be (1) to understand single-ion solvation free energies at the infinite dilution limit and then (2) to develop more accurate theories for ion-ion interactions at finite concentrations. Here we focus on the infinite dilution single-ion solvation free energy limit. There is a long history of attempts to partition easily measured pair solvation free energies into the individual contributions.^6,7,21 Two primary themes in these approaches are first to utilize a simplified molecular model to fit to measured pair thermodynamic data and then extract the single-ion quantities from the model [Marcus,²² Latimer,²³ tetraphenyl-arsonium/tetraphenyl-borate (TATB)²⁴] and second to employ a combination of bulk thermodynamic measurements and cluster formation data for many ion pairs to indirectly obtain the solvation free energy of the proton [cluster pair approximation (CPA)].^25–29 There have been many other approaches taken toward estimates of single-ion thermodynamic quantities, and this literature is thoroughly reviewed in Refs. 6 and 7.

Inserting a single ion into a solvent involves the physical process of crossing a liquid-vapor interface with a possible electrostatic potential shift due to solvent molecular features near the interface.^30–32 It is well known, however, that there is a “chicken-egg” aspect to these potentials. For example, the electromotive force (EMF) of a galvanic cell can be expressed in terms of intrinsic chemical potentials for the relevant ions (see below) with no knowledge of surface (or contact) potentials.³³ Yet the EMF is also simply the sum of all the contact potentials around the circuit. Put another way, the contact potential across an interface between two conducting liquids is determined by the intrinsic chemical potentials of the ions that partition between the bulk phases.³⁴ Nonetheless, it is physically realistic to discuss the separate ionic double layer and solvent contributions in direct calculations of interfacial potentials.

It is clear that the first single-ion strategy above (based on a simplified model and pair thermodynamic data) results in no interfacial potential contribution in terms of the net potential defined below. This is because both in the model and in the thermodynamic data, the involvement of ion pairs leads to a complete cancellation of the net interfacial potential contributions to the free energies. In a previous paper, we referred to these single-ion free energies as “bulk” chemical potentials.³⁵ Alternatively, the CPA approach has been shown to produce the “real” solvation free energy of the proton that includes interfacial potential effects.^27,28 Thus the difference between the bulk and real chemical potentials should provide an estimate of an effective electrochemical interfacial potential (solvent contribution without other interfering ions).²¹ For alternative views on whether the CPA produces the real solvation free energy, see Ref. 7.

In the context of quasichemical theory (QCT), the solvation process for a single ion consists of three steps:^36,37 (1) open an ion-sized (or slightly larger) cavity of radius λ in the solvent, (2) insert the ion into the cavity, (3) remove the cavity constraint, allowing full contact between the solvent and ion. The resulting three free energy terms are in order: packing (cavity formation), long-ranged (outer shell), and chemical (inner shell). The packing and inner shell contributions can be viewed as processes involving reversible work: the packing term is the work to grow a λ-sized cavity in the pure solvent, while the chemical term is minus the work to grow the same-sized cavity around the ion. The sum of all three terms is independent of the cavity size λ. If a continuum dielectric model is used to estimate the long-ranged term, the QCT free energy contains no net interfacial potential contribution.

Any interfacial potential contributions are thus located in the long-ranged term. We can imagine a single ion solvated at the center of a small droplet. Then there are two interfaces that can produce charge inhomogeneities that lead to potential jumps, one at the liquid-vapor interface ϕ_sp and the other one at the boundary between the cavity and the solvent ϕ_lp.³⁵ The sum of these two potentials is the net potential ϕ_np = ϕ_lp + ϕ_sp. This quantity includes significant cancellations of large molecular dipole and quadrupole contributions, with a small remainder due to the spherical geometry (with significant curvature) near the cavity.

In a previous study,³⁷ we located a consistent length scale for monovalent ions at which the packing and the long-range terms (after removing the interfacial potential contribution) cancel. This length is 6.15 Å for the SPC/E water model, and it is not expected to vary substantially with the water model. Also for the SPC/E water model and the above length scale, the local potential ϕ_lp is +9.5 kcal/mol e, while the surface potential ϕ_sp is −13.8 kcal/mol e, producing a net potential ϕ_np of −4.3 kcal/mol e. A plot of the net potential ϕ_np vs. cavity radius is only weakly dependent on the cavity radius at this and larger λ values.³⁸ 

The TATB method for determining single-ion solvation properties is attractive due to its simplicity and reasonable physical basis.²⁴ These large organic ions, decorated with 4 benzene rings tetrahedrally arranged around a central ion, exhibit a preference for lower-dielectric organic solvents.²² Even though the central ions are quite different, and thus differing charge distributions can be expected,³⁹ the aromatic rings likely shield the charge density differences substantially. In support of this notion, bulk proton hydration thermodynamic values derived from entirely different approaches such as the Marcus (Halliwell/Nyburg), Latimer, and QCT methods cited above are within 2 kcal/mol of those obtained by the TATB method. A recent QCT study by Carvalho and Pliego⁴⁰ lends further support to the above discussion; our results below are thus consistent with those recent quantum/continuum calculations. There has also been extensive experimental work examining differences in the physical behavior of the TA cations and the TB anions in the bulk, near interfaces, and proteins.^41–45 A purpose of the present work is to begin to unravel the separate contributions to the solvation from interfacial electrostatic (charge dependent) effects vs. more local specific ion effects.

Previous classical simulation studies determined the free energy change for transforming the TA cation into the TB anion in water modeled with periodic boundary conditions.^39,46–48 For the TIP3P model of water (with a charge distribution similar to the SPC/E model employed here), a significant difference of roughly 20 kcal/mol was observed for the hydration free energy of a simple spherical ion the size of the TA and TB ions, implying a significant deviation from the TATB model assumptions. On the other hand, the computed difference using the TIP5P water model was much smaller. These differences can be linked to quadrupole moment effects due to differing charge distributions.^49–51 Here we find that the main contribution to this deviation (for the SPC/E model) is twice the local potential value discussed above, as expected for a +1/−1 ion pair. In reality, this local potential effect should be subtracted out in comparing free energies since the proper comparison (for the TATB method) is with the bulk (surface-potential-free) solvation free energies. We note that, importantly, the computed deviation without subtracting out the local potential would approach a magnitude of 200 kcal/mol in more realistic quantum mechanical models with a very different intramolecular charge distribution (and thus quadrupole moment trace), compared with the simple classical models.^34,35,52,53 This fact highlights the importance of a proper handling of electrostatic effects.

The paper is organized as follows. We begin by discussing relevant theory and models that aid in analyzing our computational results. The models include thermodynamic perturbation theory (TPT), quasichemical theory (QCT), and a simple continuum model of solvation that serves as a conceptual guide to the results. Next, we present our findings from the application of QCT principles to the solvation of spherical ion pairs with sizes ranging from 3 Å to 5.5 Å in SPC/E water, OPLS/AA (Optimized Potential for Liquid Simulations/All Atom) dimethyl sulfoxide (DMSO), and OPLS/AA 1,2-dichloroethane (1,2-DCE). For clarity and brevity, data related to the non-aqueous solvents are shown in the supplementary material. These data will be discussed alongside that of aqueous systems. Finally, we summarize the conclusions of our computational results in relation to the TATB hypothesis.

In a classical model, the solvation process for a monoatomic ion can be enacted in two steps: insertion of an ion-sized van der Waals (vdW) particle followed by turning on of the ionic charge. Then, if the interaction energy of the particle with the surrounding solvent is ε_X = U(N + X) − U(N) − U(X) [where U(N + X) is the total potential energy of the solvent plus ion system, U(N) is the potential energy of the N solvent molecules, and U(X) is the potential energy of the ion], the free energy of solvation (or excess chemical potential) is

The subscript “0” means sampling with only solvent/solvent interactions. Above ε_X,es is the electrostatic part of the interaction energy. Here we focus on that dominant electrostatic portion of the free energy.

An exact rearrangement of the free energy $μX,esex$ yields mean-field and fluctuation terms,

The electrostatic interaction energy of the charge with the surrounding solvent is simply the product of the charge and the potential at the vdW particle center. The expression δε_X,es is the deviation from the mean.

The cumulant expansion to third-order is then

In Ref. 37, it was found that the third-order (skewness) term provides a nontrivial correction. The first and third order terms can thus contribute to the (charge sign-dependent) effective electrochemical surface potential result, while the second order term is closely related to a continuum Born model that is quadratic in the charge.

As mentioned above (and shown in Fig. 1 of Ref. 21), the net potential at the center of the ion is due to inhomogeneous charge distributions in two regions, one the bulk liquid-vapor interface and the other a spherical region near the ion. In between these domains, the average charge density is zero and thus the bulk water region makes no contribution to the net potential. Landau and Lifshitz³³ referred to the surface potential as the contact potential. This potential can be expressed as a simple dipole-like integral through the plane interface. In ongoing work on large water droplets, we have found that the local potential is dominated by molecular dipolar and quadrupolar contributions, as seen previously.^51,54 The distant interface producing ϕ_sp yields only a molecular dipole contribution at the cavity center since the quadrupole term from a spherical shell dies off as 1/r. For the liquid-vapor interface, the dipole contributions from the cavity region and the bulk interface nearly cancel, leaving the quadrupolar contribution near the ion as the dominant one. This is the reason for the strong model dependence of computed net potentials since the quadrupole can vary widely with the model.^49,50,54

A proper comparison between the TATB assumption and simulation then requires removal of the first and third terms above from the long-ranged part of the free energy. The work-related contributions from the inner-shell and packing terms contain no surface potential contribution. To summarize, the real electrochemical potential is

where $μX,intex=μX,bex+qϕlp$ is the intrinsic excess chemical potential that pertains to simulations in periodic boundary conditions with no explicit interfaces. As discussed above, this situation is “unbalanced” in the sense that the computed results for $μX,intex$ depend heavily on the model. Classical and quantum results for $μX,intex$ will differ by a magnitude of roughly 100 kcal/mol due to a quadrupole trace term that is very large for the realistic charge distribution of quantum matter compared with a smaller magnitude and opposite sign term in classical point charge models.^34,35,52,53 We note that the intrinsic free energy defined here differs from the corresponding free energy discussed in Ref. 7 due to the inclusion of the quadrupole trace.

A direct indication of the role of ϕ_lp in periodic boundary conditions is apparent in Ref. 55, where a QCT model employing quantum calculations for the inner-shell region and a continuum (surface potential free) model for the long-range term differ by roughly −10 kcal/mol from a model that includes a classical periodic boundary molecular dynamics (MD) model for the long-ranged term. This implies a local potential consistent with our previous result of 9.5 kcal/mol e.³⁷ 

It is interesting to ask which quantities are experimentally accessible. Based on arguments in Refs. 7, 27, and 28, it is clear that the real hydration free energy $μXex$ is experimentally accessible. In addition, the surface potential is (at least in principle) accessible via scattering experiments that can produce the average charge density through a liquid-vapor or liquid-liquid interface.^34,53 Thus, the intrinsic free energy $μX,intex$ is experimentally accessible, although it is expected to deviate substantially from other reported single-ion values due to the very large ϕ_lp contribution in realistic quantum matter. The bulk free energy $μX,bex$⁠, on the other hand, is not directly accessible to experiment due to its dependence on the chosen cavity size λ (which in turn affects the net potential ϕ_np value). There are clear indications, however, that, outside the first solvation shell, the net potential stabilizes at a near constant value.³⁸ Indirect experimental evidence hints that both the distributions of ions near interfaces⁴² and the acid-base properties of water itself^56,57 are affected by a mean potential that is around −0.4 V or −10 kcal/mol e for water.

The excess chemical potential of a monatomic ion may be evaluated using the following expression, derived from the potential distribution theorem (PDT):³⁶ 

Here M(λ) is a model potential, taken as a hard-sphere sharp-cutoff potential in the original QCT. For computational convenience, we take M(λ) to be a strongly repulsive but continuous potential, as done in Ref. 37. ΔU signifies sampling with full ion/solvent interactions.

The first (inner-shell, is) term is minus the work required to push the solvent molecules from direct contact with the ion out to the cavity radius λ. The second (packing, pk) term is correspondingly the work to create a cavity of size λ in pure water. The third (long-ranged or outer-shell, os) term is the free energy to place the ion into the cavity with the solvent held beyond the λ radius.

A reasonable goal is to choose λ large enough to observe Gaussian behavior in the long-ranged term. Such Gaussian behavior suggests that the average of two mean-field terms at the q = 0 and q = ±1 limits should be accurate. As observed long ago,⁵⁸ this is problematic for anions in water since there is a significant rearrangement of the nearby water molecules as the charge magnitude changes from 0 toward −1. This was also noted by Shi and Beck³⁷ who found that a midpoint-rule approximation was quite accurate. Thus we use that midpoint approximation in our calculations for the long-ranged term.

As mentioned above, it was also found in Ref. 37 that there is a common length scale (λ = 6.15 Å for the SPC/E water model) for small monovalent ions at which the packing and long-range contributions (with the interfacial potential contribution removed) cancel. This can be pictured as the length scale at which the cavity formation free energy balances the Born free energy. Then the bulk excess chemical potential is exactly minus the work to grow a λ-sized cavity around the ion in water. It is this bulk free energy $μX,bex$ that should be compared with the TATB model (see above). It is interesting that the net potential stabilizes at a consistent value at this length scale,³⁸ and the continuum plots in Fig. 1 (below) bear a remarkable resemblance to the simulation results of Ref. 37 for this and larger cavity sizes. Of course, the total excess free energy is independent of the cavity size λ. To summarize, the net potential is a predominantly mean-field potential that is ion-independent once the cavity size is large enough so that all near-local ion-specific effects are included within that region. Beyond this size, continuum models are more appropriate for the outer-shell term.

Finally, we reproduce the correct formula for the long-ranged term including the local potential involved in the periodic boundary simulations. The outer shell (long ranged) contribution $μosex$ can be split into vdW and electrostatic components via a cumulant expansion,^36,58

The electrostatic term has been expanded to third order. The quadratic term is the Born energy of the ion, including a self-energy correction (involving the constant ξ and the box size L) for the periodic boundary situation, while the first and third terms constitute the local potential contribution when added together. ϕ_lp is computed as described previously as the electrostatic potential at the center of a neutral cavity with the same radius used for the other terms. Convergence of the cumulant expansion can only be obtained with a sizable cavity (as discussed above), implying that higher order terms are appreciable at smaller length scales.

Here we write out a highly simplified continuum model of ion solvation as a guide to rationalize the approach taken below. We neglect dispersion effects,¹⁷ which are undoubtedly important for the total free energy of solvation but likely do not differ substantially between the TA cation and the TB anion. We include only the Born electrostatic contribution and a cavity formation term. The total bulk free energy of solvation in a Born/surface-tension model is then

where q is the ion charge, R is the ion radius, ϵ is the solvent dielectric constant, and γ is the surface tension.

The packing free energy is approximately

while the long-ranged (outer-shell) term is

Since the total is independent of λ, the inner-shell term is then

We choose the ion radius R so that the inner-shell chemical term goes to zero at the same length scale as that observed in our previous MD simulations.

This extremely simple cartoon model cannot be expected to reproduce detailed molecular scale information concerning ion solvation. We note, however, that the resulting plots for each term (as a function of cavity size, Fig. 1) do provide a decent qualitative representation of the form of the curves when compared with molecular simulations for the Cl⁻ ion, with reasonable choices for the parameters (see Fig. 4 in Ref. 37). Also we note that the apparent linearity of the inner-shell term masks its underlying nonlinear origin over the length scales displayed in the figure. Finally, there are two interesting points of cancellation: (1) the first at about 3.6 Å where the inner-shell and packing terms cancel leaving the long-range term as the total free energy and (2) the second at about 6.6 Å where the packing and long-ranged terms cancel leaving the inner-shell as the total free energy. The more accurate molecular dynamics data do not match these length scales exactly, of course, but the simple model does give a hint at the approximate distances.

Based on this model, we see that we can gain insights into expected differences between large ions by focusing on the inner-shell term, since the packing term has no dependence on the ion type, and the long-ranged term also should be very similar between the large TA cation and TB anion models. Thus we choose a large value for λ of 9 Å which includes the first solvation shell beyond the roughly 5.5 Å radius of the TA and TB ions. (Above we noted the consistent value of the local potential for cavities in the 6-9 Å size range.) We then vary the ion radius in the simulations to see how the free energy difference between the TA and TB ions varies with ion size.

All calculations were performed in double-precision with Gromacs 4.6.7⁵⁹ code installed on the Oakley cluster at the Ohio Supercomputer Center.⁶⁰ The OPLS-AA⁶¹ force field was used for DMSO and 1,2-DCE. SPC/E⁶² parameters were used for water. Bonds and angles involving hydrogen atoms were constrained with the LINCS algorithm⁶³ in all simulations. Ion parameters are listed in Table I and range from 3 Å to 5.5 Å in radius, the same size used in the Schurhammer and Wipff studies.^46–48 Dispersion energy parameters were chosen to surround the Schurhammer and Wipff value of 0.1 kcal/mol with the smaller value matching that of I⁻ in the Horinek and Netz force field.⁶⁴ The larger value was taken as one order of magnitude greater. Only spherical ions were modeled in this work, as discussed earlier.

The Good-Hope⁶⁵ and Berthelot⁶⁶ combining rules were used to generate the pair σ_ij and ε_ij parameters, respectively. The ion/cavity/solvent or cavity/solvent systems were equilibrated in the NPT ensemble for 400 ps at 1.0 bar and 298.15 K using the Berendsen barostat⁶⁷ and V-rescale algorithms⁶⁸ for pressure and temperature coupling. Real-space electrostatic and vdW interactions were truncated at 12 Å, with a 1.6 Å grid spacing and 6th-order spline for reciprocal space interactions. We used an Ewald inverse length parameter of η = 5.6/L, where L is the box length after NPT equilibration. Energies were sampled every 20 frames from 2.0 ns of a 2.2 ns production run in the NVT ensemble, also held at 298.15 K with the V-rescale algorithm.⁶⁸ For quantities obtained through thermodynamic integration, 21 unique λ values were used. Each simulation, including at each λ, was performed 3 separate times with a time step of 2.0 fs. Statistical errors were evaluated with the block-averaging method and propagated through numerical integration as well. The solvation environment around an ion with the fully grown cavity in each solvent is illustrated in Fig. 2. The cavity is kept at a 9 Å radius which did not change with ion size. The final cubic box sizes L are close to 40 Å for the water and 1,2-DCE systems and 50 Å for the DMSO system.

The solvents modeled here cover a broad range of experimental dielectric constants: water (80), DMSO (46.7), and 1,2-DCE (10.4). For additional validation of the methods, we computed the dielectric constant for each solvent. We performed 400 ps of equilibration of the pure solvent without cavities as above followed by 10 ns constant volume dynamics at 298.15 K. The first 500 ps of each trajectory was discarded as additional equilibration. The SPC/E dielectric constant was found to be 72 (71), in DMSO it was 45 (42), and in 1,2-DCE it was 12 (13). Values in parentheses were taken from the studies of Reddy and Berkowitz⁶⁹ and Caleman et al.⁷⁰ Thus the results are in reasonable agreement with experiment and previous simulations. We had also calculated the dielectric constant for the OPLS/AA acetonitrile force field, producing a dielectric constant of 20 (about half that of the experimental value).⁷⁰ Thus we did not further examine the acetonitrile case.

We begin with the discussion of the cavity-growth terms in each solvent. The packing term depends only on the solvent properties. The computed values for each solvent are (for λ = 9 Å) water (82.77 ± 0.63 kcal/mol), DMSO (83.50 ± 0.91 kcal/mol), and 1,2-DCE (49.45 ± 0.56 kcal/mol). Our result for water compares favorably with the scaled-particle theory calculation of Shi and Beck.³⁷ 

We now consider the ion specific chemical part in Fig. 3 and Figs. S1 and S6 of the supplementary material. From these figures, we see that the anion is preferentially solvated in water and 1,2-DCE relative to the cation. The opposite is observed in DMSO. Approaching the TATB size limit, the difference in the inner shell term diminishes for DMSO and water. On the other hand, the gap remains relatively consistent across all ionic radii in 1,2-DCE. This behavior is mirrored in the size of the finite-size corrected fluctuation of the electrostatic potential felt by Cl⁻- and TA⁺/TB⁻-sized ions; see Fig. 4 and Figs. S3 and S8 of the supplementary material. The increased radius leads to a significant reduction in the magnitude of the fluctuation term for water and DMSO, while very little change is observed for 1,2-DCE. Since we expect the outer shell term less the local potential to more or less follow the Born solvation model, we expect any lingering ion specificity in the inner shell term to make up the bulk of the deviation from ideal TATB-like behavior.

Figure 5 and Figs. S2 and S7 of the supplementary material show that the mean electrostatic potential at the center of large and small particles has two linear regimes across the line of zero charge. The asymmetry is noticeably greater for smaller ions, though it does not vanish entirely at the TATB limit. A best fit line through the q > 0 points intersects m_c = 0 in DMSO and water. This suggests solvent dipoles order around the neutral cavity much as they would around a cation. The opposite behavior was observed in 1,2-DCE. Deviation from m_c = 0 observed in the other trend line has been used previously as a crude estimate of ϕ_lp, though the exclusion volume is conditioned by the pair-specific Lennard Jones σ_ij term.⁷¹ 

The $μosex$ values approaching the TATB limit using the cumulant expansion formula in Eq. (6) are shown in Fig. 6 and Figs. S4 and S9 of the supplementary material. vdW contributions are left in, producing the slight curvature. Across the range of radii, these plots show a nearly constant gap between complementary sets of ions. We argue that this gap is due to the presence of a qϕ_lp factor for each ion. Thus, the figures also illustrate the trends following removal of this factor. This shifts the free energy differences into near perfect overlap for water. Some residual is present for the organic solvents due possibly to higher order effects not considered here. The local potential in each solvent [in (mean + third cumulant); total ± error format] was determined to be (8.76 + 0.20); 8.95 ± 0.14 kcal/mol e, (6.71 + −0.27); 6.44 ± 0.24 kcal/mol e, and (−1.42 + 0.16); −1.26 ± 0.20 kcal/mol e in water, DMSO, and 1,2-DCE, respectively. These potentials were sampled from configurations generated with an empty 9 Å radius cavity buried in each solvent. This method was shown by Ashbaugh to produce well-converged results at the same length scale we used here.³⁸ 

All values used to assemble the total free energy in both intrinsic and bulk thermodynamic scales are reported in the supplementary material in a tabular format. Differences between the single-ion solvation free energies listed in those tables are shown in Fig. 7 and Figs. S5 and S10 of the supplementary material. Where appropriate, we have included data points taken from the studies of Schurhammer and Wipff⁴⁶ and Carvalho and Pliego,⁴⁰ highlighting our agreement with these studies and clearly showing their results fall on different thermodynamic scales. Removal of qϕ_lp from $μosex$ nearly extinguishes the artificial stabilization of anions relative to their cation analog in both SPC/E water and OPLS/AA DMSO that produced the ∼20 kcal/mol disparity in TB versus TA solvation free energies reported by Schurhammer and Wipff.^46–48 Withdrawing qϕ_lp from $μosex$ in 1,2-DCE created a slightly greater difference in the ion solvation free energies with a deviation from the expected TATB behavior of about 6-10 kcal/mol in favor of the anion. This range is similar to that of acetonitrile, as discussed in Ref. 40.

The observation that anions below the TATB size limit have more negative solvation free energies in water and 1,2-DCE (relative to cations) is consistent with previous experiments and simulations.^{1,40,46,72–75} Our results are also consistent with previous conductometric and simulation studies that predict that small cations are better solvated in DMSO.^40,72,76,77 We find a number of qualitative points of agreement for trends in size effects and transfer free energies in the literature as well, validating our approach.^42,72,77,78

Interestingly, the authors of Ref. 48 did compute ϕ_lp = 7.4 kcal/mol as the mean potential at the center of a neutral particle using the TIP3P water model. Adding twice this to the reported Δμ^ex = −21.0 kcal/mol leaves a difference of −6.2 kcal/mol, comparing more favorably against our figure. The authors also computed ϕ_lp for TIP5P and found it to be much smaller, on the order of about 1 kcal/mol e.⁴⁸ This has been attributed to the underestimated quadrupole and octupole moments of TIP5P compared to experimental and ab initio predicted values.^49,50,54 The smaller ϕ_lp means less of a difference between the intrinsic and bulk thermodynamic scales—and so the deviation from the TATB hypothesis was found to only be −3.2 kcal/mol. The sign reflects better solvation of the anion. These results point out the strong model dependence of the classical point charge results.

With respect to the TATB hypothesis, we found a 2.1 kcal/mol deviation from cation/anion solvation parity in SPC/E water. The error bars encompass the 1.3 kcal/mol prediction of Carvalho and Pliego⁴⁰ who combined gas phase electronic structure calculations at the MP2 and MP4 levels of theory with continuum solvation using the SMD (solvation model, density) method. The authors found a similarly small deviation from the TATB scale in DMSO of 1.8 kcal/mol which fell just outside the error bars of our result. A key tenet of the TATB hypothesis is that such ions give the same solvation free energy in every solvent. Calculations in Ref. 40 for acetonitrile, however, revealed a nearly 10 kcal/mol shift from TATB behavior. We observed a comparable disparity in the solvation free energies on the bulk/Marcus scale in 1,2-dichloroethane spurred by remaining solvation asymmetry in the inner shell term favoring the anion and a negative ϕ_lp which further separated Δμ^ex(X⁺ → X⁻).

This paper has presented free energy calculations that test the TATB hypothesis of single-ion solvation.²⁴ We have discussed results clearly showing that simulations of ionic species in periodic boundary conditions (with no vapor region) include a local potential contribution to the free energy. The potential is the consequence of broken multipolar symmetry relative to the bulk solution (that has an average charge density of zero).⁵⁴ The local potential impacts the computed solvation free energy, enthalpy, and, coupled with estimates of the surface potential of a distant liquid-vapor interface, related properties such as the pK_a’s of ionizable species near interfaces.^27,28 The local potential is solvent (and model) dependent, and it applies in increments of ±q, where q is the ion charge. For example, difficulties in the development of accurate force fields for divalent cations like Mg²⁺ and Ca²⁺ may in part be due to inconsistent handling of this potential (which would impart a ∼18 kcal/mol error in μ^ex and hex in SPC/E water). Since the TATB hypothesis includes no net interfacial potential contribution in the free energy, a proper comparison for a given solvent must subtract out the local potential when simulations are performed in periodic boundary conditions.

Based on the above analysis, we have re-examined the TATB assumption which postulates that oppositely charged ions in a large, shielding ligand network have the same solvation free energy in every solvent. Our analysis focused on a spherical ion roughly the same size as the TA cation and TB anion. Thus we did not further probe expected small differences due to detailed molecular interactions in the first shell. As discussed in the Introduction, the thermodynamic evidence suggests that these inner-shell differences are of magnitude less than 2 kcal/mol.

Previous work had shown the anion of this pair to be on the order of ∼20 kcal/mol more strongly hydrated than the cation.^46–48 We showed that this difference reflected the influence of a factor of qϕ_lp for each ion and that these ions possessed nearly identical hydration free energies in the absence of the potential. We also extended our analysis to ion solvation in dimethyl sulfoxide (DMSO) and 1,2-dichloroethane (1,2-DCE) using the OPLS/AA force field. It was found that DMSO behaved similarly to water, where removing the potential allowed us to compare favorably against a recent MP4 quantum/continuum QCT model investigation,⁴⁰ while 1,2-DCE violated this assumption. We conclude that, outside of detailed molecular interactions in the first solvation shell that we expect to be relatively small based on comparisons with widely different approaches, the TATB assumption is approximately valid in water and DMSO (despite the possible classical model deficiencies). There are solvents such as 1,2-DCE that do not appear to satisfy the assumption, however, even for idealized spherical solutes. Future work should focus on quantum density functional theory (DFT) simulations that produce the free energy change to convert the TA cation into the TB anion. We reiterate that the quantum QCT calculations in Ref. 40 are consistent with the above conclusions.

We close with a brief speculation about the molecular origin of the deviation from the TATB assumption in 1,2-DCE. As discussed above, free energy differences between the TA and TB ions are largely located in the inner shell term of quasichemical theory. That contribution is equivalent to minus the work to expand a cavity of a prescribed size around the ion in the solvent (with the ion present). For the case of 1,2-DCE, the chlorine atoms (which carry partial negative charges) are significantly larger than the rest of the atoms that bear the partial positive charges. Thus, in the presence of an anion, the Cl atoms tend to be oriented away from the ion. Due to the significant curvature of the spherical space around the ion, the molecular geometry could figure prominently in the effective “surface tension” around the cavity. For the alternate case of a cation at the cavity center, configurations with the Cl atoms pointed inward will be electrostatically favored, but the larger size will create repulsive interactions that effectively reduce that “surface tension” around the cavity. Thus the work to grow a large cavity around the cation should be smaller in magnitude relative to that for an anion, consistent with our observations. This discussion implicates the geometry of the solvent molecules as a possible determinant of the accuracy of the TATB assumption. Presumably, the size disparity is not as prevalent for the water and DMSO solvent molecules.⁷⁹ Further work is required to address these interesting questions that link electrostatic effects and molecular geometry.

See supplementary material for extensive data on the quasichemical parts of the free energies of solvation, and the intrinsic and bulk solvation free energies for the ions considered in the main body of the paper. In addition, all of the figures for the DMSO and 1,2-DCE solvents that correspond to those for water in the main text are displayed.

This research was supported by NSF Grant Nos. CHE-1266105 and CHE-1565632. We are also grateful to the Ohio Supercomputer Center for a generous grant of computational resources. We thank Lawrence Pratt, Chris Mundy, Tim Duignan, and Shawn Kathmann for many helpful and informative discussions.