Describing proton transfer modes in shared proton systems with constrained nuclear–electronic orbital methods

Proton transfer is a fundamental chemical process that involves the movement of a proton from one protonation site to another. It is ubiquitous in nature and plays a vital role in a variety of chemical and biological processes, such as photosynthesis^1,2 and respiration.³ Proton transfer is often accompanied by energy transfer. In recent decades, with the ultimate goal of better harnessing energy, many experimental and theoretical efforts^4–9 have been devoted to understanding the mechanisms of proton transfer and studying proton transfer rates.

In theoretical investigations, a significant challenge for the study of proton transfer is the incorporation of nuclear quantum effects (NQEs),¹⁰ given that protons are the lightest nuclei with strong NQEs. NQEs often include zero-point effects and tunneling effects, both of which are related to the quantum delocalized nature of nuclei. NQEs are closely related to many anomalous properties of water systems^11,12 and some other hydrogen-bonded systems.^8,13 For instance, the quantum delocalization of protons in the active site of the enzyme ketosteroid isomerase (KSI) significantly stabilizes the deprotonated form of an active-site tyrosine residue, resulting in its anomalous acidity.¹⁴ Additionally, NQEs can be manifested in vibrational spectra, where a good description of NQEs will naturally describe anharmonicity, which is significant in many bond vibrations involving hydrogen atoms.^15–18 

There have been many methods developed to incorporate NQEs and anharmonicities in theoretical calculations. Based on high-dimensional potential energy surfaces (PESs), both static calculations using vibrational self-consistent field/vibrational configuration interaction (VSCF/VCI)^19,20 and dynamic simulations using multiconfiguration time-dependent Hartree (MCTDH)^21,22 can give highly accurate results. Although they are generally considered to be relatively expensive, recently both VSCF/VCI and MCTDH methods have been modified and applied to some extremely challenging systems with 10–60 atoms and successfully described the vibrational spectra.^23–25 More computationally affordable methods are vibrational perturbation theory (VPT2)²⁶ and path-integral based simulation methods, such as centroid molecular dynamics (CMD)^27,28 and ring polymer molecular dynamics (RPMD).²⁹ Compared to VSCF/VCI, VPT2 only needs local higher-order derivative information of the PES to perform calculations. Due to the relatively low computational cost, it has been widely used to compute anharmonic corrections with high accuracy. However, VPT2 faces challenges in highly anharmonic systems, for example, the shared proton systems in this Communication.^30,31 Path-integral based methods have higher costs than conventional classical molecular simulations, which depend on factors such as simulation temperature and acceleration algorithms.^32,33 They have been applied to investigate proton transfer in larger systems such as between phenol and trimethylamine in polar solvents³⁴ and through two-dimensional materials.³⁵ In general, path-integral methods are excellent for investigating static properties such as structures and free energies³⁶ and some dynamic properties such as reaction rates.³⁷ However, they have limitations in predicting some dynamic properties, particularly in calculating vibrational spectra, although recent developments have made them more reliable for these calculations.^{15–17,38,39} Recently, constrained nuclear–electronic orbital density functional theory (CNEO-DFT) has been developed in our group.^40–42 It simultaneously treats nuclei and electrons quantum mechanically but with a constraint on each quantum nuclear expectation position. The resulting CNEO-DFT energy surface is a function of quantum nuclear expectation positions as well as classical nuclear positions, and it inherently incorporates NQEs. The CNEO-DFT energy surface has many desirable properties. For instance, the harmonic vibrational frequencies obtained directly from Hessian calculations are comparable to the anharmonic frequencies obtained from DFT-VPT2.⁴² Additionally, molecular dynamics can be performed on CNEO-DFT energy surfaces (CNEO-MD), which incorporates NQEs in classical simulations with essentially the same computational cost as conventional DFT-based ab initio molecular dynamics (AIMD).^18,43 CNEO-MD has been tested and found to be capable of accurately reproducing the key features in vibrational spectra of a series of small molecules.¹⁸ Without any empirical scaling, CNEO-DFT and CNEO-MD are especially good at predicting the vibration of free X–H bonds, where X is a heavy element.^18,42

Although the low computational cost of CNEO-MD makes it directly applicable to the simulation of proton transfer processes in large-scale systems, the complexity of such systems and the interplay between different factors may make it hard to clearly assess the performance and understand the strengths and weaknesses of CNEO-MD. Fortunately, shared proton systems, in which a proton is shared between two neutral or charged moieties, have been widely investigated both experimentally and computationally.^8,44–47 As prototypical systems for the presence of protons in hydrogen-bonded networks, the shared proton systems can offer valuable insights into the nature of proton transfer processes in different chemical environment and are ideal systems for benchmarking. In aqueous solutions, several structures^9,48,49 can be involved in proton transfer processes with the Eigen form (e.g., $H9O4+$⁠)⁵⁰ and the Zundel form (e.g., $H5O2+$⁠)⁵¹ representing two limiting structures. In the Eigen structure, the proton is largely localized on the central molecule and it can be perceived as the starting and ending structure of the proton transfer process, whereas in the Zundel structure, the proton is shared equally by two adjacent water molecules, which is like a transition state for proton transfer. Out of the two, the Zundel structure is especially challenging since the shared proton lies on a highly flat and highly anharmonic PES.⁵² As a result, not only the conventional harmonic Hessian calculations but also the popular VPT2 method fail to describe the vibration spectra of the Zundel structure. In the past, more expensive methods have been applied to predict its vibrational spectrum, including diffusion Monte Carlo (DMC), MCTDH, and VSCF/VCI.^44,53

In this Communication, we will apply the CNEO method to investigate the vibrational spectra of shared proton systems or, more specifically, the challenging Zundel-like structures: $H5O2+$⁠, $H3O2−$⁠, and HCOO⁻·H⁺·⁻OOCH. We will start with a brief introduction on the theory of CNEO in Sec. II. Then, we will provide our computational details in Sec. III, followed by results and discussions in Sec. IV. We give our conclusions in the last section.

We finally note that although the CNEO theory is a multicomponent quantum chemistry theory and is closely related to the NEO theory, it is developed to tackle a completely different set of problems from the NEO theory. CNEO is designed for constructing adiabatic change to `zero-point energy corrected' effective PESs and is thus good for vibrational spectroscopy calculations as was shown in previous papers^18,42 and will be shown again in this Communication, whereas NEO is good at describing the nonadiabaticity between electrons and quantum nuclei. With quantum nuclei assumed to respond instantaneously to the motion of classical nuclei, the NEO theory is not directly applicable to vibrational frequency calculations,⁶⁰ although NEO-DFT(V)^61,62 and NEO Ehrenfest dynamics^63,64 have been developed to describe them, which are generally much more complex than the CNEO methods.^18,43

We choose three challenging Zundel and Zundel-like structures for vibrational spectra calculations: $H5O2+$⁠, $H3O2−$⁠, and HCOO⁻·H⁺·⁻OOCH. They all have the shared proton sitting at the center, and their experimental gas phase infrared (IR) spectra have been available through cold-ion vibrational action spectroscopy.^44,45,47 For comparison, we calculate their vibrational spectra using both DFT and CNEO-DFT methods. Static energy, gradient, and Hessian calculations are performed using an in-house version⁶⁵ of PySCF,^66,67 and MD simulations are performed with Atomic Simulation Environment (ASE).⁶⁸ In CNEO calculations, we only treat hydrogen nuclei quantum mechanically because heavier nuclei generally have smaller NQEs, although full-quantum calculations, which take into account NQEs for all nuclei, are also possible within the CNEO framework.^41,42 For all the calculations in this Communication, if not specified, the aug-cc-pVTZ basis set⁶⁹ is adopted for electrons and the PB4-D basis⁷⁰ is used for protons. Three electronic exchange–correlation functionals, viz., B3LYP,^71,72 PBE0,^73,74 and ωB97X,⁷⁵ are tested in both DFT and CNEO-DFT calculations. Because of the relative localized picture of quantum nuclei, we only include the nuclear–nuclear Coulomb interaction but no nuclear exchange–correlation. For electron–proton correlation,^76–78 we found that the most popular electron–proton correlation functional epc17-2⁷⁷ has negligible influence on the frequencies of the proton transfer modes (supplementary material, Table S1). Therefore, we only report results without electron–proton correlation, but we plan to perform more systematic tests in the future to find out the reason for the negligible influence. For additional comparison, DFT-VPT2 calculations are performed using Gaussian 16,⁷⁹ although we note that it has been well known in the literature that VTP2 fails badly for these shared proton systems.³¹ Since we find that DFT-ωB97X has the best agreement with CCSD(T)^80,81 in terms of harmonic frequencies for the proton transfer modes in $H5O2+$ and $H3O2−$⁠, we adopt ωB97X for MD simulations with minimal error from the choice of electronic functional. In principle, simulating IR spectra requires an ensemble average of a large number of NVE trajectories with starting configurations sampled from an equilibrated NVT trajectory. To reduce the calculation cost, we follow the practice established by existing literature^82,83 and use a single 10 ps NVT trajectory at 10 K to simulate the IR spectra. We set the initial positions and velocities of each atom following the equipartition law at 10 K with equal amount of energy added to each vibrational mode obtained from Hessian analysis.⁸⁴ In practical simulations, a time step of 0.5 fs is chosen and the Nosé–Hoover thermostat^85,86 with a large characteristic timescale being 1 ps is employed. In the end, the dipole autocorrelation function is calculated from each trajectory and then Fourier transformed to obtain the IR spectrum of each system.^82,83

$H5O2+$ and $H3O2−$ are the simplest water clusters with shared protons. They have been widely studied experimentally and theoretically.^{22,44–46,52,88–90} Their proton transfer modes are of particularly high research interest. In $H5O2+$⁠, the proton is shared by two water molecules, giving a symmetric structure.⁴⁴ Through geometry optimization, this symmetric structure can be predicted by both pure electronic structure calculations and multicomponent CNEO-DFT calculations (supplementary material, Fig. S1). In $H3O2−$⁠, the proton is shared by two OH⁻ anions, which still leads to a symmetric structure.⁴⁵ However, the geometry optimization result for $H3O2−$ strongly depends on whether the shared proton is treated classically or quantum mechanically. Most pure electronic structure calculations, except PBE0 functional, give a double-well ground state PES, and if the proton is considered to be a classical particle, it will locate closer to one of the OH⁻, essentially forming H₂O·OH⁻ (supplementary material, Fig. S2). In contrast, if path-integral simulations are performed based on the double-well potential, the quantum treatment of the proton brings the zero-point energy to the system, which is higher than the proton transfer barrier,⁹¹ and results in a proton density centered in the middle.⁴⁶ With CNEO-DFT, because the proton is treated quantum mechanically together with electrons, the proton is consistently predicted to be at the center and equally shared by the two OH⁻ anions regardless of the electronic functional used (supplementary material, Fig. S2).

Based on the individually optimized geometries, we obtain the vibrational frequencies of $H5O2+$ and $H3O2−$ by CCSD(T), DFT, DFT-VPT2, and CNEO-DFT. The results are presented in Table I. Note that the CCSD(T), DFT, and CNEO-DFT results are harmonic frequencies obtained from diagonalizing the mass-weighted Hessian matrices, whereas DFT-VPT2 performs VPT2 calculations at the DFT optimized geometries, which includes anharmonicity.

For $H5O2+$⁠, the experimental proton transfer mode frequency is 1047 cm⁻¹, while DFT harmonic calculations underestimate the frequency with all three functionals tested and the degree of underestimation greatly varies with the underlying functional. It seems that PBE0 is the best functional here, which matches the experimental frequency well. Nevertheless, we note that this seemingly great performance of PBE0 is fortuitous and primarily due to error cancellation: Based on the harmonic approximation, the CCSD(T) with the large aug-cc-pVQZ basis set predicts a frequency of 819 cm⁻¹,⁸⁷ which is more than 200 cm⁻¹ lower than the experimental reference. Since we can consider CCSD(T) to be mostly free from electronic correlation error, this huge mismatch of the gold standard method mainly comes from the classical nuclear treatment. Then, the difference between DFT-PBE0 and CCSD(T) mainly comes from the electronic error. Therefore, the good performance of DFT-PBE0 is accidently caused by error cancellation from the electronic part and the classical nuclear part. In this sense, the ωB97X functional is the best electronic functional for this system with similar performance to CCSD(T). For DFT-VPT2, although it is generally a reliable method for calculating anharmonicity, it fails badly for this system with the proton transfer mode frequency significantly overestimated by more than 300 cm⁻¹. In fact, it is known that the proton transfer modes in shared proton systems are especially challenging for VPT2 because the local higher-order derivative information is not sufficient to capture the behavior of the whole surface. In contrast, even with a harmonic approximation, CNEO-DFT is able to give reasonably good descriptions for the proton transfer modes with slightly overestimated frequencies. The frequencies by CNEO-DFT are also sensitive to the choice of electronic functional, but we find that ωB97X is the best-performing functional here, which overestimates the proton transfer mode frequency by ∼110 cm⁻¹. Considering DFT-ωB97X gives the closest results to CCSD(T) with the least error from the electronic part, this result indicates that CNEO-DFT performs better when the underlying electronic functional is more reliable, thus suggesting a good description of the quantum proton by CNEO-DFT.

As to $H3O2−$⁠, both CCSD(T) and DFT-B3LYP/ωB97X give asymmetric optimized geometries, thus there is only the bound OH stretch mode instead of the proton transfer mode. Because of this incorrect physical picture, all their vibrational frequencies are significantly higher than the experimental reference values. In order to obtain the proton transfer mode frequency, we additionally optimized the geometries requiring the proton to be in the middle, which essentially gives the transition state geometries for these methods. Based on these geometries, the real proton transfer modes can be captured, although their vibrational frequencies are imaginary. In contrast, DFT-PBE0 predicts a symmetric structure in which the shared proton sits in the middle at the optimized geometry. However, the vibrational frequency given by DFT-PBE0 is extremely low. These results show the high sensitivity to the choice of the underlying electronic structure theory but harmonic approximations based on conventional electronic structure theories all fail badly. Starting from an asymmetric geometry, DFT-VPT2 gives even more unphysical results with B3LYP and ωB97X functionals giving imaginary frequencies. Although VPT2 can be performed at a symmetric geometry using the PBE0 functional, it gives an extremely high proton transfer mode frequency. We note again that the unsatisfactory performance of VPT2 method in low-barrier double-well system is well known and it cannot cover the degenerate minima as well as the saddle point.³¹ In contrast, CNEO-DFT gives much more reliable results with a quantum mechanical treatment on nuclei. Even though the result is still sensitive to the choice of the electronic functional, we found that ωB97X, as the best functional that matches CCSD(T) most, manages to give the most accurate CNEO-DFT frequency compared to the experimental reference. All these results demonstrate that the quantum treatment on transferring proton is key to the qualitatively correct description of the proton transfer modes.

Among the three functionals tested, since ωB97X is the best one that has the least error from the electronic part, we use it for AIMD simulations. To differentiate conventional AIMD using DFT PESs from our AIMD using CNEO-DFT energy surfaces, throughout the Communication, we call them conventional AIMD and CNEO-MD, respectively. Figure 1 shows the IR predissociation spectra of the $H5O2+⋅Ne$ complex⁴⁴ and the simulated spectra for $H5O2+$ using conventional AIMD and CNEO-MD. According to the experimental spectra, the proton transfer mode has strong absorption at 1047 cm⁻¹. However, there is an additional strong absorption at 928 cm⁻¹. The capture and assignment of the 928 cm⁻¹ peak with theoretical simulations have been particularly challenging.^15,44 MCTDH is the most reliable method that has successfully reproduced this peak and characterized it as the Fermi resonance between the proton transfer mode and the combination band of the water wagging motion and O–O stretch motion.²² Conventional AIMD underestimates the proton transfer mode frequency, whereas CNEO-MD overestimates the mode frequency, which are consistent with the DFT harmonic results and CNEO-DFT harmonic results, respectively. Interestingly, for CNEO-MD, in addition to the main proton transfer mode absorption at about 1160 cm⁻¹, we additionally see a small satellite peak at 1106 cm⁻¹, which is not shown in CNEO Hessian analysis. Through normal coordinate analysis on the CNEO-MD trajectory using Trajectory Analyzer and Visualizer,^92–94 we find that both the main peak and the satellite peak have significant contributions from proton transfer motion and the wagging motion of water molecules (supplementary material, Fig. S4 and animations of the modes in GIF format). Therefore, CNEO-MD seems to have captured some of the mode coupling features, although the quantitative correspondence is not great and it remains challenging for classical MD methods to accurately describe overtones, combination bands, and Fermi resonances.⁸² Apart from the proton transfer mode, there are two peaks around 1400 cm⁻¹ corresponding to the perpendicular motions of the shared proton. Although the intensities of the two modes are low, CNEO methods can capture them with vibrational frequencies in good agreement with prior MCTDH calculations.²² Around 1600 cm⁻¹ and 1750 cm⁻¹, there are two water bending modes according to MCTDH calculations.²² The one around 1600 cm⁻¹ is a dark state while the intensity of the other mode around 1750 cm⁻¹ could vary significantly in experiments depending on the type of noble gas tag.⁴⁴ CNEO methods can accurately describe the frequencies and intensities of both modes. However, we also see a peak splitting for the bright higher energy mode, which is not expected. Using normal coordinate analysis based on the MD trajectory, we find that although the mode is dominated by water bending, it still has substantial contribution from the proton transfer motion, which can also be coupled with the water wagging motion (supplementary material, Fig. S5 and animations of the modes in GIF format). This may be the reason for the peak splitting. Because this splitting is not observed either in experimental or MCTDH spectra, it could be an artifact of our method or it could correspond to the small peak around 1850 cm⁻¹ in experiment. In the even higher frequency range (3400–4000 cm⁻¹), we can see the symmetric and asymmetric free OH stretches. CNEO-MD describes these modes excellently, which is not surprising considering we have observed the great performance of CNEO-MD in gas phase water molecules.¹⁸ 

For the $H3O2−$ system, Fig. 2 presents the IR predissociation spectra of the $H3O2−⋅Ar$ complex⁴⁵ and the simulated spectra for $H3O2−$ using conventional AIMD and CNEO-MD. The IR spectrum of $H3O2−$ is simple with most intensities on the proton transfer mode. Because DFT predicts an asymmetric placement of the shared proton, conventional AIMD gives a bound OH stretch slightly above 1800 cm⁻¹ and thus completely fails to reproduce the experimental proton transfer mode around 700 cm⁻¹. We also simulated IR spectra using conventional AIMD starting from saddle point geometries at 10 and 100 K (supplementary material, Fig. S6). When the temperature is set at 10 K, the system quickly gets trapped into one of the local wells and its spectrum resembles that starting from the optimized geometry. While at 100 K, the shared proton possesses enough energy to overcome the barrier but cannot maintain a continuous proton transfer picture due to energy fluctuation. As a result, neither of these simulations can give significant peaks around the experimental proton transfer frequency (supplementary material, Fig. S6). In contrast to the failure of conventional AIMD in H₃O₂₋, CNEO-MD correctly captures the proton transfer mode at around 700 cm⁻¹.

The good performance of CNEO-MD in $H5O2+$ and $H3O2−$ demonstrates its great potential as an accurate and efficient method for proton transfer studies, especially in large systems. Here, we extend our method to another important category of shared proton systems, proton-bound organic molecules/ions, and use the proton-bound formate dimer (AHA⁻) as an example to study.

The AHA⁻ system was previously studied by cold-ion IR action spectroscopy and discrete variable representation (DVR) analysis, both of which revealed that the bridging proton is shared equally between two formate anions.⁴⁷ It was also found that the fundamental of the proton transfer mode ν_‖ is strongly coupled to the out-of-phase carboxylate deformation mode, δ(OCO)_OOP.⁴⁷ Similar to the $H5O2+$ and $H3O2−$ cases, we perform both static calculations and molecular dynamics simulations to the AHA⁻ system.

With pure electronic DFT, the equilibrium structure of AHA⁻ is highly sensitive to the functional choice. Specifically, B3LYP and ωB97X give a shallow double-well PES and predict the proton to be on either side of the complex, whereas PBE0 gives a single-well picture and predicts the proton to be equally shared (supplementary material, Fig. S3). In contrast, the CNEO-DFT geometry is much less sensitive to functional choice with all three functionals predicting the bridging proton to be in the middle of two formate anions (supplementary material, Fig. S3). These results again demonstrate the necessity of the quantum treatment on the shared proton.

Based on the individually optimized geometries, we calculate the harmonic vibrational frequencies with DFT and CNEO-DFT. The results are shown in Table II. For comparison, we additionally include DFT-VPT2 results and literature results using the DVR analysis in two-dimensional (2D) and four-dimensional (4D) PESs.⁴⁷ Note that because pure electronic DFT-B3LYP and DFT-ωB97X calculations break the C₂ symmetry of AHA⁻, leading to completely different normal modes that cannot correspond well to the experimental results, here we do not present their poor results. In AHA⁻, there are two modes that have significant contributions from the proton transfer mode (supplementary material, Fig. S7). The reason is that the proton transfer mode ν_‖ and a hybridized carboxylate deformation mode δ(OCO)_OOP have the same B symmetry within the C₂ point group and thus can linearly combine to form two new modes, both of which involve large displacement of the shared proton along the O–O direction.⁴⁷ These two coupled modes are denoted as mode 1 and mode 2 in Table II and Fig. 3. Pure electronic DFT with the PBE0 functional can predict the symmetrically shared proton picture, but it greatly underestimates the two proton transfer modes by about 200 cm⁻¹. On top of the PBE0 reference, performing vibrational perturbative calculations with VPT2 makes the mode descriptions even worse, which, conversely, dramatically overestimates the vibrational frequencies (1000 cm⁻¹ overestimation for mode 1 and 300 cm⁻¹ overestimation for mode 2). Again, this poor performance of VPT2 is not surprising and we have observed similar performance of it in $H5O2+$⁠. In both cases, because of the high anharmonicity and high complexity of the underlying PES, local higher-order derivative information is not sufficient for an accurate description of the whole PES, leading to the failure of popular VPT2 method. In contrast, CNEO-DFT with harmonic approximations can even more accurately predict the vibrational frequencies of the two modes and give similar results across different functionals. Specifically, the errors of B3LYP and ωB97X are all within 35 cm⁻¹ for mode 1 and within 60 cm⁻¹ for mode 2 while PBE0 has larger errors when compared to experimental results. In fact, as a highly efficient method, CNEO-DFT also outperforms the DVR results based on a 2D PES and has a similar or slightly better performance than DVR results based on a 4D PES when two additional modes are included in calculations, all of which demonstrate the high competence of CNEO-DFT in describing the proton transfer modes in larger and more complex systems.

To directly compare with experimental spectra, we also simulated the vibrational spectra of AHA⁻ using both conventional AIMD and CNEO-MD. Without harmonic CCSD(T) results, we do not know which functional is the best for the electronic part description. Therefore, to be consistent with the $H5O2+$ and $H3O2−$ cases, we continue using the ωB97X functional for MD simulations. The simulated spectra for AHA⁻ along with the experimental Ar predissociation IR spectrum⁴⁷ are presented in Fig. 3. The experimental spectrum is highly complex, and many of the peaks may be contributed by combination bands or overtones, which are challenging for classical MD methods. The key features of the spectra are the two proton transfer modes with high intensities, and CNEO-MD can capture them both. Note that mode 1 matches the experimental frequency well. Although the match is not perfect for mode 2 with 60 cm⁻¹ underestimation, CNEO-MD still gives the best match to the experimental peaks among the methods that have been applied to this system. For conventional AIMD, although there is a mode at 681 cm⁻¹, which seems to match reasonably with mode 1, it is actually a HCOOH deformation mode rather than the proton transfer modes (supplementary material, Fig. S8). Therefore, the physical picture of the conventional AIMD is completely wrong. We also performed conventional AIMD simulations starting at saddle point and its performance is similar to the $H3O2−$ case. At 10 K, the proton quickly gets trapped at a local well while at 100 K, the higher temperature makes it possible for the shared proton to transfer but cannot maintain a continuous proton transfer picture due to energy fluctuation. Therefore, neither of them shows significant peaks around the experimental proton transfer frequencies (supplementary material, Fig. S10). In CNEO-MD, we observe a small satellite peak that is slightly higher in frequency than the main peak of mode 2. Although the satellite peak seems to match with CNEO Hessian peaks at 1055 cm⁻¹, we note that the Hessian 1055 cm⁻¹ peak corresponds to free OH bending and is not IR active. Therefore, the satellite peak does not correspond to the free OH bending. Instead, through normal coordinate analysis based on the MD trajectory, we find that this peak is indistinguishable in nature from the main peak of mode 2 (supplementary material, Fig. S9) and could be an artifact of our method. Another important feature of AHA⁻ is the fundamentals of $CO2−$ deformations with peaks showing up as a doublet (1685 and 1700 cm⁻¹) separated by 15 cm⁻¹. These frequencies are between the frequency of the carboxylate antisymmetric stretch (1622 cm⁻¹)⁹⁵ and the frequency of the carboxylic acid CO stretch (1818 cm⁻¹),⁹⁶ which, together with the small doublet peak splitting, is an experimental evidence for the proton being equally shared by the two formate anions. CNEO-DFT harmonic analysis predicts the doublet to be at 1772 and 1789 cm⁻¹. Although the frequencies are overestimated, the splitting of 17 cm⁻¹ is consistent with the experimental results. In CNEO-MD spectrum, the two peaks merge into one because of the small splitting between them.

In summary, with three prototypical shared proton systems, $H5O2+$⁠, $H3O2−$ and AHA⁻, we demonstrated that our CNEO methods can accurately and efficiently describe the proton transfer modes in shared proton systems. In all three systems, CNEO-DFT accurately predicted the proton to be shared equally by the two ending groups whereas pure electronic DFT with classical treatment on nuclei often cannot. Based on the physically correct geometries, both CNEO harmonic analysis and CNEO-MD can well describe the key features related to fundamental transitions in the IR spectra, which is greatly superior to conventional DFT harmonic analysis, DFT-VPT2, and conventional DFT-based AIMD. Although the descriptions of overtones, combination bands, and Fermi resonances are still challenging and high-level quantum methods, such as MCTDH and VSCF/VCI, are still needed for accurate description of these subtle features, our CNEO methods excel in their abilities to accurately and efficiently capture the key fundamental transition natures in relatively large systems. Therefore, CNEO methods are promising for future applications to larger and more complex systems like biological systems and materials systems in which proton transfer plays a crucial role.

See the supplementary material for DFT and CNEO-DFT optimized geometries of H5O2+, H3O2−, and HCOO−·H+·−OOCH (Figs. S1–S3); comparison of proton transfer mode frequencies calculated with CNEO/no epc and CNEO/epc17-2 (Table S1); normal coordinate analysis for H5O2+ and HCOO−·H+·−OOCH (Figs. S4, S5, and S9 and animations in gif format); conventional AIMD of H3O2− and HCOO−·H+·−OOCH starting with saddle point geometries (Figs. S6 and S10); proton transfer modes of HCOO−·H+·−OOCH (Fig. S7); and computational cost comparison (Table S2).

The authors thank Mark A. Johnson for providing experimental spectra for all three systems. The authors are grateful for the funding support from the National Science Foundation under Grant No. 2238473 and from the University of Wisconsin via the Wisconsin Alumni Research Foundation.

The authors have no conflicts to disclose.

Yuzhe Zhang: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Xi Xu: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (equal); Supervision (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Nan Yang: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Writing – review & editing (supporting). Zehua Chen: Methodology (supporting); Supervision (supporting); Writing – review & editing (supporting). Yang Yang: Conceptualization (lead); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Project administration (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.