# Thermal dependence of the hydrated proton and optimal proton transfer in the protonated water hexamer

Félix Mouhat, Matteo Peria, Tommaso Morresi, Rodolphe Vuilleumier, Antonino Marco Saitta, Michele Casula

► To cite this version:

Félix Mouhat, Matteo Peria, Tommaso Morresi, Rodolphe Vuilleumier, Antonino Marco Saitta, et al. Thermal dependence of the hydrated proton and optimal proton transfer in the protonated water hexamer. Nature Communications, 2023, 14 (1), pp.6930. 10.1038/s41467-023-42366-4. hal-04291382

HAL Id: hal-04291382

https://hal.science/hal-04291382v1

Submitted on 17 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Thermal dependence of the hydrated proton and optimal proton transfer in the protonated water hexamer

Received: 12 December 2022

Accepted: 25 September 2023

Published online: 30 October 2023

<figure><img src="image_3.png" /></figure>

Check for updates

Félix Mouhat Ⓔ<fnref n="1" />, Matteo Peria Ⓔ², Tommaso Morresi³, Rodolphe Vuilleumier Ⓔ⁴, Antonino Marco Saitta Ⓒ² & Michele Casula Ⓒ²✉

Water is a key ingredient for life and plays a central role as solvent in many biochemical reactions. However, the intrinsically quantum nature of the hydrogen nucleus, revealing itself in a large variety of physical manifestations, including proton transfer, gives rise to unexpected phenomena whose description is still elusive. Here we study, by a combination of state-of-the-art quantum Monte Carlo methods and path-integral molecular dynamics, the structure and hydrogen-bond dynamics of the protonated water hexamer, the fundamental unit for the hydrated proton. We report a remarkably low thermal expansion of the hydrogen bond from zero temperature up to 300 K, owing to the presence of short-Zundel configurations, characterised by proton delocalisation and favoured by the synergy of nuclear quantum effects and thermal activation. The hydrogen bond strength progressively weakens above 300 K, when localised Eigen-like configurations become relevant. Our analysis, supported by the instanton statistics of shuttling protons, reveals that the near-room-temperature range from 250 K to 300 K is optimal for proton transfer in the protonated water hexamer.

For more than 200 years and the seminal work of von Grotthus, the properties of the hydrated proton H⁺ₐₐql have intrigued the scientific community<fnref n="1" />,². Despite significant advances, the exact role of the solvated proton in proton transfer (PT) reactions in chemical and biological systems is not fully elucidated yet. The textbook picture is that the hydrated proton exists as classical hydronium cation H₃O⁺, but it looks more appropriately described as a delocalised charge defect shared by multiple molecules. The spread of this charge defect blurs the identity of the excess proton between two limiting structures, namely the Zundel³ and the Eigen⁴ ions. Indeed, the hydrated proton Infrared (IR) spectrum displays a combination of a few discrete absorption bands on top of an absorption continuum, broadly extended over the entire spectrum. Neither the symmetrically solvated hydronium H₉O₄⁺ (Eigen) ion, nor the equally shared proton in the H₅O₂⁺ (Zundel) ion, can individually rationalise this characteristic IR

fingerprint. Models involving fast inter-conversions between these two ionic species are also shown to fail⁵. Moreover, the question of whether the hydrated excess proton and the related PT mechanism should be based on an Eigen-⁶⁻⁸ or Zundel-like motif⁹,¹⁰, one or the other taken as the most dominant species, has been actively debated.

To deal with these issues, Stoyanov et al.¹¹ have introduced the stable H₁₃O₆⁺ species, the protonated water hexamer, which is Zundel-type in the sense that the excess proton is equally shared between two water molecules. The core of the cluster is characterised by a central oxygen-oxygen distance dO₁O₂ that, however, is more elongated than in the Zundel cation¹². On the other hand, recent Molecular Dynamics (MD) simulations suggest the existence of a distorted, nonsymmetric Eigen-type cation, remaining at the heart of a dynamical charge defect spanning multiple water molecules⁶,¹³. The protonated water hexamer represents one of the smallest protonated water clusters for which

<footnote n="1">¹Saint Gobain Research Paris, 39, Quai Lucien Lefranc, 93300 Aubervilliers, France. ²IMPMC, Sorbonne Université, CNRS, MNHN, UMR 7590, 4 Place Jussieu, 75252 Paris, France. ³ECT*-Fondazione Bruno Kessler*, 286 Strada delle Tabarelle, 38123 Trento, Italy. ⁴PASTEUR, Département de Chimie, École normale supérieure, PSL Research University, Sorbonne Université, CNRS, 24 Rue Lhomond, 75005 Paris, France. ✉e-mail: michele.casula@sorbonne-universite.fr</footnote>
both of these characteristic binding motifs coexist¹⁴⁻¹⁶. Snapshots of the main protonated hexamer configurations are represented in Fig. 1, for both Zundel- and Eigen-like forms, while other low-lying isomers exist¹⁷,¹⁸, but with a much lower probability to occur starting from the global minimum, due to their different topology.

The protonated hexamer Potential Energy Surface (PES) has been partially explored by IR spectroscopy¹⁹⁻²¹ and electronic structure calculations performed within Empirical Valence Bond (EVB), Density Functional Theory (DFT), experiment-directed DFT, Møller-Plesset (MP2), and Coupled Cluster (CC) approaches, also supplemented by Machine Learning (ML) techniques⁶,⁷,⁹,¹⁵,¹⁸⁻²⁰,²²⁻²⁴, confirming that the two structures introduced above are the lowest energy isomers. The quantum nature of the proton however induces a delocalised structure on this PES.

In this work, we apply MD simulations fully retaining the nuclear quantum nature of the atoms. In particular, the quantum proton, described within the Feynman path integral (PI) approach, evolves in a very accurate PES estimated by means of Quantum Monte Carlo (QMC). This stochastic technique introduces an intrinsic noise, which affects the forces driving the ion dynamics and, consequently, the simulation temperature. Relying on the generalised fluctuation-dissipation theorem, a Langevin-based approach has been developed to address this issue for classical²⁵ and quantum²⁶ ions. It allows one to sample microscopic configurations in the canonical ensemble with a high level of quantum accuracy. Although our QMC wave function ansatz leads to the dissociation limit of two water molecules ≈ 20% (≈1 kcal/mol) off with respect to CCSD(T) values, the QMC PES is sufficiently accurate to yield structural properties and quantum-thermal distributions around equilibrium in a good agreement with CCSD(T)-derived PESs for benchmark systems such as the water dimer²⁷ and the Zundel ion²⁸. Details of the method and of its accuracy are provided in the 'Methods' section and in Supplementary Notes I and II of the Supplementary Information (SI).

We find that the hydrogen bond (H-bond) mediated by the hydrated proton shows a remarkably low thermal expansion from zero temperature up to 300 K, with a nearly temperature-independent length that becomes shorter than the classical-ion counterpart in the [200-350] K temperature range. A non-trivial behaviour of the H-bond has also been found in H-rich crystals and ferroelectric materials, such as the potassium dihydrogen phosphate (KDP)²⁹, first detected by Ubbelohde in 1939 upon isotopic substitution³⁰. In the latter case, the lighter the hydrogen, the shorter the H-bond. This was interpreted as a quantum manifestation of proton delocalisation, strengthening the H-bond. In the present situation, the strength of the H-bond results from a non-trivial cooperation of nuclear quantum effects (NQEs) and thermal activation, as we will show in this work. Indeed, NQEs strongly affect the

vibrational levels of the proton shuttling mode bridging the central O₁ and O₂ oxygen atoms. These levels are then thermally occupied according to the $d_{O_1O_2}$ distance of a given configuration. We can thus distinguish three regimes (see Fig. 1): (i) “short-Zundel” configurations with the shortest $d_{O_1O_2}$, where the proton along the shuttling mode feels a quadratic potential close enough to its energy minimum and it is perfectly shared between the two central water molecules; (ii) “elongated-Zundel” configurations for intermediate $d_{O_1O_2}$, comprising the equilibrium distance, where a potential energy barrier starts to develop in between O₁ and O₂ and the proton is delocalised only due to NQEs; (iii) “distorted-Eigen” configurations at even larger $d_{O_1O_2}$, where the central barrier is large enough that the hydrated proton is localised on one of the two flanking water molecules, forming an Eigen-like complex.

Here, we show that the occurrence of short-Zundel configurations is key to understand the H-bond thermal robustness and to enhance the proton transfer dynamics. Despite being energetically disfavoured by the short $d_{O_1O_2}$ distances at the classical level, these configurations are populated thanks to the synergistic action of NQEs and temperature, yielding a sweet spot for proton transfer in the [250-300] K temperature range.

## Results

### Thermal expansion of the H-bond

To rationalise our main outcome, we first study the zero-temperature classical geometry and PES of the protonated water hexamer, and compare it with the Zundel cation. While the latter system misses a large part of water solvation effects, the former includes the full contribution of the first and second shells of the solvated proton. The zero-temperature results are reported in Supplementary Note III.1 of the SI. We simply highlight here that the Variational Monte Carlo (VMC) equilibrium $O_1-O_2$ distance is found to be $d_{O_1O_2} = d_{min} = 2.3930(5)$ Å, in good agreement with MP2 calculations, the most widely used post Hartree-Fock theory to study water clusters (see Supplementary Note I.2 of the SI for a more extended comparison between VMC and MP2). VMC has a milder scale with the system size than MP2, allowing one to perform extensive calculations of the protonated hexamer. At variance with the Zundel cation³¹,³², the protonated hexamer equilibrium geometry is asymmetric, implying that the global minimum is split into a double well, separated by an energy barrier between the two central water molecules. This barrier vanishes at $d_{O_1O_2} = d_{symm} \simeq 2.38$ Å, a distance that separates the short Zundel below from the elongated-Zundel configurations above. The height of the barrier is less than 100 K (in $k_B$ units) at $d_{min}$, rapidly increasing as a function of $d_{O_1O_2}$. We therefore expect several consequences on the hydrated proton distribution and on its mobility at finite temperature, once the NQEs are taken into account.

<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Different regimes of the protonated water hexamer $H_{13}O_4^+$, a Short-Zundel configuration with a Zundel center ($H_5O_2^+$) in colours and its first solvation shell (4 H₂O) in grey shades. **b** elongated Zundel with the quantum nature of hydrogen atoms highlighted by the full representation of its imaginary-time</figcaption></figure>

positions in a PI configuration. **c** distorted-Eigen configuration with an Eigen cation ($H_9O_4^+$) in colours accompanied by two solvating water molecules (2 H₂O) in grey shades. The O₁, O₂ and H' labels are used throughout the paper to refer to the corresponding atoms, as indicated here.
<figure><img src="image_5.png" /><figcaption>**Fig. 2 |** Classical and quantum oxygen-oxygen *g*<sub>O<sub>1</sub>O<sub>2</sub></sub> pair correlation functions as a function of temperature. **a, b** The dashed vertical lines indicate the average 〈d<sub>O<sub>1</sub>O<sub>2</sub></sub>〉 distance for each simulation, at the corresponding temperature. The dotted vertical line is located at the classical equilibrium geometry. **Panel (c)** shows</figcaption></figure>

the T-dependence of the 〈d<sub>O<sub>1</sub>O<sub>2</sub></sub>〉 average distance. The classical equilibrium geo-
metry is represented by a short-dashed horizontal black line. At 250 K and 300 K
the oxygen-oxygen distance is shortened by NQEs with respect to the classical
counterpart. Source data are provided as a Source data file.

To understand how the dynamics of the hydrated proton evolves
with temperature, QMC-driven ab initio MD simulations are relevant.
Such calculations are carried out for both classical and quantum nuclei
of the H<sub>13</sub>O<sub>6</sub><sup>+</sup> ion, within the temperature interval T ∈ [50–350] K,
thanks to the methodological developments detailed in ref. 26 and in
the ‘Methods’ section. At these conditions, the clusters are stable
during the simulated time frame (≈30 ps), allowing us to access the
thermal properties of the hydrated proton and the O<sub>1</sub>H<sup>+</sup>O<sub>2</sub> bond over
an extended temperature range. This is the advantage of performing
accurate computer simulations for this system, particularly if an
experimental investigation of the cluster will require a more extended
lifetime, with the risk of molecular evaporation at the highest tem-
peratures studied here. The details of these challenging calculations,
such as the number of iterations and their computational costs at a
given temperature, are reported in Supplementary Table 4 of the SI.

From our QMC-MD simulations, we extract the normalised Pair Correlation Functions (PCFs) g<sub>O<sub>1</sub>O<sub>2</sub></sub>, for the two oxygen atoms O<sub>1</sub> and O<sub>2</sub> of the cluster core (Fig. 2). The expected broadening of the PCFs due to nuclear quantisation is significant over the whole temperature range (Fig. 2b). Only at temperatures as high as 350 K, the classical g<sub>O<sub>1</sub>O<sub>2</sub></sub> (Fig. 2a) starts resembling the quantum distribution. This implies that the NQEs cannot be neglected for temperatures up to this value, above ambient conditions. We also notice that, when comparing to the Zundel ion results<sup>33</sup>, the peak position is shifted up by at least ~0.01 Å. Thus, it appears that the H<sub>13</sub>O<sub>6</sub><sup>+</sup> cluster frequently adopts elongated-Zundel configurations<sup>8,34,35</sup> at the lowest temperatures considered here. This is at variance with the protonated water dimer, where the hydrated proton lives in a single minimum symmetrically located between the two water molecules.

Focusing our attention to 〈d<sub>O<sub>1</sub>O<sub>2</sub></sub>〉 (Fig. 2c), its classical and quantum behaviours are remarkably different as a function of temperature. On the one hand, the classical d<sub>O<sub>1</sub>O<sub>2</sub></sub> keeps increasing with temperature, as more energy is given to the intermolecular vibration modes. On the other hand, the quantum d<sub>O<sub>1</sub>O<sub>2</sub></sub> displays a nearly flat behaviour with the cluster temperature, up to 300 K. This very low thermal expansion extended over a wide temperature range leads to a temperature regime where d<sub>O<sub>1</sub>O<sub>2</sub></sub> for the quantum system become shorter than the classical values at the same temperatures. This is clearly seen in Fig. 2c. We will come back to this point later.

Finally, as the temperature further increases, the NQEs reduction
weaken the central H-bond strength. Consequently, d<sub>O<sub>1</sub>O<sub>2</sub></sub> spreads
out, due to stochastic fluctuations of the core and the solvent, and a

more classical regime is reached, when the averaged d<sub>O<sub>1</sub>O<sub>2</sub></sub> values for
classical and quantum nuclei meet again. The PCF distributions display
longer tails, with more configurations covering regions with d<sub>O<sub>1</sub>O<sub>2</sub></sub> ∈
[2.5 – 2.7] Å, and the peak position rapidly shifts to larger values.
Configurations with such a large 〈d<sub>O<sub>1</sub>O<sub>2</sub></sub>〉 are of distorted-Eigen type<sup>8,36</sup>.

A cooperative thermal-quantum species: the short-Zundel ion

To refine our structural analysis, we compute the bidimensional dis-
tribution function ρ<sub>2D</sub>, which correlates the oxygen-oxygen (O<sub>1</sub>O<sub>2</sub>)
and the oxygen-proton (O<sub>1</sub>/2H<sup>+</sup>) distances, and study its temperature
dependence ρ<sub>2D</sub> = ρ<sub>2D</sub>(T). In Fig. 3, we show the contour plot of the
temperature-driven ρ<sub>2D</sub> variation (see also Supplementary Note III.2 of
the SI). By taking ρ<sub>2D</sub>(250 K) as reference, four temperature variations
are explored: 100 K, 200 K, room temperature (RT), and 350 K (from
the top to the bottom of Fig. 3).

In the classical protonated hexamer (Fig. 3, left column), rising the temperature from 250 K up to 350 K tends to stretch 〈d<sub>O<sub>1</sub>O<sub>2</sub></sub〉, by promoting configurations from the elongated Zundel (blue central distribution with d<sub>O<sub>1</sub>O<sub>2</sub></sub> ∈ [2.38,2.5] Å in Fig. 3) to an Eigen-like arrangement with larger d<sub>O<sub>1</sub>O<sub>2</sub></sub> and a proton much more localised on one of the two central oxygen atoms (red wings). The situation is reversed at lower temperatures (100 K and 200 K) if compared to the 250 K reference, with positive (red) variations in the elongated Zundel and negative (blue) variations in the wings. Thus, for classical nuclei, there is a progressive depletion of the elongated Zundel and a corresponding population of the distorted-Eigen wings upon temperature rise. Short-Zundel configurations, highlighted in Fig. 3 by a grey background, seem to play a very marginal role in the temperature-driven density distribution shift.

The scenario is strikingly different with quantum nuclei (right column), particularly at the lowest temperatures (100 K and 200 K). In this regime, distorted-Eigen configurations are barely populated or depleted, and the density shift upon rising temperature takes place between the elongated-Zundel region and the short-Zundel sector. The latter is significantly more populated at 250 K than at lower temperatures at the expense of the elongated Zundel, which instead loses density with respect to the classical counterpart at the same temperature.

In the higher-temperature limit, at 350 K, NQEs are less relevant
and, by consequence, the classical and quantum variations have a
qualitatively similar behaviour. In both classical and quantum case, we
notice the presence of red wings at large oxygen-oxygen distances
Classical particles

<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Bidimensional oxygen-oxygen/oxygen-proton distributions. Difference between bidimensional oxygen-oxygen/oxygen-proton distributions $\rho_{2D}$ obtained by QMC-driven LD simulations for classical (left panels) and quantum (right panels) particles, computed at different temperatures. The bidimensional distribution computed at 250 K is taken as reference. Positive (negative) regions are</figcaption></figure>

Quantum particles

<figure><img src="image_7.png" /></figure>

in red (blue) colour. The black filled circles correspond to the zero-temperature equilibrium geometries of the $H_{13}O_6^+$ ion at a fixed $d_{O_1O_2}$ distance. The coloured background highlights the three different regimes explained in the paper: the short Zundel (grey), the elongated Zundel (yellow), and the distorted Eigen (green) species. Source data are provided as a Source Data file.

($d_{O_1O_2} \in [2.5,2.7]$ Å), which are the signature of thermally activated Eigen-like states, with a strongly localised proton. This is related to less frequent elongated-Zundel configurations, indicated by the depleted distribution for $d_{O_1O_2} < 2.5$ Å, confirming that the distorted-Eigen configurations are indeed promoted by high temperature. For quantum nuclei, the corresponding depletion goes well below the elongated-Zundel region, by touching also short-Zundel configurations, down to

$d_{O_1O_2} \sim 2.3$ Å, at variance with the classical case, where the short-Zundel configurations are not involved.

To interpret these results, we first construct an accurate effective potential by projecting the full PES, computed during QMC-driven classical MD calculations, onto the degrees of freedom mostly relevant to understand the dynamics of the hydrated proton. These are the $d_{O_1O_2}$ distance and the proton sharing coordinate $\delta$, referenced to the
midpoint of the O₁H⁺O₂ complex: $\delta \equiv \bar{d}_{O_{1/2}H^+} - d_{O_1O_2}/2$, with $\bar{d}_{O_{1/2}H^+}$ the O₁/₂-H⁺ distance projected onto the O₁O₂ direction. The resulting two-dimensional (2D) potential is $V_{2D} = V_{2D}(d_{O_1O_2}, \delta)$. We refer the reader to Supplementary Notes IV and V of the SI for technical details about the PES projection. We highlight that the potential $V_{2D}$ is derived here at VMC quality. We also notice that $\delta$ is the vibrational coordinate of the proton shuttling mode, while $d_{O_1O_2}$ is related to the stretching mode of the two water molecules in the cluster core.

Given $V_{2D}(d_{O_1O_2}, \delta)$, we then proceed to quantize the variable $\delta$. Indeed, while $d_{O_1O_2}$ can be taken as classical, for it is related to the motion of heavier oxygen atoms of mass $m_O$, the $\delta$ coordinate must be quantised, owing to the light mass ($m_H$) of the hydrated proton. At the leading order in $2m_H/(m_O + m_H)$, we separate the stretching mode from the shuttling one, by invoking an adiabatic Born-Oppenheimer type of approximation for the two species³⁷. We finally solve quantum-mechanically the Hamiltonian of a proton in the potential $V_\delta = V_{2D}(\alpha, \delta)|_{\alpha=d_{O_1O_2}}$ at fixed $d_{O_1O_2}$ value. In Fig. 4a–c we plot the ground state distribution and eigenvalues obtained for three distances, i.e. at $d_{O_1O_2} = 2.375$ Å, in the short-Zundel region close to the boundary between the short and the elongated Zundel, at $d_{O_1O_2} = 2.495$ Å, in the elongated-Zundel region close to the frontier between the elongated Zundel and the distorted Eigen, and finally at $d_{O_1O_2} = 2.585$ Å, deep into the distorted Eigen regime.

One can notice three different quantum behaviours of the vibrational shuttling mode, that provide a more quantitative ground to the three-regime distinction made at the beginning. In the short Zundel, $V_\delta$ is indeed a quadratic potential with a single minimum at the core center, which widens as $d_{O_1O_2}$ gets close to $d_{symm} \approx 2.38$ Å, a distance where it becomes quartic because its curvature falls to zero before

changing sign. The ground state energy, i.e. the zero point energy (ZPE) of the shuttling mode, decreases as the potential widens, as reported in Fig. 4d. In the elongated Zundel, a central barrier starts to develop, with a ground-state proton distribution that stays uni-modal thanks to a ZPE larger than its height, till $d_{O_1O_2} \simeq 2.5$ Å, where the ZPE equals the barrier height. In this regime, for $d_{O_1O_2} \in [d_{symm}, 2.5$ Å], the ZPE is particularly small, due to the quartic nature of $V_\delta$, and weakly $d_{O_1O_2}$-dependent, as shown in Fig. 4d. Finally, for $d_{O_1O_2} > 2.5$ Å, we enter the distorted-Eigen regime, with an even larger central barrier (>1000 K), such that the quantum proton is instantaneously localised in one of the two wells, and its distribution is then bimodal. The ZPE starts to rise again as $d_{O_1O_2}$ is stretched, with a slope steeper–in absolute value–than the ZPE decrease in the short Zundel, because it is now set by the much deeper lateral minima of the double-well potential. This can be seen again in Fig. 4d.

We can now correct the classical O₁-O₂ potential, defined as $V_{O_1O_2} \equiv V_{2D}(d_{O_1O_2}, \delta)|_{\delta=\delta_{min}}$, where $\delta_{min}$ is the V₂D minimum at fixed $d_{O_1O_2}$ value, by adding the ZPE $\forall d_{O_1O_2}$, obtained from the quantisation of the shuttling mode $\delta$. The resulting potential is plotted in Fig. 4e. Remarkably, the anharmonic classical $V_{O_1O_2}$ potential becomes harmonic after ZPE-correction. It is a consequence of the much larger ZPE in the distorted-Eigen configurations than in the short Zundel, which compensates for the underlying $V_{O_1O_2}$ anharmonicity. This rationalises two main features. On the one hand, it explains the very low thermal expansion of $\langle d_{O_1O_2} \rangle$, being the average position in a harmonic potential temperature-independent. On the other hand, it proves that NQEs enhance the occurrence of short-Zundel configurations upon heating, while the distorted Eigen is penalised by its large ZPE with respect to the classical counterpart. The enhancement of the

<figure><img src="image_8.png" /><figcaption>**Fig. 4** | NQEs on the shuttling mode, and their impact on the O₁-O₂ interatomic potential $V_{O_1O_2}$. We quantize the proton shuttling mode $\delta$, defined as the displacement along the segment connecting the two oxygen atoms in the core of the cluster from its midpoint position. We study the ground-state wave function and the first 5 eigenvalues for the confining potential $V_\delta$, as a function of $d_{O_1O_2}$. Panels (a), (b) and (c) report the ground state wave function and the lowest 5 energy levels for $d_{O_1O_2} = 2.375$, 2.495 and 2.585 Å, respectively. In panel (d), the variation of the</figcaption></figure>

zero-point (ground-state) energy (ZPE) as a function of $d_{O_1O_2}$ is explicitly plotted. While the ZPE dependence is very flat in the elongated-Zundel region (depicted by the yellow shaded area), the ZPE increases in both short-Zundel (grey shaded area) and distorted-Eigen (green shaded area) regions, with a much steeper slope in the latter. In panel (e), the ZPE is added to the classical interatomic potential $V_{O_1O_2}$ (solid blue line) to yield the quantum-corrected effective interatomic potential (solid dark-pink line) between the two inner oxygen atoms.
occurrence of Zundel configurations by NQEs is also revealed by the population analysis presented in Supplementary Note VI and Supplementary Fig. 17 of SI. This is also in agreement with a similar analysis carried out in ref. 6 in bulk water.

Above RT, the distorted Eigen configurations will eventually become dominant again. This can be understood within this framework as well. Indeed, thermal excitations are energetically more available in the distorted Eigen, where the spacing between the ZPE and the first-excited state shrinks, and higher excited states are piled up more densely than in the short and elongated Zundel (see Fig. 4a-c). The full thermal dependence of the short Zundel, elongated Zundel and distorted Eigen populations is reported in Supplementary Fig. 16 of Supplementary Note VI in SI.

## Optimal proton transfer from instantons statistics

The analysis made so far highlights the paramount importance of the NQEs to set the non-trivial temperature behaviour of the $H_{13}O_6^+$ cluster. At this stage, direct information about the excess proton dynamics along the QMC-PIMD trajectory is necessary to estimate more quantitatively its impact on the PT processes occurring in the system.

One way to achieve this goal is by analysing the statistics of selected transition-state (TS) configurations, defined by means of instanton theory. Within the PI formalism, the instanton path seamlessly connects the reactants and products minima, along the minimal action trajectory, periodic in the quantum imaginary time $\tau = \beta\hbar^{38}$. It provides a generalisation of the TS theory for anharmonic quantum systems$^{39}$, and it has been very recently applied in a QMC framework$^{40,41}$, by efficiently recovering the proper scaling of ground-state tunnelling rates. TS configurations are therefore identified as those where each half of the instanton path is located on either side of the central $O_1O_2$ midpoint, sampled during the QMC-PIMD dynamics.

With the aim at resolving the contribution of the three different regimes to the PT dynamics, we collect the instanton events and compute their statistical distribution as a function of $d_{O_1O_2}$. We plot the instanton density distribution function in Fig. 5a at various temperatures. To deepen our analysis, we compute also the cumulative density distribution function in Fig. 5b, after normalising it based on the algorithmic frequency of the instanton occurrences, as counted during our QMC-PIMD simulations. Although this does not give direct access to real-time quantities, the ring-polymer MD with Langevin

thermostat has been shown to yield physically reliable information on frequencies and frequency variations$^{42}$. Note that the coupling with the Langevin thermostat is kept constant across the full temperature range analysed here$^{42}$. The fully integrated frequency distribution gives the total proton hopping frequency, plotted in Fig. 5c as a function of temperature. This shows a clear maximum located in the [250–300] K temperature range. Consequently, we expect the hydrated proton mobility to be optimal in a near-RT window, with a maximised Grotthus diffusion. To understand the source of this temperature sweet spot, in the same panel (c) we plot the contribution to the total frequency of instanton events occurring in the short-Zundel region. This is yielded by the cumulative frequency distribution of panel (b) evaluated at the boundaries between short and elongated Zundel, i.e. at $d_{O_1O_2} = d_{symm}$. The short-Zundel contribution to the total frequency shows a peak of the same intensity as the total one in the same temperature range, clearly pointing to the key role played by thermally activated short-Zundel configurations to the PT dynamics. The short-Zundel arrangement enables instantaneous proton jumps between the two sides of the cation, since there is no barrier to cross. Thus, the sweet spot constitutes the best compromise between acquiring enough thermal energy to access short-distance configurations, boosted by NQEs, and controlling the amplitude of the chemical (covalent or H-) bonds fluctuations, that might trap the proton into an asymmetric well. Indeed, at larger temperatures (>300 K), the onset of distorted-Eigen and the corresponding fall of short-Zundel configurations localize the hydrated proton around its closest oxygen atom, thus reducing its shuttling probability. A similar non-monotonous PT behaviour has experimentally been found in bulk water by assessing the limiting conductivities of the $H_3O^+$ and $D_3O^+$ species$^{43}$. Thanks to these measures, performed at 20 MPa, the excess molar conductivities due to PT have been estimated. They show a peak located at a temperature in between 420 K and 430 K. In this temperature range and at the pressure conditions of the experiment, the water density is only 7–8% smaller than the standard conditions$^{44}$, a regime comparable to the one of our cluster.

Beside this PT mechanism, which is adiabatic in nature and driven by the synergy of ZPE and thermal effects, NQEs could also contribute to the proton diffusion by means of instantaneous tunnelling, which can further accelerate the PT dynamics. By computing the root-mean-square (RMS) displacement correlation

<figure><img src="image_9.png" /><figcaption>**Fig. 5** | Instanton statistics and proton hopping frequency. a Instanton distribution resolved as a function of the $d_{O_1O_2}$ distance for different temperatures. b Cumulative distribution of (a) normalised by the occurrence frequency of the instanton (proton hopping) events during the PIMD simulations. c Proton hopping frequency as a function of temperature, together with the contribution coming</figcaption></figure>

from the short Zundel configurations, with $d_{O_1O_2} < d_{symm} = 2.38$ Å. The $d_{symm}$ value is reported as vertical dashed line in (a) and (b). Here, we report simulations performed also at 400 K, a temperature at which the cluster is still stable or metastable. Source data are provided as a Source data file.
functions⁴⁵ over the instantons population, we verified that tunneling events could take place only in the distorted Eigen and in the intermediate temperature range (see Supplementary Note VII of SI for a detailed analysis). This additional PT channel has however a marginal effect with respect to the main mechanism unveiled here. Indeed, Fig. 5c shows that the sweet spot is mainly due to PT events originating in short-Zundel configurations, where quantum tunneling is not relevant.

## Discussion

Using highly accurate QMC-PIMD simulations of the $H_{13}O_6^+$ cation at finite temperature, we found a remarkably low thermal expansion of the protonated water hexamer core. It stems from a cooperative action of both NQEs and thermal effects, which leads to the emergent behaviour of short-Zundel species as PT booster, where the excess proton is perfectly shared between two neighbouring water molecules. The relevance of short-Zundel configurations is enhanced by NQEs, which instead penalize the distorted-Eigen states, having a larger ZPE. In the intermediate temperature range, comprising RT, the occurrence of short-Zundel events is maximised by thermal population, leading to a sweet spot in the PT dynamics. Around these temperatures, distorted-Eigen states can still contribute to PT with quantum tunnelling processes, although occurring at much lower rates. The cluster core spreads out again at larger temperatures, as soon as stronger thermal fluctuations favour the formation of more classical distorted-Eigen structures, where the proton gets strongly localised in one of the flanking molecules.

The short-Zundel quantum species is crucial for an efficient proton diffusion, as the shortness of its structure enables a fast charge redistribution during the adiabatic PT process. Recent progress in ultrafast broadband two-dimensional (2D) IR spectroscopy⁴⁶,⁴⁷ allowed to probe the vibrational properties of protonated water at vibrational frequencies around the hydrated proton stretching mode, by measuring the lowest-lying excitations in the mid-infrared continuum⁴⁷. These state-of-the-art experiments revealed a strongly inhomogeneous behaviour of the pump-probe spectra, implying large structural distributions in proton asymmetry and $O_1O_2$ distance. Therefore, the traditional "Zundel limit"³ needs to be revisited and extended, in order to cover the broad range of structures detected experimentally⁹,⁴⁸. In particular, the occurrence of qualitatively different short hydrogen-bond configurations, straightforwardly connected with the short-Zundel species described here, has been detected and highlighted in a recent fully solvated ($HF_2$)($H_2O_6$) experiment through femtosecond 2D IR spectroscopy in ref. 49. The present work crucially extends those findings by providing a temperature resolved analysis of the short hydrogen-bond events and by revealing their fundamental relation with the PT dynamics.

While proton transfer and proton transport occur in a variety of environments, from solutions to membrane proteins and fuel-cell membranes, the protonated water hexamer is one of the smallest clusters to incorporate most of the PT experimental features and solvation effects at the leading order. According to ref. 50, one more hydration layer is needed to reach the water bulk limit. From this viewpoint, the hexamer is close to that limit, and some relevant effects, emerging already at this size, can be transferred to larger systems. Our findings thus call for further efforts to explore the temperature behaviour of the proton dynamics and transport both in aqueous systems and in other extended environments, by keeping the same accuracy as the one delivered by our QMC-driven PIMD approach in the protonated water hexamer. That goal could be achieved by training efficient atomistic ML potentials on QMC⁵¹,⁵² or other high-quality datasets²⁴,⁵³, which could allow one to generalise quantitatively the PT behaviour unveiled here to a wider class of aqueous systems.

## Methods

### Zero-temperature electronic structure calculations

Before running finite-temperature calculations, we build a quantum Monte Carlo (QMC) variational wave function. The molecular dynamics (MD) will develop on the potential energy surface (PES) generated by the variational energy of this wave function. All zero- and finite-temperature calculations have been carried out with the TurboRVB code⁵⁴. We choose the variational ansatz such that the chemical accuracy (1 kcal/mol) in the binding energy of the Zundel ion and water dimer is reached. We highlight the fact that benchmarking the binding energy is much stricter than taking energy differences of geometries around the minimum, because the configurations involved in the binding energy are very different.

The variational wave function ($\Psi_q$) we used in our work is written as a Jastrow Antisymmetrised Geminal Power (JAGP) product⁵⁵

$$ \Psi_q(\mathbf{x}_1, \dots, \mathbf{x}_{N_{\text{el}}}) = J_q(\mathbf{r}_1, \dots, \mathbf{r}_{N_{\text{el}}}) \Psi_{\text{AGP},q}(\mathbf{x}_1, \dots, \mathbf{x}_{N_{\text{el}}}). \quad (1) $$

The set {$\mathbf{x}_i = (\mathbf{r}_i, \sigma_i)$}$_{i=1,...,N_{\text{el}}}$ represents spatial and spin coordinates of the $N_{\text{el}}$ electrons, and $\mathbf{q}$ is the vector of nuclear coordinates.

We report here the main ingredients of the JAGP wave function⁵⁶⁻⁵⁹. The Bosonic Jastrow factor is written as a product of one-body, two-body and three/four-body terms $J_q = J_{1,q/2,q/3,q}$. The one-body term reads

$$ J_{1,q} = \exp \left( - \sum_{i}^{N_{\text{el}}} \sum_{j}^{N} (2Z_j)^{3/4} u \left( \left( (2Z_j)^{1/4} \right) |\mathbf{r}_i - \mathbf{q}_j| \right) \right) \quad (2) $$

with $u(|\mathbf{r} - \mathbf{q}|) = \frac{1 - e^{-b|\mathbf{r} - \mathbf{q}|}}{2b}$ and $b$ is a variational parameter, which satisfies the electron-ion Kato cusp conditions. $N$ the number of atoms and $Z_j$ the electric charge of the j-th atom. In the protonated hexamer Hamiltonian, the hydrogen keeps the bare Coulomb potential, while the oxygen atoms are replaced by the Burkatzki-Filippi-Dolg (BFD) pseudopotential⁶⁰, smooth at the electron-ion coalescence points. Thus, $J_{1,q}$ is applied only to the hydrogen atom. The two-body correlations are dealt with by the higher-order Jastrow factors. The two-body Jastrow factor is defined as

$$ J_{2,q} = \exp \left( \sum_{i<j}^{N_{\text{el}}} u \left( |\mathbf{r}_i - \mathbf{r}_j| \right) \right), \quad (3) $$

where $u$ is a function of the same form as in Eq. (2), but with a different variational parameter. The three-four body Jastrow factor is

$$ J_{3,q} = \exp \left( \sum_{i<j}^{N_{\text{el}}} \Phi_{j_q}(\mathbf{r}_i, \mathbf{r}_j) \right), \quad (4) $$

with

$$ \Phi_{j_q}(\mathbf{r}_i, \mathbf{r}_j) = \sum_{a,b} \sum_{\mu,\nu}^{N_{\text{basis}}^j} g_{\mu,\nu}^{a,b} \Psi_{a,\mu}^j(\mathbf{r}_i - \mathbf{q}_a) \Psi_{b,\nu}^j(\mathbf{r}_j - \mathbf{q}_b), \quad (5) $$

where $N_{\text{basis}}^j$ is the number of the basis set functions $\Psi_{a,\mu}^j$. We used optimally contracted geminal embedded orbitals (GEOs)⁶¹ as basis set, expanded over a primitive O(3s,2p,1d) H(2s,1p) Gaussian basis. The convergence study of the water dimer binding energy as a function of the Jastrow GEO expansion is reported in Supplementary Note I.1 of the SI.
The Fermionic part of the wave function is expressed as an anti-symmetrised product of the spin singlet geminals or pairing (AGP) functions $\Phi_q(\mathbf{x}_i, \mathbf{x}_j)$:

$$ \Psi_{AGP,q}(\mathbf{x}_1, \dots, \mathbf{x}_{N_{\text{el}}}) = \hat{A}[\Phi_q(\mathbf{x}_1, \mathbf{x}_2), \dots, \Phi_q(\mathbf{x}_{N_{\text{el}}-1}, \mathbf{x}_{N_{\text{el}}})]. \quad (6) $$

The spatial part $\phi_q(\mathbf{r}_i, \mathbf{r}_j)$ of the spin singlets $\Phi_q(\mathbf{x}_i, \mathbf{x}_j)$ is expanded over $N_{\text{basis}}^{\text{AGP}}$ optimally contracted GEOs, such that

$$ \phi_q(\mathbf{r}_i, \mathbf{r}_j) = \sum_{a,b}^{N} \sum_{\mu,\nu}^{N_{\text{basis}}^{\text{AGP}}} \lambda_{\mu,\nu}^{a,b} \Psi_{a,\mu}^{\text{AGP}}(\mathbf{r}_i - \mathbf{q}_a) \bar{\Psi}_{b,\nu}^{\text{AGP}}(\mathbf{r}_j - \mathbf{q}_b). \quad (7) $$

The AGP GEOs are linear combination of primitive O(5s5p2d) H(4s2p) Gaussian basis functions. We highlight that the Fermionic part is fully optimised at the QMC level for each MD step, and not generated by on-the-fly DFT calculations.

The GEOs are very effective in reducing the total number of variational parameters, by keeping a high level of accuracy. This makes the wave function optimisation$^{5,62,63}$ much more efficient. Dealing with a compact wave function is very important, if one wants to use it as variational ansatz in a MD simulation, because in a MD framework the wave function needs to be optimised at every MD iteration. Indeed, the wave function optimisation is by far the most time-consuming step of our QMC-driven MD approach.

Previous works on the Zundel ion$^{26,59}$ found that the optimal balance between accuracy and computational cost for the determinantal part is reached by the O[8]H[2] contracted GEO basis, in self-explaining notations. As the protonated water hexamer is a very similar system, in this work we used the same O[8]H[2] GEO contraction for the AGP part. Moreover, we further simplified the variational wave function previously developed for the Zundel ion, by contracting also the Jastrow basis set, using the same GEO embedding scheme. We found that the O[6]H[2] GEO basis set for the Jastrow factor is a very good compromise between accuracy and number of parameters, as we checked in the water dimer (see Supplementary Note 1.1 of the SI). The final accuracy is about 1 kcal/mol in the dissociation energy of the water dimer, and supposedly it is much higher around the stable geometries of water clusters. Thus, we used the O[6]H[2] GEO basis set for the Jastrow factor, and the O[8]H[2] GEO basis for the AGP part in all our subsequent MD simulations. For the protonated water hexamer, this results into a total number of 6418 variational parameters, comprising $g_{\mu,\nu}^{a,b}$ in Eq. (5), $\lambda_{\mu,\nu}^{a,b}$ in Eq. (7), the parameters of the homogeneous one-body (Eq. (2)) and two-body (Eq. (3)) Jastrow factors, and the linear coefficients of the Jastrow and determinantal basis sets.

## Finite-temperature calculations

**Path integral Langevin dynamics.** At finite temperature, we carried out path integral Langevin dynamics simulations to include NQEs. To do so, we used the recently developed algorithm published in ref. 26, which combines a path integral approach with very accurate Born-Oppenheimer (BO) forces computed by QMC. It is an efficient approach, alternative to the coupled electron ion Monte Carlo (CEIMC) developed by Ceperley and coworkers$^{64–66}$. The intrinsic noise of the QMC force estimator is treated by the noise correction scheme developed in refs. 25,26,67, which is based on the fulfilment of the fluctuation-dissipation theorem. This implies that the friction matrix $\gamma$ governing the dumped dynamics is related to the random force $\eta$ via the $\alpha$ matrix:

$$ \alpha(\mathbf{q}) = 2k_B T \gamma(\mathbf{q}), \quad (8) $$

$$ \langle \eta_i(t) \eta_j(t') \rangle = \alpha_{i,j}(\mathbf{q}) \delta(t - t'), \quad (9) $$

with $\mathbf{q}$ the vector of nuclear coordinates. The $\mathbf{q}$-dependence comes from the QMC noise correction, implemented through the relations:

$$ \alpha_{i,j}(\mathbf{q}) = \gamma_{BO}/(2k_B T) \delta_{i,j} + \Delta_0 \alpha_{i,j}^{\text{QMC}}(\mathbf{q}) \quad (10) $$

$$ \alpha_{i,j}^{\text{QMC}}(\mathbf{q}) = \langle \langle (\mathbf{f}_i(\mathbf{q}) - \langle \langle \mathbf{f}_i(\mathbf{q}) \rangle \rangle) \rangle \langle (\mathbf{f}_j(\mathbf{q}) - \langle \langle \mathbf{f}_j(\mathbf{q}) \rangle \rangle) \rangle, \quad (11) $$

where $\alpha^{\text{QMC}}$ in Eq. (11) is the covariance matrix of QMC ionic forces, measuring their stochastic fluctuations, and $\gamma_{BO}$ and $\Delta_0$ parameters taking values for an optimal Langevin dynamics. The random force $\eta$ is then used to thermalise the system to a target temperature, according to the Langevin thermostat of Eq. (8).

The quantum particles are described by necklaces extended in the imaginary time interval $[0, \hbar\beta]$, with $\beta = 1/(k_B T)$, following the quantum-to-classical isomorphism. This imaginary time interval is divided into M slices, the so-called “beads”, leading to an effective classical system of NM particles. The path integral Langevin dynamics we developed in ref. 26 is very efficient, because the quantum harmonic forces and the Langevin thermostat—thermalising the quantum degrees of freedom—are evolved together by means of an exact propagator, without Trotter breakup. The evolution of quantum particles in a thermal bath à la Langevin represents a quantum Ornstein-Uhlenbeck dynamics. Therefore, we dubbed our integration scheme as “path integral Ornstein Uhlenbeck dynamics (PIOUD)”. The algorithm is detailed in ref. 26.

In order to evolve the system and to sample the thermal quantum partition function, nuclear forces must be estimated at each iteration by computing the gradients $\mathbf{f} = -\nabla_q V(\mathbf{q})$. The potential energy surface $V(\mathbf{q})$ is evaluated by VMC, namely:

$$ V(\mathbf{q}) = \frac{\langle \Psi_{\mathbf{q}} | H(\mathbf{q}) | \Psi_{\mathbf{q}} \rangle}{\langle \Psi_{\mathbf{q}} | \Psi_{\mathbf{q}} \rangle}, \quad (12) $$

where $|\Psi_{\mathbf{q}}\rangle$ is the QMC wave function, which minimizes the expectation value of $H(\mathbf{q})$ for each bead configuration $\mathbf{q}^{(k)}$, according to the Born-Oppenheimer approximation.

The electronic variational wave function depends on the coordinates $\mathbf{q}^{(k)}$ of the k-th bead in two ways. Directly, through the explicit dependence on the ion positions provided by the localised basis set, and in an indirect way, through the wave function parameters optimised at each $\mathbf{q}^{(k)}$, in compliance with the Born-Oppenheimer approximation. While the former dependence is of leading order, the latter can be neglected in a first approximation. A clever way to do so is to average the optimal parameters across different beads, by gaining a significant fraction of statistics in the QMC energy minimisation. More precisely, the beads are gathered in groups of $N_{\text{groups}}$ members each, to share the same set of wave function parameters. This is called “bead grouping approximation”, and it has been introduced in ref. 26.

Once the number of groups $N_{\text{groups}}$ is set for the beads, the electronic wave function is optimised at each new ionic position generated by the dynamics. This is done by energy minimisation via the most advanced optimisation techniques$^{57,63}$. Between two consecutive steps of ion dynamics, one needs to perform $N_{\text{opt}}$ steps of energy minimisation, in an iterative fashion. $N_{\text{opt}}$ must be large enough to converge the wave function for each new ionic configuration. Thus, this parameter is tuned such that the BO approximation is fulfilled, and the dynamics follows the correct PES along the PIOUD trajectories.

During the dynamics, the GTO exponents in both the Jastrow and the AGP part of the wave function are kept frozen to make the simulation stable. Due to the continuity of the nuclear trajectories, the number of energy minimisation steps is significantly smaller than the one required for a wave function optimisation from scratch.
There is a set of sensitive parameters one needs to tune to have stable and unbiased simulations. They are:

(i) Convergence for quantum effects set by M. In our calculations, we used $M=32$ for $T=250-400$ K, $M=64$ for $T=150-200$ K, and $M=128$ for $T=50-100$ K. The bead grouping approximation is made with $N_{group}=1$. Therefore the whole ring shares the same wave function parameters, except for the nuclear coordinates;

(ii) Time step $\delta t$ for the integration of the equations of motion. We used $\delta t=1$ fs for a controlled time integration error, yielding a difference between the virial and primitive estimators of the quantum kinetic energy below a 25 mHa threshold along all the trajectories;

(iii) A stable target temperature is reached despite the QMC noise by setting the parameters $\gamma_{BO}$ and $\Delta_0$, defining the $\alpha$ matrix in Eq. (10). We used $\gamma_{BO}=0$ and $\Delta_0=0.5\delta t$, optimal values for other similar systems, such as the Zundel ion$^{26}$. Indeed, the simulation is efficient when the damping in the BO sector is minimised. The thermalisation of the system is guaranteed thanks to the optimal damping condition$^{68}$ applied to the internal modes of the ring-polymer, with the damping parameter of the center-of-mass translational modes set to $0.231 \text{ fs}^{-1}$;

(iv) To enforce the fulfilment of the BO approximation, we allowed for $N_{opt}=5$ iterative wave function optimisation steps at each MD iteration, such that the electronic energy minimum is reached within the statistical accuracy for every ionic configuration, and the PES is sampled without stochastic bias.

**Classical Langevin dynamics.** For classical MD calculations, we used an improved variant of the original algorithm developed by Attaccalite and Sorella$^{25}$. This variant has been detailed in refs. 26,67. It includes a better integration scheme, involving a Langevin noise touching both coordinates and momenta, which are therefore correlated.

As in the PIOUD calculations, relevant parameters are the time step $\delta t=1$ fs, the QMC noise correction parameters $\gamma_{BO} = 0.2 \text{ fs}^{-1}$ and $\Delta_0 = \delta t$, and the number of QMC optimisation steps per nuclear iteration $N_{opt} = 5$.

## Error estimates

Datapoints obtained from stochastic sampling, based either on the electron or on the joint electron-nuclear distribution, are provided with the associated stochastic error bars, which are obtained via the reblocking technique$^{69}$. The error bars on the fitting parameters plotted in Supplementary Fig. 14 of the SI are yielded by the error propagation in the Levenberg-Marquardt solution of the nonlinear least squares curve-fitting problem. In all cases, they represent the standard deviation of the data.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

Source data are provided with this paper.

## Code availability

The TurboRVB code$^{54}$ used to carry out our QMC-driven PIMD calculations is available under the GPLv3 license. It can be downloaded from the following open source repository: https://github.com/sissaschool/turborvb. In this work we have used Sandro Sorella's legacy version (v1.0.0), which can be found at this link: https://github.com/sissaschool/turborvb/releases/tag/v1.0.0.

## References

1. Marx, D. Proton transfer 200 years after von Grotthuss: Insights from ab initio simulations. *ChemPhysChem* **7**, 1848–1870 (2006).

2. Cukierman, S. Et tu, Grotthuss! and other unfinished stories. *Biochim. Biophys. Acta (BBA) - Bioenerg.* **1757**, 876–885 (2006).

3. Zundel, G. & Metzger, H. Energiebänder der tunnelnden überschußprotonen in flüssigen säuren. Eine IR-spektroskopische untersuchung der natur der gruppierungen. *Z. Physik. Chem.* **58**, 225–245 (1968).

4. Eigen, M. Proton transfer, acid-base catalysis, and enzymatic hydrolysis. Part I: elementary processes. *Angew. Chem. Int. Ed.* **3**, 1–19 (1964).

5. Vuilleumier, R. & Borgis, D. Transport and spectroscopy of the hydrated proton: a molecular dynamics study. *J. Chem. Phys.* **111**, 4251–4266 (1999).

6. Calio, P. B., Li, C. & Voth, G. A. Resolving the structural debate for the hydrated excess proton in water. *J. Am. Chem. Soc.* **143**, 18672–18683 (2021).

7. Roy, S. et al. Resolving heterogeneous dynamics of excess protons in aqueous solution with rate theory. *J. Phys. Chem. B* **124**, 5665–5675 (2020).

8. Markovitch, O. et al. Special pair dance and partner selection: elementary steps in proton transport in liquid water. *J. Phys. Chem. B* **112**, 9456–9466 (2008).

9. Daly Jr, C. A. et al. Decomposition of the experimental Raman and infrared spectra of acidic water into proton, special pair, and counterion contributions. *J. Phys. Chem. Lett.* **8**, 5246–5252 (2017).

10. Kundu, A. et al. Hydrated excess protons in acetonitrile/water mixtures: solvation species and ultrafast proton motions. *J. Phys. Chem. Lett.* **10**, 2287–2294 (2019).

11. Stoyanov, E. S., Stoyanova, I. V. & Reed, C. A. The structure of the hydrogen ion ($H_{aq}^+$) in water. *J. Am. Chem. Soc.* **132**, 1484–1485 (2010).

12. Bell, R. A., Christoph, G. G., Fronczek, F. R. & Marsh, R. E. The cation $H_{13}O_6^+$: a short, symmetric hydrogen bond. *Science* **190**, 151–152 (1975).

13. Knight, C. & Voth, G. A. The curious case of the hydrated proton. *Acc. Chem. Res.* **45**, 101–109 (2012).

14. Yagi, K. & Thomsen, B. Infrared spectra of protonated water clusters, $h^+$ ($h2o$) 4, in eigen and Zundel forms studied by vibrational quasi-degenerate perturbation theory. *J. Phys. Chem. A* **121**, 2386–2398 (2017).

15. Jiang, J.-C. et al. Infrared spectra of H. *J. Am. Chem. Soc.* **122**, 1398–1410 (2000).

16. Mizuse, K. & Fujii, A. Infrared photodissociation spectroscopy of $H^+(H_2O)_6 \cdot M_m$ ($M$ = Ne, Ar, Kr, Xe, $H_2$, $N_2$, and $CH_4$): messenger-dependent balance between $H_3O^+$ and $H_5O_2^+$ core isomers. *Phys. Chem. Chem. Phys.* **13**, 7129 (2011).

17. Finney, J. M. et al. Isotope effects in the Zundel–Eigen isomerization of $H^+(H_2O)_6$. *J. Phys. Chem. Lett.* **14**, 4666–4672 (2023).

18. Heindel, J. P., Yu, Q., Bowman, J. M. & Xantheas, S. S. Benchmark electronic structure calculations for $H_3O^+(H_2O)_n$, $n=0-5$, clusters and tests of an existing 1, 2, 3-body potential energy surface with a new 4-body correction. *J. Chem. Theory Comput.* **14**, 4553–4566 (2018).

19. Wang, Y.-S. et al. Investigations of protonated and deprotonated water clusters using a low-temperature 22-pole ion trap. *J. Phys. Chem. A* **107**, 4217–4225 (2003).

20. Headrick, J. M. et al. Spectral signatures of hydrated proton vibrations in water clusters. *Science* **308**, 1765–1769 (2005).

21. Heine, N. et al. Isomer-selective detection of hydrogen-bond vibrations in the protonated water hexamer. *J. Am. Chem. Soc.* **135**, 8266–8273 (2013).

22. Wei, D. & Salahub, D. R. Hydrated proton clusters and solvent effects on the proton transfer barrier: a density functional study. *J. Chem. Phys.* **101**, 7633–7642 (1994).
23. Napoli, J. A., Marsalek, O. & Markland, T. E. Decoding the spectroscopic features and time scales of aqueous proton defects. *J. Chem. Phys.* **148**, 222833 (2018).

24. Schran, C., Brieuc, F. & Marx, D. Transferability of machine learning potentials: protonated water neural network potential applied to the protonated water hexamer. *J. Chem. Phys.* **154**, 051101 (2021).

25. Attaccalite, C. & Sorella, S. Stable liquid hydrogen at high pressure by a novel Ab Initio molecular-dynamics calculation. *Phys. Rev. Lett.* **100**, 114501 (2008).

26. Mouhat, F., Sorella, S., Vuilleumier, R., Saitta, A. M. & Casula, M. Fully quantum description of the Zundel ion: combining variational quantum Monte Carlo with path integral Langevin dynamics. *J. Chem. Theory Comput.* **13**, 2400–2417 (2017).

27. Yu, Q. et al. q-AQUA: a many-body CCSD(T) water potential, including four-body interactions, demonstrates the quantum nature of water from clusters to the liquid phase. *J. Phys. Chem. Lett.* **13**, 5068–5074 (2022).

28. Huang, X., Braams, B. J. & Bowman, J. M. Ab initio potential energy and dipole moment surfaces for $H_5O_2^+$. *J. Chem. Phys.* **122**, 044308 (2005).

29. Koval, S., Kohanoff, J., Migoni, R. & Tosatti, E. Ferroelectricity and isotope effects in hydrogen-bonded KDP crystals. *Phys. Rev. Lett.* **89**, 187602 (2002).

30. Robertson, J. M. & Ubbelohde, A. R. J. P. Structure and thermal properties associated with some hydrogen bonds in crystals I. The isotope effect. *Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci.* **170**, 222–240 (1939).

31. Jones, D. J., Rozière, J., Penfold, J. & Tomkinson, J. Incoherent inelastic neutron scattering studies of proton conducting materials trivalent metal acid sulphate hydrates: Part I. The vibrational spectrum of $H_5O_2^+$. *J. Mol. Struct.* **195**, 283–291 (1989).

32. Stoyanov, E. S. & Reed, C. A. IR spectrum of the $H_5O_2^+$ cation in the context of proton disolvates $L-H^+-L$. *J. Phys. Chem. A* **110**, 12992–13002 (2006).

33. Suzuki, K., Tachikawa, M. & Shiga, M. Temperature dependence on the structure of Zundel cation and its isotopomers. *J. Chem. Phys.* **138**, 184307 (2013).

34. Kreuer, K.-D., Paddison, S. J., Spohr, E. & Schuster, M. Transport in proton conductors for fuel-cell applications: simulations, elementary reactions, and phenomenology. *Chem. Rev.* **104**, 4637–4678 (2004).

35. Asthagiri, D., Pratt, L. & Kress, J. Ab initio molecular dynamics and quasichemical study of $H^+(aq)$. *Proc. Natl Acad. Sci. USA* **102**, 6704–6708 (2005).

36. Schmitt, U. W. & Voth, G. A. The computer simulation of proton transport in water. *J. Chem. Phys.* **111**, 9361–9381 (1999).

37. Huchmala, R. M. & McCoy, A. B. Exploring the origins of spectral signatures of strong hydrogen bonding in protonated water clusters. *J. Phys. Chem. A* **126**, 1360–1368 (2022).

38. Richardson, J. O. & Althorpe, S. C. Ring-polymer molecular dynamics rate-theory in the deep-tunneling regime: Connection with semiclassical instanton theory. *J. Chem. Phys.* **131**, 214106 (2009).

39. Väinhtein, A. I., Zakharov, V. I., Novikov, V. A. & Shifman, M. A. ABC of instantons. *Sov. Phys. Usp.* **25**, 195–215 (1982).

40. Jiang, Z. et al. Scaling analysis and instantons for thermally assisted tunneling and quantum Monte Carlo simulations. *Phys. Rev. A* **95**, 382 (2017).

41. Mazzola, G., Smelyanskiy, V. N. & Troyer, M. Quantum Monte Carlo tunneling from quantum chemistry to quantum annealing. *Phys. Rev. B* **96**, 483 (2017).

42. Hele, T. J. On the relation between thermostatted ring-polymer molecular dynamics and exact quantum dynamics. *Mol. Phys.* **114**, 1461–1471 (2016).

43. Arcis, H., Plumridge, J. & Tremaine, P. R. Limiting conductivities of strong acids and bases in $D_2O$ and $H_2O$: Deuterium isotope effects on proton hopping over a wide temperature range. *J. Phys. Chem. B* **126**, 8791–8803 (2022).

44. Linstrom, P. J. & Mallard, W. G. National Institute of Standards and Technology (NIST) chemistry webbook, NIST standard reference database **69**, published by NIST, (Washington, D.C., 1997).

45. Chandler, D. & Leung, K. Excess electrons in liquids: geometrical perspectives. *Annu. Rev. Phys. Chem.* **45**, 557–591 (1994).

46. Dahms, F., Fingerhut, B. P., Nibbering, E. T. J., Pines, E. & Elsaesser, T. Large-amplitude transfer motion of hydrated excess protons mapped by ultrafast 2D IR spectroscopy. *Science* **357**, 491–495 (2017).

47. Fournier, J. A., Carpenter, W. B., Lewis, N. H. C. & Tokmakoff, A. Broadband 2D IR spectroscopy reveals dominant asymmetric $H_5O_2^+$ proton hydration structures in acid solutions. *Nat. Chem.* **10**, 932–937 (2018).

48. Carpenter, W. B. et al. Decoding the 2D IR spectrum of the aqueous proton with high-level VSCF/VCI calculations. *J. Chem. Phys.* **153**, 124506 (2020).

49. Dereka, B. et al. Crossover from hydrogen to chemical bonding. *Science* **371**, 160–164 (2021).

50. Lapid, H., Agmon, N., Petersen, M. K. & Voth, G. A. A bond-order analysis of the mechanism for hydrated proton mobility in liquid water. *J. Chem. Phys.* **122**, 014506 (2005).

51. Tirelli, A., Tenti, G., Nakano, K. & Sorella, S. High-pressure hydrogen by machine learning and quantum Monte Carlo. *Phys. Rev. B* **106**, 041105 (2022).

52. Niu, H. et al. Stable solid molecular hydrogen above 900 K from a machine-learned potential trained with diffusion quantum Monte Carlo. *Phys. Rev. Lett.* **130**, 076102 (2023).

53. Zhang, L., Wang, H., Car, R. & Weinan, E. Phase diagram of a deep potential water model. *Phys. Rev. Lett.* **126**, 236001 (2021).

54. Nakano, K. et al. Turborvb: a many-body toolkit for ab initio electronic simulations by quantum Monte Carlo. *J. Chem. Phys.* **152**, 204121 (2020).

55. Casula, M., Attaccalite, C. & Sorella, S. Correlated geminal wave function for molecules: an efficient resonating valence bond approach. *J. Chem. Phys.* **121**, 7110 (2004).

56. Casula, M., Filippi, C. & Sorella, S. Diffusion Monte Carlo method with lattice regularization. *Phys. Rev. Lett.* **95**, 1431 (2005).

57. Sorella, S., Casula, M. & Rocca, D. Weak binding between two aromatic rings: feeling the van der Waals attraction by quantum Monte Carlo methods. *J. Chem. Phys.* **127**, 014105 (2007).

58. Marchi, M., Azadi, S., Casula, M. & Sorella, S. Resonating valence bond wave function with molecular orbitals: application to first-row molecules. *J. Chem. Phys.* **131**, 154116 (2009).

59. Dagrada, M., Casula, M., Saitta, A. M., Sorella, S. & Mauri, F. Quantum Monte Carlo study of the protonated water dimer. *J. Chem. Theory Comput.* **10**, 1980–1993 (2014).

60. Burkatzki, M., Filippi, C. & Dolg, M. Energy-consistent pseudopotentials for quantum Monte Carlo calculations. *J. Chem. Phys.* **126**, 234105 (2007).

61. Sorella, S., Devaux, N., Dagrada, M., Mazzola, G. & Casula, M. Geminal embedding scheme for optimal atomic basis set construction in correlated calculations. *J. Chem. Phys.* **143**, 244112 (2015).

62. Sorella, S. Generalized Lanczos algorithm for variational quantum Monte Carlo. *Phys. Rev. B* **64**, 024512 (2001).

63. Umrigar, C. J., Toulouse, J., Filippi, C., Sorella, S. & Hennig, R. G. Alleviation of the fermion-sign problem by optimization of many-body wave functions. *Phys. Rev. Lett.* **98**, 110201 (2007).

64. Morales, M. A., Pierleoni, C., Schwegler, E. & Ceperley, D. M. Evidence for a first-order liquid-liquid transition in high-pressure
hydrogen from ab initio simulations. *Proc. Natl Acad. Sci. USA* **107**,
12799–12803 (2010).

65. Liberatore, E., Morales, M. A., Ceperley, D. M. & Pierleoni, C. Free energy methods in coupled electron ion Monte Carlo. *Mol. Phys.* **109**, 3029–3036 (2011).

66. Rillo, G., Morales, M. A., Ceperley, D. M. & Pierleoni, C. Coupled electron-ion Monte Carlo simulation of hydrogen molecular crystals. *J. Chem. Phys.* **148**, 102314 (2018).

67. Mazzola, G., Yunoki, S. & Sorella, S. Unexpectedly high pressure for molecular dissociation in liquid hydrogen by electronic simulation. *Nat. Commun.* **5**, 1–6 (2014).

68. Ceriotti, M., Parrinello, M., Markland, T. E. & Manolopoulos, D. E. Efficient stochastic thermostatting of path integral molecular dynamics. *J. Chem. Phys.* **133**, 124104 (2010).

69. Flyvbjerg, H. & Petersen, H. G. Error estimates on averages of correlated data. *J. Chem. Phys.* **91**, 461–466 (1989).

## Acknowledgements

F.M., R.V., A.M.S. and M.C. acknowledge Dominik Marx for useful discussions. M.P., R.V., A.M.S. and M.C. thank the CNRS for the 80PRIME doctoral grant allocation. F.M. and M.C. acknowledge computational resources provided by the PRACE project number 2016163936. M.C. thanks GENCI for providing computational resources at TGCC and IDRIS under the grant number 0906493, and the Grands Challenges DARI for allowing calculations on the Joliot-Curie Rome HPC cluster under the project number gch0420. T.M. and M.C. are grateful to the European Centre of Excellence in Exascale Computing TREX-Targeting Real Chemical Accuracy at the Exascale, which partially supported this work. This project has received funding from the European Unions Horizon 2020 Research and Innovation programme under Grant Agreement No. 952165.

## Author contributions

R.V., A.M.S. and M.C. conceived the project. F.M., M.P., and M.C. carried out the PIMD simulations. F.M. and M.C. prepared the variational wave function used in the QMC calculations. All the authors analysed and interpreted the data. F.M., R.V., A.M.S. and M.C. wrote the manuscript, with the help from M.P. and T.M.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-42366-4.

**Correspondence** and requests for materials should be addressed to Michele Casula.

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023