molecules

Article
Interaction Energy Analysis of Monovalent Inorganic Anions in
Bulk Water Versus Air/Water Interface
John M. Herbert * and Suranjan K. Paul

Department of Chemistry & Biochemistry, The Ohio State University, Columbus, OH 43210, USA;
paul.926@osu.edu
* Correspondence: herbert@chemistry.ohio-state.edu

Abstract: 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 ClO− n
(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
Citation: Herbert, J.M.; Paul, S.K. surface activity. However, CT stabilization amounts to <20% of the total induction energy for each of
Interaction Energy Analysis of the ions considered here, and CT-free polarization energies are systematically larger in bulk water in
Monovalent Inorganic Anions in
all cases. As such, the role of these effects in the surface activity of soft anions remains unclear. This
Bulk Water Versus Air/Water
analysis complements our recent work suggesting that the short-range solvation structure around
Interface. Molecules 2022, 26, 6719.
these ions is scarcely different at the air/water interface from what it is in bulk water. Together, these
https://doi.org/10.3390/
observations suggest that changes in first-shell hydration structure around soft anions cannot explain
molecules26216719
observed surface activities.
Academic Editor: Eric Glendening
Keywords: air–water interface; Hofmeister series; hydrogen bonding; charge transfer; symmetry-
Received: 13 October 2021 adapted perturbation theory; noncovalent interactions
Accepted: 3 November 2021
Published: 6 November 2021

Publisher’s Note: MDPI stays neutral 1. Introduction

with regard to jurisdictional claims in One of the earliest results of surface-sensitive vibrational sum-frequency generation
published maps and institutional affil-
(VSFG) experiments [1,2] was the observation that soft anions impact the vibrational
iations.
lineshape in the O–H stretching region, but that hard anions do not [3–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
Copyright: © 2021 by the authors. anions is often discussed in terms of polarizability [9–17], it should be noted that polyvalent
Licensee MDPI, Basel, Switzerland. anions such as SO24− (aq) are quite polarizable [18] but the presence of polyvalent anions in
This article is an open access article solution does not affect the O–H lineshape measured in VSFG experiments [19]. Molecular
distributed under the terms and dynamics (MD) simulations suggest that hard anions, including polyvalent species but
conditions of the Creative Commons
also fluoride, are repelled from the air/water interface [20,21].
Attribution (CC BY) license (https://
The nature of the surface activity exhibited by soft anions remains a matter of debate.
creativecommons.org/licenses/by/
Whereas continuum electrostatics predicts that all ions are repelled from the air/water
4.0/).

Molecules 2022, 26, 6719. https://doi.org/10.3390/molecules26216719 https://www.mdpi.com/journal/molecules

Molecules 2022, 26, 6719 2 of 20

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–29], which aligns with the interpreta-
tion 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 in-
terface as compared to that in bulk water [7]. This observation comes from MD simula-
tions 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. Io-
dide 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, per-
haps 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 femtosec-
ond vibrational pump–probe experiments [37–40], although the effects on the long-range
hydrogen-bonding structure of water are less clear. Neutron diffraction experiments [41–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–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− , NO2− , NO3− , and ClO− n (n = 1, 2, 3, 4). Some of these
are typical surface-active ions (e.g., Br− , I− , SCN− , and ClO4− ), whereas others (such as
CN− , OCN− , and NO2− ) 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
NO3− [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–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 environ-
ment engenders any discernible changes in the ion–water interactions, relative to what is
observed for the same ion in bulk water.

2. Computational Methods
2.1. 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/cm3 . The size of the periodic simulation cell
Molecules 2022, 26, 6719 3 of 20

(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 h Nw i
depending on both the size of the ion and how tightly hydrated it is. In bulk water,
these averages range from h Nw i ≈ 28 for Cl− (aq) up to h Nw i ≈ 43–44 for Br− (aq) and
I− (aq), with h Nw i = 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 (Eint ) between the
ion and its first two hydration shells. The quantity Eint 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].

2.2. Symmetry-Adapted Perturbation Theory

Quantum-mechanical values of Eint 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 Eint 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]

Eint = Eelst + Eexch + Eind + Edisp . (1)

The energy components [51,65] include electrostatics (Eelst ), meaning the Coulomb inter-
action between isolated-monomer charge densities; exchange or Pauli repulsion (Eexch ),
which is the penalty to antisymmetrize the isolated-monomer wave functions; induction
Molecules 2022, 26, 6719 4 of 20

(Eind ), which includes both polarization and charge transfer (CT); and, finally, dispersion
(Edisp ). In our approach,
(1)
Eelst ≡ Eelst (2)
and
(1)
Eexch ≡ Eexch (3)
are the first-order SAPT electrostatic and exchange energies, while Edisp is the dispersion
energy computed using the MBD model [62]. The induction energy comes from second-
order 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 Eind (by an average of 1.6 kcal/mol) and
overestimates Eelst (by an average of 4.2 kcal/mol), whereas Eexch and Edisp are nearly
identical in both basis sets.

–80
6-311+G(d,p)
–85
def2-TZVPD
–90

–95
Eint (kcal/mol)

–100

–105

–110

–115

–120

–125
0 10 20 30 40 50
snapshot index

Figure 1. Total interaction energies for snapshots of Cl− (aq) in bulk water, computed at the
SAPT0 + MBD level using two different basis sets. Solid horizontal lines show the ensemble-averaged
values obtained using either basis set. These averages are h Eint i = −106.6 ± 8.4 kcal/mol for
6-311+G(d,p) and h Eint i = −104.6 ± 8.2 kcal/mol for def2-TZVPD, where the uncertainties represent
one standard deviation.

More important than these relatively small differences is the fact that instantaneous
values of Eint fluctuate from snapshot to snapshot in a similar way in either basis set. For
these calculations, which involve Cl− (H2 O)n with an average of n = 28 water molecules,
SAPT0 + MBD/6-311+G(d,p) calculations are 17× faster than the corresponding calcula-
tions 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.
Molecules 2022, 26, 6719 5 of 20

Interaction energies defined in Equation (1) do not include relaxation of the monomer
geometries, so Eint is an interaction energy in the “vertical” sense, not a binding energy or
a solvation energy. In considering the ion–water clusters X− (H2 O)n extracted from MD
simulations, we treat the entire water cluster (H2 O)n as a single monomer for the purpose of
computing Eint and its components, then average over the ensemble of snapshots. Even so,
the ensemble-averaged value h Eint i 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 range-
separation 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,
IE = −ε HOMO , (4)
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.277a0−1 , which represents an average over
several cluster geometries. For the ions, we tune ω individually at the optimized gas-
phase geometry of each, resulting in a range of values from ω = 0.248a0−1 for iodide and
ω = 0.261a0−1 for bromide, where the tails of the anion’s density are most diffuse, up
to ω = 0.398a0−1 for cyanate and ω = 0.405a0−1 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–75]. However, the present study does
not make use of any charge embedding, and instead the X− (H2 O)n system is treated as
a dimer with (H2 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.

2.3. Polarization and Charge Transfer

In our calculations, the induction term in Equation (1) is defined as

(2) (2)
Eind = Eind + Eexch-ind + δEHF . (5)

The first two terms are the second-order (SAPT0) induction and exchange-induction ener-
gies, and
HF (1) (1) (2) (2)
δEHF = ∆Eint


− Eelst + Eexch + Eind,resp + Eexch-ind,resp (6)

is the so-called “δHF” correction [51]. It uses a counterpoise-corrected, supramolecular

HF ) to correct the SAPT0 interaction energy for induction effects
HF interaction energy (∆Eint
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).
Molecules 2022, 26, 6719 6 of 20

As defined in SAPT, the induction energy in Equation (5) contains both polarization
and CT,
Eind = Epol + ECT , (7)
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 ECT . The CT energy thus obtained is then
subtracted from the SAPT induction energy to obtain the CT-free polarization energy,
Epol = Eind − ECT [77,79–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 ECT 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 pur-
pose [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 Z
NA = w A (r) ρ(r) dr, (8)

where ρ(r) is the supramolecular electron density, which is integrated subject to a weighting
function w A (r). That function is defined as

ρ0A (r)
w A (r) = , (9)
∑ B ρ0B (r)

where ρ0X (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 ECT 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
Molecules 2022, 26, 6719 7 of 20

different depending on which method is used to count electrons, with the FBH approach
compressing the CT energy into the interval 0 > ECT > −2 kcal/mol, whereas the
Becke procedure affords values of | ECT | as large as 20 kcal/mol. The latter value is
comparable to the the average magnitude of the total SAPT0 induction energy, which
is h Eind i = −22.3 kcal/mol for I− (aq) in bulk water. (Note that energy components
corresponding to attractive interactions are negative.)

0.0 –2
FBH
–0.2 Becke –4

–0.4
–6
–0.6

ECT, Becke (kcal/mol)

–8
ECT, FBH (kcal/mol)

–0.8
–10
–1.0
–12
–1.2
–14
–1.4
–16
–1.6

–1.8 –18

–2.0 –20
0 10 20 30 40 50
snapshot index

Figure 2. 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 (Epol = Eind − ECT ) that is obtained using
either the Becke or the FBH weighting function to define the charge constraint. (Both
definitions of Epol start from the same SAPT0 induction energy, Eind .) 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:
| ECT (Becke)| > | ECT (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 | ECT | amongst the
anions in our data set, at least when the FBH definition is used. Figure 4 considers both
definitions and examines how ECT fluctuates from snapshot to snapshot. Becke’s partition
predicts very little CT for ClO− (aq) in bulk water (h ECT i = −1.2 kcal/mol), whereas the
FBH definition results in an average value h ECT i = −6.2 kcal/mol. In either case, ECT is
consistently larger for the interfacial snapshots.
Molecules 2022, 26, 6719 8 of 20

–5

–10
Epol (kcal/mol)
–15

–20

–25
FBH
Becke
–30
0 10 20 30 40 50
snapshot index

Figure 3. CT-free polarization energy (Epol = Eind − ECT ) for snapshots of I− (aq) in bulk water. This
quantity is obtained by removing ECT from the total SAPT0 induction energy using either of two
schemes (FBH weighting or Becke partition with atomic size adjustments) to integrate the charge
constraint that defines the CT-free reference state. Solid horizontal lines show the ensemble-averaged
values, which are h Epol i = −21.0 ± 3.9 kcal/mol (FBH) and h Epol i = −12.6 ± 6.1 kcal/mol (Becke).
Uncertainties represent one standard deviation.

–5
ECT (kcal/mol)

–10

–15

–20
FBH, bulk
FBH, interface
–25
Becke, bulk
Becke, interface
–30
0 10 20 30 40 50
snapshot index

Figure 4. CT energies for ClO− (aq), computed using either Becke partition or else FBH weighting to
define the CT-free reference state, and considering both the bulk and interfacial data sets.

We will use the FBH definition of ECT 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 ECT 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 charge-
constrained 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
Molecules 2022, 26, 6719 9 of 20

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 correc-
tions, 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.

3. Results and Discussion

Figure 5 presents ensemble-averaged interaction energies for the sets of X− (H2 O)n
structures that are considered here, where X− is one of 12 monovalent inorganic anions.
Two solvation shells of surrounding water are treated as a single monomer for the pur-
pose of the SAPT calculations. Results are presented both at the quantum-mechanical
SAPT0 + MBD/6-311+G(d,p) level (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− (H2 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.

(a) SAPT0 + MBD (b) AMOEBA

–40 –40
bulk bulk
interface interface
–60 –60

–80 –80
Eint (kcal/mol)

Eint (kcal/mol)
–100 –100

–120 –120

–140 –140

–160 –160
4 —
2 —

3 —

N—

2 —

3 —
—
r—
4 —

N—

N—
2 —

l—
3 —

N—

2 —

3 —
—
r—

I—
N—

N—
l—

I—

lO
lO

lO
lO

lO
C

O
lO

B
lO

lO
C

C
B

SC
C
C

SC

N
C

C
N

O
C
C
C

O
C

Figure 5. Total ion–water interaction energies between various anions and their first two hydration shells, computed using
(a) quantum chemistry at the level of SAPT0 + MBD/6-311+G(d,p), versus (b) the AMOEBA force field that was used to
obtain the structures. Side-by-side box and whisker plots are shown for each ion in bulk and interfacial environments. Box
plots show the mean value h Eint i (yellow bar) and extend for one standard deviation in both directions. Whiskers indicate
minimum and maximum values of Eint .

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 h Eint i for a given ion are small compared to the fluctuations in the instantaneous
value of Eint . Bulk values of h Eint i are systematically (slightly) larger in magnitude than
interfacial values, except for CN− , OCN− , and NO3− 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 h Eint i 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
Molecules 2022, 26, 6719 10 of 20

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.

Table 1. Ensemble-averaged interaction energies h Eint i and induction energies h Eind i, computed at
the SAPT0 + MBD/6-311+G(d,p) level including the δHF correction for induction a .

h Eint i (kcal/mol) h Eind i (kcal/mol)

Ion
Bulk Interface Bulk Interface
Cl− −106.6 ± 8.4 −104.4 ± 9.0 −26.4 ± 2.4 −24.6 ± 2.7
Br− −103.4 ± 9.6 −95.5 ± 9.9 −24.6 ± 2.7 −22.4 ± 3.3
I− −94.1 ± 7.9 −87.5 ± 7.8 −22.3 ± 2.7 −20.3 ± 2.8
ClO− −126.7 ± 9.1 −118.9 ± 8.7 −47.0 ± 5.1 −44.5 ± 4.3
ClO2− −136.2 ± 8.8 −132.2 ± 9.7 −52.3 ± 4.7 −50.2 ± 5.3
ClO3− −85.8 ± 10.0 −83.6 ± 12.5 −42.1 ± 4.1 −41.7 ± 4.4
ClO4− −82.5 ± 7.1 −77.0 ± 9.1 −18.8 ± 1.8 −16.2 ± 2.1
CN− −105.6 ± 8.9 −105.4 ± 8.8 −25.0 ± 2.6 −23.8 ± 2.2
OCN− −110.4 ± 8.2 −110.4 ± 8.0 −29.9 ± 3.0 −29.3 ± 2.9
SCN− −92.2 ± 6.9 −83.5 ± 9.5 −24.6 ± 2.2 −22.4 ± 3.5
NO2− −112.7 ± 8.3 −110.5 ± 10.3 −34.5 ± 3.2 −33.0 ± 3.5
NO3− −106.2 ± 7.6 −106.2 ± 7.9 −29.9 ± 3.2 −29.5 ± 3.3
a Uncertainties represent one standard deviation.

A second interesting observation is the generally strong correlation between classical

(AMOEBA) and quantum-mechanical (SAPT) values of Eint , even if the former are sys-
tematically 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 ClO3− . 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− (H2 O) complexes [55]. This underscores
the importance of considering larger ion–water clusters, given the many-body nature of
polarization in aqueous systems [87–92]. Simulation of the hydration free energy of chlo-
ride 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–96]. In view of this, the
systematic difference of 17 kcal/mol between AMOEBA and SAPT0 + MBD values of h Eint i
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–99].
The chlorate (ClO3− ) 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 (h Eint i = −126.6 ± 9.9 kcal/mol in
bulk water) than it is by SAPT0 + MBD (h Eint i = −85.8 ± 10.0 kcal/mol). Considering
the chlorine oxyanions as a group, the trend amongst the AMOEBA values of |h Eint i| is
ClO2− > ClO3− & ClO−  ClO4− . The fact that perchlorate (ClO4− ) 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 h Eint i, h Eelst i, and h Eind i all follow the same trend exhibited by
the gas-phase dipole moments of the ions in question. However, this means that the trend
Molecules 2022, 26, 6719 11 of 20

amongst total interaction energies is different from that predicted by AMOEBA. Instead,
from the quantum-mechanical calculations the trend (from strongly to weakly interacting)
is ClO2− > ClO−  ClO3− & ClO4− .

Table 2. Ensemble-averaged energy components for ClO−

n (aq) in bulk water computed at the
SAPT0 + MBD/6-311+G(d,p) level of theory.

Dipole Energy Components (kcal/mol) a

Ion
b
Moment (D) h Eint i h Eelst i h Eexch i h Eind i h Edisp i
ClO− 3.04 −126.7 −137.6 90.2 −47.0 −32.2
ClO2− 3.20 −136.2 −149.6 103.6 −52.3 −37.8
ClO3− 2.46 −85.8 −120.6 125.1 −42.1 −48.2
ClO4− 0.00 −82.5 −78.1 41.1 −18.8 −26.7
a hE i = hE b
int elst i + h Eexch i + h Eind i + h Edisp i, up to roundoff error in the averaging. ωB97X-V/6-311+G(d) level
of theory at the optimized gas-phase geometry, with the center of nuclear charge as the origin.

In contrast to the AMOEBA results, the SAPT0 + MBD calculations afford similar
ensemble-averaged interaction energies for both ClO3− and ClO4− . Given that all of the
chlorine oxyanions except for ClO4− have sizable dipole moments, the ClO3− 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 interac-
tion 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.

(a) total interaction (b) elst + exch

–40 60
bulk bulk
interface 40 interface
–60
Eelst + Eexch (kcal/mol)

20
–80
Eint (kcal/mol)

0
–100
–20
–120
–40

–140
–60

–160 –80
4 —

4 —
2 —

2 —
3 —

3 —
N—

N—
2 —

3 —

2 —

3 —
—

—
r—

r—
N—

N—
N—

N—
l—

l—
I—

I—
lO

lO
lO

lO
lO

lO
lO

lO
C

C
O

O
B

B
C

C
C

C
SC

SC
C

C
N

N
C

C
C

C
O

O
C

Figure 6. Box-and-whisker plots showing mean values, standard deviations, and extremal values of (a) the total interaction
energy (Eint ) and (b) first-order electrostatics plus exchange (Eelst+exch ) for ion–water clusters. These calculations were
performed at the SAPT0(KS) + MBD/6-311+G(d,p) level of theory and therefore include intramolecular electron correlation
effects. Panel (a) should be compared to Figure 5a, as the difference lies solely in whether HF or KS molecular orbitals
are used within the SAPT0 formalism, and the vertical scales are the same in both figures. Similarly, panel (b) should be
compared to Figure 7a although the vertical scales are slightly different.

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 environ-
ments and ensemble-averaging over either data set. In considering the energy decom-
position in Equation (1), we have opted to group first-order electrostatics and exchange
together,
Eelst+exch = Eelst + Eexch , (10)
because their sum approximates the electrostatic interaction between antisymmetrized
monomer wave functions. This combination of “primitive” electrostatics (Eelst , which is
the Coulomb interaction between isolated-monomer charge densities) and Pauli repulsion
Molecules 2022, 26, 6719 12 of 20

(Eexch ) has proven to be easier to interpret for halide–water systems as compared to

electrostatics alone [65,81], in part because Eelst and Eexch 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 h Eexch i 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 h Eelst i, suggesting a hydration sphere that is not as
tight around the ion as it is in the smaller (but electrostatically much more attractive) ClO− ,
ClO2− , and ClO3− ions.

(a) electrostatics + exchange (b) induction (polarization + CT)

40 –10
bulk
interface
20 –20
Eelst + Eexch (kcal/mol)

Eind (kcal/mol)
0 –30

–20 –40

–40 –50

–60 –60 bulk

interface
–80 –70
4 —
2 —

3 —

N—

2 —

3 —
—
r—

N—

N—
l—

4 —
2 —

3 —

N—
I—

2 —

3 —
—
r—

N—

N—
l—

I—
lO

lO
lO
lO

lO
C

O
B

lO
lO

lO
C

O
C

B
C

SC

C
C

SC
C

C
C
C

O
C

N
C
C

O
(c) dispersion C
–10

–15

–20
Edisp (kcal/mol)

–25

–30

–35

–40

–45

–50
bulk
interface
–55
4 —
2 —

3 —

N—

2 —

3 —
—
r—

N—

N—
l—

I—

lO

lO
lO

lO
C

O
B

C
C

SC
C

N
C
C

O
C

Figure 7. Box-and-whisker plots showing mean values, standard deviations, and extremal values of ion–water interaction
energy components computed at the SAPT0 + MBD/6-311+G(d,p) level: (a) first-order electrostatics and exchange, Eelst +
Eexch ; (b) Eind from Equation (5); and (c) Edisp from the MBD model. The sum of the energy components in panels (a–c)
equals Eint in Figure 5a.

Statistical distributions of Eelst+exch are shown in Figure 7a for all of the ions, and
ClO3− immediately stands out as the only one for which h Eelst+exch i > 0. In other words,
the sum of first-order interactions is net repulsive for ClO3− 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
Eelst+exch are tiny in comparison to the instantaneous fluctuations, as was the case for Eint .
Furthermore, this anomalous prediction regarding ClO3− 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 Eelst+exch at that level of theory
that are plotted in Figure 6b. We note that the largest values of Eexch often correspond
to the largest (most attractive) total interaction energies, as is seen for example in SAPT
−
calculations of ClO− n · · · C6 H6 complexes (n = 1, 2, 3, 4) [100]. In the present case, ClO3
bucks this trend, according to the energy components listed in Table 2.
Molecules 2022, 26, 6719 13 of 20

A possible explanation for the apparently anomalous behavior of ClO3− 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 ClO3− is the
appearance of two distinct peaks in the RDF for Cl· · · Ow (where Ow 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 ClO3−
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 ClO3− is an indication of an extraneous water molecule
present at short range, then this could explain the anomalously repulsive values of Eelst+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· · · Ow RDF
is also observed in our previous simulations of SO23− (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
ClO3− (aq) and SO23− (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 com-
ponents 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 | Edisp | decreases in the order ClO3− > ClO2− > ClO− > ClO4− .
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 Edisp using the
third-generation ab initio dispersion potential aiD3 [50], which consists of atom–atom C6
and C8 potentials fitted to dispersion-only data from high-level SAPT calculations. Dis-
persion 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 disper-
sion 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
C6 coefficients [104]. For small molecules, three-body dispersion effects are typically quite
small [90].) Notably, in the aiD3 model the C6 and C8 coefficients depend only on atomic
number and do not respond to the electronic structure of the monomers.
The sharp drop in dispersion between chlorate (ClO3− ) 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 ClO3− and ClO4− . As compared to
Molecules 2022, 26, 6719 14 of 20

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.

Table 3. Dispersion energies for ClO−

n (aq) computed using different models and averaged over the
bulk and interfacial data sets.

h Edisp i (kcal/mol) a
System
aiD3 MBD b
ClO− bulk −28.0 ± 2.3 −32.2 ± 2.6
ClO− interface −24.4 ± 2.2 −27.9 ± 2.5
ClO2− bulk −35.3 ± 2.8 −37.8 ± 2.9
ClO2− interface −34.8 ± 3.2 −37.2 ± 3.4
ClO3− bulk −49.7 ± 3.8 −48.2 ± 3.3
ClO3− interface −47.5 ± 4.4 −46.0 ± 4.1
ClO4− bulk −28.4 ± 2.6 −26.7 ± 2.5
ClO4− interface −22.7 ± 4.5 −20.9 ± 4.2
a Uncertainties represent one standard deviation. b Based on HF monomer wave functions.

Polarization is often invoked in discussions of ions at the air/water interface [9–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 h Eind i 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 CT-
free polarization energy (Epol ) is much larger than the CT energy (ECT ), the result is that
Epol follows essentially the same trend from ion to ion as does the total induction energy,
Eind . 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− , ClO4− , and SCN− [19,34,35,105], the polarization energy is
significantly smaller in the interfacial environment, by at least the standard deviation of
Epol in bulk water; see Figure 8a.

(a) polarization (b) charge transfer

–5 0
–10
–2
–15
–20 –4
ECT (kcal/mol)
Epol (kcal/mol)

–25
–6
–30
–35
–8
–40
–45 –10

–50 bulk bulk

–12
–55
interface interface

–60 –14
4 —

4 —
2 —

2 —
3 —

N—

3 —

N—
2 —

3 —
—

2 —

3 —
r—

—
r—
N—

N—

N—

N—
l—

l—
I—

I—
lO

lO
lO

lO
lO

lO
lO

lO
C

O
B

B
C

C
C

SC

SC
C

C
N

N
C

C
C

C
O
C

O
C

Figure 8. Mean values, standard deviations, and extremal values of (a) the CT-free polarization energy (Epol ) and (b) the CT
energy (ECT ), computed at the SAPT0 + MBD/6-311+G(d,p) level. Together, Epol + ECT = Eind , and the sum (induction
energy) is plotted in Figure 7b.

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 ECT
for all of the ions are plotted in Figure 8b. In contrast to other energy components, only for
Molecules 2022, 26, 6719 15 of 20

ECT 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 | ECT | . 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− , ClO4− , 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 h ECT i 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 ECT .
In the context of the Hofmeister series [106,107], the anions I− , ClO4− , 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, NO3− 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 ECT
that we obtain for NO3− 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 | ECT | 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.

4. 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− , ClO4− , and SCN− ), as well as for hard anions that are not considered to be surface
active (CN− , OCN− , and NO2− ). That said, systematic differences in the mean values h Eint i
and h Epol i in bulk versus interfacial environments are rather small in comparison to the
magnitude of the instantaneous fluctuations in Eint and Epol .
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− , ClO4− , and SCN− . However, NO3− can also
be detected in surface-sensitive vibrational spectroscopy [19], yet for that species h ECT i
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–113].
Molecules 2022, 26, 6719 16 of 20

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 (first-
shell) solvation structure is extremely similar in the bulk and interfacial environments [7].
The detailed mechanism of soft anion surface activity remains an open question.

Supplementary Materials: The following are available online: AMOEBA parameter file and coordi-
nates for each of the structures used.
Author Contributions: Conceptualization, J.M.H.; methodology, J.M.H.; software, J.M.H. and S.K.P.;
validation, J.M.H.; formal analysis, J.M.H.; investigation, J.M.H. and S.K.P.; resources, J.M.H.; data
curation, J.M.H.; writing—original draft preparation, J.M.H.; writing—review and editing, J.M.H. and
S.KP.; visualization, J.M.H.; supervision, J.M.H.; project administration, J.M.H.; funding acquisition,
J.M.H. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the U.S. National Science Foundation under grant number
CHE-1955282.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: The data that support this study are available from the corresponding
author upon reasonable request.
Acknowledgments: It is a pleasure to wish Frank Weinhold a happy 80th birthday. J.M.H. thanks
Frank for his excellent teaching (in both thermodynamics and quantum chemistry) when the author
was a graduate student at the University of Wisconsin.
Conflicts of Interest: J.M.H. serves on the board of directors of Q-Chem Inc.
Sample Availability: Samples of the compounds are available from the authors.

References
1. Verreault, D.; Hua, W.; Allen, H.C. From conventional to phase-sensitive vibrational sum frequency generation spectroscopy:
Probing water organization at aqueous interfaces. J. Phys. Chem. Lett. 2012, 3, 3012–3028. [CrossRef]
2. Tang, F.; Ohto, T.; Sun, S.; Rouxel, J.R.; Imoto, S.; Backus, E.H.G.; Mukamel, S.; Bonn, M.; Nagata, Y. Molecular structure and
modeling of water–air and ice–air interfaces monitored by sum-frequency generation. Chem. Rev. 2020, 120, 3633–3667. [CrossRef]
3. Raymond, E.A.; Richmond, G.L. Probing the molecular structure and bonding of the surface of aqueous salt solutions. J. Phys.
Chem. B 2004, 108, 5051–5059. [CrossRef]
4. Liu, D.; Ma, G.; Levering, L.M.; Allen, H.C. Vibrational spectroscopy of aqueous sodium halide solutions and air–liquid interfaces:
Observation of increased interfacial depth. J. Phys. Chem. B 2004, 108, 2252–2260. [CrossRef]
5. Mucha, M.; Frigato, T.; Levering, L.M.; Allen, H.C.; Tobias, D.J.; Dang, L.X.; Jungwirth, P. Unified molecular picture of the
surfaces of aqueous acid, base, and salt solutions. J. Phys. Chem. B 2005, 109, 7617–7623. [CrossRef]
6. Garrett, B. Ions at the air/water interface. Science 2004, 303, 1146–1147. [CrossRef] [PubMed]
7. Paul, S.K.; Herbert, J.M. Probing interfacial effects on ionization energies: The surprising banality of anion–water hydrogen
bonding at the air/water interface. J. Am. Chem. Soc. 2021, 143, 10189–10202. [CrossRef]
8. Collins, K.D. Charge density-dependent strength of hydration and biological structure. Biophys. J. 1997, 72, 65–76. [CrossRef]
9. Jungwirth, P.; Tobias, D.J. Ions at the air/water interface. J. Phys. Chem. B 2002, 106, 6361–6373. [CrossRef]
10. Dang, L.X. Computational study of ion binding to the liquid interface of water. J. Phys. Chem. B 2002, 106, 10388–10394. [CrossRef]
11. Herce, D.H.; Perera, L.; Darden, T.A.; Sagui, C. Surface solvation for an ion in a water cluster. J. Chem. Phys. 2005, 122, 024513.
[CrossRef] [PubMed]
12. Archontis, G.; Leontidis, E.; Andreou, G. Attraction of iodide ions by the free water surface, revealed by simulations with a
polarizable force field based on Drude oscillators. J. Phys. Chem. B 2005, 109, 17957–17966. [CrossRef] [PubMed]
13. Jungwirth, P.; Tobias, D.J. Specific ion effects at the air/water interface. Chem. Rev. 2006, 106, 1259–1281. [CrossRef] [PubMed]
Molecules 2022, 26, 6719 17 of 20

14. Chang, T.M.; Dang, L.X. Recent advances in molecular simulations of ion solvation at liquid interfaces. Chem. Rev. 2006, 106,
1305–1322. [CrossRef] [PubMed]
15. Ishiyama, T.; Morita, A. Molecular dynamics study of gas–liquid aqueous sodium halide interfaces. I. Flexible and polarizable
molecular modeling and interfacial properties. J. Phys. Chem. C 2007, 111, 721–737. [CrossRef]
16. Dweik, J.; Srour, M.; Karaky, K.; Kobeissi, M.; Joumaa, W.; Abou-Saleh, K. Molecular simulation of ion transport at the water/
vapor interface. Open J. Phys. Chem. 2012, 2, 147–155. [CrossRef]
17. Sun, L.; Li, X.; Tu, Y.; Ågren, H. Origin of ion selectivity at the air/water interface. Phys. Chem. Chem. Phys. 2015, 17, 4311–4318.
[CrossRef] [PubMed]
18. Jungwirth, P.; Curtis, J.E.; Tobias, D.J. Polarizability and aqueous solvation of the sulfate dianion. Chem. Phys. Lett. 2003, 367,
704–710. [CrossRef]
19. Jubb, A.M.; Hua, W.; Allen, H.C. Environmental chemistry at vapor/water interfaces: Insights from vibrational sum frequency
generation spectroscopy. Annu. Rev. Phys. Chem. 2012, 63, 107–130. [CrossRef]
20. Jungwirth, P.; Tobias, D.J. Molecular structure of salt solutions: A new view of the interface with implications for heterogeneous
atmospheric chemistry. J. Phys. Chem. B 2001, 105, 10468–10472. [CrossRef]
21. Caleman, C.; Hub, J.S.; van Maaren, P.J.; van der Spoel, D. Atomistic simulation of ion solvation in water explains surface
preference of halides. Proc. Natl. Acad. Sci. USA 2011, 108, 6838–6842. [CrossRef]
22. Petersen, P.B.; Saykally, R.J.; Mucha, M.; Jungwirth, P. Enhanced concentration of polarizable anions at the liquid water surface:
SHG spectroscopy and MD simulations of sodium thiocyanide. J. Phys. Chem. B 2005, 109, 10915–10921. Erratum: ibid. 2015, 109,
13402. [CrossRef]
23. Horinek, D.; Herz, A.; Vrbka, L.; Sedlmeier, F.; Mamatkulov, S.I.; Netz, R.R. Specific ion adsorption at the air/water interface: The
role of hydrophobic solvation. Chem. Phys. Lett. 2009, 479, 173–183. [CrossRef]
24. Netz, R.R.; Horinek, D. Progress in modeling of ion effects at the vapor/water interface. Annu. Rev. Phys. Chem. 2012, 63, 401–418.
[CrossRef]
25. Baer, M.D.; Mundy, C.J. Toward an understanding of the specific ion effect using density functional theory. J. Phys. Chem. Lett.
2011, 2, 1088–1093. [CrossRef]
26. Baer, M.D.; Stern, A.C.; Levin, Y.; Tobias, D.J.; Mundy, C.J. Electrochemical surface potential due to classical point charge models
drives anion adsorption to the air-water interface. J. Phys. Chem. Lett. 2012, 3, 1565–1570. [CrossRef] [PubMed]
27. Ou, S.; Hu, Y.; Patel, S.; Wan, H. Spherical monovalent ions at aqueous liquid–vapor interfaces: Interfacial stability and induced
interface fluctuations. J. Phys. Chem. B 2013, 117, 11732–11742. [CrossRef] [PubMed]
28. Ou, S.C.; Cui, D.; Patel, S. Molecular modeling of ions at interfaces: Exploring similarities to hydrophobic solvation through the
lens of induced aqueous interfacial fluctuations. Phys. Chem. Chem. Phys. 2016, 18, 30357–30365. [CrossRef]
29. Ishiyama, T.; Imamura, T.; Morita, A. Theoretical studies of structures and vibrational sum frequency generation spectra at
aqueous interfaces. Chem. Rev. 2014, 114, 8447–8470. [CrossRef] [PubMed]
30. Levin, Y.; dos Santos, A.P.; Diehl, A. Ions at the air-water interface: An end to a hundred-year-old mystery? Phys. Rev. Lett. 2009,
103, 257802. [CrossRef]
31. Levin, Y.; dos Santos, A.P. Ions at hydrophobic interfaces. J. Phys. Condens. Matt. 2014, 26, 203101. [CrossRef] [PubMed]
32. Wick, C.D.; Xantheas, S.S. Computational investigation of the first solvation shell structure of interfacial and bulk aqueous
chloride and iodide ions. J. Phys. Chem. B 2009, 113, 4141–4146. [CrossRef] [PubMed]
33. Baer, M.D.; Mundy, C.J. An ab initio approach to understanding the specific ion effect. Faraday Discuss. 2013, 160, 89–101. [CrossRef]
34. Viswanath, P.; Motschmann, H. Oriented thiocyanate anions at the air–electrolyte interface and its implications on interfacial
water—A vibrational sum frequency spectroscopy study. J. Phys. Chem. C 2007, 111, 4484–4486. [CrossRef]
35. Viswanath, P.; Motschmann, M. Effect of interfacial presence of oriented thiocyanate on water structure. J. Phys. Chem. C 2008,
112, 2099–2103. [CrossRef]
36. Cooper, R.J.; O’Brien, J.T.; Chang, T.M.; Williams, E.R. Structural and electrostatic effects at the surfaces of size- and charge-selected
aqueous nanodrops. Chem. Sci. 2017, 8, 5201–5213. [CrossRef]
37. Bakker, H.J. Structural dynamics of aqueous salt solutions. Chem. Rev. 2008, 108, 1456–1473. [CrossRef]
38. Omta, A.W.; Kropman, M.F.; Woutersen, S.; Bakker, H.J. Negligible effect of ions on the hydrogen-bond structure in liquid water.
Science 2003, 301, 347–349. [CrossRef] [PubMed]
39. Park, S.; Fayer, M.D. Hydrogen bond dynamics in aqueous NaBr solutions. Proc. Natl. Acad. Sci. USA 2007, 104, 16731–16738.
[CrossRef]
40. Moilanen, D.E.; Wong, D.; Rosenfeld, D.E.; Fenn, E.E.; Fayer, M.D. Ion–water hydrogen-bond switching observed with 2D IR
vibrational echo chemical exchange spectroscopy. Proc. Natl. Acad. Sci. USA 2009, 106, 375–380. [CrossRef]
41. Botti, A.; Bruni, F.; Imberti, S.; Ricci, M.A.; Soper, A.K. Ions in water: The microscopic structure of concentrated NaOH solutions.
J. Chem. Phys. 2004, 120, 10154–10162. [CrossRef]
42. Botti, A.; Bruni, F.; Imberti, S.; Ricci, M.A.; Soper, A.K. Ions in water: The microscopic structure of a concentrated HCl solution. J.
Chem. Phys. 2004, 121, 7840–7848. [CrossRef] [PubMed]
43. Imberti, S.; Botti, A.; Bruni, F.; Cappa, G.; Ricci, M.A.; Soper, A.K. Ions in water: The microscopic structure of concentrated
hydroxide solutions. J. Chem. Phys. 2005, 122, 194509. [CrossRef]
Molecules 2022, 26, 6719 18 of 20

44. Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M.A.; Soper, A.K. Perturbation of water structure due to monovalent ions in solution.
Phys. Chem. Chem. Phys. 2007, 9, 2959–2967. [CrossRef] [PubMed]
45. Chandra, A. Effects of ion atmosphere on hydrogen-bond dynamics in aqueous electrolyte solutions. Phys. Rev. Lett. 2000, 85,
768–771. [CrossRef]
46. Chialvo, A.A.; Simonson, J.M. The effect of salt concentration on the structure of water in CaCl2 aqueous solutions. J. Mol. Liq.
2004, 112, 99–105. [CrossRef]
47. Tielrooij, K.J.; Garcia-Araez, N.; Bonn, M.; Bakker, H.J. Cooperativity in ion hydration. Science 2010, 328, 1006–1009. [CrossRef]
48. Migliorati, V.; Mancini, G.; Chillemi, G.; Zitolo, A.; D’Angelo, P. Effect of the Zn2+ and Hg2+ ions on the structure of liquid water.
J. Phys. Chem. A 2011, 115, 4798–4803. [CrossRef] [PubMed]
49. Seidel, R.; Winter, B.; Bradforth, S.E. Valence electronic structure of aqueous solutions: Insights from photoelectron spectroscopy.
Annu. Rev. Phys. Chem. 2016, 67, 283–305. [CrossRef] [PubMed]
50. Lao, K.U.; Herbert, J.M. Accurate and efficient quantum chemistry calculations of noncovalent interactions in many-body systems:
The XSAPT family of methods. J. Phys. Chem. A 2015, 119, 235–253. [CrossRef]
51. Carter-Fenk, K.; Lao, K.U.; Herbert, J.M. Predicting and understanding non-covalent interactions using novel forms of symmetry-
adapted perturbation theory. Acc. Chem. Res. 2021, 54, 3679–3690. [CrossRef] [PubMed]
52. Francisco, E.; Pendás, A.M. Energy partition analyses: Symmetry-adapted perturbation theory and other techniques. In Non-
Covalent Interactions in Quantum Chemistry and Physics; de la Roza, A.O., DiLabio, G.A., Eds.; Elsevier: Amsterdam, The Netherlan-
das, 2017; Chapter 2, pp. 27–64. [CrossRef]
53. Patkowski, K. Recent developments in symmetry-adapted perturbation theory. WIREs Comput. Mol. Sci. 2020, 10, e1452. [CrossRef]
54. Ren, P.; Ponder, J.W. Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B 2003, 107,
5933–5947. [CrossRef]
55. Grossfield, A.; Ren, P.; Ponder, J.W. Ion solvation thermodynamics from simulation with a polarizable force field. J. Am. Chem.
Soc. 2003, 125, 15671–15682. [CrossRef]
56. Ren, P.; Wu, C.; Ponder, J.W. Polarizable atomic multipole-based molecular mechanics for organic molecules. J. Chem. Theory
Comput. 2011, 7, 3143–3161. [CrossRef]
57. Rackers, J.A.; Wang, Z.; Lu, C.; Laury, M.L.; Lagardère, L.; Schnieders, M.J.; Piquemal, J.P.; Ren, P.; Ponder, J.W. Tinker 8: Software
tools for molecular design. J. Chem. Theory Comput. 2018, 14, 5273–5289. [CrossRef]
58. Parker, T.M.; Burns, L.A.; Parrish, R.M.; Ryno, A.G.; Sherrill, C.D. Levels of symmetry adapted perturbation theory (SAPT). I.
Efficiency and performance for interaction energies. J. Chem. Phys. 2014, 140, 094106. [CrossRef]
59. Gray, M.; Herbert, J.M. Simplified tuning of long-range corrected density functionals for symmetry-adapted perturbation theory.
J. Chem. Phys. 2021, 155, 034103. [CrossRef]
60. Lao, K.U.; Herbert, J.M. Atomic orbital implementation of extended symmetry-adapted perturbation theory (XSAPT) and
benchmark calculations for large supramolecular complexes. J. Chem. Theory Comput. 2018, 14, 2955–2978. [CrossRef] [PubMed]
61. Ambrosetti, A.; Reilly, A.M.; DiStasio, R.A., Jr.; Tkatchenko, A. Long-range correlation energy calculated from coupled atomic
response functions. J. Chem. Phys. 2014, 140, 18A508. [CrossRef] [PubMed]
62. Carter-Fenk, K.; Lao, K.U.; Liu, K.Y.; Herbert, J.M. Accurate and efficient ab initio calculations for supramolecular complexes:
Symmetry-adapted perturbation theory with many-body dispersion. J. Phys. Chem. Lett. 2019, 10, 2706–2714. [CrossRef]
63. Gray, M.; Herbert, J.M. Comprehensive basis-set testing of extended symmetry-adapted perturbation theory. 2021. [CrossRef]
64. Epifanovsky, E.; Gilbert, A.T.B.; Feng, X.; Lee, J.; Mao, Y.; Mardirossian, N.; Pokhilko, P.; White, A.F.; Coons, M.P.; Dempwolff,
A.L.; et al. Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package. J. Chem.
Phys. 2021, 155, 084801. [CrossRef] [PubMed]
65. Herbert, J.M. Neat, simple, and wrong: Debunking electrostatic fallacies regarding noncovalent interactions. J. Phys. Chem. A
2021, 125, 7125–7137. [CrossRef]
66. Rohrdanz, M.A.; Martins, K.M.; Herbert, J.M. A long-range-corrected density functional that performs well for both ground-state
properties and time-dependent density functional theory excitation energies, including charge-transfer excited states. J. Chem.
Phys. 2009, 130, 054112. [CrossRef] [PubMed]
67. Lao, K.U.; Herbert, J.M. Symmetry-adapted perturbation theory with Kohn-Sham orbitals using non-empirically tuned, long-
range-corrected density functionals. J. Chem. Phys. 2014, 140, 044108. [CrossRef] [PubMed]
68. Hapka, M.; Rajchel, L.; Modrzejewski, M.; Chałasiǹski, G.; Szczȩśniak, M.M. Tuned range-separated hybrid functionals in the
symmetry-adapted perturbation theory. J. Chem. Phys. 2014, 141, 134120. [CrossRef] [PubMed]
69. Baer, R.; Livshits, E.; Salzner, U. Tuned range-separated hybrids in density functional theory. Annu. Rev. Phys. Chem. 2010, 61,
85–109. [CrossRef] [PubMed]
70. Alam, B.; Morrison, A.F.; Herbert, J.M. Charge separation and charge transfer in the low-lying excited states of pentacene. J. Phys.
Chem. C 2020, 124, 24653–24666. [CrossRef]
71. Uhlig, F.; Herbert, J.M.; Coons, M.P.; Jungwirth, P. Optical spectroscopy of the bulk and interfacial hydrated electron from ab
initio calculations. J. Phys. Chem. A 2014, 118, 7507–7515. [CrossRef]
72. Jacobson, L.D.; Herbert, J.M. An efficient, fragment-based electronic structure method for molecular systems: Self-consistent
polarization with perturbative two-body exchange and dispersion. J. Chem. Phys. 2011, 134, 094118. [CrossRef]
Molecules 2022, 26, 6719 19 of 20

73. Herbert, J.M.; Jacobson, L.D.; Lao, K.U.; Rohrdanz, M.A. Rapid computation of intermolecular interactions in molecular and ionic
clusters: Self-consistent polarization plus symmetry-adapted perturbation theory. Phys. Chem. Chem. Phys. 2012, 14, 7679–7699.
[CrossRef]
74. Jacobson, L.D.; Richard, R.M.; Lao, K.U.; Herbert, J.M. Efficient monomer-based quantum chemistry methods for molecular and
ionic clusters. Annu. Rep. Comput. Chem. 2013, 9, 25–58. [CrossRef]
75. Liu, K.Y.; Carter-Fenk, K.; Herbert, J.M. Self-consistent charge embedding at very low cost, with application to symmetry-adapted
perturbation theory. J. Chem. Phys. 2019, 151, 031102. [CrossRef] [PubMed]
76. Moszyński, R.; Cybulski, S.M.; Chałasiński, G. Many-body theory of intermolecular induction interactions. J. Chem. Phys. 1994,
100, 4998–5010. [CrossRef]
77. Lao, K.U.; Herbert, J.M. Energy decomposition analysis with a stable charge-transfer term for interpreting intermolecular
interactions. J. Chem. Theory Comput. 2016, 12, 2569–2582. [CrossRef]
78. Kaduk, B.; Kowalczyk, T.; Van Voorhis, T. Constrained density functional theory. Chem. Rev. 2012, 112, 321–370. [CrossRef]
[PubMed]
79. Řezáč, J.; de la Lande, A. Robust, basis-set independent method for the evaluation of charge-transfer energy in noncovalent
complexes. J. Chem. Theory Comput. 2015, 11, 528–537. [CrossRef]
80. Řezác, J.; de la Lande, A. On the role of charge transfer in halogen bonding. Phys. Chem. Chem. Phys. 2017, 19, 791–803. [CrossRef]
[PubMed]
81. Herbert, J.M.; Carter-Fenk, K. Electrostatics, charge transfer, and the nature of the halide–water hydrogen bond. J. Phys. Chem. A
2021, 125, 1243–1256. [CrossRef] [PubMed]
82. Becke, A.D. A multicenter numerical integration scheme for polyatomic molecules. J. Chem. Phys. 1988, 88, 2547–2553. [CrossRef]
83. Medvedev, N.N. The algorithm for three-dimensional Voronoi polyhedra. J. Comput. Phys. 1986, 67, 223–229. [CrossRef]
84. Slater, J.C. Atomic radii in crystals. J. Chem. Phys. 1964, 41, 3199–3204. [CrossRef]
85. Dasgupta, S.; Herbert, J.M. Standard grids for high-precision integration of modern density functionals: SG-2 and SG-3. J.
Comput. Chem. 2017, 38, 869–882. [CrossRef] [PubMed]
86. Wu, Q.; Chen, C.L.; Van Voorhis, T. Configuration interaction based on constrained density functional theory: A multireference
method. J. Chem. Phys. 2007, 127, 164119. [CrossRef] [PubMed]
87. Richard, R.M.; Lao, K.U.; Herbert, J.M. Aiming for benchmark accuracy with the many-body expansion. Acc. Chem. Res. 2014, 47,
2828–2836. [CrossRef]
88. Lao, K.U.; Liu, K.Y.; Richard, R.M.; Herbert, J.M. Understanding the many-body expansion for large systems. II. Accuracy
considerations. J. Chem. Phys. 2016, 144, 164105. [CrossRef]
89. Liu, K.Y.; Herbert, J.M. Understanding the many-body expansion for large systems. III. Critical role of four-body terms,
counterpoise corrections, and cutoffs. J. Chem. Phys. 2017, 147, 161729. [CrossRef]
90. Herbert, J.M. Fantasy versus reality in fragment-based quantum chemistry. J. Chem. Phys. 2019, 151, 170901. [CrossRef]
91. Heindel, J.P.; Xantheas, S.S. The many-body expansion for aqueous systems revisited: I. Water–water interactions. J. Chem. Theory
Comput. 2020, 16, 6843–6855. [CrossRef]
92. Heindel, J.P.; Xantheas, S.S. The many-body expansion for aqueous systems revisited: II. Alkali metal and halide ion–water
interactions. J. Chem. Theory Comput. 2021, 17, 2200–2216. [CrossRef]
93. Tissandier, M.D.; Cowen, K.A.; Feng, W.Y.; Gundlach, E.; Cohen, M.H.; Earhart, A.D.; Coe, J.V.; Tuttle, T.R., Jr. The proton’s absolute
aqueous enthalpy and Gibbs free energy of solvation from cluster-ion solvation data. J. Phys. Chem. A 1998, 102, 7787–7794.
Erratum: ibid. 1998, 102, 9308. [CrossRef]
94. Prasetyo, N.; Hünenberger, P.H.; Hofer, T.S. Single-ion thermodynamics from first principles: Calculation of the absolute
hydration free energy and single-electrode potential of aqueous Li+ using ab initio quantum mechanical/molecular mechanical
dynamics simulations. J. Chem. Theory Comput. 2018, 14, 6443–6459. [CrossRef]
95. Malloum, A.; Fifen, J.J.; Conradie, J. Determination of the absolute solvation free energy and enthalpy of the proton in solutions.
J. Mol. Liq. 2021, 322, 114919. [CrossRef]
96. Herbert, J.M. Dielectric continuum methods for quantum chemistry. WIREs Comput. Mol. Sci. 2021, 11, e1519. [CrossRef]
97. Rackers, J.A.; Wang, Q.; Liu, C.; Piquemal, J.P.; Ren, P.; Ponder, J.W. An optimized charge penetration model for use with the
AMOEBA force field. Phys. Chem. Chem. Phys. 2017, 19, 276–291. [CrossRef]
98. Deng, S.; Wang, Q.; Ren, P. Estimating and modeling charge transfer from the SAPT induction energy. J. Comput. Chem. 2017, 38,
2222–2231. [CrossRef]
99. Jing, Z.; Liu, C.; Ren, P. Advanced electrostatic model for monovalent ions based on ab initio energy decompoosition. J. Chem. Inf.
Model. 2021, 61, 2806–2817. [CrossRef]
100. Elgengehi, S.M.; El-Taher, S.; Ibrahim, M.A.A.; El-Kelany, K.E. Unexpected favourable noncovalent interaction between chlorine
oxyanions (ClO− x ; x = 1–4) and benzene: Benchmarking DFT and SAPT methods with respect to CCSD(T). Comput. Theor. Chem.
2021, 1199, 113214. [CrossRef]
101. Eklund, L.; Hofer, T.S.; Persson, I. Structure and water exchange dynamics of hydrated oxo halo ions in aqueous solution using
QMCF MD simulations, large angle X-ray scattering and EXAFS. Dalton Trans. 2015, 44, 1816–1828. [CrossRef] [PubMed]
102. Eklund, L.; Hofer, T.S.; Pribil, A.B.; Rode, B.M.; Persson, I. On the structure and dynamics of the hydrated sulfite ion in aqueous
solution—An ab initio QMCF MD simulation and large angle x-ray scattering study. Dalton Trans. 2012, 41, 5209–5216. [CrossRef]
Molecules 2022, 26, 6719 20 of 20

103. Lao, K.U.; Herbert, J.M. A simple correction for nonadditive dispersion within extended symmetry-adapted perturbation theory
(XSAPT). J. Chem. Theory Comput. 2018, 14, 5128–5142. [CrossRef]
104. Dobson, J.F. Beyond pairwise additivity in London dispersion interactions. Int. J. Quantum Chem. 2014, 114, 1157–1161. [CrossRef]
105. Tong, Y.; Zhang, I.Y.; Campen, R.K. Experimentally quantifying anion polarizability at the air/water interface. Nat. Commun.
2018, 9, 1313. [CrossRef] [PubMed]
106. Zhang, Y.; Cremer, P.S. Interactions between macromolecules and ions: The Hofmeister series. Curr. Opin. Struc. Biol. 2006, 10,
659–663. [CrossRef] [PubMed]
107. Jungwirth, P.; Cremer, P.S. Beyond Hofmeister. Nat. Chem. 2014, 6, 261–263. [CrossRef] [PubMed]
108. Rembert, K.B.; Paterová, J.; Heyda, J.; Hilty, C.; Jungwirth, P.; Cremer, P.S. Molecular mechanisms of ion-specific effects on
proteins. J. Am. Chem. Soc. 2012, 134, 10039–10046. [CrossRef] [PubMed]
109. Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M.A.; Soper, A.K. Hydration of sodium, potassium, and chloride ions in solution and the
concept of structure maker/breaker. J. Phys. Chem. B 2007, 111, 13570–13577. [CrossRef] [PubMed]
110. Wick, C.D.; Lee, A.J.; Rick, S.W. How intermolecular charge transfer influences the air-water interface. J. Chem. Phys. 2012, 137,
154701. [CrossRef]
111. Vácha, R.; Marsalek, O.; Willard, A.P.; Bonthuis, D.J.; Netz, R.R.; Jungwirth, P. Charge transfer between water molecules as the
possible origin of the observed charging at the surface of pure water. J. Phys. Chem. Lett. 2012, 3, 107–111. [CrossRef]
112. Samson, J.S.; Scheu, R.; Smolentsev, N.; Rick, S.W.; Roke, S. Sum frequency spectroscopy of the hydrophobic nanodroplet/water
interface: Absence of hydroxyl ion and dangling OH bond signatures. Chem. Phys. Lett. 2014, 615, 124–131. [CrossRef]
113. Poli, E.; Jong, K.H.; Hassanali, A. Charge transfer as a ubiquitous mechanism in determining the negative charge at hydrophobic
interfaces. Nat. Commun. 2020, 11, 901. [CrossRef] [PubMed]