Quasi-Chemical Theory for Anion Hydration and Specific Ion Effects: Cl$^-$(aq) \emph{vs.} F$^-$(aq)

TL;DR

This paper applies quasi-chemical theory and dynamics calculations to compare hydration structures and free energies of chloride and fluoride ions, providing insights into specific ion effects relevant to biological and chemical systems.

Contribution

It introduces a combined approach of electronic structure and dynamics to accurately model anharmonicity and competing hydration structures of Cl$^-$ and F$^-$ ions.

Findings

Cl$^-$ and F$^-$ hydration structures differ significantly.

Predicted hydration free energy difference matches experimental data.

Cl$^-$ exhibits more flexible hydration interactions than F$^-$.

Abstract

Anion hydration is complicated by H-bond donation between neighboring water molecules in addition to H-bond donation to the anion. This situation can lead to competing structures for chemically simple clusters like (H$_2$O)$_n$Cl$^-$ and to anharmonic vibrational motions. Quasi-chemical theory builds from electronic structure treatment of isolated ion-water clusters, partitions the hydration free energy into inner-shell and outer-shell contributions, and provides a general statistical mechanical framework to study complications of anion hydration. The present study exploits dynamics calculations on isolated (H$_2$O)$_n$Cl$^-$ clusters to account for anharmonicity, utilizing ADMP (atom-centered basis sets and density-matrix propagation) tools. Comparing singly hydrated F$^-$ and Cl$^-$ clusters, classic OH-bond donation to the anion occurs for F$^-$, while Cl$^-$ clusters exhibit more flexible but dipole-dominated interactions between ligand and ion. The predicted Cl$^-$ -- F$^-$ hydration free energy difference agrees well with experiment, a significant theoretical step for addressing issues like Hofmeister ranking and selectivity in ion channels.

Equations

$$\begin{gathered} {\mu }_{{\mathrm{X}}^{-}}^{(\mathrm{ex})}=-RT\ln K_n^{(0)}{\rho }_{{\mathrm{H}}_2\mathrm{O}}{}^n+RT\ln p_{{\mathrm{X}}^{-}}(n) \\ +\left({\mu }_{{\left({\mathrm{H}}_2\mathrm{O}\right)}_n{\mathrm{X}}^{-}}^{(\mathrm{ex})}-n{\mu }_{{\mathrm{H}}_2\mathrm{O}}^{(\mathrm{ex})}\right), \end{gathered}$$

$$n{\mathrm{H}}_2\mathrm{O}+\mathrm{X}{}^{-}\stackrel{K_n^{(0)}}{\rightleftharpoons }{\left({\mathrm{H}}_2\mathrm{O}\right)}_n\mathrm{X}{}^{-}$$

$$K_n^{(0)}$$

$$\begin{gathered} \frac{K_n{}^{(0)}}{K_{n-1}{}^{(0)}{K}_1{}^{(0)}}= \\ \frac{\mathcal{Q}\left({\gamma }_n\sigma \right)/\mathcal{Q}(\sigma )}{n\left[\mathcal{Q}\left({\gamma }_{n-1}\sigma \right)/\mathcal{Q}(\sigma )\right]\left[\mathcal{Q}\left(\gamma \sigma \right)/\mathcal{Q}(\sigma )\right]}. \end{gathered}$$

$$\begin{gathered} \Delta U=E\left({\gamma }_n\sigma \right)-E\left({\gamma }_{n-1}\sigma \right) \\ -E\left(\gamma \sigma \right)+E(\sigma ), \end{gathered}$$

$$n{K}_n{}^{(0)}=\frac{K_1{}^{(0)}{K}_{n-1}{}^{(0)}}{{⟨e^{\beta \Delta U}⟩}_n}.$$

$${\mu }_{{({\mathrm{H}}_2\mathrm{O})}_n{\mathrm{Cl}}^{-}}^{(\mathrm{ex})}=-RT\ln \left[\left(\frac{1}{N}\right)\sum _{j=1}^N{\mathrm{e}}^{-{\bar{\epsilon }}_j/RT}\right],$$

Full text

1 Introduction

Here, we build a fundamental molecular statistical mechanical theory for hydration of chloride ion (Cl-). Theory of the type targeted here, specifically quasi-chemical theory (QCT), permits transparent comparison of Cl- to analogous cases such as fluoride hydration (F-). Comparisons of that sort underpin specific ion effects in areas of broad current interest, including the Hofmeister ranking1, 2, 3, 4, 5, 6 of ions and the mechanism of specific ion permeation in membrane ion channels. An example of membrane ion channels is the CLC family, which regulate membrane transport of Cl-.7 Though we comment further below on such long-term goals, we focus primarily on extending QCT to hydrated anions. One motivation for comparisons such as Cl- with F- is that they permit rigorous experimental thermodynamic testing of single-ion free energies without discussion of potentials of the phase,8, 9 or surface potentials.10, 11, 12

Adequate molecular simulation is necessary for compelling molecular statistical thermodynamic theory of liquid solutions, and extensive simulation studies have been carried out.13, 14, 15, 16, 17, 18, 19, 20, 21, 22 Available computational work has nicely delineated distinctions in hydration structure of anions compared to simple metal ions, and distinctions between different halide ions.23, 24, 25, 26 In this crowded landscape of computational work, molecular-scale statistical mechanical theory for these systems has been meager, with striking exceptions.27

Quasi-chemical theory (QCT), in the present incarnation based on clusters,28, 29 provides a natural starting point for the theory sought.30, 6 QCT31, 32, 33, 34, 35, 28, 29 naturally divides the free energy determination into inner-shell and outer-shell considerations, as is apparent from Eq. (1) below. In this approach, the necessary contributions are physically meaningful,36 and can then be evaluated to the desired accuracy — even judiciously using experimental results for natural intermediate results.37 QCT, considered broadly over nearly two decades, has been applied successfully to a wide range of systems, including the hard sphere fluid,38 liquid water,35, 39, 36, small molecule solvation40, 41, 42, 43, 44, 45 also involving ion channels46, 47, 48, 49 and other macromolecules,50, 51, 52 cations in water and non-aqueous solvents,32, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62 and biomolecule hydration.63, 64

The QCT approach provides a complete statistical mechanical framework for evaluating solvation free energies of ions by building up from carefully defined isolated ion-solvent clusters.28 Experimental measurements addressing the energetics of such clusters, based on mass spectroscopy,65, 66, 67 have existed for some decades. For hydrous clusters, a compilation of those experimental values are available from Tissandier, et al.68

QCT has been applied less to the challenging cases of anion hydration than to metal cations, comparatively. A primary difficulty is hydrogen bonding between water molecules of clusters such as (H2O)${}_{n}\mathrm{Cl}^{-}$, leading to structures that are less susceptible to simple approximation in the implementation of QCT. Recent efforts have been directed toward addressing those initial approximations in the application of QCT. Some of the refinements include quantification of anharmonic effects41, 44, 48, 37 on free energies of solute-water clusters and the sufficiency of the polarizable continuum model69, 41, 60 (PCM) for the hydration free energy of those clusters.

Anticipating results below (Section 3.2), a harmonic approximation to the potential energy surface of (H2O)${}_{n}\mathrm{Cl}^{-}$ clusters is satisfactory for cluster formation enthalpies, but not for free energies. To investigate such issues, we use a molecular dynamics approach that explicitly includes electron coordinates, as we now discuss.

Classical trajectory calculations that include electron coordinates fall into two major categories: Born-Oppenheimer (BO) molecular dynamics (MD) and extended Lagrangian (EL) dynamics. In the BO approach, the atomic dynamics develop on a potential energy surface obtained from an optimized-in-detail electronic structure calculation. In EL approaches such as Car-Parrinello molecular dynamics (CPMD),70 electron coordinates are introduced as fictitious dynamic variables that couple electronic and atomic motion. The dynamics should keep the fictitious electronic degrees of freedom near the ground state BO potential energy surface, yielding sufficiently accurate forces for the atomic trajectory.

The atom-centered basis sets and density matrix propagation (ADMP)71 tool is an extended Lagrangian approach for molecular dynamics. In contrast with CPMD, ADMP emphasizes finite non-periodic systems, and thus is convenient for the isolated clusters of interest here (Figure 1). ADMP has been shown to provide similar functionality to BO molecular dynamics.71

1.1 Perspective from interest in Cl- and F- ion channels

The regulation of ion concentration is essential for several physiological functions in cells.72 For anions such as Cl-, this control is achieved through membrane channels and transport proteins from the CLC family, which facilitate the passage of Cl- through electro-chemical potential gradients.7 Numerous efforts73, 74, 75, 76 have been directed toward deciphering the structures of those proteins and understanding their functional characteristics, such as conductance, gating, and ion selectivity.

While CLC proteins select for Cl-, other channel proteins discriminate against Cl-. An interesting example is FLUC, a family of fluoride-specific ion channels with dual-topology architecture.77 These channels display an astonishingly high selectivity of $10^{4}$ for F*-* over Cl- despite close similarity in size and identical charge of the ions. FLUC channels provide a ladder of hydrogen-bond donating residues that apparently create a polar track for F-. How that polar track leads to the unusually high level of discrimination between F- and Cl- remains an open question. Understanding such mechanisms for selectivity in ion transport has been a primary target for many recent modeling and simulation studies, emphasizing K+/Na+ selectivity of potassium ion channels.78, 46, 79, 80, 50, 81, 82, 83, 47, 84, 85, 86, 87, 88, 89, 90, 91, 92 But computational studies of anion transport mechanisms demonstrated by CLC channels93, 94, 95, 96, 97, 98, 99 that address comparison to alternatives such as FLUC are less mature.

To address ion selectivity between different media, a rigorous treatment would require comparison of the total free energy of ion solvation between those environments. However, ion channels are thought to function primarily by manipulation of the local environment of the ion. That manipulation may arise from constraints on local structure imposed by the structural and chemical properties of the surrounding environment or the binding sites on the channel itself.46, 100, 79, 80, 50, 57, 82, 47, 48, 29, 52, 49 Furthermore, contributions to solvation free energy from long-ranged electrostatic interactions between an ion and the distant environment might cancel between ions of the same charge. Therefore, analysis of selectivity might not require a detailed description of the long-ranged interactions in the molecular theory. Instead, rigorous treatment of local interactions are likely important for accurate treatment of ion selectivity.

Here, we compare the hydration of Cl- with F-, focusing on rigorous evaluation of local interactions. Broader issues such as the Hofmeister ranking of ions and selectivity in membrane channels are reserved for the future. In the next two sections we lay out the theory and results, followed by a concluding discussion. The details of the numerical implementation are specified in Section 5.

2 Theory: Step-wise evaluation of the

isolated cluster free energy

QCT formulates the net hydration free energy as

[TABLE]

collecting the individual quasi-chemical contributions.28, 29 The first and last term represent the inner-shell and outer-shell contributions defined based on assignment of a clustering radius ($\lambda$). $p_{\mathrm{X}^{-}}(n)$ then identifies the thermal probability of observing $n$ waters within that prescribed $\lambda$. The several features of this expression will be identified in the following discussion.

To that end, consider the chemical process

[TABLE]

for formation of the isolated cluster, as in an ideal gas phase (X = Cl here). We seek the free energy (inner-shell) contribution $-RT\ln K_{n}^{(0)}\rho_{\mathrm{H_{2}O}}{}^{n}$, with

[TABLE]

a classic aspect of statistical thermodynamics.34, 37 Here $V$ is the system volume, $\sigma$ and $\gamma$ represent the solute (Cl-) and solvent (H2O) molecules, respectively, and the $\mathcal{Q}\left(\gamma_{n}\sigma\right)$ are configurational integrals associated with single molecule/cluster canonical partition functions.37

Our scheme for evaluating $K_{n}{}^{(0)}$ proceeds step-wise in $n$ according to

[TABLE]

The numerator on the right of Eq. (4) involves integration carried over canonical distributions of configurations of $n$ solvent molecules, while the denominator treats configurations of $n-1$ and 1, respectively, treated independently. Our scheme is then based on evaluation of

[TABLE]

so that

[TABLE]

The brackets, $\langle\ldots\rangle_{n}$, indicate the thermal average utilizing configurations from the canonical simulation stream for the $\gamma_{n}\sigma$ cluster.

For $n=1$, Eq. (6) reduces to the trivial case of $K_{0}{}^{(0)}=1$. In the evaluation of $K_{n}{}^{(0)}$ for $n\geq 2$, the value of $K_{1}{}^{(0)}$ can be supplied from experiment68 or theory. This term describes the interaction between Cl- and one water molecule. We note that for (H2O)Cl-, a harmonic approximation gives a sufficiently good estimate for the free energy, but not for bigger cluster sizes (Figure 5). Carrying out subsequent steps in this scheme then addresses the issues that make anion hydration more challenging; that is, the sampling of cluster configurations, 37 including competing H-bonding interactions between neighboring water molecules in those clusters.

Here, we take advantage of the QCT partitioning to evaluate the inner-shell contribution based on accurate quantum mechanical treatment of dynamics (using ADMP) in ion-water clusters. The evaluation of $\left(\mu^{\mathrm{(ex)}}_{\mathrm{(H_{2}O)}_{n}\mathrm{Cl}^{-}}-n\mu^{\mathrm{(ex)}}_{\mathrm{H_{2}O}}\right)$ utilizes the same trajectory, but approximates the outer-shell waters using the polarizable continuum model69 (PCM), as described in section 5. Experimental comparisons (section 7) can then indicate the sufficiency of such an approximation.

3 Results

3.1 Solution structure

Neighborship is a fundamental concept in QCT.46, 28, 59, 37, 101 The H-atom of water is the closest neighbor to the anions, F- and Cl-. Proximity distinctions between those cases are sharper when H atoms are utilized rather than O atoms in the radial distributions. Distributions of water H atoms relative to the ions (Figure 3) then guide the assignment of the clustering radius, $\lambda$. Further, the neighborship-ordered contributions help identify cluster sizes in which waters are always in direct contact with the ion. Notice that the distribution of the 6${}^{\text{th}}$ nearest water around Cl- extends beyond the first minimum of $g_{\mathrm{H|Cl}}\left(r\right)$. For instance, a choice of $\lambda_{\mathrm{Cl}^{-}}\leq 0.26$ nm eliminates such split-shell clusters.46, 59, 102 The thermal probability ($p_{\mathrm{Cl}^{-}}(n)$), contributing to the second term of Eq. (1), is then evaluated based on this constraint.

Compared to Cl-(aq), the radial distribution of F-(aq) displays a larger maximum and a nearly zero first minimum, leading to a flat behavior of the running coordination number above the minimum (Figure 3). For F-(aq), this structure suggests tightly bound inner-shell waters, with small deviations from orientations that take the nearest H-atom away from the ion. By comparison, inner-shell waters of Cl-(aq) are more flexibly structured. The differences highlighted here also feature in the behavior of waters in the gas phase clusters of $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{Cl}^{-}$ and $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{F}^{-}$ (Figures 1 and 2).

3.2 Isolated cluster energetics: $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{Cl}^{-}$

To evaluate cluster energetics, a harmonic approximation was applied initially to the lowest energy geometry-optimized structures (SI, Figure S1). The enthalpies of $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{Cl}^{-}$ cluster calculated this way agree well with experiment (Figure 4). The UPBE1PBE functional with the aug-cc-pvdz basis set showed the best agreement among the two tested functionals. However, the free energy evaluated using the harmonic approximation deviates from experiment, and that deviation increases with $n$ (left panel, Figure 5). We conclude that the entropy of the $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{Cl}^{-}$ cluster is not satisfactorily predicted in our harmonic approximation.

Notice further that the harmonic approximation with the UPBE1PBE functional gives an accurate estimate of free energy for $\left(\mathrm{H}_{2}\mathrm{O}\right)\mathrm{Cl}^{-}$. Therefore, we expect the discrepancy in the free energy of $\left(\mathrm{H}_{2}\mathrm{O}\right)_{2}\mathrm{Cl}^{-}$ to emerge from interactions between water molecules, roughly with an average strength of an H-bond. Hence, we perform molecular dynamics of the cluster with ADMP for better sampling of relative configurations of water in the cluster. Those configurations are then utilized in the scheme above (Section 2) for evaluation of the free energy. The free energies evaluated this way are in excellent agreement with experimental gas phase cluster data (Figure 5). Our scheme (Section 2) is highly effective here because it is applied to the isolated cluster and interactions with only a single surface water molecule are manipulated.

3.3 Comparison with $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{F}^{-}$

In light of the successful treatment of chloride clusters with ADMP, we next take up the case of fluoride. For comparison, the same steps are followed as before (Section 2) to treat the $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{F}^{-}$ cluster. Beyond the complication of ligand-ligand H-bonding, we initially note that the hydration structures of $\mathrm{H_{2}O}-\mathrm{F}^{-}$ and $\mathrm{H_{2}O}-\mathrm{Cl}^{-}$ are qualitatively different. In the former case, classic OH-bond donation characterizes the clusters, while in the latter case, more flexible but ‘dipole-dominated’ configurations occur (Figure 1). Another interesting observation is that, in contrast to $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{Cl}^{-}$, the harmonic approximation satisfactorily predicts free energy in the case of $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{F}^{-}$ for $n\leq 4$ (right panel, Figure 5).

These observations on the differences between $\mathrm{Cl}^{-}$ and $\mathrm{F}^{-}$ can be explained as follows. Firstly, spectroscopic studies103 of halide-water clusters based on Ar predissociation have identified vibrational bands that correspond to inter-water H-bonding. The spectra revealed that those interactions gradually weaken when going from $\left(\mathrm{H}_{2}\mathrm{O}\right)_{2}\mathrm{I}^{-}$ to $\left(\mathrm{H}_{2}\mathrm{O}\right)_{2}\mathrm{F}^{-}$ clusters. For $\left(\mathrm{H}_{2}\mathrm{O}\right)_{2}\mathrm{F}^{-}$, that vibrational band disappears, with water molecules separating entirely. This result is consistent with our geometry optimization calculations for $\left(\mathrm{H}_{2}\mathrm{O}\right)_{2}\mathrm{F}^{-}$ (Figure 6). The fact that a harmonic approximation reproduces the experimental free energy at 298 K when applied to that minimum energy structure indicates that inter-water H-bonding between waters is not significant in $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{F}^{-}$. That argument extends up to $n=4$ because no H-bond exists in those energy-optimized structures.

Secondly, previous ab initio MD work 104 compared the computed IR absorption spectrum for $\left(\mathrm{H}_{2}\mathrm{O}\right)_{6}\mathrm{Cl}^{-}$ with experimental spectra obtained for $\left(\mathrm{H}_{2}\mathrm{O}\right)_{5}\mathrm{Cl}^{-}$. That work suggested that accounting for anharmonicity and coupling between modes should be important for the treatment of inter-water H-bond dynamics in aqueous clusters. Here, our step-wise scheme (Section 2) based on the ADMP approach naturally includes those dynamical features. This quantitative treatment addresses the competing inter-water H-bond dynamics in $\left(\mathrm{H}_{2}\mathrm{O}\right)_{n}\mathrm{Cl}^{-}$ ($n\geq 2$), leading to cluster free energies in excellent agreement with experiments. The earlier suggestion also explains why the harmonic approximation works well for the case of $\left(\mathrm{H}_{2}\mathrm{O}\right)\mathrm{Cl}^{-}$, where inter-water interactions are absent.

3.4 Free Energy of Hydration

Summing the quasi-chemical components of Eq.(1), we obtain an estimate of the net hydration free energy ($\mu^{\mathrm{(ex)}}_{\mathrm{X}^{-}}$, Figure 7). Note that the smallest contribution to the hydration free energy comes from $RT\ln p_{\mathrm{X}^{-}}(n)$, which accounts for heterogeneity in occupancy of the inner solvation shell. For the most probable occupation cases ($n=4,5$), this contribution is inconsequential ($\approx-0.5~{}$kcal/mol) compared to the others. Next, the inner-shell interactions are built up by augmenting the isolated cluster free energy (evaluated at 1 atm, Figure 5) with $-nRT\ln 1354$, the ligand replacement contribution that accounts for the actual density of ligands available in liquid phase at standard conditions, 1 g/cm3. The full inner shell free energy contribution $\left(-RT\ln K^{(0)}_{n}\rho_{\mathrm{H_{2}O}}{}^{n}\right)$ increases in magnitude with $n$. The outer-shell contribution $\left(\mu^{\mathrm{(ex)}}_{\mathrm{(H_{2}O)_{n}X^{-}}}-n\mu^{\mathrm{(ex)}}_{\mathrm{H_{2}O}}\right)$ is treated using the polarizable continuum model (PCM) and decreases in magnitude with $n$. Finally, these contributions balance to produce a hydration free energy that is independent of $n$ (Figure 7). These evaluations agree reasonably well with experimental hydration free energies105 for Cl- ($-81.3$ kcal/mol) and F- ($-111.1$ kcal/mol).

As an indication of the numerical sensitivity of such calculations, we note a previous QCT treatment that arrived at a value of $-228$ kcal/mol for the LiF (aq) pair.37 For comparison, an AIMD calculation based on an electron density functional obtained $-240$ kcal/mol.106 These two results bracket experimental values, nearly equally. When the latter result was corrected ex post facto by MP2 calculations on a finite solution fragment for the F- case, the corrected value splits the indicated difference. This comparison suggests the importance of detailed treatment of dispersion effects for anions, though not necessarily for cations.107

Since QCT explicitly discriminates between inner- and outer-shell effects, comparison of the different free energy contributions to hydration of Cl- and F- leads to some interesting observations (Figure 7). The cluster hydration contributions, $\left(\mu^{\mathrm{(ex)}}_{\mathrm{(H_{2}O)_{n}X^{-}}}-n\mu^{\mathrm{(ex)}}_{\mathrm{H_{2}O}}\right)$, are nearly identical. This quantitatively supports the intuition that long-ranged interactions, predominantly electrostatic, do not discriminate between Cl- and F-. In contrast, the inner-shell interactions, $-RT\ln K^{(0)}_{n}\rho_{\mathrm{H_{2}O}}{}^{n}$, differ by 24 kcal/mol, essentially making up the difference (26.3 kcal/mol) in the hydration free energies. Clearly, disparities in the size and electronic structure of these ions contribute to the differences reported here, but those effects are localized to the inner solvation shell. Thus, the treatment of specific ion effects can focus primarily on accurate evaluation of inner-shell contributions to free energy through rigorous treatment of dynamics in small clusters, as demonstrated here.

4 Discussion and Conclusions

These results confirm again that a foremost difficulty in treating hydrated anions such as Cl-(aq) is significant water-water H-bonding within the inner solvation shell (Figure 6). Nevertheless, the hydration structures of one-water molecule cases, $\mathrm{H_{2}O}-\mathrm{F}^{-}$ and $\mathrm{H_{2}O}-\mathrm{Cl}^{-}$, differ qualitatively (Figures 1 and 2), with the latter complex being qualitatively a dipole-dominated structure. Of course, this feature facilitates inter-ligand H-bonding for the Cl-(aq) case. The ADMP computational tool works satisfactorily for dynamics of the isolated inner-shell clusters that are treated in QCT. But details of the trajectory thermostat algorithm deserve further elucidation (Supporting Information). A one-ligand stepwise evaluation of the cluster free energies works surprisingly well, and those computational results agree nicely with experimental cluster free energies (Figure 5).

Collection of the several factors contributing to the QCT free energy (Figure 7) produces net results in good agreement with experimental tabulations.108, 105 More decisively, the difference $\mu^{\mathrm{(ex)}}_{\mathrm{Cl}^{-}}$ - $\mu^{\mathrm{(ex)}}_{\mathrm{F}^{-}}$ is not subject to an ambiguity of a potential of the phase (or ‘surface potential’). Thermodynamic results then provide an unambiguous test for the difference. Here, the QCT result for the difference is about 2-3 kcal/mol smaller than the experimental value for that difference, about -29 kcal/mol. This discrepancy is likely due to use of the PCM treatment for the cluster-hydration contribution, which treats long-range interactions. Thus, differential outer-shell structuring not considered by PCM is likely a principal contribution to the difference. The AIMD results presented here (Figure 3) characterize that outer-shell structure. Note that the QCT approach uses the PCM approximation in a context where it is insensitive to adjustment of dielectric radii, so the discrepancy identified here is likely to be intrinsic to the PCM model.

The QCT analysis assigns (Figure 7) 90% of the hydration free energy difference $\mu^{\mathrm{(ex)}}_{\mathrm{Cl}^{-}}$

• $\mu^{\mathrm{(ex)}}_{\mathrm{F}^{-}}$ to inner-shell contributions deriving from the isolated hydrous cluster. This result is an important idea to keep in mind in addressing Hofmeister effects and selectivity of ion channels.109, 110 Precision in analyzing single-ion free energies should assist in understanding selectivity in ion channels by establishing end-point free energies bracketing a path for ion transit.

5 Methods

5.1 ADMP

We evaluated molecular dynamics of the isolated (H2O)${}_{n}\mathrm{Cl}^{-}$ clusters for $1\leq n\leq 5$ using atom-centered basis sets and the density matrix propagation (ADMP) approach71 available through the Gaussian09 software package.111 The UPBE1PBE density functional was utilized with the aug-cc-pvdz basis set. The fictitious masses that couple electron motion with the nuclei was set to 0.1 amu. The initial velocities of the individual atoms were sampled randomly from a Boltzmann distribution, and the initial density matrix velocity was chosen to be zero. With a time step of 0.1 fs, dynamic simulations were carried out for 1 ps at 300 K for all cluster sizes, and the last 0.5 ps were used for analysis.

One reservation in using Gaussian09 for dynamics is that description of the only implemented thermostat option has been sketchy.71 Operationally, this algorithm constrains the total nuclear kinetic energy of the system. Intuitive, non-canonical thermostats are known to cause problems in some cases, as in the flying ice cube effect.112 Nevertheless, due to lack of alternatives in Gaussian09, we used that thermostat to maintain the system at 300 K.

Though reservations with that thermostat should be borne in mind, they do not appear to be serious here. First, the atomic velocities are verified to be Maxwell-Boltzmann distributed (SI, Figure S2). In addition, we performed 10 ps of ab initio molecular dynamics at 300 K with the Nose-Hoover thermostat113 for $\left(\mathrm{H}_{2}\mathrm{O}\right)_{2}\mathrm{Cl}^{-}$ using CP2K.114 The PBE density functional115 was utilized with GTH116 pseudopotentials in the Gaussian and plane wave schemes.117 Molecularly optimized DZVP-MOLOPT-GTH118 basis sets were obtained from the CP2K website. The structures sampled by this CP2K trajectory were then analyzed with Gaussian09 using the UPBE1PBE density functional and the aug-cc-pvdz basis set. The isolated cluster free energy obtained this way is within 0.1 kcal/mol of the ADMP trajectory result. A detailed comparison of the cluster simulation methods discussed above is provided with the Supporting Information (SI, Figure S3).

5.2 PCM for ADMP sampled clusters

The outer-shell cluster contribution $\left(\mu^{\mathrm{(ex)}}_{\mathrm{(H_{2}O)}_{n}\mathrm{Cl}^{-}}-n\mu^{\mathrm{(ex)}}_{\mathrm{H_{2}O}}\right)$ to the hydration free energy is treated using the polarizable continuum model69 (PCM) in Gaussian09.111 Configurations that obey the clustering constraint were extracted from the last 1 ps of the ADMP trajectory. The sampled geometrical structures were subjected to two (2) single point electronic calculations, separately, one for the isolated cluster and a second with the external (dielectric) medium described by the PCM tool. The difference, $\bar{\varepsilon}_{j}$, is employed in computing

[TABLE]

where $N$ is the number of configurations from the simulation stream that obey the clustering constraint. Eq. (7) corresponds to the potential distribution theorem (PDT) approach,34, 28 recognizing that thermal fluctuations implicit in the PCM69 approach complete the PDT averaging.

5.3 AIMD

A system consisting of a single halide ion (Cl- or F-) and 64 waters was simulated using the VASP AIMD simulation package.119, 120 A wide range of alternative electron density functionals are available, in principal,121 for such calculations. Here we chose the PW91 generalized gradient approximation principally for consistency with previous calculations.122 The plane wave basis had a high kinetic energy cutoff of 400 eV.

A cubic cell of 1.24 nm was used to set a satisfactory density of the water. The calculation adopted the Nosé thermostat procedure113 at 350 K. With respect to the ADMP calculations with $T$ = 300 K and with the water density set, that temperature adjustment amounts to a scaling of the potential energy by 300/350. That scaling at constant density avoids glassy behavior that can result otherwise.123

A time step of 0.5 fs was chosen to realize a 100 ps trajectory. The last 50 ps were used for analysis. The radial distribution functions, g${}_{\mathrm{H}|\mathrm{X}}(r)$, and the inner-shell occupancy probabilities, $p_{\mathrm{X}^{-}}(n)$, were evaluated on that basis. The radial distribution functions (Figure 3) provide qualitative guidance to the present statistical theory, but not numerical input. The occupancy probabilities, $p_{\mathrm{X}^{-}}(n)$, contribute to Eq. (1) but at the boundary of significance for these studies (Figure 7).

6 Acknowledgements

This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. The work was supported by Sandia National Laboratories (SNL) LDRD program. SNL is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE’s National Nuclear Security Administration under contract DE-NA-0003525. The views expressed in the article do not necessarily represent the views of the U.S. DOE or the United States Government. LRP thanks Dilip Asthagiri and Diego Gomez for helpful discussions.

Supporting Information

Geometry optimizations

Thermostat testing

Comparison of cluster simulation methods

Gaussian09 trajectory:

The ADMP approach utilized the UPBE1PBE density functional with the aug-cc-pvdz basis set. The reservation in using Gaussian09 for dynamics is that the only available thermostat option is not fully described. Nevertheless, the velocity distribution obtained here (Figure 9) matches Maxwell-Boltzmann distribution. For further testing, we use CP2K, as described below.

CP2K trajectory:

CP2K allows molecular dynamics with the Nose-Hoover thermostat.113 The PBE density functional115 was utilized with pseudopotentials proposed by Goedecker, Teter and Hutter (GTH116) in the hybrid Gaussian and plane wave scheme.117 Molecularly optimized basis sets denoted as DZVP-MOLOPT-GTH118 were obtained from the CP2K website. The trajectory obtained is then analysed in Gaussian09 using the UPBE1PBE density functional with the aug-cc-pvdz basis set.