Interaction Energy Analysis of Monovalent Inorganic Anions in Bulk Water Versus Air/Water Interface

Soft anions exhibit surface activity at the air/water interface that can be probed using surface-sensitive vibrational spectroscopy, but the structural implications of this surface activity remain a matter of debate. Here, we examine the nature of anion–water interactions at the air/water interface using a combination of molecular dynamics simulations and quantum-mechanical energy decomposition analysis based on symmetry-adapted perturbation theory. Results are presented for a set of monovalent anions, including Cl−, Br−, I−, CN−, OCN−, SCN−, NO2−, NO3−, and ClOn− (n=1,2,3,4), several of which are archetypal examples of surface-active species. In all cases, we find that average anion–water interaction energies are systematically larger in bulk water although the difference (with respect to the same quantity computed in the interfacial environment) is well within the magnitude of the instantaneous fluctuations. Specifically for the surface-active species Br−(aq), I−(aq), ClO4−(aq), and SCN−(aq), and also for ClO−(aq), the charge-transfer (CT) energy is found to be larger at the interface than it is in bulk water, by an amount that is greater than the standard deviation of the fluctuations. The Cl−(aq) ion has a slightly larger CT energy at the interface, but NO3−(aq) does not; these two species are borderline cases where consensus is lacking regarding their surface activity. However, CT stabilization amounts to <20% of the total induction energy for each of the ions considered here, and CT-free polarization energies are systematically larger in bulk water in all cases. As such, the role of these effects in the surface activity of soft anions remains unclear. This analysis complements our recent work suggesting that the short-range solvation structure around these ions is scarcely different at the air/water interface from what it is in bulk water. Together, these observations suggest that changes in first-shell hydration structure around soft anions cannot explain observed surface activities.


Introduction
One of the earliest results of surface-sensitive vibrational sum-frequency generation (VSFG) experiments [1,2] was the observation that soft anions impact the vibrational lineshape in the O-H stretching region, but that hard anions do not [3][4][5][6]. The term "soft" is chosen carefully here, as an alternative to "polarizable"; it can be roughly interpreted as monovalent and polarizable, equivalent to having a low surface charge density [7], and such ions are sometimes called "chaotropic" [8]. Although the surface activity of certain anions is often discussed in terms of polarizability [9][10][11][12][13][14][15][16][17], it should be noted that polyvalent anions such as SO 2− 4 (aq) are quite polarizable [18] but the presence of polyvalent anions in solution does not affect the O-H lineshape measured in VSFG experiments [19]. Molecular dynamics (MD) simulations suggest that hard anions, including polyvalent species but also fluoride, are repelled from the air/water interface [20,21].
The nature of the surface activity exhibited by soft anions remains a matter of debate. Whereas continuum electrostatics predicts that all ions are repelled from the air/water interface [13], a first wave of MD simulations using polarizable force fields suggested that soft anions are not only present at the interface but in fact partition preferentially there [9,13,20,22]. More recent work suggests that these concentration enhancements were exaggerated by the force fields in use at the time [23][24][25][26][27][28][29], which aligns with the interpretation of some of the early experiments [3]. According to this point of view, surface activity may simply reflect the absence of depletion of soft anions at the interface [30,31], rather than a concentration enhancement.
To this debate, the present authors have recently added the observation that the first-shell hydration structure around soft anions is hardly different at the air/water interface as compared to that in bulk water [7]. This observation comes from MD simulations using polarizable force fields, and such similarities had been noted previously in simulations of I − (aq) [32] and SCN − (aq) [33], in the latter case using ab initio MD. Iodide and thiocyanate are archetypal examples of ions that perturb the O-H lineshape in VSFG experiments [3,4,19,34,35]. Our work [7] considered a larger set of anions, and the structural similarities that we observe, including the number and orientation of the ion-water hydrogen bonds, suggest that the origins of anion-induced changes in the O-H vibrational lineshape must be rather subtle effects on water-water hydrogen bonds, perhaps due to ion-induced changes in local electric fields [36]. These observations need to be reconciled with the prevailing modern view that monovalent ions have little effect on the long-range hydrogen-bonding dynamics of liquid water [37], as measured by femtosecond vibrational pump-probe experiments [37][38][39][40], although the effects on the long-range hydrogen-bonding structure of water are less clear. Neutron diffraction experiments [41][42][43][44] and some MD simulations [45] do suggest that even monovalent ions alter the tetrahedral ordering of water beyond the first solvation shell, for solutions of NaOH(aq), HCl(aq), NaCl(aq), and KCl(aq). Pronounced structural changes have been documented in some cases involving polyvalent ions [46][47][48].
Our previous work [7] was limited to structural characterization of the ions in question, along with a detailed examination of their ionization energies in order to make contact with liquid microjet photoelectron spectroscopy [49]. The present work adds another dimension to this analysis as we compute anion-water interaction energies for the same set of anions: Cl − , Br − , I − , CN − , OCN − , SCN − , NO − 2 , NO − 3 , and ClO − n (n = 1, 2, 3, 4). Some of these are typical surface-active ions (e.g., Br − , I − , SCN − , and ClO − 4 ), whereas others (such as CN − , OCN − , and NO − 2 ) visit the interface much less frequently, according to the MD simulations [7], and are not classified as surface-active. Intermediate cases where the surface activity is weak, or where experimental consensus is lacking, include Cl − and NO − 3 [19]. Amongst these ions, our simulations indicate that even the ones that are not considered surface active nevertheless spend enough time near the air/water interface that it is possible to assemble an interfacial data set for them. These cases offer a useful comparison to the canonical surface-active anions.
We present a detailed analysis of the (ensemble-averaged) interaction between each of these ions and its short-range hydration sphere, in both bulk water and at the air/water interface, using the quantum-chemical methods of symmetry-adapted perturbation theory (SAPT) [50][51][52][53]. The SAPT family of methods [50,51] is designed for accurate calculation of noncovalent interaction energies, as well as a physically-motivated energy decomposition analysis of those energies [51,52]. Of key interest will be whether the interfacial environment engenders any discernible changes in the ion-water interactions, relative to what is observed for the same ion in bulk water.

Classical MD Simulations
MD simulations of the aforementioned ions in a periodic slab configuration were reported previously [7] and the same set of simulations is used here to obtain snapshots for interaction energy analysis. These simulations were performed under NVT conditions at T = 298 K and a bulk density of 0.997 g/cm 3 . The size of the periodic simulation cell (31.3 Å × 31.3 Å × 156.7 Å) was previously shown to afford converged results [7]. The simulation data were subsequently partitioned into bulk and interfacial parts depending on the position of the ion relative to the Gibbs dividing surface (GDS) that we take to define the air/water interface. For the snapshots classified as "interfacial", the ion's center of mass lies no more than 3 Å below the GDS. Anything beyond this cutoff is considered to be a bulk water environment, as this interior region of the periodic slab affords properties that are essentially indistinguishable from results performed in an isotropic simulation that has no interface [7]. Simulations were performed using the AMOEBA force field for water [54], whose parameterization includes some of the ions in question, such as the halides [55]. Parameters for the remaining ions were developed along similar lines [7], following an established protocol [56], and are included in the Supplementary Material. Energetic analyses with the AMOEBA force field were performed using Tinker software, v. 8 [57].
Following an equilibration period, snapshots were extracted that include two solvation shells around the ion, according to distance criteria described previously [7]. The number of water molecules varies from one snapshot to the next, with the average number N w depending on both the size of the ion and how tightly hydrated it is. In bulk water, these averages range from N w ≈ 28 for Cl − (aq) up to N w ≈ 43-44 for Br − (aq) and I − (aq), with N w = 35-37 for the remaining ions. The interfacial snapshots contain fewer water molecules, on average, as the water density is reduced in the interfacial region. In the analysis that follows, we consider interaction energies (E int ) between the ion and its first two hydration shells. The quantity E int is intensive with respect to system size and this insulates our analysis against the step-to-step fluctuations in the number of water molecules that are included in these calculations. Ensemble averages reported below represent 51 snapshots for each ion in bulk water, as well as 51 snapshots for each ion at the air/water interface, with individual snapshots separated by 10 ps in time. These ensembles were taken from our previous work [7], and coordinate files for these data sets are provided in the Supplementary Material. The bulk ensemble represented by these 51 snapshots affords statistical distributions that are indistinguishable from results obtained from an isotropic simulation, and the interfacial ensemble affords similar distributions regardless of whether the interface is defined by GDS − 3 Å (as in the present work) versus GDS − 1 Å or GDS − 5 Å [7].

Symmetry-Adapted Perturbation Theory
Quantum-mechanical values of E int were computed using SAPT based on Hartree-Fock (HF) wave functions for the monomers and second-order perturbation theory for the intermolecular Coulomb operators. This approach is usually called SAPT0 [51,58,59] and is closely related to second-order Møller-Plesset perturbation theory (MP2). However, because second-order dispersion is far from quantitative [50,51,60], we replace it in these calculations with a many-body dispersion (MBD) model [51,61,62], in what we have termed a "hybrid" or "extended" form of SAPT [51]. This hybrid method will be designated as SAPT0 + MBD. At this level of theory, results for small-molecule data sets suggest that errors in E int are within ∼1 kcal/mol of the best-available benchmarks [59,62], provided that adequate basis sets are employed [59,63]. All electronic structure calculations were performed using Q-Chem software, v. 5.4 [64].
The interaction energy computed using SAPT0 + MBD is naturally partitioned as [50,51] The energy components [51,65] include electrostatics (E elst ), meaning the Coulomb interaction between isolated-monomer charge densities; exchange or Pauli repulsion (E exch ), which is the penalty to antisymmetrize the isolated-monomer wave functions; induction (E ind ), which includes both polarization and charge transfer (CT); and, finally, dispersion (E disp ). In our approach, and are the first-order SAPT electrostatic and exchange energies, while E disp is the dispersion energy computed using the MBD model [62]. The induction energy comes from secondorder SAPT but warrants additional discussion, which we defer until Section 2.3. Previous basis-set testing of SAPT0 + MBD reveals that polarized triple-ζ basis sets, augmented with diffuse functions, are both necessary and sufficient to obtain converged energetics [59,63]. This is a unique feature of our hybrid approach to SAPT [51], which replaces the very slow basis-set convergence of perturbative dispersion with a model (namely, MBD) that converges quickly [63]. Tests for Cl − (aq) in Figure 1 demonstrate that interaction energies computed using the 6-311+G(d,p) basis set agree with SAPT0 + MBD/ def2-TZVPD values to within an average of 2.0 kcal/mol, in a total interaction energy that averages −106 kcal/mol. Relative to the more complete def2-TZVPD basis set, the Pople basis set systematically underestimates E ind (by an average of 1.6 kcal/mol) and overestimates E elst (by an average of 4.2 kcal/mol), whereas E exch and E disp are nearly identical in both basis sets. More important than these relatively small differences is the fact that instantaneous values of E int fluctuate from snapshot to snapshot in a similar way in either basis set. For these calculations, which involve Cl − (H 2 O) n with an average of n = 28 water molecules, SAPT0 + MBD/6-311+G(d,p) calculations are 17× faster than the corresponding calculations with def2-TZVPD. (This speedup results largely from the absence of diffuse functions on hydrogen but also benefits from Q-Chem's very efficient handling of sp shells in Pople basis sets.) In the present work, we are concerned with comparisons between bulk and interfacial behavior rather than absolute interaction energies, and the need for ensemble averaging requires high throughput. As such, 6-311+G(d,p) is used for all subsequent SAPT calculations.
Interaction energies defined in Equation (1) do not include relaxation of the monomer geometries, so E int is an interaction energy in the "vertical" sense, not a binding energy or a solvation energy. In considering the ion-water clusters X − (H 2 O) n extracted from MD simulations, we treat the entire water cluster (H 2 O) n as a single monomer for the purpose of computing E int and its components, then average over the ensemble of snapshots. Even so, the ensemble-averaged value E int corresponds to vertical removal of the ion. It includes the change in electronic polarization of the water molecules upon removal of the ion but does not include the orientational reorganization energy of the water to fill the void left behind by the ion.
Unless otherwise specified, all of the SAPT0 calculations reported herein use HF wave functions for the monomers. However, we will also report a few SAPT0(KS) calculations [51,59], in which Kohn-Sham (KS) molecular orbitals from density functional theory (DFT) are used in place of HF orbitals. These SAPT0(KS) calculations employ the long-range corrected (LRC) density functional LRC-ωPBE [66]. Previous work has emphasized the importance of using an asymptotically correct exchange potential in SAPT calculations [59,60,67,68], and this condition can be achieved in practice via monomer-specific tuning of the rangeseparation parameter (ω) in the LRC-ωPBE functional. Although "optimal tuning" of LRC functionals [69,70] is sometimes accomplished using the ionization energy (IE) theorem of DFT, a more robust procedure in the present context is the "global density-dependent" (GDD) or "ω GDD " procedure [59,60,68]. This approach, which adjusts ω based on the size of the exchange hole, mitigates the strong dependence on system size that is observed when using IE tuning [59], which might otherwise be a problem when studying water clusters of varying size [71]. For water, we use ω = 0.277a −1 0 , which represents an average over several cluster geometries. For the ions, we tune ω individually at the optimized gasphase geometry of each, resulting in a range of values from ω = 0.248a −1 0 for iodide and ω = 0.261a −1 0 for bromide, where the tails of the anion's density are most diffuse, up to ω = 0.398a −1 0 for cyanate and ω = 0.405a −1 0 for cyanide, where the density is most compact. (Note that LRC functionals switch from semilocal exchange to HF exchange on a length scale of ∼1/ω.) In previous work, we have often used self-consistent charge embedding of the SCF monomer wave functions as a means to incorporate many-body polarization effects into a pairwise SAPT calculation, albeit implicitly [50,[72][73][74][75]. However, the present study does not make use of any charge embedding, and instead the X − (H 2 O) n system is treated as a dimer with (H 2 O) n as one monomer. In principle, charge embedding could be used to describe these clusters more efficiently as (n + 1)-body systems, but we have chosen not to do so here. The dimer approach makes the SAPT interaction energies more directly comparable to those obtained using the AMOEBA force field.

Polarization and Charge Transfer
In our calculations, the induction term in Equation (1) is defined as The first two terms are the second-order (SAPT0) induction and exchange-induction energies, and is the so-called "δHF" correction [51]. It uses a counterpoise-corrected, supramolecular HF interaction energy (∆E HF int ) to correct the SAPT0 interaction energy for induction effects beyond second order in perturbation theory, which is crucial for the accurate description of hydrogen bonds [51,59]. See reference [76] for a definition of the second-order response ("resp") energies that appear in Equation (6).
As defined in SAPT, the induction energy in Equation (5) contains both polarization and CT, for reasons that are discussed in reference [77]. In the analysis of hydrogen bonding it is often of interest to separate these effects but that separation has historically been considered problematic. The dilemma is not limited to SAPT and many schemes for separating polarization from CT exhibit strong dependence on the choice of basis set [77].
To accomplish the separation in Equation (7) in a robust way that converges rapidly with respect to basis set, we use the machinery of a charge-constrained self-consistent field (SCF) calculation [78] to define a CT-free reference state. Here, the monomers are allowed to polarize one another but their charge densities are constrained to integrate to integer numbers of electrons. Because the SCF procedure is variational, lifting of this constraint necessarily lowers the energy (to that of the fully-relaxed SCF solution), and this energy lowering is taken to define E CT . The CT energy thus obtained is then subtracted from the SAPT induction energy to obtain the CT-free polarization energy, [77,[79][80][81]. CT energies defined in this way are very nearly converged already in double-ζ basis sets [77]. This approach has previously been used to demonstrate that E CT furnishes a driving force for formation of quasi-linear hydrogen bonds in binary halide-water complexes [65,81]. Implementation of the charge-constrained SCF procedure requires a method to count electrons, and Becke's multicenter partition scheme [82] is commonly used for this purpose [78]. Becke's approach first divides space into Voronoi cells [83], which are regions of space that are closest to a particular nucleus, and then applies a smoothing function at the boundaries of these polyhedra. Alternatively, and specifically for the purpose of defining a CT-free reference state in order to affect the partition suggested in Equation (7), a counting procedure based on fragment-based Hirshfeld (FBH) weighting has also been suggested [79,81]. In the latter approach, the number of electrons contained in fragment A is defined as where ρ(r) is the supramolecular electron density, which is integrated subject to a weighting function w A (r). That function is defined as where ρ 0 X (r) is the charge density of isolated fragment X. The denominator in Equation (9) is thus a superposition of isolated-fragment densities.
The Becke scheme can also be conceptualized as a form of Equation (8) in which w A (r) is a smoothed version of a Heaviside step function, which switches rapidly between w A (r) = 0 and w A (r) = 1 at the boundaries of the Voronoi polyhedra. In practice, our implementation of Becke's procedure uses "atomic size adjustments" [82], in which a set of empirical atomic radii [84] are used to adjust the boundaries of the Voronoi cells away from the midpoints of the internuclear vectors. As discussed below, this adjustment is crucial for systems with substantial size mismatch between nearby atoms.
Even so, the FBH approach strikes us as the more reasonable one, especially where anions are involved, because Becke's method depends only on the positions of the atoms (along with the empirical atomic radii), whereas the weight function defined in Equation (9) respects the diffuseness of the isolated anion's wave function. In the present context, this almost inevitably means that the extent of anion → water CT is smaller when the FBH approach is used, because the tails of the X − wave function cause a larger region of space to contribute to that fragment's integrated number of electrons. As an example, Figure 2 presents E CT computed using both Becke partition (with atomic size adjustments) and FBH weighting, for each snapshot of I − (aq) in bulk water. The results are considerably different depending on which method is used to count electrons, with the FBH approach compressing the CT energy into the interval 0 > E CT > −2 kcal/mol, whereas the Becke procedure affords values of |E CT | as large as 20 kcal/mol. The latter value is comparable to the the average magnitude of the total SAPT0 induction energy, which is E ind = −22.3 kcal/mol for I − (aq) in bulk water. (Note that energy components corresponding to attractive interactions are negative.) CT energies for snapshots of I − (aq) in bulk water, computed using a charge-constrained SCF procedure with the charge constraint defined using either fragment-based Hirshfeld (FBH) weighting (scale at left), or else Becke's multicenter partitioning procedure (scale at right). Results using the Becke scheme include the "atomic size adjustments" that are described in reference [82], wherein Slater's set of atomic radii [84] are used to adjust the boundaries of the Voronoi cells based on atomic size. Figure 3 shows the polarization energy (E pol = E ind − E CT ) that is obtained using either the Becke or the FBH weighting function to define the charge constraint. (Both definitions of E pol start from the same SAPT0 induction energy, E ind .) It is apparent that the two definitions afford step-to-step fluctuations that do not seem to correlate with one another. In the Becke definition, the size and shape of the Voronoi cell that contains iodide is sensitive to the instantaneous values of all iodide-water distances in the first solvation shell, whereas the FBH definition uses a spherically-symmetric charge density for the isolated anion in order to define the charge constraint. The latter definition proves to be less sensitive to fluctuations in the atomic coordinates, although it remains sensitive to the presence of hydrogen bonds [65,81].
For I − (aq), the CT-free reference state defined using Becke partition consistently results in CT energies that are larger in magnitude as compared to the FBH scheme: |E CT (Becke)| > |E CT (FBH)|. This is evident from the rather different energy scales in Figure 2, but the situation is not the same for all of the anions examined here. As a second example we consider ClO − (aq), which exhibits the largest values of |E CT | amongst the anions in our data set, at least when the FBH definition is used. Figure 4 considers both definitions and examines how E CT fluctuates from snapshot to snapshot. Becke's partition predicts very little CT for ClO − (aq) in bulk water ( E CT = −1.2 kcal/mol), whereas the FBH definition results in an average value E CT = −6.2 kcal/mol. In either case, E CT is consistently larger for the interfacial snapshots.    We will use the FBH definition of E CT for the remainder of this work. Our main interest lies in understanding how various energy components compare when the ion is in a bulk versus an interfacial environment, but the magnitude of E CT can depend strongly on the method that is used to count electrons, as noted above. This observation suggests that in other applications of constrained DFT [78], which is the more common form of chargeconstrained SCF calculation (in contrast to the constrained HF calculations employed here), the results should be checked carefully to ensure that conclusions are robust with respect to the details of how the constraints are implemented.
The SG-3 quadrature grid [85] is used to integrate the SCF constraint equations as well as Equation (8). As a technical aside, we note that the atomic size adjustments mentioned above are crucial in order to obtain results that are even remotely sensible when Becke partition is used to implement the charge constraint. However, the original implementation of constrained DFT in the Q-Chem program did not include these corrections [86], which were added recently for the purpose of SAPT-based CT analysis [81]. Absent these corrections, the Voronoi cell boundaries are placed at midpoints of the internuclear vectors, which affords unreasonable results in cases where neighboring atoms have very different size. This includes covalent bonds to hydrogen, where the midpoint definition causes too much density to be assigned to the smaller hydrogen atom, often leading to a negative charge assigned to hydrogen [81]. In the present work, neglecting the atomic size corrections leads to a significant fraction of the iodide's charge being assigned to first-shell water molecules, resulting in completely unrealistic CT energies whose magnitudes exceed the total SAPT0 induction energy. In our view, constrained DFT based on Becke partition should not be used without the atomic size corrections.  Figure 5a) and also using the AMOEBA force field (Figure 5b), where the latter is the same force field that was used for the simulations from which these X − (H 2 O) n structures were extracted. Bulk and interfacial data are averaged separately, with the criterion GDS − 3 Å used to decide whether a particular snapshot represents a bulk or an interfacial solvation environment. There are two interesting observations to be made from the interaction energy data in Figure 5. Foremost is the fact that differences between the bulk and interfacial mean values E int for a given ion are small compared to the fluctuations in the instantaneous value of E int . Bulk values of E int are systematically (slightly) larger in magnitude than interfacial values, except for CN − , OCN − , and NO − 3 where the averages are essentially identical in both environments. In all cases, the difference between bulk and interfacial average values is well within the standard deviation in either quantity; see the numerical values that are provided in Table 1. For the halides, modest reductions in E int at the interface (up to 7-8 kcal/mol for bromide and iodide) are consistent with results from classical MD simulations indicating that the average ion-water interaction is reduced, for all of the halides, as the ion moves towards the interface [21]. It should be noted that the simulations reported in reference [21] indicate that the enthalpic portion of the potential of mean force is more favorable for the heavier halides at the interface, as compared to its value in bulk water. As such, the rather subtle differences between ion-water interactions that are documented in our quantum-mechanical calculations are more than compensated by ion-induced changes in the water-water interactions [21]. This is consistent with our detailed structural analysis of the ions [7], which indicates very little change in the first-shell structure at the interface as compared to that in bulk water. A second interesting observation is the generally strong correlation between classical (AMOEBA) and quantum-mechanical (SAPT) values of E int , even if the former are systematically smaller than the latter, e.g., by 15-19 kcal/mol for the halide ions. (Systematic differences are smaller for the other ions, except in the case of ClO − 3 . The latter is discussed in detail below.) For the halide ions, we use AMOEBA parameters that were originally developed by Ponder and co-workers [55], and we note that the discrepancies between the force field and the quantum chemistry that are documented in Figure 5 are much larger than those reported previously for binary X − (H 2 O) complexes [55]. This underscores the importance of considering larger ion-water clusters, given the many-body nature of polarization in aqueous systems [87][88][89][90][91][92]. Simulation of the hydration free energy of chloride using AMOEBA results in an error of 11.9 kcal/mol with respect to experiment [55], assuming that the reference value is defined using the proton solvation energy of Tissandier et al. [93], which has since emerged as the consensus value [94][95][96]. In view of this, the systematic difference of 17 kcal/mol between AMOEBA and SAPT0 + MBD values of E int in bulk water (see Table 1) is not so dissimilar from previous results. Improvements to the AMOEBA force field for ions, using SAPT energy components as benchmark data, is a topic of contemporary interest [97][98][99].

Results and Discussion
The chlorate (ClO − 3 ) ion represents the lone exception to an otherwise systematic correlation between classical and quantum-chemical interaction energies. This particular species is much more strongly solvated by AMOEBA ( E int = −126.6 ± 9.9 kcal/mol in bulk water) than it is by SAPT0 + MBD ( E int = −85.8 ± 10.0 kcal/mol). Considering the chlorine oxyanions as a group, the trend amongst the AMOEBA values of | E int | is The fact that perchlorate (ClO − 4 ) is an outlier is easy to rationalize in terms of its tetrahedral symmetry and vanishing dipole moment, but the trend amongst the other three chlorine oxyanions is more puzzling. Ensemble-averaged SAPT0 + MBD energy components for the four species ClO − n (aq) are listed in Table 2, from which it may be seen that E int , E elst , and E ind all follow the same trend exhibited by the gas-phase dipole moments of the ions in question. However, this means that the trend amongst total interaction energies is different from that predicted by AMOEBA. Instead, from the quantum-mechanical calculations the trend (from strongly to weakly interacting)  In contrast to the AMOEBA results, the SAPT0 + MBD calculations afford similar ensemble-averaged interaction energies for both ClO − 3 and ClO − 4 . Given that all of the chlorine oxyanions except for ClO − 4 have sizable dipole moments, the ClO − 3 interaction energy seems anomalously small when computed using SAPT0 + MBD. As a sanity check, we recomputed interaction energies for all of the ions using SAPT0(KS) + MBD, which includes intramolecular electron correlation. These results are plotted in Figure 6a, which should be compared to the corresponding SAPT0 + MBD results in Figure 5a. Total interaction energies at either level of theory are quite comparable, and in particular both methods exhibit the same (seemingly anomalous) trend amongst the ClO − n ions, which differs from the trend predicted by AMOEBA. To investigate this further, we examine the SAPT0 + MBD energy components. These are plotted for each of the ions in Figure 7, again separating bulk and interfacial environments and ensemble-averaging over either data set. In considering the energy decomposition in Equation (1), we have opted to group first-order electrostatics and exchange together, because their sum approximates the electrostatic interaction between antisymmetrized monomer wave functions. This combination of "primitive" electrostatics (E elst , which is the Coulomb interaction between isolated-monomer charge densities) and Pauli repulsion (E exch ) has proven to be easier to interpret for halide-water systems as compared to electrostatics alone [65,81], in part because E elst and E exch are the largest energy components (in magnitude) but often have opposite signs, such that their sum is more comparable in magnitude to the remaining energy components. An example can be found in the ensemble-averaged energy components for the ClO − n (aq) species (Table 2), where the much less repulsive value of E exch for perchlorate at first seems at odds with the larger size of this ion. However, the reduced Pauli repulsion in this case is actually commensurate with a much less attractive value of E elst , suggesting a hydration sphere that is not as tight around the ion as it is in the smaller (but electrostatically much more attractive) ClO − , ClO − 2 , and ClO − 3 ions. Statistical distributions of E elst+exch are shown in Figure 7a for all of the ions, and ClO − 3 immediately stands out as the only one for which E elst+exch > 0. In other words, the sum of first-order interactions is net repulsive for ClO − 3 but is net attractive for each of the other anions. This observation is independent of whether one considers the bulk or interfacial data set because differences between bulk and interfacial mean values of E elst+exch are tiny in comparison to the instantaneous fluctuations, as was the case for E int . Furthermore, this anomalous prediction regarding ClO − 3 is not unique to the SAPT0 level of theory that is used in Figure 7. A similar anomaly is evident in the SAPT0(KS) results, which can been seen from the statistical distributions of E elst+exch at that level of theory that are plotted in Figure 6b. We note that the largest values of E exch often correspond to the largest (most attractive) total interaction energies, as is seen for example in SAPT calculations of ClO − n · · · C 6 H 6 complexes (n = 1, 2, 3, 4) [100]. In the present case, ClO − 3 bucks this trend, according to the energy components listed in Table 2.
A possible explanation for the apparently anomalous behavior of ClO − 3 can be found by examining radial distribution functions (RDFs), g(r), obtained from the MD simulations. (These can be found in the Supporting Information for reference [7] but the salient details are described here.) Amongst the chlorine oxyanions, a unique feature of ClO − 3 is the appearance of two distinct peaks in the RDF for Cl· · · O w (where O w denotes water oxygen), one at r ≈ 3.5 Å and another at r ≈ 4.1 Å. For each of the other ClO − n species, the RDF consists of a single well-resolved feature at r ≈ 3.5-3.7 Å. The shorter-r feature for ClO − 3 does not appear to be present in simulations based on a hybrid quantum mechanics/ molecular mechanics (QM/MM) formalism, which were used to interpret x-ray scattering results [101]. If the small-r feature for ClO − 3 is an indication of an extraneous water molecule present at short range, then this could explain the anomalously repulsive values of E elst+exch that we then compute using quantum mechanics applied to snapshots extracted from the classical MD simulations. The presence of such a water molecule in those simulations, however, suggests that something in AMOEBA's ion-water interaction is compensating for the short-range repulsion, or perhaps that the latter is simply not repulsive enough. Although polyvalent anions are not considered in the present work (because they are excluded from the air/water interface), it is notable that a short-r peak in the S· · · O w RDF is also observed in our previous simulations of SO 2− 3 (aq) [7]. That feature is absent from QM/MM simulations and x-ray scattering experiments [102]. In view of this, AMOEBA parameterizations for both of these ions ought to be revisited. This is beyond the scope of the present work, though it is interesting to note the way that SAPT analysis of ion-water clusters was able to detect an anomaly. Notably, vertical ionization energies computed for ClO − 3 (aq) and SO 2− 3 (aq) based on these simulations are no less accurate, as compared to experimental values [49], than what we obtain for other inorganic anions, including other ClO − n ions [7]. The typical accuracy is ∼0.2 eV [7], considerably smaller than the widths of the corresponding photoelectron spectra.
Returning exclusively to the monovalent ions and examining the other energy components whose statistics are summarized in Figure 7, another curiosity arises in regard to dispersion energies for the chlorine oxyanions. Dispersion is size-extensive, so that all else being equal it should scale in proportion to the number of electrons. For the ClO − n species, however, we observe that |E disp | decreases in the order ClO − 3 > ClO − 2 > ClO − > ClO − 4 . This time, perchlorate is the apparent anomaly. Dispersion energies in Figure 7c were computed using the MBD model [62], so as a sanity check, we recomputed E disp using the third-generation ab initio dispersion potential aiD3 [50], which consists of atom-atom C 6 and C 8 potentials fitted to dispersion-only data from high-level SAPT calculations. Dispersion energies obtained for the ClO − n (aq) species with either MBD or aiD3 are provided in Table 3 in the form of ensemble averages. Both models afford rather similar dispersion energies, consistent with previous tests for cases where many-body effects are not significant [62]. (In the context of dispersion, "many-body" implies an effect that cannot be described by pairwise atom-atom potentials [60,103]. Many-body dispersion effects typically arise in conjugated molecules where screening significantly modifies the effective C 6 coefficients [104]. For small molecules, three-body dispersion effects are typically quite small [90].) Notably, in the aiD3 model the C 6 and C 8 coefficients depend only on atomic number and do not respond to the electronic structure of the monomers.
The sharp drop in dispersion between chlorate (ClO − 3 ) and perchlorate is a feature of both the MBD and aiD3 dispersion models, suggesting that this is not an artifact. A likely explanation is that, in perchlorate, the addition of a fourth oxygen atom around the central (and more polarizable) chlorine atom screens the water molecules from this polarizable center and thus significantly attenuates chlorine's contribution to the dispersion energy. In contrast, for the other ClO − n ions the chlorine atom remains solvent-exposed, and the dispersion is much larger. This mechanism would be reflected in both dispersion models, if only as a function of increased chlorine-water distance in the aiD3 case. Also in support of this hypothesis are the data in Figure 7b for SAPT0 + MBD induction energies, which also exhibit a pronounced drop in magnitude between ClO − 3 and ClO − 4 . As compared to dispersion interactions, polarization effects decay more slowly with distance, e.g., as r −4 for charge-dipole polarization, but this dependence is still rather steep. Polarization is often invoked in discussions of ions at the air/water interface [9][10][11][12][13][14][15][16][17], so it is interesting to note that induction energies are systematically smaller in the interfacial environment; see Figure 7b. As with the total interaction energies, however, the difference between bulk and interfacial mean values E ind is small in comparison to the instantaneous fluctuations as measured by the standard deviation. (Numerical data corresponding to Figure 7b are provided in Table 1.) Note that "polarization" as it is typically understood means strictly intramolecular redistribution of charge, with CT considered as a separate effect; these two parts of the induction energy are separated in Figure 8. Because the CTfree polarization energy (E pol ) is much larger than the CT energy (E CT ), the result is that E pol follows essentially the same trend from ion to ion as does the total induction energy, E ind . In particular, this means that the polarization energy is systematically smaller in the interfacial environment, for each of the ions considered here. Indeed, for the canonical surface-active anions Br − , I − , ClO − 4 , and SCN − [19,34,35,105], the polarization energy is significantly smaller in the interfacial environment, by at least the standard deviation of E pol in bulk water; see Figure 8a. That observation, in turn, is a direct result of CT energies that are systematically larger at the interface for precisely those four surface-active anions. Statistical distributions of E CT for all of the ions are plotted in Figure 8b. In contrast to other energy components, only for E CT do we observe a pronounced difference between averages computed for the bulk and interfacial data sets. That said, the overall scale of the CT energies is a rather small part of either the total induction energy or the total interaction energy, with |E CT | 10 kcal/mol except in the case of interfacial ClO − . (Although CT energies smaller than 10 kcal/mol do play a pivotal role in establishing the directionality of hydrogen bonds [65,81], that kind of detailed analysis of a potential energy surface is not attempted in the present work, where we are interested in ensemble-averaged properties.) For Br − , I − , ClO − 4 , and SCN − , the average CT energy at the air/water interface is larger than its mean value in bulk water by at least one standard deviation in the bulk value. For Cl − (aq), the interfacial average value of E CT is larger in magnitude than the bulk value, though not quite by a full standard deviation. It is perhaps noteworthy that outliers for the CT energies tend to be larger at the interface, particularly towards negative (more stabilizing) values of E CT .
In the context of the Hofmeister series [106,107], the anions I − , ClO − 4 , and SCN − have especially large binding constants to protein [107,108], which is historically associated with the definition of chaotropes or "structure breakers" [8]. (Note that Ricci and co-workers [109] point out that the kosmotrope and chaotrope or "structure maker" and "structure breaker" labels are largely thermodynamic in origin and should not be taken too seriously in terms of their implications for microscopic hydrogen-bonding structure). In comparison to the aforementioned ions, Cl − binds to proteins more weakly [108]. That said, NO − 3 is usually categorized as a structure-breaker on par with Br − in the Hofmeister series [106], and as weakly surface-active on the basis of VSFG measurements [19], yet the mean values of E CT that we obtain for NO − 3 are essentially identical in the bulk and interfacial environments, albeit with larger outliers in the interfacial case. The hypochlorite ion (ClO − ) stands out in our analysis, with a significantly larger mean value of |E CT | in the interfacial environment. This species is not typically discussed in the context of the Hofmeister series or in VSFG studies of the air/water interface, due to its limited stability in aqueous solution.

Conclusions
Detailed analysis of anion-water clusters extracted from MD simulations reveals that the total ion-water interaction energy (considering two solvation shells around the ion) is systematically larger for a given ion in bulk water than it is for the same ion near the air/water interface. The same is true for the CT-free polarization component of the total interaction energy, which is interesting given that polarization is often assumed to play a central role in surface activity [13], although this contention is disputed [23,24]. In any case, we observe systematically larger polarization energies in bulk water for both the "soft" anions with low surface charge density that are usually considered to be surface active (Br − , I − , ClO − 4 , and SCN − ), as well as for hard anions that are not considered to be surface active (CN − , OCN − , and NO − 2 ). That said, systematic differences in the mean values E int and E pol in bulk versus interfacial environments are rather small in comparison to the magnitude of the instantaneous fluctuations in E int and E pol .
Anion-to-water CT stands out as the only energy component whose magnitude is larger at the air/water interface for some of the ions. In fact, it is larger specifically for the traditional surface-active anions: Br − , I − , ClO − 4 , and SCN − . However, NO − 3 can also be detected in surface-sensitive vibrational spectroscopy [19], yet for that species E CT is essentially the same at the interface as it is in bulk water. The Cl − ion is a borderline case whose average CT energy is slightly more stabilizing at the interface, albeit by less than one standard deviation in the fluctuations. In all cases, the CT energy constitutes less than 20% of the total induction energy, meaning that it is at least 5× smaller than the CT-free polarization energy, the latter of which does not exhibit a surface preference and is in fact larger in bulk water. Nevertheless, the consequences of this "excess" CT for soft anions at the air/water interface seem worth considering in future work, especially in the context of VSFG experiments. Intermolecular CT mechanisms have been invoked in the past to explain the surface charge of liquid water that is inferred from electrophoretic measurements [110][111][112][113].
Considering the halide ions as a series that ranges from kosmotropic to chaotropic [8], or equivalently whose surface activities decrease in the order I − > Br − > Cl − F − , it has previously been noted that no single mechanistic explanation for this ordering can be gleaned from atomistic simulations [21,24]. Changes in the water-water interactions as the the anion approaches the interface appear to play a role [21]. The present analysis, based on accurate quantum-mechanical calculations of ion-water interaction energies, supports the notion that ion-water interactions alone do not readily afford any kind of a diagnostic (let alone a mechanism) to determine whether an ion resides in a bulk or interfacial environment. This null result complements our recent conclusion that short-range (firstshell) solvation structure is extremely similar in the bulk and interfacial environments [7]. The detailed mechanism of soft anion surface activity remains an open question.

Data Availability Statement:
The data that support this study are available from the corresponding author upon reasonable request.