Dynamics of protonated oxalate from machine-learned simulations and experiment: infrared signatures, proton transfer dynamics and tunneling splittings

Abstract

The infrared spectroscopy and proton transfer dynamics together with the associated tunneling splittings for H/D-transfer in oxalate are investigated using a machine learning-based potential energy surface (PES) of CCSD(T) quality, calibrated against the results of new spectroscopic measurements. Second order vibrational perturbation calculations (VPT2) very successfully describe both the framework and H-transfer modes compared with the experiments. In particular, a newly observed low-intensity signature at 1666 cm⁻¹ was correctly predicted from the VPT2 calculations. An unstructured band centered at 2940 cm⁻¹ superimposed on a broad background extending from 2600 to 3200 cm⁻¹ is assigned to the H-transfer motion. The broad background involves a multitude of combination bands but a major role is played by the COH-bend. For the deuterated species, VPT2 and molecular dynamics simulations provide equally convincing assignments, in particular for the framework modes. Finally, based on the new PES the tunneling splitting for H-transfer is predicted as Δ_H = 35.0 cm⁻¹ from ring polymer instanton calculations using higher-order corrections. This provides an experimentally accessible benchmark to validate the computations, in particular the quality of the machine-learned PES.

---

1 Introduction

Vibrational spectroscopy is a nondestructive technique to characterize the structure and dynamics of chemical and biological materials in the gas- and condensed-phase.^1–3 One of the main applications in chemistry is the identification and classification of compounds which is essential in analytical, materials, and synthetic chemistry. Of prime importance in the practical application of vibrational spectroscopy to concrete problems is the ability to assign spectral responses to the underlying motions in the molecule or material of interest. This can sometimes be achieved through analogy with related compounds for which assignments of spectra to motions have already been made.⁴ Alternatively, isotopic substitutions can be used to assist with this assignment which, however, can be tedious from a preparatory perspective.

One of the most direct ways to arrive at assignment and identification of spectral patterns is through computation. Vibrational spectroscopy is an area where “experiment” and “theory/computation” meet in a most natural and direct fashion. However, fruitful interplay between these two approaches requires high-level computer methods which typically need to go beyond the harmonic approximation. It is generally believed that calculations at the coupled cluster with singles, doubles and perturbative triples (CCSD(T)) are needed for reliable computational vibrational spectroscopic work.⁵ Also, the nuclear dynamics needs to be described at a sufficiently advanced level. Standard normal mode analysis is often sufficient to obtain an overview of the spectral features. If quantitative agreement with experiment is sought, more advanced techniques such as vibrational perturbation theory (VPT2)⁶ or spectra from explicit molecular dynamics (MD) simulations are needed, though.^2,7

A chemically particularly interesting motif involves hydrogen bonds and shared hydrogen atoms/protons between an acceptor and a donor. H-bonds and hydrogen/proton transfer (HT/PT) play important roles in governing molecular structure, stability, and dynamics and their energetics and dynamics are of fundamental importance in biology and chemistry.^8–10 HT/PT is primarily influenced by the height of the barrier which is, however, difficult to determine reliably through direct experimentation. Possibilities include high resolution spectroscopy where the splitting of spectral lines can provide information about the barrier height,¹¹ or nuclear magnetic resonance (NMR) experiments^12,13 whereas kinetic isotope effects or shifts of bands in the vibrational spectrum alone can not be used directly to determine the energetics for PT. Proton transfer in systems containing X–H*⋯Y motifs – where X and Y are the donor and acceptor atoms, respectively, and H* is the transferring hydrogen – can lead to characteristically broadened features in vibrational spectra.¹⁴ This broadening reflects strong coupling between the X–H stretch and other framework modes of the environment and structural heterogeneity.¹⁵ The broadening can even persist down to low temperatures and cooling the species does not necessarily lead to sharper bands.^16,17

In experiments, the broadening of spectral lines in X−H*⋯Y species has been extensively observed and reported in both the liquid and the gas phase.¹ For example, in liquid water, the maximum in the OH-stretching spectrum shifts to the red as the temperature decreases from 47 °C to −6 °C, with concomitant increase of the intensity by 16% even without hydrogen or proton transfer to take place.¹⁸ In the high-density liquid, the line shapes in the OH stretching spectra of supercritical methanol was found to be sensitive to the hydrogen bonding network which depends on temperature.¹

Spectral shifts and their broadening are even more pronounced in systems where a transferring proton H* is shared by a donor and acceptor moiety and potentially provides information about the proton transfer energetics. An empirical relationship between the position of the infrared (IR) absorption and the height of the proton transfer barrier has been found in combined computational/experimental investigations of acetylacetone¹⁹ and formic acid dimer.²⁰ Earlier studies of the infrared spectroscopy of protonated ammonia dimer have also established that the IR-signatures are broad and correlated with the barrier height for proton transfer.²¹

Here, a high-level theoretical study is combined with improved spectroscopic measurements to elucidate the vibrational dynamics at play in the protonated oxalate (HO₂CCO₂⁻) anion, hereafter referred to as OxH and OxD for the hydrogenated and deuterated forms of the protonated dianion, see Fig. 1. OxH is a compelling system to study the dynamics and ensuing spectroscopic signatures of a shared proton coupled to framework modes of the surrounding molecular scaffold.^15,22 Previous computational work concerned the kinetic isotope effects from electronic structure calculations,²³ the solution dynamics of hydrated p-Oxa,²⁴ and the IR spectroscopy in the gas phase using a reactive empirical energy function.²² The experimental IR spectrum features a rather diffuse/broad signal between 2600 and 2900 cm⁻¹. One potential source is strong mechanical coupling between the high-frequency O–H stretch and (potentially) the low-frequency COH bends and/or skeletal deformations.

Fig. 1 Atom labels for OxH and OxD molecules used for internal coordinates in power spectra.

Studying the dynamics and spectroscopy of a transferring proton from MD simulations requires energy functions capable of describing bond-formation and bond-breaking. One such method is molecular mechanics with proton transfer (MMPT) which combines an accurate representation of the energetics for the transferring proton and a lower-level empirical force field for the remaining degrees of freedom which allows long simulation times.^25–27 Using “PES morphing”²⁸ provides flexibility to adapt characteristics such as the barrier height and position of the minima. More recently, alternative ways to represent high-dimensional potential energy surfaces have been explored. One of them is based on neural networks and was applied to bulk silicon or the H + HBr reaction.^29,30 The evaluation of the neural network, once trained, is orders of magnitude faster than the underlying reference quantum chemical calculations the model was trained on. To reach the highest levels of quantum chemical theory, such as coupled cluster-level, which is required for spectroscopic studies, Transfer learning (TL) can be used. Such approaches have been shown to be data efficient for obtaining validated high quality PESs.^31–40 TL is based on the notion that the topography of the energy landscape of a molecular system is very similar even at different levels of theory. This allows fine-tuning the learnable parameters of a neural network trained on a lower level of theory using small amounts of data calculated at a high level. Here, the MP2 level of theory provides a computationally efficient route to a broad coverage of the energy landscape, while the smaller high-accuracy dataset allows more precise modeling.

With such “vetted” PESs, it should be even possible to predict experimentally accessible observables. Similar to the shared proton situation in malonaldehyde or tropolone, OxH should exhibit tunneling splitting for the (quantum) motion between the two symmetrical wells. Quantum mechanical tunnelling is a key concept in chemistry and physics and is well established,⁴¹ and includes, inter alia, electron transfer in solution,^42–45 enzymatic reactions,^46–48 and the structure and dynamics of water.^49–52 In contrast to the often challenging measurement of chemical rates with high precision, tunnelling splittings can be accurately measured through spectroscopic methods, e.g., using microwave⁵³ or infrared spectroscopy.⁵⁴ Tunneling splittings provide valuable benchmarks for validating state-of-the-art theoretical methodologies, which require both an accurate PES and an accurate treatment of the resulting nuclear Hamiltonian.⁵⁵ The calculation of accurate tunneling splittings for multidimensional systems using fully quantum mechanical methods is a challenging task. The ring-polymer instanton (RPI) approach is a semiclassical approximation method for tunnelling calculations which is capable of scaling well with system size.^50,56,57 However, instanton theory makes a local harmonic approximation around the optimal tunnelling pathway, resulting in deviations from quantum–mechanical benchmarks that can reach 20%.⁵⁸ This error can be mitigated by perturbatively correcting the results using information about the third and fourth derivatives along the path.⁵⁹ It has been demonstrated that correction of a low level potential energy surface to a higher level using TL, in combination with ring-polymer instanton, can be used to quantitatively characterise tunnelling splittings in such systems as malonaldehyde³⁴ and larger systems consisting of 15-atom tropolone and 12-atom propiolic acid–formic acid dimer (PFD).⁴⁰

The present work is structured as follows. First, the computational and experimental methods are described, followed by a validation of the machine learned PESs. Then, results on the IR spectroscopy of OxH and OxD, the H-transfer dynamics and HT tunneling splittings are presented. Finally, conclusions are drawn.

2 Methods

2.1 Machine learning-based PESs

Machine learning is used to learn a full-dimensional, reactive potential energy surface (PES) for OxH based on ab initio reference data of MP2/aug-cc-pVTZ quality. The data, available from previous work,⁶⁰ was used to train a message-passing neural network (PhysNet).⁶¹ Structures for OxH were sampled by running NVT MD simulations at multiple temperatures ([100, 300, 500, 1000, 1500] K) using the semiempirical GFN2-xTB method⁶² (2500 structures each except for 1500 K for which only 1000 were used). The region around the proton transfer transition state was sampled with a constraining harmonic potential (1500 structures) and at T = 500 K. Additionally, normal mode sampling⁶³ at progressively higher temperatures ([100, 300, 500, 1000, 1500, 2000] K) was carried out for both the optimized and transition state structure of OxH (800 for each T). This yields a total of 22 100 structures, for which reference energies (E^ref), forces (F^ref_i,α), and dipole moments (p^ref_α) (α = x, y, z) were determined at the MP2/aug-cc-pVTZ level of theory using MOLPRO⁶⁴ and served as the reference data set. PhysNet was then trained by minimizing the loss functionusing the Adam optimizer.^65,66 The hyperparameters^61,67 w_i, i ∈ {E, F, Q, p} differentially weigh the contributions to the loss function and were w_E = 1 [1/energy], w_F ∼ 52.92 [length/energy], w_Q ∼ 14.39 [1/charge] and w_p ∼ 27.21 [1/charge/length], respectively, and the term is a “nonhierarchical penalty” that regularizes the loss function.⁶¹ For training, a 80/10/10% split of the data as training/validation/test sets was used. Two independent PhysNet models were trained, which were subsequently employed for active learning.⁶⁸

Adaptive sampling⁶⁹ and diffusion Monte Carlo simulations (DMC)⁷⁰ were run to ensure the robustness of the PES and to locate so-called “holes”. Adaptive sampling was done at 1000 K in the NVT ensemble and structures (to be included in training the next models) were detected if the predictions of the two independent models differed by more than 0.5 kcal mol⁻¹. All DMC simulations employed a GPU implementation of the DMC algorithm³³ for PhysNet and were run with 300 walkers, a step size of 5 a.u. and for a total of 60 000 steps. Unphysical structures were saved if the predicted energy was lower than the energy of the minimum energy structure. During the first active learning cycle, adaptive sampling provided 80 structures and 20 defective structures were located from DMC simulations. Ab initio MP2 data for the 100 structures was determined and added to the MP2 data set before retraining PhysNet. A second active learning cycle was performed, however, no additional deficient structures were found. The final and robust MP2-based PES constitutes the “base model”.

To improve the performance of the base model for molecular spectroscopy, TL was used to refine and elevate the PES from the MP2/aug-cc-pVTZ to the gold-standard CCSD(T)/aug-cc-pVTZ level of theory.^71–73 The TL data set contains a total of 2688 structures selected semi-randomly from the original data set. It comprises 1067 structures that were sampled using MD simulations at 1000 K, 960 and 661 structures obtained from normal mode sampling around the minimum and transition state, respectively, for which energies, forces and dipole moments were determined at the CCSD(T)/aug-cc-pVTZ level of theory using MOLPRO.⁶⁴ These structures and corresponding quantum chemical information were then used to fine-tune the learnable parameters of the MP2 PES with the same hyperparameters as for the low-level training except for the learning rate and the weighting hyperparameters in the loss function. The learning rate was reduced to 10⁻⁴ and w_E = w_F = w_Q = w_p = 1.

2.2 Molecular dynamics simulations

The molecular dynamics simulations were run using the atomic simulation environment (ASE)⁷⁴ with energies and forces obtained from PhysNet. The structure of OxH was first optimized. Next, random momenta were drawn from a Maxwell–Boltzmann distribution corresponding to T = [300, 600, 900] K, which were assigned to the atoms. The MD simulations were carried out in the NVE ensemble using the Velocity Verlet algorithm⁷⁵ and a time step of Δt = 0.2 fs to conserve energy as all bonds including those involving hydrogen atoms were flexible. The molecular system was equilibrated for 2 ps, followed by production simulations of 200 ps. For each temperature and system (OxH and OxD), 500 independent trajectories were run using different initial momenta. This accumulated to a total of 100 ns simulation time per temperature T.

2.3 Infrared spectroscopy

IR spectra I(ω) were calculated from the Fourier transform of the dipole–dipole auto-correlation function^76,77 according toHere, μ_i(t) is the molecular dipole moment along direction i at time t and Q(ω) is a quantum correction factor⁷⁸

2.4 Ring polymer instanton calculations

The RPI method is a semiclassical approximation for computing tunneling splittings in molecular systems.^79,80 RPI locates the optimal tunneling pathway, known as the instanton, and is defined as an imaginary-time τ → ∞ path connecting two degenerate wells which minimizes the total action, S. The path is constructed by optimizing a path discretized into N ring-polymer beads and taking the limit N → ∞ (typically N ∼ 1000 is sufficient for convergence). The potential U_N of a ring polymer is given bywith x = (x₁,…,x_N) being the mass-scaled coordinates of the beads, and β_N = β/N. The first term in eqn (3) corresponds to the sum over all single bead potentials and the second term represents the harmonic springs with frequency 1/(β_Nħ) that connect adjacent beads. As the potential U_N and the action S are related by division with (β_Nħ), the action S contains information along the instanton path (IP). A detailed description of the method is provided, e.g., in ref. 79 and 80.

Within standard RPI (sRPI) theory, fluctuations around the instanton path are included up to second order and the information is combined into a contribution Φ. More specifically, fluctuations around the path are determined from the Hessian matrix at each of the beads.^34,57 With this, the leading-order tunneling splitting in a symmetric double-well system is

For malonaldehyde and formic acid dimer sRPI tunneling splittings on high quality PESs were within ∼20% of the experimental values.^34,40,81 This error arises due to approximations underlying sRPI theory⁵⁸ and remaining deficiencies in the PES because Δ_RPI depends exponentially on the action S, see eqn (4). One limitation of the sRPI method for determining tunneling splittings is that fluctuations around the instanton are only included up to second order.⁸⁰ To a good approximation, it is expected that this captures the dominant tunneling contribution, except for cases in which anharmonic effects perpendicular to the instanton are significant or where the barrier is low. For this reason, a perturbatively corrected RPI (pcRPI) theory was recently developed,⁵⁹ that accounts for anharmonicity by including information from the third and fourth order derivatives of the potential along the instanton. The tunneling splitting obtained from pcRPI theory is denoted as Δ_PC in the following which, in practice, is obtained from scaling Δ_RPI with a correction factor c_PC according to Δ_PC = c_PC·Δ_RPI.

2.5 Experimental methods

The vibrational predissociation spectrum of oxalate with two H₂-tags was recorded on the Yale triple focusing photofragmentation time-of-flight (TOF) mass spectrometer that has been described previously.^82–84 1 mM solutions of oxalic acid in pure acetonitrile were prepared and the ionic species were generated via electrospray ionization (ESI) of the oxalic acid solutions in negative ion mode. ESI results in deprotonation of one of the carboxylic acid head groups, generating the oxalate anion. Ions are then transferred via radiofrequency guides to a temperature-controlled (5–40 K), 3D (Paul) ion trap. A few ms before the ions arrive in the trap, He buffer gas doped with 20% H₂ was pulsed into the trap where collisions with He thermalize the ions to the trap temperature, allowing H₂ tags to condense onto the ions.^85–89 After ∼90 ms, the tagged ions are ejected from the trap into the TOF region of the instrument. At a transient focus of the TOF, the ions are excited with an IR pulse generated from a LaserVision OPO/OPA (5 ns, 2–20 mJ per pulse) in the 800–4200 cm⁻¹ range. Photofragments were generated when the ions absorb a single resonant IR photon, evaporating the H₂ tags. The photofragment yield was measured as a function of photon energy to generate an infrared photodissociation (IRPD) spectrum in a linear action mode. The reported spectra are normalized to the laser fluence to account for the wavelength dependence of the output power.

3 Results

3.1 Validation of the potential energy surfaces

3.1.1 The MP2 PES. The out-of-sample errors of the MP2 base model are reported in Table 1 and their correlation and signed errors are shown in Fig. S1. The energy barrier for hydrogen transfer on the PhysNet PES is 2.355 kcal mol⁻¹ and differs from the ab initio MP2/aug-cc-pVTZ barrier by only ΔE_a = 8 × 10⁻⁴ kcal mol⁻¹. Harmonic frequencies, obtained from diagonalizing the Hessian matrix characterize the shape of the PES around stationary points. PhysNet achieves a mean absolute error for the harmonic frequencies, MAE(ω), for the harmonic frequencies at the global minimum and the TS for H-transfer of 0.09 and 0.11 cm⁻¹ with respect to the ab initio harmonic frequencies. All harmonic frequencies are given in Table S1.

Table 1 Mean absolute and root mean squared errors (MAE and RMSE) for energies E and forces F on a hold-out test set containing 2220 and 270 random OxH structures for MP2 and CCSD(T) level of theory, respectively. In addition, the deviation 1 − R² of the correlation coefficient from 1, the barrier height E_a of the PhysNet model and the deviation ΔE_a from the corresponding electronic structure calculation is reported. MAE(ω) is the difference between the normal mode frequencies calculated from the ML-PES and at the respective level of theory. For normal modes see also Tables S1 and S2. Energies are reported in kcal mol⁻¹, forces in kcal mol⁻¹ Å⁻¹, and harmonic frequencies in cm⁻¹

	PhysNet (MP2)	PhysNet (CCSD(T))
MAE(E)	0.009	0.005
RMSE(E)	0.047	0.012
MAE(F)	0.065	0.023
RMSE(F)	0.459	0.065
1 − R^2	5.9 × 10^−06	7.8 × 10^−07
Ea	2.355	3.351
ΔEa	8 × 10^−04	7 × 10^−04
MAE(ω)	0.09	0.34

---

3.1.2 The transfer-learned CCSD(T) PES. Transfer-learning to the CCSD(T) level of theory used only 10% of the data underlying the MP2-base model. The out-of-sample errors for the TL-PES are given in Table 1. Their correlation and signed errors are shown in Fig. S2. The energy barrier for hydrogen transfer on the PhysNet PES is 3.351 kcal mol⁻¹ and differs from the ab initio CCSD(T)/aug-cc-pVTZ barrier by only ΔE_a = 7 × 10⁻⁴ kcal mol⁻¹. The MAE(ω) for the minimum and H-transfer TS between harmonic frequencies from the TL-PES and normal modes from explicit CCSD(T)/aVTZ calculations are 0.34 and 0.57 cm⁻¹ with respect to the ab initio harmonic frequencies. All harmonic frequencies are given in Table S2.

3.2 Infrared spectroscopy

The infrared spectrum of OxH-(D₂)₂ recorded for this study, hereafter denoted IR2025, is displayed in Fig. 2. This spectrum is characterized by a number of sharper bands for the region between 1200 and 2000 cm⁻¹ and a considerably broader feature between 2500 and 3200 cm⁻¹. This feature can be viewed as a band centered around 2950 cm⁻¹, superimposed on a diffuse background absorption spanning over 700 cm⁻¹. Because there is only one hydrogen atom in OxH, the band in the high-frequency region must be associated with fundamentals and combination modes involving primarily the OH-stretching and COH-bending modes.

Fig. 2 Experimental predissociation vibrational spectrum of doubly-H₂ tagged oxalate. For clarity, the spectrum in the range 2400–3400 cm⁻¹ is expanded by a factor of 5 in the inset. The sharper feature at 4090 cm⁻¹ is due to the nominally IR forbidden H₂ fundamental, which appears −72 cm⁻¹ below the vibrational quantum in isolated (bare) H₂.⁹⁰

One specific reason to revisit the previously published IR spectrum¹⁵ was to reassess the broad background between 2500 and 3200 cm⁻¹ and establish the relative intensities of the bands across the entire spectrum. The complete spectrum also includes the nominally forbidden fundamental of the H₂ stretch, (observed at 4090 cm⁻¹, labeled ν_H2) which is red-shifted by −72 cm⁻¹ when attached to the molecular anion. The IR2015 (grey dashed) and IR2025 (black solid) spectra are reported in Fig. 3A. The plateau-like shape between 2550 to 2820 cm⁻¹ in the IR2015 spectrum features a more gradual increase in IR2025 but overall, the two spectra are remarkably similar, see Fig. S3. The spectrum from the present work, however, reveals a better defined shoulder near 1660 cm⁻¹ on the red side of the strong band at 1696 cm⁻¹ than was evident in the IR2015 spectrum (Fig. S3). All other features agree between the two spectra. One reason for the differences between IR2015 and IR2025 concerned normalization with respect to laser power which was much improved in IR2025.

Fig. 3 Infrared spectra of OxH. Panel A: Measured (H₂-tagged) gas phase spectra IR2015¹⁵ (dashed grey) and IR2025 from the present work (solid black). Panel B: Spectrum from second-order vibrational perturbation theory, VPT2,⁶ using Lorentzian Broadening with HWHM = 4 cm⁻¹ and bands labelled from a to g for identification. Asterisks indicate fundamental vibrations, see also Table S3. Panels C to F: IR spectra at different simulation temperatures as indicated from the Fourier transform of the dipole–dipole autocorrelation function. Grey shaded areas highlight the position and width of the measured framework modes for direct comparison between computations and experiment.

Second order vibrational perturbation theory⁶ provides a useful way to go beyond the harmonic approximation underlying normal mode analysis, also for H-bonded systems.⁹¹ The VPT2 spectrum from using the TL-PES is reported in Fig. 3B. Focusing first on the framework modes below 2000 cm⁻¹ it is found that all features in the measured spectrum line up well with the VPT2 spectrum (bands a to g; fundamentals a, b₁, c, e₂, f indicated by asterisks). This includes the multiplett-structures around 1360 cm⁻¹ (bands b₁, b₂, c, and d) and centered at 1750 cm⁻¹ (bands e₁, e₂, f, and g). Most notably, a low-intensity signature e₁ at 1666 cm⁻¹ which was not clearly present in the IR2015 spectrum, is manifest in the VPT2 calculations and the IR2025 spectrum, see Fig. S3, and assigned to the ν₅ + ν₁₀ combination band, see Table S3. Furthermore, the low-frequency ν₉ mode (band a) at 928 cm⁻¹ in the IR2025 spectrum is at 937 cm⁻¹ in the VPT2 calculations and a small feature at 1091 cm⁻¹ from the measurements and at 1090 cm⁻¹ from the VPT2 calculations (see entry with asterisk in Table S3, not labelled). The VPT2 spectrum using the MP2 ML-PES is only qualitatively correct, see Fig. S4, and is not suitable for quantitative analysis.

The high-frequency part of the measured spectrum shows a wide absorption (h₁) between 2500 and 3200 cm⁻¹ onto which a less broad signature h₂ centered at ∼2900 cm⁻¹ is superimposed. VPT2 calculations feature a considerable number of discrete peaks covering the frequency range consistent with the experiments. As the experiment does not probe a single conformation but rather a distribution of structures, it is of interest to obtain VPT2 frequencies for perturbed minimum energy structures. For this, the optimized oxalate structure was distorted by perturbing all bond lengths with a standard deviation of 0.01 Å followed by a VPT2 calculation for each structure. The Boltzmann-averaged spectrum, shown in Fig. S5, demonstrates that the computed signal is more representative of the measured frequency distribution than the VPT2 calculation for a single minimum energy structure.

Because the entire molecular structure was perturbed, some of the framework modes also shift away from the frequencies corresponding to the minimum energy structure, see Fig. S6. However, if only the transferring hydrogen atom is displaced away from its equilibrium position by 0.002 Å along the OH-bond vector, the VPT2 calculations still find the pattern in the H-transfer region to broaden and shift by several 10 cm⁻¹ to the red or blue, respectively, whereas the framework modes respond primarily with respect to their relative intensities but the peak positions remain at the measured frequencies, see Fig. S7. It should be noted that VPT2 calculations use the harmonic approximation as its zeroth-order approximation to which anharmonic corrections are added. Hence, for strongly anharmonic modes or modes involving Fermi resonances a VPT2 treatment may be unreliable.^32,92

H-transfer during the dynamics will also be reflected in the infrared spectroscopy. To include finite-temperature effects, MD simulations using the TL-PES were carried out for temperatures between 50 K and 900 K, see Fig. 3C–F. The OH-stretch is found at 3200 cm⁻¹ and broadens considerably as temperature increases. The blue shift away from the maximum of H-stretch mode at 2940 cm⁻¹ from experiments is well-understood and occurs because MD simulations only sample a narrow range around the bottom of the well without accessing the anharmonic parts of the PES, in particular for the OH-stretch mode.⁹³ As temperature increases, broader signatures emerge which overlap with the peak positions from the VPT2 calculations and the measured spectrum.

For partial assignment of the measured peaks, power spectra for different internal coordinates were determined, see Fig. 4. Due to coupling between different degrees of freedom, straightforward and unambiguous assignments are not particularly meaningful. Two infrared signatures at 1700 and 1770 cm⁻¹ are clearly assigned in the IR spectrum and involve the C_A–O_A, C_A–O_B, C_B–O_C, and C_B–O_D distances. The peak to the red, at 1666 cm⁻¹ and part of the high-energy triplet, does not appear in the IR and power spectra, but is clearly visible in the VPT2 calculations which also report combination modes and hot bands, see Table S3. The triplet centered around 1360 cm⁻¹ is also reproduced by the computed IR spectra and involves the C_AO_AH angle and all four CO-stretches.

Fig. 4 Power spectra for OxH molecule at 300 K for different internal coordinates (low-frequency range).

In order to further investigate the nature of the broad peak in the 2600–3200 cm⁻¹ region, MD trajectories run at 600 K were analyzed. For this, separate IR and power spectra were generated for trajectories that did or did not feature H-transfer, see Fig. 5 and Fig. S8, S9. The internal coordinates considered were the C_A–O_A–H angle and the C_A–O_B, C_B–O_C, C_B–O_D, C_A–O_A, and O_A–H separations. As already mentioned, the blue shift of the OH-stretching vibration is well-understood.^93,94 However, the broad [2600, 3200] cm⁻¹ background from the measurements is rather convincingly associated with signatures in the C_A–O_A–H angle especially for trajectories featuring H-transfer. Interestingly, shifting the OH-power spectrum to the red by −260 cm⁻¹ not only aligns the experimentally observed peak at 2940 cm⁻¹ but also the faint measured feature around 2450 cm⁻¹. For formic acid monomer,⁹⁵ the shift required to align the computed with the measured infrared signature arising from the OH-stretch was −290 cm⁻¹, which is consistent with the observation for the OH-stretch in oxalate.

Fig. 5 Power spectra for the COH-bend (green) and OH-stretch (orange) from simulations at 600 K that with (solid line, 362 trajectories) and without (dashed line, 68 trajectories) H-transfer during the dynamics. The indigo trace is the IR spectrum (300 K; no H-transfer). The blue dashed line is the shifted OH-stretch power spectrum to best overlap with the feature in the measured spectrum. For spectra at higher T that feature H-transfer, see Fig. 3.

To further corroborate this assignment, simulations were carried out whereby the mass of the hydrogen atom ranged from 1m_H to 2m_H (i.e. m_D) in increments of 0.2 mass units, see Fig. S10. As expected, the spectral response shifts to the red as the mass of the hydrogen atom increases. It is interesting to note that with increasing temperature in Fig. 3 the region above 2600 cm⁻¹ gains intensity (note also the different y-scaling factors). Analysis of the MD trajectories shows that with increasing temperature the probability for H-transfer between the oxygen atoms O_A and O_D increases. Hence, the power spectra for different internal coordinates were separately analyzed for trajectories that did (dashed lines) and did not (solid) feature hydrogen transfer, see Fig. 5. In the region of interest (above 2600 cm⁻¹) the power spectrum of the COH bending motion (green trace) features clear signals which are totally absent from simulations that do not feature H-transfer. The OH-stretch motion (orange) is also clearly visible and coupled to the COH-bend.

3.3 Spectroscopy of OxD

For the deuterated species, OxD, the IR2015 spectrum¹⁵ was used to further validate the TL-PES and to provide spectroscopic assignments. First, the measurements were compared with the VPT2 calculations, see Fig. 6. Again, for the framework modes (bands a to f below 2000 cm⁻¹) excellent agreement is found, see Table S4 for numerical values and assignments. The spectral patterns related to the deuterium motion, which are located above 2000 cm⁻¹ are also very well captured, in particular finer and low-intensity peaks around 3000 cm⁻¹ and above. However, there is a prominent pattern between 1860 and 1980 cm⁻¹ (blue shaded area) which is largely empty for the VPT2 calculations. To further explore the spectroscopy in this frequency range a perturbed minimum energy structure was generated with a bond length standard deviation of 0.01 Å. The result of the VPT2 calculations were compared with the VPT2 calculation for the minimum energy structure (red) and experimental data (brown), see Fig. S11. It is found that perturbation away from the minimum energy geometry leads to a shift in the OD-stretch frequency, consistent with the findings for OxH.

Fig. 6 Infrared spectra for OxD. The top trace reports the experimentally determined (H₂-tagged) gas phase spectrum¹⁵ and the bottom trace (light coral) is the computed spectrum from second-order vibrational perturbation theory⁶ using the CCSD(T) ML-PES. Note that data within the dashed box have been scaled by a factor of 10 for the two traces.

Next, the infrared spectra were obtained from MD simulations at 600 K together with power spectra for important internal degrees of freedom. Again, trajectories without and with D-transfer were analyzed separately, see Fig. 7A and B. The OD-stretch vibration likely corresponds to the peak observed at 2110 cm⁻¹ in the experimental spectrum, whereas in the spectrum computed from MD simulations, this feature appears at 2320 cm⁻¹. The shift of approximately −200 cm⁻¹ for OxD is consistent with the isotope-corrected shift observed for OxH . The maximum peak position at 2110 cm⁻¹ in the experimental spectrum of OxD, similar to that in OxH, is associated with motions along the O_A D-stretch and C_AO_A D-bending motions, see olive and blue traces Fig. 7B.

Fig. 7 Power Spectra and IR spectra for OxD (high-frequency region). Panel A: Without deuterium transfer, Panel B: with deuterium transfer.

3.4 H-Transfer and tunneling splittings

First, the classical “over-the-barrier” H-transfer process was considered from equilibrium MD simulations. Running classical dynamics simulations at different temperatures, the number of barrier crossings was determined from which the H-transfer rate can be obtained. For this 500 MD simulations, 200 ps each (aggregate of 100 ns) were carried out at temperatures T ∈ [50, 900] using ASE. These simulations yield proton transfer (PT) rates 0, 4.8, 131.4 and 339.0 ns⁻¹ at T = 50, 300, 600, and 900 K, respectively. Using an Arrhenius expression one obtains E_a = 3.84 kcal mol⁻¹, and A = 3057.3 ns⁻¹. This compares and is consistent with a barrier height of 3.35 kcal mol⁻¹ from the CCSD(T)/aug-cc-pVTZ calculations.

For a first estimation of the tunneling splittings, sRPI theory was employed together with the single-precision PhysNet representation of the PES. Such a treatment is considered appropriate for cases in which anharmonic effects perpendicular to the instanton path are less pronounced and the barrier for H-transfer is not particularly low. All leading-order splitting calculations in the present work were carried out with up to N = 4096 beads at an effective temperature of T_eff = ħ/k_Bτ = 25 K. Using the MP2-level PES, featuring a barrier height E_B = 2.355 kcal mol⁻¹, Δ_H = 116.5 cm⁻¹ was obtained. As for the higher-level TL-PES E_B increases to 3.351 kcal mol⁻¹, the tunneling splitting decreases to Δ_H = 36.0 cm⁻¹. The numerical convergence of the leading-order tunneling splittings as a function of N is reported in Table S5.

One limitation of sRPI for determining tunnelling splittings is that fluctuations around the instanton are treated within a harmonic approximation (i.e. the Hessians are used). As a remedy, a perturbatively corrected approach was developed which is based on third and fourth order derivatives of the PES along the instanton path.⁵⁹ For numerical stability, pcRPI requires highly accurate higher-order derivatives. Within the context of machine learning-based PESs it is common to use single-precision arithmetics for training and evaluation, which was found to be insufficient for pcRPI.⁹⁶ Hence, a double-precision TL-PES was obtained from a retrained MP2 PhysNet model also using double precision arithmetics throughout. The higher order instanton calculations were carried out at T_eff = 25 K and with N = 1024 beads as pcRPI is considerably more memory intensive than sRPI, and the perturbative correction was found to converge significantly more rapidly than the leading-order term.⁴⁰ The correction factor for H-transfer was c_PC = 0.98 within pcRPI which leads to Δ^pcRPI_H = 35 cm⁻¹. This is 1 cm⁻¹ smaller than the tunneling splitting of 36.0 cm⁻¹ determined on the single precision TL-PES. Assuming a conservative error estimate of 20% for sRPI theory applied to the computed tunneling splitting using the TL-PES brackets the prediction for H-transfer splitting to around Δ_H ∈ [28, 42] cm⁻¹. More accurate predictions will likely require further refinement of the TL-PES, with a specific focus on tunneling splittings^34,40 which are very sensitive to the shape of the PES, in particular around the global minimum and the barrier to H-transfer.

Earlier, recent work using transfer-learned PESs to the CCSD(T) level of theory for malonaldehyde and tropolone showed excellent agreement with measured tunneling splittings.^34,40,96 For malonaldehyde the measured⁹⁷ splitting Δ^expt._H = 21.6 cm⁻¹ compares with a computed sRPI (pcRPI) Δ^TL_H = 25.3 (22.1) cm⁻¹ on a PES transfer learned to the CCSD(T) level of theory and featuring a barrier height of E_B = 3.89 kcal mol⁻¹. For tropolone,⁹⁸ on the other hand, Δ^expt._H = 0.97 cm⁻¹ is considerably smaller because E_B = 6.62 kcal mol⁻¹ at the CCSD(T) level of theory and pcRPI Δ^theo._H = 0.94 cm⁻¹. Hence increasing the barrier by ca. 2 kcal mol⁻¹ reduces Δ_H by approximately 1 order of magnitude. This is consistent with an increase to Δ_H = 35.0 cm⁻¹ for oxalate because the barrier to HT decreases by 0.5 kcal mol⁻¹ (to E_B = 3.35 kcal mol⁻¹) compared with malonaldehyde.

4 Summary and conclusion

The present work provides a molecularly resolved picture for the dynamics, spectroscopy and H/D-transfer dynamics in the oxalate anion. VPT2 calculations and finite-temperature MD simulations using a CCSD(T)-quality PES provide very realistic spectral patterns for the framework modes (below 2000 cm⁻¹), including the prediction of the ν₅ + ν₁₀ combination band (feature e₁) at 1652 cm⁻¹ from VPT2 calculations which is clearly defined (observed at 1666 cm⁻¹) in the new measurements (IR2025, see Fig. S3).

A particular focus was on exploring the spectral signatures between 2400 and 3200 cm⁻¹ for OxH and 1900 to 2500 cm⁻¹ for OxD, respectively. The VPT2 calculations find a high density of states in these frequency ranges and slightly perturbing the structure away from the minimum energy configuration further increases the density of vibrational signatures. These calculations together with the finite-temperature MD simulations identify a center-frequency at ∼2940 and ∼2110 cm⁻¹ as the main H-/D-transfer mode. However, unequivocal assignment of one feature to a single motional pattern is not meaningful because the fundamentals, combination and hot bands are heavily mixed.

Consistent with earlier work,¹⁵ it is found here that the COH angle strongly couples to the OH stretching frequencies. This is particularly apparent for trajectories that feature H-transfer. The VPT2 calculations combined with MD simulations at the highest level of theory reported in the literature provided assignments of the vibrational modes and correct previous mode descriptions. The mismatch in assignment with the one made for IR2015 is mostly related to using a comparatively low level of quantum chemical theory (B3LYP/6-311++G(d,p)). Using a reactive force field, the effective barrier height for H-transfer was also inferred from matching the measured IR2015 spectrum with that computed from the dipole-autocorrelation function from finite-temperature MD simulations. An effective barrier height of 4.2 kcal mol⁻¹ was found to best reproduce the IR2015 in the region of the H-transfer mode. This compares with an electronic barrier height of 3.35 kcal mol⁻¹ from the present work and 3.4 kcal mol⁻¹ from previous CCSD(T)/aug-cc-pVTZ calculations.⁹⁹ However, the electronic and effective barrier heights describe different physical quantities and thus can not be compared directly. The electronic E_B represents the pure potential energy difference between reactant and the transition state whereas the effective barrier height scales the PES along a single (reaction) coordinate (here the HT-motion), depends on temperature, and accounts for zero-point effects.

Tunneling splittings for OxH have been determined with RPI theory. A rather conservative estimate for the tunneling splitting is Δ^sRPI_H ∈ [28, 42] cm⁻¹ which, however, is likely to require improvements of the TL-PES for quantitative predictions of potential future experiments. Nevertheless, the predicted range offers a reference for future spectroscopic efforts and experimental confirmation would not only validate the RPI framework but also aid in the iterative improvement of the PES guided by experiments as has been recently done for small ionic complexes through morphing.^28,100,101

In summary, an updated measured infrared spectrum for the oxalate anion was discussed in the context of state-of-the art VPT2 calculations and MD simulations using a transfer-learned CCSD(T)-level PES. Assignments were made based on inspecting the VPT2 vibrations and the power spectra. In the region of the framework modes considerable coupling between the different degrees of freedom was found. The agreement between theory and experiment and the fact that a specific framework mode was correctly located by the VPT2 calculations highlights the predictive power of the theoretical approaches, whereby vibrational motion can be reliably related to observed spectral features. Examples such as OxH and OxD are particularly valuable as they allow computer-based methods to approach “prediction mode”, finding the “right answer for the right reason”.

Conflicts of interest

The authors declare no conflict of interest.

Data availability

The PhysNet codes and input files for the MP2 and CCSD(T) PESs are available from https://github.com/MMunibas/oxalate/.

Supplementary information: Fig. S1 and S2 present the performance of PhysNet models trained at different levels of theory. Fig. S3–S11 provide supplementary analyses and spectral data that support and extend the results discussed in the main manuscript. Tables S1–S4 report a detailed comparison between the vibrational frequencies predicted by the models and those obtained from reference quantum chemical methods and experimental measurements, including peak assignments. Table S5 summarizes the calculated tunneling splittings together with their convergence. See DOI: https://doi.org/10.1039/d5cp03085d.

Acknowledgements

The authors gratefully acknowledge financial support from the Swiss National Science Foundation through grants 200020_219779 (MM), 200021_215088 (MM), the University of Basel (MM). MAJ gratefully acknowledges the Department of Energy through the condensed phase and interfacial molecular science (CPIMS) program under grants DE-SC0021012. Research reported in this publication was supported by the National Institute Of General Medical Sciences of the National Institutes of Health under Award Number T32GM149438 (to ELB). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

References

1. H. J. Bakker and J. L. Skinner, Vibrational spectroscopy as a probe of structure and dynamics in liquid water, Chem. Rev., 2010, 110, 1498–1517.
2. M.-P. Gaigeot, Theoretical spectroscopy of floppy peptides at room temperature. A DFTMD perspective: gas and aqueous phase, Phys. Chem. Chem. Phys., 2010, 12, 3336–3359.
3. E. Kraka, W. Zou and Y. Tao, Decoding chemical information from vibrational spectroscopy data: local vibrational mode theory, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2020, 10, e1480.
4. S. Lammers and M. Meuwly, Investigating the Relationship between Infrared Spectra of Shared Protons in Different Chemical Environments: A Comparison of Protonated Diglyme and Protonated Water Dimer, J. Phys. Chem. A, 2007, 111, 1638–1647.
5. C. Qu and J. M. Bowman, Quantum approaches to vibrational dynamics and spectroscopy: is ease of interpretation sacrificed as rigor increases?, Phys. Chem. Chem. Phys., 2019, 21, 3397–3413.
6. V. Barone, J. Bloino, C. A. Guido and F. Lipparini, A fully automated implementation of VPT2 infrared intensities, Chem. Phys. Lett., 2010, 496, 157–161.
7. D. Koner, S. M. Salehi, P. Mondal and M. Meuwly, Perspective: Non-conventional Force Fields for Applications in Spectroscopy and Chemical Reaction Dynamics, J. Chem. Phys., 2020, 153, 010901.
8. W. Cleland and M. Low-barrier Kreevoy, Hydrogen-bonds and Enzymatic Catalysis, Science, 1994, 264, 1887–1890.
9. A. Warshel, A. Papazyan and P. Kollman, On Low-barrier Hydrogen-bonds and Enzyme Catalysis, Science, 1995, 269, 102–104.
10. S. Shan, S. Loh and D. Herschlag, The Energetics of Hydrogen Bonds in Model Systems: Implications for Enzymatic Catalysis, Science, 1996, 272, 97–101.
11. D. Borst, J. Roscioli, D. Pratt, G. Florio, T. Zwier, A. Muller and S. Leutwyler, Hydrogen Bonding and Tunneling in the 2-pyridone·2-hydroxypyridine Dimer. Effect of Electronic Excitation, Chem. Phys., 2002, 283, 341–354.
12. P. Frey, S. Whitt and J. Tobin, A Low-barrier Hydrogen-bond in the Catalytic Triad of Serine Proteases, Science, 1994, 264, 1927–1930.
13. J. Braun, M. Schlabach, B. Wehrle, M. Kocher, E. Vogel and H. Limbach, NMR-study of the Tautomerism of Porphyrin Including the Kinetic HH/HD/DD Isotope Effects in the Liquid and the Solid State, J. Am. Chem. Soc., 1994, 116, 6593–6604.
14. L. Bondesson, K. V. Mikkelsen, Y. Luo, P. Garberg and H. Ågren, Hydrogen bonding effects on infrared and Raman spectra of drug molecules, Spectrochim. Acta, Part A, 2007, 66, 213–224.
15. C. T. Wolke, A. F. DeBlase, C. M. Leavitt, A. B. McCoy and M. A. Johnson, Diffuse vibrational signature of a single proton embedded in the oxalate scaffold, HO₂CCO₂⁻, J. Phys. Chem. A, 2015, 119, 13018–13024.
16. C. M. Leavitt, A. F. DeBlase, C. J. Johnson, M. van Stipdonk, A. B. McCoy and M. A. Johnson, Hiding in Plain Sight: Unmasking the Diffuse Spectral Signatures of the Protonated N-terminus in Isolated Dipeptides Cooled in a Cryogenic Ion Trap, J. Phys. Chem. Lett., 2013, 4, 3450–3457.
17. C. J. Johnson, L. C. Dzugan, A. B. Wolk, C. M. Leavitt, J. A. Fournier, A. B. McCoy and M. A. Johnson, Microhydration of Contact Ion Pairs in M²⁺OH⁻(H₂O)_n=1–5 (M = Mg, Ca) Clusters: Spectral Manifestations of a Mobile Proton Defect in the First Hydration Shell, J. Phys. Chem. A, 2014, 118, 7590–7597.
18. J.-B. Brubach, A. Mermet, A. Filabozzi, A. Gerschel and P. Roy, Signatures of the hydrogen bonding in the infrared bands of water, J. Chem. Phys., 2005, 122, 184509.
19. D. L. Howard, H. G. Kjaergaard, J. Huang and M. Meuwly, Infrared and near-infrared spectroscopy of acetylacetone and hexafluoroacetylacetone, J. Phys. Chem. A, 2015, 119, 7980–7990.
20. K. Mackeprang, Z.-H. Xu, Z. Maroun, M. Meuwly and H. G. Kjaergaard, Spectroscopy and dynamics of double proton transfer in formic acid dimer, Phys. Chem. Chem. Phys., 2016, 18, 24654–24662.
21. M. Meuwly and M. Karplus, Simulation of Proton Transfer Along Ammonia Wires: An ab Initio and Semiempirical Density Functional Comparison of Potentials and Classical Molecular Dynamics, J. Chem. Phys., 2002, 116, 2572–2585.
22. Z.-H. Xu and M. Meuwly, Vibrational spectroscopy and proton transfer dynamics in protonated oxalate, J. Phys. Chem. A, 2017, 121, 5389–5398.
23. T. N. Truong and J. A. McCammon, Direct dynamics study of intramolecular proton transfer in hydrogenoxalate anion, J. Am. Chem. Soc., 1991, 113, 7504–7508.
24. O. Kroutil, B. Minofar and M. Kabeláč, Structure and dynamics of solvated hydrogenoxalate and oxalate anions: a theoretical study, J. Mol. Model., 2016, 22, 1–10.
25. S. Lammers, S. Lutz and M. Meuwly, Reactive Force Fields for Proton Transfer Dynamics, J. Comput. Chem., 2008, 29, 1048–1063.
26. S. Lutz, I. Tubert-Brohman, Y. Yang and M. Meuwly, Water-assisted Proton Transfer in Ferredoxin I, J. Biol. Chem., 2011, 286, 23679–23687.
27. Y. Yang and M. Meuwly, A generalized reactive force field for nonlinear hydrogen bonds: hydrogen dynamics and transfer in malonaldehyde, J. Chem. Phys., 2010, 133, 064503.
28. M. Meuwly and J. M. Hutson, Morphing ab initio potentials: a systematic study of Ne-HF, J. Chem. Phys., 1999, 110, 8338–8347.
29. J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Phys. Rev. Lett., 2007, 98, 146401.
30. A. Pukrittayakamee, M. Malshe, M. Hagan, L. Raff, R. Narulkar, S. Bukkapatnum and R. Komanduri, Simultaneous fitting of a potential-energy surface and its corresponding force fields using feedforward neural networks, J. Chem. Phys., 2009, 130, 134101.
31. P. O. Dral, A. Owens, A. Dral and G. Csányi, Hierarchical machine learning of potential energy surfaces, J. Chem. Phys., 2020, 152, 204110.
32. S. Käser, E. D. Boittier, M. Upadhyay and M. Meuwly, Transfer Learning to CCSD(T): Accurate Anharmonic Frequencies from Machine Learning Models, J. Chem. Theory Comput., 2021, 17, 3687–3699.
33. S. Käser and M. Meuwly, Transfer learned potential energy surfaces: accurate anharmonic vibrational dynamics and dissociation energies for the formic acid monomer and dimer, Phys. Chem. Chem. Phys., 2022, 24, 5269–5281.
34. S. Käser, J. O. Richardson and M. Meuwly, Transfer Learning for Affordable and High-Quality Tunneling Splittings from Instanton Calculations, J. Chem. Theory Comput., 2022, 18, 6840–6850.
35. V. Zaverkin, D. Holzmüller, L. Bonfirraro and J. Kästner, Transfer learning for chemically accurate interatomic neural network potentials, Phys. Chem. Chem. Phys., 2023, 25, 5383–5396.
36. M. S. Chen, J. Lee, H.-Z. Ye, T. C. Berkelbach, D. R. Reichman and T. E. Markland, Data-efficient machine learning potentials from transfer learning of periodic correlated electronic structure methods: liquid water at AFQMC, CCSD, and CCSD(T) accuracy, J. Chem. Theory Comput., 2023, 19, 4510–4519.
37. A. Nandi, G. Laude, S. S. Khire, N. D. Gurav, C. Qu, R. Conte, Q. Yu, S. Li, P. L. Houston and S. R. Gadre, et al., Ring-polymer instanton tunneling splittings of tropolone and isotopomers using a Δ-machine learned CCSD(T) potential: theory and experiment shake hands, J. Am. Chem. Soc., 2023, 145, 9655–9664.
38. S. Käser and M. Meuwly, Transfer-learned potential energy surfaces: toward microsecond-scale molecular dynamics simulations in the gas phase at CCSD(T) quality, J. Chem. Phys., 2023, 158, 214301.
39. A. Nandi, P. Pandey, P. L. Houston, C. Qu, Q. Yu, R. Conte, A. Tkatchenko and J. M. Bowman, Δ-Machine Learning to Elevate DFT-Based Potentials and a Force Field to the CCSD(T) Level Illustrated for Ethanol, J. Chem. Theory Comput., 2024, 20, 8807–8819.
40. S. Käser, J. O. Richardson and M. Meuwly, Transfer Learning for Predictive Molecular Simulations: Data-Efficient Potential Energy Surfaces at CCSD(T) Accuracy, J. Chem. Theory Comput., 2025, 21, 6633–6643.
41. V. I. Gol'danskii, L. I. Trakhtenberg and V. N. Fleurov, Tunneling phenomena in chemical physics, Routledge, 2021.
42. D. Chandler, Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific, 1998, pp. 25–49.
43. A. R. Menzeleev, N. Ananth and T. F. Miller, Direct simulation of electron transfer using ring polymer molecular dynamics: comparison with semiclassical instanton theory and exact quantum methods, J. Chem. Phys., 2011, 135, 074106.
44. W. Fang, R. A. Zarotiadis and J. O. Richardson, Revisiting nuclear tunnelling in the aqueous ferrous-ferric electron transfer, Phys. Chem. Chem. Phys., 2020, 22, 10687–10698.
45. J. E. Lawrence and D. E. Manolopoulos, Confirming the role of nuclear tunneling in aqueous ferrous-ferric electron transfer, J. Chem. Phys., 2020, 153, 154114.
46. S. Hammes-Schiffer, Hydrogen tunneling and protein motion in enzyme reactions, Acc. Chem. Res., 2006, 39, 93–100.
47. S. Hu, S. C. Sharma, A. D. Scouras, A. V. Soudackov, C. A. M. Carr, S. Hammes-Schiffer, T. Alber and J. P. Klinman, Extremely elevated room-temperature kinetic isotope effects quantify the critical role of barrier width in enzymatic C–H activation, J. Am. Chem. Soc., 2014, 136, 8157–8160.
48. J. B. Rommel, Y. Liu, H.-J. Werner and J. Kästner, Role of tunneling in the enzyme glutamate mutase, J. Phys. Chem. B, 2012, 116, 13682–13689.
49. F. N. Keutsch and R. J. Saykally, Water clusters: untangling the mysteries of the liquid, one molecule at a time, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 10533–10540.
50. J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso, G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate and S. C. Althorpe, Concerted hydrogen-bond breaking by quantum tunneling in the water hexamer prism, Science, 2016, 351, 1310–1313.
51. M. T. Cvitaš and J. O. Richardson, Molecular Spectroscopy and Quantum Dynamics, Elsevier, 2021, pp. 301–326.
52. M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales and T. E. Markland, Nuclear quantum effects in water and aqueous systems: experiment, theory, and current challenges, Chem. Rev., 2016, 116, 7529–7550.
53. W. F. Rowe Jr, R. W. Duerst and E. B. Wilson, The intramolecular hydrogen bond in malonaldehyde, J. Am. Chem. Soc., 1976, 98, 4021–4023.
54. M. Ortlieb and M. Havenith, Proton transfer in (HCOOH)₂: an IR high-resolution spectroscopic study of the antisymmetric C–O stretch, J. Phys. Chem. A, 2007, 111, 7355–7363.
55. J. Cruzan, L. Braly, K. Liu, M. Brown, J. Loeser and R. Saykally, Quantifying hydrogen bond cooperativity in water: VRT spectroscopy of the water tetramer, Science, 1996, 271, 59–62.
56. J. O. Richardson and S. C. Althorpe, Ring-polymer instanton method for calculating tunneling splittings, J. Chem. Phys., 2011, 134, 054109.
57. J. O. Richardson, Ring-polymer instanton theory, Intern. Rev. Phys. Chem., 2018, 37, 171–216.
58. J. O. Richardson, Perspective: ring-polymer instanton theory, J. Chem. Phys., 2018, 148, 200901.
59. J. E. Lawrence, J. Dušek and J. O. Richardson, Perturbatively corrected ring-polymer instanton theory for accurate tunneling splittings, J. Chem. Phys., 2023, 159, 014111.
60. S. Käser, D. Koner and M. Meuwly, The Bigger the Better? Accurate Molecular Potential Energy Surfaces from Minimalist Neural Networks, arXiv, 2024, preprint, arXiv:2411.18121 DOI:10.48550/arXiv.2411.18121.
61. O. T. Unke and M. Meuwly, PhysNet: a neural network for predicting energies, forces, dipole moments, and partial charges, J. Chem. Theory Comput., 2019, 15, 3678–3693.
62. C. Bannwarth, S. Ehlert and S. Grimme, GFN2-xTB—An accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions, J. Chem. Theory Comput., 2019, 15, 1652–1671.
63. J. S. Smith, O. Isayev and A. E. Roitberg, ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost, Chem. Sci., 2017, 8, 3192–3203.
64. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona and R. Lindh, et al., MOLPRO, version 2019, a package of ab initio programs, 2019.
65. D. Kingma and J. Ba, Adam: a method for stochastic optimization. arXiv, 2014, preprint, arXiv:1412.6980 DOI:10.48550/arXiv.1412.6980.
66. S. J. Reddi, S. Kale and S. Kumar, On the Convergence of Adam and Beyond, arXiv, 2019, preprint, arXiv:1904.09237; https://arxiv.org/abs/1904.09237.
67. K. Song, S. Käser, K. Töpfer, L. I. Vazquez-Salazar and M. Meuwly, PhysNet Meets CHARMM: A Framework for Routine Machine Learning/Molecular Mechanics Simulations, arXiv, 2023, preprint, arXiv:2304.12973 DOI:10.48550/arXiv.2304.12973.
68. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Machine learning force fields, Chem. Rev., 2021, 121, 10142–10186.
69. G. Csányi, T. Albaret, M. Payne and A. De Vita, “Learn on the fly”: a hybrid classical and quantum-mechanical molecular dynamics simulation, Phys. Rev. Lett., 2004, 93, 175503.
70. I. Kosztin, B. Faber and K. Schulten, Introduction to the diffusion Monte Carlo method, Am. J. Phys., 1996, 64, 633–644.
71. M. E. Taylor and P. Stone, Transfer learning for reinforcement learning domains: a survey, J. Mach. Learn. Res., 2009, 10, 1633–1685.
72. S. J. Pan and Q. Yang, A survey on transfer learning, IEEE Trans. Knowl. Data Eng., 2009, 22, 1345–1359.
73. J. S. Smith, B. T. Nebgen, R. Zubatyuk, N. Lubbers, C. Devereux, K. Barros, S. Tretiak, O. Isayev and A. E. Roitberg, Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning, Nat. Commun., 2019, 10, 1–8.
74. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dulak, J. Friis, M. N. Groves, B. Hammer and C. Hargus, et al., The atomic simulation environment—a Python library for working with atoms, J. Phys.: Condens. Matter, 2017, 29, 273002.
75. L. Verlet, Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev., 1967, 159, 98.
76. R. Gordon, Adv. Magn. Opt. Res., Elsevier, 1968, vol. 3, pp. 1–42.
77. B. J. Berne and R. Pecora, Dynamic light scattering: with applications to chemistry, biology, and physics, Courier Corporation, 2000.
78. C. P. Lawrence and J. L. Skinner, Quantum corrections in vibrational and electronic condensed-phase spectroscopy: line shapes and echoes, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6720–6725.
79. J. O. Richardson and S. C. Althorpe, Ring-polymer instanton method for calculating tunneling splittings, J. Chem. Phys., 2011, 134, 054109.
80. J. O. Richardson, Ring-polymer instanton theory, Int. Rev. Phys. Chem., 2018, 37, 171–216.
81. J. O. Richardson, Full- and reduced-dimensionality instanton calculations of the tunnelling splitting in the formic acid dimer, Phys. Chem. Chem. Phys., 2017, 19, 966–970.
82. Techniques for the Study of Ion–Molecule Reactions, ed. J. M. Farrar, J. Saunders and H. William, Techniques of Chemistry, John Wiley & Sons, New York, 1988, vol. 20, p. 652.
83. F. S. Menges, E. H. Perez, S. C. Edington, C. H. Duong, N. Yang and M. A. Johnson, Integration of high-resolution mass spectrometry with cryogenic ion vibrational spectroscopy, J. Am. Soc. Mass Spectrom., 2019, 30, 1551–1557.
84. N. Yang, C. H. Duong, P. J. Kelleher and M. A. Johnson, Unmasking Rare, Large-Amplitude Motions in D₂-Tagged I-˙(H₂O) 2 Isotopomers with Two-Color, Infrared-Infrared Vibrational Predissociation Spectroscopy, J. Phys. Chem. Lett., 2018, 9, 3744–3750.
85. A. B. Wolk, C. M. Leavitt, E. Garand and M. A. Johnson, Cryogenic ion chemistry and spectroscopy, Acc. Chem. Res., 2014, 47, 202–210.
86. J. Mikosch, H. Kreckel, R. Wester, R. Plašil, J. Glosk, D. Gerlich, D. Schwalm and A. Wolf, Action spectroscopy and temperature diagnostics of H₃⁺ by chemical probing, J. Chem. Phys., 2004, 121, 11030–11037.
87. O. V. Boyarkin and V. Kopysov, Cryogenically cooled octupole ion trap for spectroscopy of biomolecular ions, Rev. Sci. Instr., 2014, 85, 033105.
88. R. Wester, Radiofrequency multipole traps: tools for spectroscopy and dynamics of cold molecular ions, J. Phys. B: At., Mol. Opt. Phys., 2009, 42, 154001.
89. O. Asvany and S. Schlemmer, Numerical simulations of kinetic ion temperature in a cryogenic linear multipole trap, Int. J. Mass Spectrom., 2009, 279, 147–155.
90. O. C. Moss, T. Schleif, J. P. Messinger, A. G. Rullán Buxó, K. Greis, E. H. Perez and M. A. Johnson, Hydrogen Tag Shifts as Vibrational Reporters for Positional Isomers of Formylphenides: Surprising Mobility of the Carbanion Center Upon Collisional Decarboxylation of the Parent Benzoates, J. Phys. Chem. Lett., 2024, 15, 1969–1974.
91. J. C. Howard, J. L. Gray, A. J. Hardwick, L. T. Nguyen and G. S. Tschumper, Getting down to the Fundamentals of Hydrogen Bonding: Anharmonic Vibrational Frequencies of (HF)₂ and (H₂O)₂ from Ab Initio Electronic Structure Computations, J. Chem. Theory Comput., 2014, 10, 5426–5435.
92. A. A. Kananenka and J. L. Skinner, Fermi resonance in OH-stretch vibrational spectroscopy of liquid water and the water hexamer, J. Chem. Phys., 2018, 148, 244107.
93. A. Nejad and M. A. Suhm, Concerted Pair Motion Due to Double Hydrogen Bonding: The Formic Acid Dimer Case, J. Indian Inst. Sci., 2020, 100, 5–19.
94. K. Töpfer, S. Käser and M. Meuwly, Double proton transfer in hydrated formic acid dimer: interplay of spatial symmetry and solvent-generated force on reactivity, Phys. Chem. Chem. Phys., 2022, 24, 13869–13882.
95. S. Käser and M. Meuwly, Transfer Learned Potential Energy Surfaces: Accurate Anharmonic Vibrational Dynamics and Dissociation Energies for the Formic Acid Monomer and Dimer, Phys. Chem. Chem. Phys., 2022, 24, 5269–5281.
96. S. Kaeser and M. Meuwly, Numerical Accuracy Matters: Applications of Machine Learned Potential Energy Surfaces, J. Phys. Chem. Lett., 2024, 15, 3419–3424.
97. T. Baba, T. Tanaka, I. Morino, K. M. Yamada and K. Tanaka, Detection of the tunneling-rotation transitions of malonaldehyde in the submillimeter-wave region, J. Chem. Phys., 1999, 110, 4131–4133.
98. K. Tanaka, H. Honjo, T. Tanaka, H. Kohguchi, Y. Ohshima and Y. Endo, Determination of the proton tunneling splitting of tropolone in the ground state by microwave spectroscopy, J. Chem. Phys., 1999, 110, 1969–1978.
99. D. Boutwell, D. Pierre-Jacques, O. Cochran, J. Dyke, D. Salazar, C. Tyler and M. Kaledin, Intramolecular proton transfer in the hydrogen oxalate anion and the cooperativity effects of the low-frequency vibrations: a driven molecular dynamics study, J. Phys. Chem. A, 2022, 126, 583–592.
100. K. P. Horn, L. I. Vazquez-Salazar, C. P. Koch and M. Meuwly, Improving potential energy surfaces using measured Feshbach resonance states, Sci. Adv., 2024, 10, eadi6462.
101. B. Gazdy and J. M. Bowman, An adjusted global potential surface for HCN based on rigorous vibrational calculations, J. Chem. Phys., 1991, 95, 6309–6316.

---

Footnote

† Present Address: Roche Pharma Research and Early Development, Pharmaceutical Sciences, Roche Innovation Center Basel, F. Hoffmann-La Roche Ltd, Basel, Switzerland.

---