# Nuclear Quantum Effects on the Electronic Structure of Water and Ice

Margaret Berrens  Department of Chemistry, University of California Davis, One
Shields Ave. Davis, CA, 95616.  Arpan Kundu  Pritzker School of Molecular
Engineering, University of Chicago, Chicago, IL 60637  Marcos F. Calegari
Andrade  Quantum Simulations Group, Materials Science Division, Lawrence
Livermore National Laboratory, Livermore, California, 94550-5507  Tuan Anh
Pham  Quantum Simulations Group, Materials Science Division, Lawrence
Livermore National Laboratory, Livermore, California, 94550-5507  Giulia Galli
Pritzker School of Molecular Engineering, University of Chicago, Chicago, IL
60637  Department of Chemistry, University of Chicago, Chicago, Illinois
60637, United States  Materials Science Division and Center for Molecular
Engineering, Argonne National Laboratory, Lemont, Illinois 60439, United
States  Davide Donadio<EMAIL_ADDRESS>Department of Chemistry,
University of California Davis, One Shields Ave. Davis, CA, 95616.

###### Abstract

The electronic properties and optical response of ice and water are
intricately shaped by their molecular structure, including the quantum
mechanical nature of hydrogen atoms. In spite of numerous studies appeared
over decades, a comprehensive understanding of the effect of the nuclear
quantum motion on the electronic structure of water and ice at finite
temperatures remains elusive. Here, we utilize molecular simulations that
harness the efficiency of machine-learning potentials and many-body
perturbation theory to assess the impact of nuclear quantum effects on the
electronic structure of water and hexagonal ice. By comparing the results of
path-integral and classical simulations, we find that including nuclear
quantum effects leads to a larger renormalization of the fundamental gap of
ice, compared to that of water, eventually leading to a comparable gap in the
two systems, consistent with experimental estimates. Our calculations suggest
that the quantum fluctuations responsible for an increased delocalization of
protons in ice, relative to water, are a key factor leading to the enhancement
of nuclear quantum effects on the electronic structure of ice.

Water and ice provide a natural solvation environment for most chemical
processes of atmospheric and biological importance. The electronic structure
of water affects the chemistry of biological, geochemical, and environmental
reactions, as well as the properties of energy conversion devices, such as
photoelectrochemical cells.Zhong  _et al._ (2019); Sverjensky  _et al._
(2014); Walter  _et al._ (2010); Grätzel (2001) Understanding the relation
between the structural and electronic properties of water will help shed light
on the chemistry of natural processes and aid the optimization of
technological applications. There has been extensive research on the
electronic properties of liquid waterPham  _et al._ (2014); Chen  _et al._
(2016); Gaiduk  _et al._ (2018); Bischoff  _et al._ (2021); Tal  _et al._
(2024), but significantly less attention has been directed towards the
electronic structure of hexagonal ice (Ice Ih), water’s most common solid
form. A thorough understanding of its electronic structure is essential to
elucidate the role of ice as a substrate or a catalyst in environmental
chemistry,Grannas  _et al._ (2014); Klanova  _et al._ (2003); Takenaka and
Bandow (2007); Grannas  _et al._ (2007) and astrochemistry.Mifsud  _et al._
(2021); Jenniskens and Blake (1994)

Data inferred from experiments and interpreted by Chen et al. point at water
and ice having similar fundamental gaps ($E_{g}$), despite considerable
structural differences: $9.0\pm 0.2$ for waterBernas  _et al._ (1997); Winter
_et al._ (2004); Ambrosio  _et al._ (2017); Gaiduk  _et al._ (2018) at room
temperature and $9.4\pm 0.3$ for iceCampbell  _et al._ (1979); Shibaguchi  _et
al._ (1977); Baron  _et al._ (1978) at 77 K. Interpreting experiments at the
molecular level and connecting the structural and dynamical properties to the
electronic structure of water and ice remain challenging tasks, especially due
to the disordered and dynamic nature of these systems. Molecular dynamics (MD)
simulations combined with electronic structure calculations offer a promising
means to address this knowledge gap.Gaiduk  _et al._ (2018) Whereas nuclei are
treated as classical particles in first-principles MD (FPMD), the inclusion of
nuclear quantum effects (NQE) is necessary to describe the electron-phonon
interactions affecting the properties of materials with light
elements.Ceriotti  _et al._ (2016); Markland and Ceriotti (2018); Kundu  _et
al._ (2021, 2022); Alvertis and Engel (2022); Kundu and Galli (2023, 2024)
Nuclear quantum dynamics can be modeled using Feynman’s path integrals (PI) or
related methods,Berne and Thirumalai (1986); Marx and Parrinello (1996);
Herrero and Ramírez (2014); Huppert  _et al._ (2022) although at a
substantially higher computational cost.

Using trajectories generated with the Vydrov and Van Voorhis (rVV10) van der
Waals density functional,Sabatini  _et al._ (2013) path-integral MD (PIMD)
simulations with 6 beads, and supercells with 32 water molecules, Chen et al.
estimated the bandgap of water at 300 K using self-consistent GW calculations
with a two-point exchange-correlation kernel.Chen  _et al._ (2016) They
obtained a band gap of 9.8 and 8.9 eV without and with vertex corrections,
respectively, and in both cases a renormalization ($\Delta E_{g}^{NQE}$) of
-0.7 eV. A similar computational protocol was used by Tal et al. Tal  _et al._
(2024) who however increased the number of frames over which they conducted
MBPT calculations and slightly modified the treatment of the exchange-
correlation Kernel, obtaining a band gap of 9.2 eV when using vertex
corrections. Gaiduk et al. reported a bandgap of 10 eV with a $\Delta
E_{g}^{NQE}$ of -0.5 eV, using G0W0 calculations starting from hybrid DFT
(they used both range-separated (RSH) and self-consistent hybrids)Gaiduk  _et
al._ (2018). They conducted simulations with the MB-pol potential, PIMD with
32 beads and used 64 water molecule cells. The major difference between the
results obtained by Tal et al. Tal  _et al._ (2024) and Gaiduk et al.
originates from the different positions of the conduction band predicted in
the two papers, while the position of the valence band is similar. For ice Ih,
Engel et al. and Monserrat et al. estimated a large energy gap
renormalization, -1.5 eV, using stochastic methods based on the quasi-harmonic
approximationEngel  _et al._ (2015); Monserrat  _et al._ (2015), though the
accuracy of these methods in a strong anharmonic regime is not
guaranteed.Alvertis and Engel (2022); Kundu and Galli (2023). Bischoff et al.
reported a band gap of ice of 9.8 and 9.3 eV (and that of water of 9.6 and 9.1
eV), depending on the type of self-consistent GW method used. In the case of
ice, they applied just the zero-point renormalization (ZPR) obtained in
Ref.Engel  _et al._ (2015); Monserrat  _et al._ (2015) to electronic structure
calculations at zero T, carried out for a Bernal-Fowler unit cell with 12
molecules.Bischoff  _et al._ (2021) The origin of the difference between the
impact of NQE on the electronic properties of water and ice has not yet been
clarified, nor has $\Delta E_{g}^{NQE}$ been computed for ice by MD at finite
temperatures. In addition, comparisons have been made for calculations carried
out at different levels of electronic structure theory and with different
structural models (different force fields and/or density functionals).
Clearly, a fair and robust comparison calls for modeling the dynamical
properties of water and ice at finite temperatures on an equal footing, using
quantum simulations and computing the electronic structure of both systems
within the same approximations. Such a comparison is still amiss.

In this letter, we investigate the effect of NQEs on the electronic structure
of water and ice Ih by carrying out MBPT calculations at the G0W0 level of
theory for an ensemble of configurations obtained from classical and PIMD
simulations, to provide a consistent estimate of the bandgap renormalization
of both systems. Our calculations confirm that NQEs have a significantly
larger effect on the electronic properties of ice $I_{h}$ compared to those of
liquid water. A comparative analysis of hydrogen-bonding configurations in
water and ice suggests that the degree of transient proton transfer from
hydrogen bond donors to acceptors is the key structural feature affecting the
band-gap renormalization.

We use machine learning potentials (MLP) to accelerate both classical and PIMD
simulations. We specifically employ two different MLPs derived from density
functional theory (DFT) with two functionals: a van der Waals corrected hybrid
functional (revPBE0-D3),Perdew  _et al._ (1996a); Adamo and Barone (1999);
Grimme  _et al._ (2010) and the strongly constrained and appropriately normed
(SCAN) exchange-correlation functional.Sun  _et al._ (2015) Using revPBE0-D3,
we obtained a neural network potential fitted with an evolutionary algorithm;
this potential has been shown to reproduce accurately the thermodynamic
properties of water and ice, including the density anomaly of waterFan  _et
al._ (2021); Chen  _et al._ (2024), and it is of similar quality as previously
reported neural network models fitted on the same density functional. However
it is more computationally efficient.Cheng  _et al._ (2019); Schran  _et al._
(2020) In the case of the SCAN functional, we used a previously trained and
validated deep neural network potential (DNNP) that accurately reproduces
structural Zhang  _et al._ (2021) and vibrationalM. Sommers  _et al._ (2020)
properties of water and ice. We compared electronic structure results obtained
for trajectories with two different MLPs to rule out systematic biases
originating from a specific MLP scheme or the underlying density functional
approach. The differences and similarities of the two underlying density
functionals used in this work have been the focus of several ab intio MD
studies and are briefly discussed in the Supporting Information.Ohto  _et al._
(2019); Palos  _et al._ (2022) Our results are further validated using
trajectories generated with a DNNP fitted to the first-principles many-body
force field MB-pol.Bore and Paesani (2023)

We performed quantum simulations with a generalized Langevin equation
thermostat (PIGLET),Ceriotti and Manolopoulos (2012); Kapil  _et al._ (2019)
which provides well-converged results, compared to PIMD, using a small number
of beads ($p=8$) in the ring-polymer.Kapil  _et al._ (2016) MD simulations
were carried out at 300 K and 230 K for water and ice respectively. A
temperature of 230 K was chosen for ice to investigate the NQE on electronic
properties at environmentally relevant conditions. It is expected that as the
temperature increases the gap will only slightly decrease for Ice
$I_{h}$.Engel  _et al._ (2016) For the PIGLET simulations, at 230 K, the
average quantum kinetic and potential energy of ice is converged within 3
meV/atom with only 8 beads (Figure S2). Size effects on the calculated bandgap
of liquid water were carried out by Gaiduk  _et al._ (2018) Here, we tested
size effects on the bandgap of bulk ice using supercells with 96 and 192
molecules, which showed that a simulation cell of 96 water molecules is
sufficient to obtain an estimate of the bandgap within the statistical
uncertainty of our calculations (Table S1). The choice of system size in our
MD simulations was mainly constrained by the cost of performing several tens
of MBPT calculations over classical MD and PIMD trajectories.

Electronic structure calculations are performed using the G0W0 method on
electronic structure calculations at the PBE Perdew  _et al._ (1996b) level of
theory, for 50 equally spaced frames selected from MD and PIMD trajectories of
100 ps for water and ice, respectively (convergence as a function of the
number of frames is shown in Figure S3). The electronic structure calculations
from the PIMD runs were performed for a randomly chosen ring-polymer bead. All
G0W0 calculations were carried out with the WEST code.Govoni and Galli (2015)
We note that G0W0@PBE underestimates the $E_{g}$ of water, i.e. the energy
difference between the conduction band minimum (CBM) and valence band maximum
(VBM), of liquid water by approximately 1 eV lower, compared to higher level
theories using hybrid functionals as starting points for MBPT and compared to
experimental estimatesGaiduk  _et al._ (2018); Pham  _et al._ (2014). The use
of the quasiparticle self-consistent GW (QSGW) method that includes vertex
corrections also yields larger gapsChen  _et al._ (2016); Bischoff  _et al._
(2021). However, for liquid water one-shot G0W0@PBE calculations yield an
estimate of the effect of NQE on the electronic structure which is similar to
that of the QSGW method.Chen  _et al._ (2016). Since the focus of this work is
to estimate NQEs on the gap rather than the absolute values of the $E_{g}$, we
have utilized G0W0@PBE, which allowed us to sample a well-converged ensemble
of configurations for relatively large systems.

Figure 1: The difference between the fundamental gap obtained with classical
and quantum trajectories ($\Delta E_{g}^{NQE}$), upper panel, and the
difference between the position of the VBM (middle panel) and CBM (lower
panel) obtained with classical and quantum trajectories, for liquid water and
Ice $I_{h}$. The results on the left panels were obtained at the DFT level
using the revPBE0 functional for 100 frames (left) and those on the right
panels at the G0W0 level for 50 frames (right). The x-axis indicates the Ml
potentials used in the calculations. $\Delta$ represents the difference in
$\Delta E_{g}^{NQE}$ between water and ice for each ML potential.

The difference between the $E_{g}$ of water and ice obtained with classical
and quantum simulations, as well as the difference between the positions of
the VBM, and CBM, obtained from G0W0@PBE and hybrid calculations, are shown in
Figure 1. Regardless of the underlying potential used for sampling and the
electronic structure method, Figure 1 shows that NQEs considerably reduce the
bandgap of both water and ice by raising the energy of the VBM and lowering,
to a much smaller extent, the energy of the CBM. Consistent with previous
findings, we find that NQEs have a much larger impact in ice, for which the
$\Delta E_{g}^{NQE}$ amounts to 1.1 $\pm$ 0.32 – 1.5 $\pm$ 0.43 eV on the band
gap, which instead decreases by 0.70 $\pm$ 42 – 0.81 eV $\pm$ 0.33 in water.
Hence, the larger NQEs in ice overall lead to a similar $E_{g}$ for ice and
liquid water, consistent with data inferred from experiments ($9.0\pm 0.2$ eV
for water and $9.4\pm 0.3$ eV for ice). Our results for $\Delta E_{g}^{NQE}$
are consistent with previous calculations for both water (-0.7 eV in Chen  _et
al._ (2016) and -0.5 in Gaiduk  _et al._ (2018)) and ice (-1.52 eV in Engel
_et al._ (2015)). We note that the value for ice Ih in Ref. Engel  _et al._
(2015) was obtained for a configuration optimized at zero temperature. We find
a similar bandgap renormalization of 1.65 eV due to NQEs using a geometry-
optimized structure.

Overall, sampling with NEP enhances the difference in the NQEs-induced bandgap
renormalization between ice and water, although it remains within the same
uncertainty as the DNNP results. The error bars on the $\Delta E^{NQE}$ for
the VBM are much larger than for the CBM. The main effect of quantum dynamics
is not only to raise the energy of the VBM state, especially in ice, but also
to broaden its energy distribution. In turn, the CBM is only slightly affected
and does not undergo significant broadening.

Using the hybrid revPBE0 functional Perdew  _et al._ (1996a); Adamo and Barone
(1999); Grimme  _et al._ (2010), we verified that the trends reported above
for $\Delta E^{NQE}$ for water and ice hold for trajectories generated with
MB-pol data.Bore and Paesani (2023); Medders and Paesani (2014); Babin  _et
al._ (2014, 2013) This model provides an accurate description of the
structural, thermodynamic, dynamical, and spectroscopic properties of
waterCole  _et al._ (2016); Moberg  _et al._ (2017); Medders and Paesani
(2016)Bore and Paesani (2023). Finally, the corresponding $\Delta
E_{g}^{NQE}$’s computed from the Kohn-Sham orbitals using the SCAN functional
are reported in Figure S5. The values of the VBM and CBM energies from the
MBPT and the Kohn-Sham DFT calculations are also reported in Tables S2, S3,
and S4, respectively. We find that inexpensive DFT calculations (such as those
using the SCAN functional) reproduce the trends for NQEs obtained from
accurate MBPT, but, unsurprisingly, MBPT is necessary to achieve bandgaps in
reasonable agreement with experiments. As expected, the use of a hybrid
exchange and correlation functional brings DFT calculations much closer to
MBPT results and results for $\Delta E_{g}^{NQE}$ are within statistical
uncertainty. The underestimation of the bandgap with revPBE0 can be in part
attributed to the fact that the inverse of the appropriate fraction of exact
exchange ($\alpha$) in the hybrid functional is more than two times larger
than that of the electronic dielectric constant ($\varepsilon_{\infty}$) of
water.Pham  _et al._ (2014)

Figure 2: Probability density of the inverse participation ratio (IPR), fit
with a Gaussian kernel density estimation, for the VB calculated with SCAN
eigenstates, for liquid water (left) and ice (right). For each water system,
we report the IPR for the DNNP with (purple) and without (turquoise) NQEs, and
for the NEP with (light purple) and without (green) NQEs. A higher IPR
indicates stronger localization. The IPR values are scaled by 10000.

To understand the significant effect of quantum dynamics on the position of
the VBM we calculated the inverse participation ratio (IPR) for an ensemble of
snapshots from either PIMD or classical MD simulations (Figure 2). The IPR of
the $i$th Kohn-Sham orbital is calculated as
$\int|\psi_{i}|^{4}d^{3}r/(\int|\psi_{i}|^{2}d^{3}r)^{2}$ where a higher IPR
value corresponds to a more localized single-particle wave function. We
calculated IPRs for the Kohn-Sham orbitals using the SCAN functional including
states within 50 meV of the VBM, which we refer to as the valence band (VB).
For each water system and potential, we calculated the Gaussian kernel density
estimation of the IPR for the classical and quantum simulations. The histogram
of IPR values for both the VB and the conduction band (CB) is reported in
Figure S6. The IPR distribution for the VB spans values over two orders of
magnitude larger than those for the CB, consistent with the known localized
nature of occupied levels as opposed to delocalized empty levels. For both
water and ice, the NQEs shift the VB distribution toward higher IPR values,
indicating that quantum delocalization of the nuclei enhances the localization
of the electronic levels. This phenomenon, previously observed in liquid water
by Chen et al.,Chen  _et al._ (2016) was attributed to Anderson localization
resulting from an increased disorder in the proton distributions, due to
quantum fluctuations. The increased localization of electronic levels, akin to
surface effects observed in water,Gaiduk  _et al._ (2018) may also stem from
the weakening of hydrogen bonds due to quantum delocalization. Since water and
ice have different cell sizes, we cannot compare the absolute shift of the IPR
as IPR depends on cell size. Therefore we report relative IPR shift with
respect to the position of the peak of IPR distribution in the classical
simulation. For water and ice, there is an average shift in the relative peak
position of the IPR distribution of the VB of 18% and 66%, respectively. The
more pronounced shift for the VB of ice Ih is responsible for the larger
bandgap renormalization.

Figure 3: Joint probability distribution of the proton-transfer coordinate
$\nu$ and the distance d(O–O’) between the covalently bound and acceptor
oxygen atoms for liquid water (left) and Ice $I_{h}$ (right). The top two rows
are a comparison of the sampling with DNNP with and without NQEs, and the
bottom two rows are a comparison of the sampling with NEP with and without
NQEs.

Previous studies investigating the NQEs on the structural properties of liquid
water found that the inclusion of NQEs softens the local structure of water
due to proton delocalization.Markland and Ceriotti (2018) In particular, the
oxygen-hydrogen and hydrogen-hydrogen radial distribution functions (gOH and
gHH) exhibit the most apparent structural impatc from NQEs. Specifically, the
NQEs cause a broadening of the first two peaks of gOH, which suggests a higher
probability of proton transfer. The calculated radial and angular distribution
functions in Figure S7 confirm that both MLPs reproduce these structural
features for both water and ice. Similarly, NQEs cause a broadening of the
distribution function of the H-O-H angle for water and ice. However, the
causes of the difference in the NQE-induced band-gap renormalization between
water and ice, cannot be inferred from differences in radial distribution
functions. To investigate why $\Delta E_{g}^{NQE}$ is much larger in ice Ih
than in liquid water, we calculated the distributions of proton transfer
coordinates and oxygen-oxygen distances for each hydrogen-bonded pair of water
molecules (Figure 3) in our simulation cells. The proton transfer coordinate
is defined as $\nu=d_{OH}-d_{O^{\prime}H}$, where $d_{OH}$ is the distance
between O and H in the covalent bond and $d_{O^{\prime}H}$ is the distance
between the same hydrogen atom and the acceptor hydrogen-bonded oxygen
$O^{\prime}$.Ceriotti  _et al._ (2013); Fritsch  _et al._ (2014) Incorporating
NQEs in the simulations enhances hydrogen delocalization, yielding a small
fraction of configurations with $\nu>0$ that corresponds to short-lived
autoprotolysis events. Here, for liquid water, we see that NQEs broaden the
joint probability distribution $P(d_{OO^{\prime}},\nu)$ with a tail extending
into the $\nu>0$ region. Comparing the distributions obtained with different
MLPs and PIGLET, we observe a two times larger fraction of hydrogen bond
configurations with $\nu>0$ with the DNNP. In classical MD simulations, the
probability distribution $P(d_{OO^{\prime}},\nu)$ for ice Ih is much narrower
compared to that of water. However, the inclusion of NQEs results in a
broadening of the $P(d{OO^{\prime}},\nu)$ distribution for ice, causing it to
resemble that of water. Notably, this broadened distribution exhibits a higher
likelihood of transient autoprotolysis in ice than in water after accounting
for the NQEs.

Figure 4: Scatter plot of the standard deviation in the proton transfer
coordinate ($\sigma_{\nu}$) for all hydrogen bonds in a given water system
versus the corresponding VBM-CBM bandgap for a given frame calculated at the
DFT level (revPBE0 functional) for liquid water (left) and Ice $I_{h}$
(right).

To verify whether the observed change in the structure of the hydrogen bonds
is responsible for the NQEs on Eg, we plot Eg (calculated at the DFT level
with the revPBE0 functional) against the standard deviation of the proton
transfer coordinate $(\sigma_{\nu})$ for each snapshot; this coordinate is a
measure of the spatial fluctuations in the proton transfer coordinates. Our
results for classical and quantum simulations for both water and ice show that
there is a correlation between the value of Eg and the proton fluctuations in
the hydrogen bonds of the system, with a larger $\sigma_{\nu}$ associated with
a lower band-gap.

When comparing water and ice, it is clear that for each system, a larger value
of $\sigma_{\nu}$ correlates with lowered bandgaps. Significant distinctions
in the distributions of $\sigma_{\nu}$ emerge between classical and quantum
simulations for ice, with pronounced variations accentuated when employing the
NEP potential, which yields larger NQEs on the bandgap. In contrast, the
distributions for liquid water exhibit not only considerably smaller
disparities between classical and quantum simulations but also similar trends
for both MLPs potential energy landscapes. These findings confirm that there
exists a correlation between the fluctuations of protons in the hydrogen bond
and the influence of NQEs on the bandgap. Previous investigations on liquid
water Ceriotti  _et al._ (2013) have also reported that NQEs induce
substantial rearrangements in the positions of Wannier centers along a proton
transfer coordinate, underscoring the pivotal role of proton delocalization in
shaping the electronic structure of liquid water. Figure 4 highlights that
NQEs exert an even more substantial impact on hydrogen bond fluctuations in
Ice $I_{h}$, consistent with their greater influence on the bandgap
renormalization.

In summary, we have calculated the finite-temperature bandgap renormalization
of water and ice Ih using many-body perturbation theory on ensembles of
configurations generated using either classical MD or PIMD simulations, and
employing two different machine-learned potentials fitted to DFT forces
obtained with meta-GGA or hybrid functionals. Consistent with earlier
studies,Bischoff  _et al._ (2021) we find that nuclear quantum effects induce
a significant decrease in the electronic bandgap of water, and in ice, the
reduction is even more pronounced, reaching twice the magnitude observed in
water. As a consequence water and Ice $I_{h}$ turn out to have similar
fundamental gaps, consistent with data inferred from experiments. The two sets
of simulations with different machine-learned potentials exhibit the same
trends on nuclear quantum effects, albeit with a few quantitative differences,
stemming from slight distinctions in the structural features of hydrogen bonds
in classical simulations with the two potentials. A statistical analysis of
the differences in the local structure of hydrogen bonds in quantum and
classical simulations suggests that proton delocalization, and the resulting
enhanced transient autoprotolysis, observed in PIMD simulations are the main
cause for the much larger nuclear quantum effects on the bandgap
renormalization of ice Ih, compared to water. These results highlight the
critical importance of accounting for nuclear quantum effects when modeling
the electronic properties of water, and even more so of ice Ih. The observed
influence of quantum proton delocalization and the ensuing bandgap
renormalization are poised to have a considerable impact on molecular
simulation studies examining ice as a solvation medium for environmental
pollutants and as a catalyst for atmospheric chemistry reactions within
snowpacks or clouds, a topic that has garnered increasing interest.Bartels-
Rausch  _et al._ (2014); Bononi  _et al._ (2020); Berrens  _et al._ (2023);
Loerting  _et al._ (2006)

## Conflicts of interest

There are no conflicts to declare.

## Supporting Information

Computational Methods, discussion of differences between MLPs, vibrational
density of states, convergence of number of beads and number of frames used
for electronic structure calculations, discussion of sizes effects on the
proton transfer coordinate and bandgap values, a summary of calculated VBM CBM
gap values, radial and angular distribution functions, histograms of IPR, and
electronic density of states.

## Data availability

Data are available in the Materials Cloud Archive (www.materialscloud.org)
with ID materialscloud:xxxx

## Acknowledgements

We are grateful to Jiawei Zhan and Wenzhe Yu for their assistance with the
G0W0 calculations, and to Yifan Li for his assistance in the proton transfer
coordinate calculations. We thank Francesco Paesani for useful suggestions
about the choice of ML potentials. is This work partly supported by the
National Science Foundation under Grants No. 2305164 and No. 2053235, and by
the Advanced Materials for Energy-Water Systems Center, an Energy Frontier
Research Center funded by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences. Many body perturbation theory calculations were carried
out with the WEST code, whose development is supported by MICCoM, as part of
the Computational Materials Sciences Program funded by the U.S. Department of
Energy, Office of Science, Basic Energy Sciences, Materials Sciences, and
Engineering Division through Argonne National Laboratory, under Contract No.
DE-AC02-06CH11357. G0W0 calculations have been performed using resources of
the Argonne Leadership Computing Facility, which is a DOE Office of Science
User Facility, supported by an ALCC grant. The work at the Lawrence Livermore
National Laboratory was performed under the auspices of the U.S. Department of
Energy under contract DE-AC52-07NA27344. M.F.C.A. and T.A.P. were supported as
part of the Center for Enhanced Nanofluidic Transport, an Energy Frontier
Research Center funded by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences under award DE-SC0019112.

## References

  * Zhong  _et al._ (2019) J. Zhong, M. Kumar, J. Anglada, M. Martins-Costa, M. Ruiz-Lopez, X. Zeng, and J. S. Francisco, Annu. Rev. of Phys. Chem. 70, 45 (2019).
  * Sverjensky  _et al._ (2014) D. A. Sverjensky, B. Harrison, and D. Azzolini, Geochimica et Cosmochimica Acta 129, 125 (2014).
  * Walter  _et al._ (2010) M. G. Walter, E. L. Warren, J. R. McKone, S. W. Boettcher, Q. Mi, E. A. Santori, and N. S. Lewis, Chem. Rev. 110, 6446 (2010).
  * Grätzel (2001) M. Grätzel, Nature 414, 338 (2001).
  * Pham  _et al._ (2014) T. A. Pham, C. Zhang, E. Schwegler, and G. Galli, Phys. Rev. B 89, 060202 (2014).
  * Chen  _et al._ (2016) W. Chen, F. Ambrosio, G. Miceli, and A. Pasquarello, Phys. Rev. Lett. 117, 186401 (2016).
  * Gaiduk  _et al._ (2018) A. P. Gaiduk, T. A. Pham, M. Govoni, F. Paesani, and G. Galli, Nat Commun 9, 247 (2018), number: 1 Publisher: Nature Publishing Group.
  * Bischoff  _et al._ (2021) T. Bischoff, I. Reshetnyak, and A. Pasquarello, Phys. Rev. Res. 3, 023182 (2021).
  * Tal  _et al._ (2024) A. Tal, T. Bischoff, and A. Pasquarello, Proc. Nat. Acad. Sci. 12, e2311472121 (2024).
  * Grannas  _et al._ (2014) A. M. Grannas, L. P. Pagano, B. C. Pierce, R. Bobby, and A. Fede, Environ. Sci. Technol. 48, 10725 (2014).
  * Klanova  _et al._ (2003) J. Klanova, P. Klan, J. Nosek, and I. Holoubek, Environ. Sci. Technol. 37, 1568 (2003).
  * Takenaka and Bandow (2007) N. Takenaka and H. Bandow, J. Phys. Chem. A 111, 8780 (2007).
  * Grannas  _et al._ (2007) A. M. Grannas, A. E. Jones, J. Dibb, M. Ammann, C. Anastasio, H. J. Beine, M. Bergin, J. Bottenheim, C. S. Boxe, G. Carver, G. Chen, J. H. Crawford, F. Dominé, M. M. Frey, M. I. Guzman, D. E. Heard, D. Helmig, M. R. Hoffmann, R. E. Honrath, L. G. Huey, M. Hutterli, H. W. Jacobi, P. Klan, B. Lefer, J. McConnell, J. Plane, R. Sander, J. Savarino, P. B. Shepson, W. R. Simpson, J. R. Sodeau, R. von Glasow, R. Weller, E. W. Wolff, and T. Zhu, Atmos. Chem. Phys. 7, 4329 (2007).
  * Mifsud  _et al._ (2021) D. V. Mifsud, Z. Kaňuchova, P. Herczku, S. Ioppolo, Z. Juhasz, S. T. S. Kovacs, N. J. Mason, R. W. McCullough, and B. Sulik, Space Sci Rev 217, 14 (2021).
  * Jenniskens and Blake (1994) P. Jenniskens and D. F. Blake, Science 265, 753 (1994).
  * Bernas  _et al._ (1997) A. Bernas, C. Ferradini, and J.-P. Jay-Gerin, Chemical Physics 222, 151 (1997).
  * Winter  _et al._ (2004) B. Winter, R. Weber, W. Widdra, M. Dittmar, M. Faubel, and I. V. Hertel, J. Phys. Chem. A 108, 2625 (2004).
  * Ambrosio  _et al._ (2017) F. Ambrosio, G. Miceli, and A. Pasquarello, J. Phys. Chem. Lett. 8, 2055 (2017).
  * Campbell  _et al._ (1979) M. J. Campbell, J. Liesegang, J. D. Riley, R. C. G. Leckey, J. G. Jenkin, and R. T. Poole, J. Phys. Chem. A 15, 83 (1979).
  * Shibaguchi  _et al._ (1977) T. Shibaguchi, H. Onuki, and R. Onaka, J. Phys. Soc. Jpn. 42, 152 (1977).
  * Baron  _et al._ (1978) B. Baron, D. Hoover, and F. Williams, J. Chem. Phys. 68, 1997 (1978).
  * Ceriotti  _et al._ (2016) M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales, and T. E. Markland, Chem. Rev. 116, 7529 (2016).
  * Markland and Ceriotti (2018) T. Markland and M. Ceriotti, Nat. Rev. Chem. 2 (2018).
  * Kundu  _et al._ (2021) A. Kundu, M. Govoni, H. Yang, M. Ceriotti, F. Gygi, and G. Galli, Phys. Rev. Materials 5, L070801 (2021).
  * Kundu  _et al._ (2022) A. Kundu, Y. Song, and G. Galli, Proc. Natl. Acad. Sci. U.S.A. 119, e2203083119 (2022).
  * Alvertis and Engel (2022) A. M. Alvertis and E. A. Engel, Phys. Rev. B 105, L180301 (2022).
  * Kundu and Galli (2023) A. Kundu and G. Galli, J. Chem. Theory Comput. 19, 4011 (2023), https://doi.org/10.1021/acs.jctc.3c00424 .
  * Kundu and Galli (2024) A. Kundu and G. Galli, J. Phys. Chem. Lett. 15, 802 (2024), https://doi.org/10.1021/acs.jpclett.3c03269 .
  * Berne and Thirumalai (1986) B. J. Berne and D. Thirumalai, Annu. Rev. Phys. Chem. 37, 401 (1986).
  * Marx and Parrinello (1996) D. Marx and M. Parrinello, J. Chem. Phys. 104, 4077 (1996).
  * Herrero and Ramírez (2014) C. P. Herrero and R. Ramírez, J. Phys. Condens. Matter 26, 233201 (2014).
  * Huppert  _et al._ (2022) S. Huppert, T. Plé, S. Bonella, P. Depondt, and F. Finocchi, Applied Sciences 12 (2022), 10.3390/app12094756.
  * Sabatini  _et al._ (2013) R. Sabatini, T. Gorni, and S. de Gironcoli, Phys. Rev. B 87, 041108 (2013).
  * Engel  _et al._ (2015) E. A. Engel, B. Monserrat, and R. J. Needs, J. Chem. Phys. 143, 244708 (2015).
  * Monserrat  _et al._ (2015) B. Monserrat, E. A. Engel, and R. J. Needs, Phys. Rev. B 92, 140302 (2015).
  * Perdew  _et al._ (1996a) J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105, 9982 (1996a).
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J Chem Phys 132, 154104 (2010).
  * Sun  _et al._ (2015) J. Sun, A. Ruzsinszky, and J. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Fan  _et al._ (2021) Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen, and T. Ala-Nissila, Phys. Rev. B 104, 104309 (2021).
  * Chen  _et al._ (2024) Z. Chen, M. L. Berrens, K.-T. Chan, Z. Fan, and D. Donadio, J. Chem. Eng. Data 69, 128 (2024).
  * Cheng  _et al._ (2019) B. Cheng, E. Engel A., J. Behler, and M. Ceriotti, Proc. Natl Acad. Sci. U.S.A. 116, 1110 (2019).
  * Schran  _et al._ (2020) C. Schran, K. Brezina, and O. Marsalek, J. Chem. Phys. 153, 104105 (2020).
  * Zhang  _et al._ (2021) L. Zhang, H. Wang, R. Car, and W. E, Phys. Rev. Lett. 126, 236001 (2021).
  * M. Sommers  _et al._ (2020) G. M. Sommers, M. F. C. Andrade, L. Zhang, H. Wang, and R. Car, Phys. Chem. Chem. Phys. 22, 10592 (2020).
  * Ohto  _et al._ (2019) T. Ohto, M. Dodia, J. Xu, S. Imoto, F. Tang, F. Zysk, T. D. Kühne, Y. Shigeta, M. Bonn, X. Wu, and Y. Nagata, J. Phys. Chem. Lett. 10, 4914 (2019).
  * Palos  _et al._ (2022) E. Palos, E. Lambros, S. Swee, J. Hu, S. Dasgupta, and F. Paesani, J. Chem. Theory Comput. 18, 3410 (2022).
  * Bore and Paesani (2023) S. Bore and F. Paesani, Nat. Commun. 14 (2023), https://doi.org/10.1038/s41467-023-38855-1.
  * Ceriotti and Manolopoulos (2012) M. Ceriotti and D. E. Manolopoulos, Phys. Rev. Lett. 109, 100604 (2012).
  * Kapil  _et al._ (2019) V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, B. A. Helfrecht, P. Juda, S. P. Bienvenue, W. Fang, J. Kessler, I. Poltavsky, S. Vandenbrande, J. Wieme, C. Corminboeuf, T. D. Kühne, D. E. Manolopoulos, T. E. Markland, J. O. Richardson, A. Tkatchenko, G. A. Tribello, V. Van Speybroeck, and M. Ceriotti, Comput. Phys. Commun. 236, 214 (2019).
  * Kapil  _et al._ (2016) V. Kapil, J. Behler, and M. Ceriotti, J. Chem. Phys. 145, 234103 (2016), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4971438/15520102/234103_1_online.pdf .
  * Engel  _et al._ (2016) E. A. Engel, B. Monserrat, and R. J. Needs, J. Chem. Phys. 145, 044703 (2016).
  * Perdew  _et al._ (1996b) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996b).
  * Govoni and Galli (2015) M. Govoni and G. Galli, J. Chem. Theory Comput. 11, 2680 (2015).
  * Medders and Paesani (2014) V. Medders, G.R. Babin and F. Paesani, J. Chem. Theory Comput. 10, 2906–2910 (2014).
  * Babin  _et al._ (2014) V. Babin, G. R. Medders, and F. Paesani, J. Chem. Theory Comput. 10, 1599–1607 (2014).
  * Babin  _et al._ (2013) V. Babin, C. Leforestier, and F. Paesani, J. Chem. Theory Comput. 9, 5395–5403 (2013).
  * Cole  _et al._ (2016) W. T. S. Cole, J. D. Farrell, D. J. Wales, and R. J. Saykally, Science 352, 1194 (2016).
  * Moberg  _et al._ (2017) D. R. Moberg, S. C. Straight, C. Knight, and F. Paesani, J. Phys. Chem. Lett. 8, 2579–2583 (2017).
  * Medders and Paesani (2016) G. R. Medders and F. Paesani, J. Am. Chem. Soc. 138, 3912–3919 (2016).
  * Ceriotti  _et al._ (2013) M. Ceriotti, J. Cuny, M. Parrinello, and D. E. Manolopoulos, Proc. Natl. Acad. Sci. U.S.A. 110, 15591 (2013).
  * Fritsch  _et al._ (2014) S. Fritsch, R. Potestio, D. Donadio, and K. Kremer, J. Chem. Theory Comput. 10, 816 (2014).
  * Bartels-Rausch  _et al._ (2014) T. Bartels-Rausch, H.-W. Jacobi, T. F. Kahan, J. L. Thomas, E. S. Thomson, J. P. D. Abbatt, M. Ammann, J. R. Blackford, H. Bluhm, C. Boxe, F. Domine, M. M. Frey, I. Gladich, M. I. Guzman, D. Heger, T. Huthwelker, P. Klan, W. F. Kuhs, M. H. Kuo, S. Maus, S. G. Moussa, V. F. McNeill, J. T. Newberg, J. B. C. Pettersson, M. Roeselova, and J. R. Sodeau, Atmos. Chem. Phys. 14, 1587 (2014).
  * Bononi  _et al._ (2020) F. C. Bononi, Z. Chen, D. Rocca, O. Andreussi, T. Hullar, C. Anastasio, and D. Donadio, J. Phys. Chem. A 124, 9288 (2020).
  * Berrens  _et al._ (2023) M. L. Berrens, Z. Chen, K.-T. Chan, C. Anastasio, and D. Donadio, ACS Earth Space Chem. 7, 1761 (2023).
  * Loerting  _et al._ (2006) T. Loerting, A. F. Voegele, C. S. Tautermann, K. R. Liedl, L. T. Molina, and M. J. Molina, J. Geophys. Res. 111 (2006), 10.1029/2006JD007065.