The role of exact exchange on the structure of water dimer radical cation: Hydrogen bond vs hemibond

Electronic and structural properties based on density-functional theory (DFT) have been used successfully in molecular and condensed-matter studies. The generalized gradient approximation (GGA) has frequently been employed in DFT applications for exchange and correlation.¹ Although this approximation is numerically efficient, it suffers from serious drawbacks. A poor description of chemical reaction barriers and dissociating energies of molecular ions may result from the delocalization error of approximate functionals.^2,3 This shortcoming can be partly cured by mixing the exchange–correlation functional with some percentage of exact exchange, which makes it typically superior to GGA at describing structural and electronic features.^4–6 

One important example where most DFT functionals fail seriously is water cluster cations $[(H2O)n+]$ in iodized water.⁷ These clusters are obtained in many chemical reactions and have a lifetime of several tens of femtoseconds⁸ after which water dimer radical cations $(H2O)2+$ are formed.⁹ The ionization of water was experimentally produced under special conditions in the gas phase.^9–11 Furthermore, oxidized water in the form of HO^· radicals has also been observed in oxidation processes over supposedly inert electrode materials, such as glassy carbon^12,13 or boron doped diamond.^14–16 These oxidized water clusters may also be produced through Fenton’s reaction where hydrogen peroxide is decomposed catalytically through the addition of Fe(II) salts.^17,18 Some theoretical investigations for water cluster cations suggested the hemi-bonded structure where the dimer is stabilized through a weak bond between the oxygen atoms in a peroxide-like fashion.^19,20 Other work proposed a proton-transfer between the water molecules to form a H₃O⁺ oxonium ion and a HO^· radical.²¹ In this complex, the radical is somewhat stabilized through a strong hydrogen bond between both molecules. Most DFT functionals, such as B97D, M06L, BLYP, and BPW91, predict that the hemi-bonded structure is more stable than or at least comparable to the proton-transfer structure.^22–24 However, experimental results and high-level theories, such as CCSD(T), show that the proton-transfer structure is much more stable by 0.51 eV than the hemi-bonded structure.^23,25,26 Thus, the diffusion mechanism via a hydrogen exchange reaction might not be fully captured with the conventional exchange–correlation density functionals due to the overestimation of the hemibond strength.²⁷ 

In 1988, for example, the proton-transfer (C_s symmetry) and the hemi-bonded (C_2h symmetry) structures were characterized as two minima at the MP4/6-311G(MC)**//MP2/6-31G* level of theory by Gill and Radom.²⁸ They reported the hemi-bonded structure as metastable and predicted it to be less stable than the proton-transfer structure by 0.38 eV with a barrier of 0.34 eV between these two minima. Based on these results, they suggested that the hemi-bonded structure may also be sufficiently stable to be detected experimentally. A few years later, the ionized water dimer $(H2O)2+$⁠, including two C₁ minima (one hydrogen-bonded structure and one hemi-bonded structure) and three transition states (C_s symmetry), was studied by Sodupe et al. at the SCF and MP2 level of theory.²⁹ They found a transition state that was structurally similar to the proton-transfer structure but with out of the plane rotation of OH. It was estimated that there would be a negligible energy difference between these two configurations. By releasing the symmetry constraints, they found the minimum three-electron hemibonded complex with a covalent bond and a complete delocalization of the positive charge.

In 2009, ab initio coupled-cluster electronic structure theory with single and double excitations (CCSD) and CCSD with perturbative triple excitations CCSD(T) were employed in order to investigate fourteen stationary points for the water dimer radical cation.²⁶ They showed that stationary point geometries and energetics depend on basis sets, and at the highest level of theory, aug-cc-pVQZ CCSD(T), the global minimum of the water dimer radical cation is a hydrogen-bonded system, which lies 0.31 eV lower in energy than the hemi-bonded isomer. The energy barrier between these two structures was reported to be 0.65 eV. Lee and Kim²³ also investigated the structures, energetics, spectra, and dynamics of the water dimer cations using CCSD(T) calculations compared with the DFT results, conveying the closeness of the DFT results with MPW1K and BH&HLYP functionals to the CCSD(T)/CBS results. Although in all theoretical studies the water dimer cation dissociation is likely to occur via the oxonium channel, the hemibond, as evidenced by delocalization of the spin density, is present but obscured by the hydrogen-bonded feature.³⁰ Simulations also provide evidence that the substantially altered ultraviolet absorption spectra of the aqueous radical indicate the presence of hemibonded OH_(aq).³¹ The recent study³² offers comprehensive and high-level computational findings on the structures and relative energies of diverse types of water dimer cations, hemibonded form included. The investigation indicates that the detection of hemibonded structures in trajectory analyses of ionized water clusters or liquid is straightforward and that hemibonded cations may emerge as a competing channel with proton transfer during the radiolysis of regular water. Additionally, the study provides prognostications of the optical absorption and relative resonance Raman spectra intensities of the hemibonded dimer cation, offering crucial insights for experimental examinations.

In the present research, a systematic theoretical study employing the DFT and high-level ab initio theory is carried out to investigate structural and electronic properties of the water dimer radical cations. We will characterize the structures of the two local minima and examine the pathway between them. These results are then used to identify the most suitable DFT functionals. Furthermore, the effect of exact exchange in exchange–correlation functionals is addressed in order to compare DFT with CCSD(T) levels of theory. Our results will serve as a fundamental criterion in computing the ionized water clusters with approximate models and offer insights into the likely channel for dissociation of the water dimer cation.

Computations were performed using Gaussian 16 (Revision C.01).³³ All DFT computations were performed using the Def2-TZVPD basis set in combination with the SMD³⁴ implicit solvation model. The functional dependence was tested by comparing results obtained with B97D,³⁵ M06L,³⁶ B3LYP-D3,^37–39 PBE0-D3,^40,41 B3PW91-D3,³⁸ M06-2X,⁴² BHandH,⁴³ M06-HF,^44,45 HSE06,⁴⁶ and ω-B97XD.⁴⁷ For functionals with less than 25% exact exchange, a Grimme D3 correction⁴⁸ for dispersion interactions was added. Transition states were computed through following the appropriate normal mode and identified by the presence of a single imaginary mode, which connects reactant and product states. These data were compared to single-point CCSD(T) computations. Here, a quadruple-ζ correlation consistent aug-cc-pvqz basis set was used. Computations were again performed with the SMD³⁴ implicit solvation model.⁴⁹ Structures, entropy, and zero-point energy corrections were taken from B3LYP-D3/Def2-TZVPD computations. Electrochemical potentials are computed using the effective absolute potential method in which an effective absolute potential is computed for each computational setup and used to map the absolute potentials to the standard hydrogen electrode (SHE) scale in water.⁵ This method was found to allow for accurate prediction of electrochemical potentials due to favorable error cancellation.

The stability of neutral water dimers was evaluated through the above-described methods in Gaussian 16. However, owing to their instability, only single point computations at the structures obtained with PBE0r were performed. Changes in the zero-point energy and entropy contributions were neglected since we were unable to perform a vibrational analysis at the respective converged structures. This is not expected to change the results significantly since no major changes of these properties are expected for the three structures.

Note that the dispersion correction should be paired with all functionals to take long-range interaction into account, otherwise missing physics. However, an empirically parameterized meta-GGA functional improves dispersion energies and reduces mean errors for both kinds of weak bonding. Thus, this functional class captures more of the correlation relevant for weak and nonlocal interactions.^42,50 For the BHandH functional, there is no implementation of D3 correction in Gaussian 16.

DFT-PBE0r local hybrid functional⁵¹ calculations were performed using the Car–Parrinello projector augmented-wave (CP-PAW) code,⁵² employing the projector augmented-wave (PAW) method.⁵³ The PBE0r functional, which has been recently published, is a type of local hybrid exchange–correlation functional that is based on PBE0. The Kohn–Sham orbitals are transformed into a small set of tight-binding orbitals that are centered on individual atoms. These tight-binding orbitals are then used to determine the onsite HF exchange terms, which include the exchange interaction between the core and valence electrons. However, all other exchange contributions that involve tight-binding orbitals centered on different atoms are disregarded. Therefore, PBE0r can be classified as a range-separated hybrid functional, with the range of the exchange interaction being determined by the localized tight-binding orbitals. The “r” in PBE0r indicates that it uses a range-separated functional. The Fock admixture of the hybrid functional a_x was set to 0.1 for O and H. The augmentation of the PAW method included the 1s orbital of H and the 2s, 2p, and 3d orbitals of O. The auxiliary wave functions were constructed as nodeless partial waves. The plane wave cutoff for the auxiliary wave function is 544 eV and 1088 eV for the density. A k-point grid was set to 2 × 2 × 2 points. We specified a face-centered cubic unit cell with a lattice constant of 10 Å to avoid overlapping the wave functions with those from periodic images. In addition, a technique for removing the electrostatic interaction of periodic images has been applied to the calculations.⁵⁴ This approach calculates the electrostatic interaction between the periodic images and subtracts its energy from the Hamiltonian. Calculations are performed in the vacuum without any implicit solvation model. Spin-polarized calculations were carried out for the ionized system, and all atomic positions were optimized without symmetry constraints.

The water molecule contains ten electrons, of which eight are from the oxygen atom and one from each of the hydrogen atoms. Two of these electrons occupy the O-1s orbital tightly bound to the oxygen nucleus ([core]). The remaining eight electrons fill the states corresponding to the O-2s and O-2p orbitals interacting with H-1s orbitals. Figure 1 shows the calculated molecular orbitals in the ground state of the isolated water molecule. The electronic ground state has the electron configuration of [core]$(2a1)2(1b2)2(3a1)2(1b1)2(4a1)0$⁠. The lowest energy orbital $(2a1)2$ is formed by the 1s orbitals of the hydrogen atoms and the 2s orbital of the oxygen atom. On top of that, the $(1b2)2$ orbital is formed from O-2p_x and H-1s. At higher energy, $(3a1)2$ has the contribution of O-2p_y and H-1s. The highest occupied molecular orbital (HOMO), $(1b1)2$⁠, is predominantly O-2p_z with no contribution of the H-1s orbital. The HOMO is non-bonding and highly localized on the oxygen atom without sp³ hybridization characteristic. The three orbitals 2a₁, 1b₂, and 3a₁ are O–H bonding orbitals. The lowest unoccupied molecular orbital (LUMO) 4a₁ contributes to the O–H antibonds.

For the isolated molecule, the equilibrium bond length is 0.987 Å, and the bond angle of H–O–H is 104.8° extracted from our PBE0r calculations. The experimental values for gaseous water molecule are 0.957 Å for O–H bond length and 104.5° for H–O–H bond angle.⁵⁵ Our calculated O–H bond length is about 3% larger than the actual value, and the bond angle agrees with the experimental value. The calculated HOMO–LUMO gap is about 5.61 eV. At the CCSD(T) level, the O–H bond length and H–O–H bond angle of the water monomer are 0.962 Å and 104.2°, becoming more accurate and closer to the experimental values.²⁶ A systematic series of simulations for the O–H bond length was reported by Leung and Rempe.⁵⁶ Santra et al.⁵⁷ investigated the deformation energies calculated with PBE, BLYP, and PBE0 functionals, showing that the symmetric stretching of the O–H bonds of a water monomer is underestimated by PBE and BLYP, while it is accurately given by PBE0.

A distinct type of chemical bonding responsible for binding the water molecules is the hydrogen bond, which arises through the interaction of a hydrogen atom lying between two electronegative O atoms. In the hydrogen bond, both molecule fragments are connected through a 3 center 4 electron (3c4e) bond.⁵⁸ The strength of the hydrogen bond is between a covalent bond and the van der Waals interaction. The simple example of the hydrogen bond is that in water dimer (H₂O)₂. The optimized configurations of the neutral water dimer are shown in Fig. 2. The three configurations are classified as a linear water dimer with a C_s symmetry, a cyclic water dimer with a C_2h symmetry, and a bifurcated water dimer with a C_2ν symmetry molecular geometry. The lowest energy configuration of these dimers has a mirror plane of C_s symmetry in the equilibrium, in which the dipole moments of the two water molecules point as far as possible in opposite directions. At a short distance between O atoms (⁠$<2.66$ Å), the cyclic water dimer with a doubly hydrogen bonded structure has been found to be more stable than the C_s dimer.^26,55 The C_2ν dimer is typically about 31 meV less stable than the C_s dimer.^59,60 It is a transition structure for the interchange between the molecules being the hydrogen donor or acceptor.

The relative energies for the three different water dimer configurations (H₂O)₂ are listed in Table I. Our PBE0r calculations indicate that all three dimers should have almost identical binding energies. Interestingly, the calculated bandgap varies by 0.8 eV between the different structures despite their energetic similarities with the linear water dimer displaying the lowest bandgap. These results are, however, not reproduced by our Gaussian 16/DFT and CCSD(T) computations. These indicate that the linear water dimer is roughly 0.1 eV more stable than the cyclic or bifurcated configurations. The latter two configurations are energetically similar. The BHandH, ω-B97XD, and HSE06 functionals predict the highest energy for the bifurcated structure due to the overestimation of the hydrogen bond strength.⁶¹ It is noteworthy that CCSD(T) and all considered functional irrespective of the amount of exact exchange make identical predictions (Table I). This clearly indicates that this structure is well described by any electronic structure method, provided that hydrogen bonding is treated sufficiently well. Indeed, the lack of a suitable correction for the dispersion interaction in PBE0r computations may be the origin of the qualitative differences between both sets of computations. The calculated HOMO–LUMO gap on the other hand varies significantly with respect to the choice of the functional. In line with previous studies, functionals with 50% of exact exchange are in excellent agreement with an experimental value of 11.2 eV,^62,63 estimated from the experimental ionization potential and vertical attachment energy and other theoretical approaches.⁶⁴ The HF and CCSD(T) severely overestimate this energy by almost 3 eV and give identical orbital energies. These quantitative differences between the functionals do, however, not affect the qualitative trends between the three water structures. For example, the HOMO–LUMO gap is almost identical for all three structures with the linear configuration showing the lowest gap. Among the three structures, the linear isomer is found to have a much stronger hydrogen bond of about 1.95 Å, while the cyclic and bifurcated isomers show two weak hydrogen bonds of about 2.23 and 2.52 Å, respectively. The relative stabilities and bandgap correlate with the strength of the hydrogen bond.

The electronic ground state of the linear neutral water dimer (C_s) has the electron configuration of [core](3a′)²(4a′)²(1a″)²(5a′)²(6a′)²(7a′)²(8a′)²(2a″)². The 2a″ (HOMO) orbital represents a free electron pair formed by the O-2p_z orbital perpendicular to the molecular C_s plane at the O site. The HOMO is located at the water molecule whose oxygen atom does not receive a hydrogen bond. The 8a′ orbital mainly consists of the O-2p orbital in the C_s plane, with the lower orbital energy than that of the 2a″ orbital. The other orbitals are formed from the 2a₁ orbitals (3a′ and 4a′) and the overlap between the 3a₁ orbital and the 1b₁ orbital. Of particular importance is the molecular orbitals responsible for the hydrogen bond. The two orbitals overlapping across the hydrogen bond support the partial covalent nature of the hydrogen bond. The bond distance of the O⋯H hydrogen bond is 1.984 Å in DFT-PBE0r calculations, close to 1.976 Å of CCSD/aug-cc-pVTZ calculation³² and 1.951 Å of the CCSD(T)/aug-cc-pVQZ level of theory.²⁶ 

The geometry of an ionized system may be different from the neutral one. The ionization energy can refer to the vertical ionization energy or to the adiabatic ionization energy. Vertical ionization energy is defined in the same geometry as the neutral system, while adiabatic ionization refers to the ionization energy after relaxation of the geometry. The vertically ionized water dimer cation displays an O–O distance of 2.91 Å, and its two hydrogen-bonded monomers lie in perpendicular planes. Dissociating the vertical cation into H₂O⁺ and H₂O requires an energy of 0.99 eV. For the linear water dimer, the vertical ionization energy is 13.77 eV and the bandgap increases to 16.70 eV in the HF/CCSD(T) level of theory. Upon transferring a proton to the other oxygen atom, the proton-transferred structure becomes stable without a barrier, reducing the O–O distance to ∼2.50 Å. The EOM-IP-CCSD/aug-cc-pVTZ³² and CCSD(T)/aug-cc-pVQZ²⁶ methods give the O–O distance 2.47 and 2.51 Å, respectively. Stabilizing this new structure requires about 1.04 eV of energy. Furthermore, dissociating the proton-transferred structure into H₃O⁺ and OH necessitates 0.95 eV of energy. The obtained radical cation is in a doublet state, which has a total spin of 0.5 ℏ.

Two minima for the water dimer radical cation have been identified: a hydrogen bonded $OH⋅⋯OH3+$ adduct [Fig. 3(a)] and the hemibonded structure [Fig. 3(c)]. Many exchange–correlation functionals in DFT calculations fail to predict the presence of hydrogen bonding interaction²⁶ owing to the serious self-interaction errors in density functional approximations.

The hydrogen-bonded equilibrium structure of the water dimer radical cation has the electron configuration of [core](3a)²(4a)²(5a)²(6a)²(7a)²(8a)²(9a)²(10a). The single electron occupying the 10a orbital is predominantly O-2p character from the HO^· radical monomer. The hemibonded structure of the water dimer radical cation has the electron configuration of [core](2a)²(2b)²(3a)²(3b)²(4a)²(4b)²(5a)²(5b), with the O-2p character on both oxygens in the single electron occupying the (5b) orbital. The rotation of a water molecule in the hemibonded structure delocalizes the positive charge between both oxygens, while in the hydrogen bonded structure, the positive charge is located at one of the oxygen atoms and attracts a proton to make a hydronium H₃O⁺. The orbital interaction diagrams of hemibonded and hydrogen bonded configurations of a water dimer radical cation are shown in Figs. 4(a) and 4(b), respectively.

In the hemibonded structure, the oxygen atoms point directly toward each other, while in the hydrogen bonded structure, a hydronium points to an OH radical through a hydrogen bond. As given in Table II, the O–O distances R_OO are about 2.1 and 2.5 to 2.6 Å for the hemibonded and hydrogen bonded structures, respectively. Energetically, the hydrogen bonded structure is predicted to be 0.48 eV lower in energy than the hemibonded structure at our highest level of theory CCSD(T). The EOM-IP-CC(2, 3)/6-311++G**³² and CCSD(T)/aug-cc-pVQZ²⁶ methods predicted that the conversion of the hemibonded structure to the hydrogen bonded structure would give about 0.39 eV, including zero point vibrational energies. However, DFT-based calculations, even with the PBE0r functionals, including 10% of the exact exchange, predict the ground state of the water dimer radical cation to be a hemibonded structure by 0.29 eV. Thus, to find an accurate hybrid mixing factor, the energy difference between the hemibonded and hydrogen bonded structures is plotted as a function of the exact exchange percentages, as shown in Fig. 5.

The energy difference between two minima correlates linearly with the exact exchange percentage (see Table III). B97D, such as PBE functionals, predicts the hemibonded structure to be most stable. M06-L on the other hand predicts, despite not possessing any exact exchange, both structures to be equally stable. This is due to the more accurate level of theory in the meta-GGA functional. With 20% of the exact exchange, B3LYP-D3 and B3PW91-D3 stabilize the hydrogen bonded structure by ≈0.1 eV. A comparable result is also observed for PBE0-D3 with 25% of the exact exchange. Following this trend, BHandH and M06-2X, which both possess 50%–60% exact exchange, finally clearly favor the hydrogen bonded structure by 0.4 eV. This prediction is comparable to what is observed at the CCSD(T) level of theory. M06-HF pushes the energy difference even further and predicts an energy difference of more than 0.6 eV. Thus, it is less accurate for predicting both energetics and bandgaps. Overall, hybrid functionals with a larger fraction of the exact Hartree–Fock exchange can more accurately reproduce the results obtained with CCSD(T), while functionals with lower amounts of exact exchange (such as PBE0-D3 or B3LYP-D3) are insufficient to eliminate delocalization of unpaired spins. Various systems containing hydrogen bonds were also investigated using a range of functionals, including dispersion corrections, which indicated that none of the studied functionals are more accurate than MP2 for the investigated species, and the M06 + D3 functional performed well and gave minor errors when long-range dispersion was included.⁶⁵ 

The GGA functionals tend to overestimate the energy of the hemibond due to the reduction of (artificial) self-interaction, exploiting low-lying vacant orbitals, namely, delocalization error. Consequently, the hemibond isomer wrongly becomes the ground-state structure of $(H2O)2+$ with the M06L and B97D functionals. Correcting the self-interaction through the admixture of exact exchange rigorously reduces this error. Our calculations show that a 50% mixing in BHandH and M06-2X provides better agreement in determining the proton-transferred isomer to be the ground-state structure consistent with the CCSD(T) level of theory. The energy difference between the two isomers is about 0.48 eV. Thus, delocalized three-electrons systems need more careful density functional calculations to achieve accurate predictions of the ground-state structures.⁷ However, the functionals with 50% of exact exchange still underestimate the bandgap by several eV. To correctly predict this property, functionals with a higher percentage of exact exchange are required. Note that using a functional with a higher percentage of exact exchange leads to a severe overestimation of the stability of the hydrogen bonded structure. This makes it unsuitable for predicting reaction energies.

The potential energy surface has been calculated for both minima taking into account the R_OO bond distance, as shown for B97D, B3LYP-D3, and M06-2X calculations in Fig. 6. All calculations obtain two minima as expected for the hemibonded and hydrogen bonded structures, given in the supplementary material. The system remains in the same configuration changing the R_OO bond distance. The approach coordinate translates into electron interaction, while far away (stabilizing hole), there is just repulsion, which reduces by having a hole in between. Repulsion grows again as the molecules approach together more than the ground-state configurations. Interestingly, hybrid methods not only predict the right ground-state structure but also give an accurate R_OO for both minima.

The two minima of the water dimer radical cation are connected with the isomerization transition state, as shown in Fig. 3(b). The transition from ground-state configuration to the local minimum, where the intrinsic reaction coordinate method was used to track the proton transfer reaction, primarily involves hydrogen transfer from one O atom to another. The reaction path shows that the hydrogen donor water rotates its hydrogen atom toward the hydrogen acceptor water. This connects to the transition state. Then, the rotated hydrogen moves toward the oxygen of the hydrogen acceptor water. An energy barrier of 0.58 eV is calculated for the PT transfer process involving a hemibonded structure, which is close to 0.57 eV of the EOM-IP-CCSD/aug-cc-pVTZ³² level of theory and somewhat lower than 0.66 eV of the CCSD(T)/aug-cc-pVQZ²⁶ method. It should be noted that the O–H⋯O angle changes are used to describe the transition state properly.

Interestingly, all functionals give the same energy for moving from hemibonded configurations to a transition state, as given in Table IV. However, the energy difference between hydrogen-bonded configuration and transition state enormously depends on the exact exchange, and higher admixture produces a larger energy difference. Furthermore, the calculated binding energies demonstrate a negligible influence of the exact exchange on the hydrogen bond strength, while the hemibond strength increases by 1 eV, as listed in Table IV.

The most reliable CCSD(T) results show that the transition state is about 0.58 eV, while DFT functionals without exact exchange admixture, such as B97D and M06L, mispredict the stable structure and energy difference between the transition state and proton transferred structure. Even with 25% of the exact exchange, the energy difference between the transition state and proton transferred structure is underestimated by 0.2 eV. However, the M06-2X and BHandH calculations are closer to the CCSD(T) results for stabilizing the proton transferred configuration and getting reasonable transition state energy, as shown in Fig. 7. Therefore, the DFT functionals with 50% exact exchange are an alternative to expensive ab initio method for calculating water cluster energetics.

The two minima of the water dimer radical cation are investigated by adding an extra hole into the systems. In the proton-transferred structure, the hydrogen bond between the HO^· radical and H₃O⁺ breaks down, and the R_OO distance increases caused by repulsion. In the hemibonded structure, the second hole is inserted between the two oxygen atoms. This corresponds to the removal of one electron from the σ* O–O antibonding orbital [Fig. 8(a)]. Accordingly, a doubly protonated hydrogen peroxide with a true O–O bond (bond distance: ≈1.5 Å) is formed rather than a water dimer. The $(H2O)22+$ configuration is shown in Fig. 8(b). Under experimental conditions, this structure can be expected to rapidly release its two additional protons to form H₂O₂.

DFT calculations based on GGA functionals, extensively used to describe the water clusters, are, owing to large errors in the energetics of the oxidized water dimer and its electronic structure, clearly unable to properly describe even this simple system. A detailed analysis of the two minima for the water dimer radical cation, namely, the hydrogen-bonded and hemi-bonded geometries, shows that the DFT-GGA wrongly predicts the hemi-bonded configuration as a stable one, overestimating the stability of the symmetric H₂O⋯OH₂. The origin of the error is the low-energy unoccupied orbital, which lies close to the highest occupied orbital energy. Then, the unpaired electron can be in different orbitals with different spatial distributions and energies. Accordingly, it is an energy-driven error with the even charge-split between the two fragments due to symmetry. Using these orbitals, thus, the system alleviates the artificial self-interaction error, which is the appearance of a delocalization error. The error is always to be expected for situations such as the presence of a two-center three-electron bond where a non-integer number of electrons are involved. Even the DFT-B3LYP hybrid functional fails to describe appropriately these two minima. This is the typical indicator of systems with significant delocalization error. Thus, DFT-PBE and also hybrid functionals, such as B3LYP, overestimate the stability of the hemibonded configuration of the $(H2O)2+$ dimer, causing significant discrepancies between the DFT and CCSD(T) level of theory.

The trend between the exact exchange as an admixture and the calculated binding energy for the two low-lying structural isomers of the $(H2O)2+$ dimer identified here demonstrates that the percentage of exact exchange has a negligible influence on the hydrogen bond strength, while increasing the hemibond strength by 1 eV. Furthermore, the dispersion interaction (i.e., the interactions caused by fluctuations in the charge density) is much less prominent than the singly shared electron and has a small impact on this system. Therefore, the failure of most DFT functionals to describe hydrogen-bonded configuration is mainly due to the manifesting the delocalization error in cation systems, leading to large energy error and qualitatively incorrect structures.

Having said this, it is unsurprising that functionals with ∼50% exact exchange display a superior performance as indicated by the excellent agreement with CCSD(T). These functionals are also able to accurately predict transition state energies (Table IV). Identical trends were, indeed, also observed for dispersion interactions between simple organic molecules.⁶⁶ However, despite their good performance in predicting reaction energetics, these functionals still underestimate the binding energy of the water dimer by 0.1 eV and tend to significantly underestimate the bandgap of the water dimer by several eV. Indeed, functionals with higher percentage of exact exchange or range-separated hybrids with a full long-range exact exchange⁶⁷ are needed to predict this property correctly. The functional with a higher percentage of exact exchange, however, severely overestimates the stability of the hydrogen bonded structure and is, therefore, not suited for predicting reaction energies. On the other hand, the range-separated functionals are highly system dependent, and a single parameterization cannot be used across all chemistry groups. This clearly highlights that a self-interaction correction optimized for predicting binding and for optimizing bandgap predictions is not necessarily suited for the prediction of reaction energies.

Ab initio electronic structure methods have been employed to investigate the water dimer radical cation. Two minima have been located on the electronic ground-state potential energy surface: hydrogen-bonded and hemibonded configurations. At our highest level of theory, CCSD(T), the global minimum of the water dimer radical cation has been identified to be a hydrogen-bonded system, which lies 0.4 eV lower in energy than the hemibonded isomer. On the other hand, most DFT calculations with various functionals favor hemi-bonded configuration as a stable structure. We found that geometries and energetics are sensitive to basis sets and functionals, and nevertheless, the DFT results with 50% of exact exchange, such as M06-2X and BHandH functionals, are very close to the CCSD(T) level of theory. Thus, the energy comparison of the two minima serves as a criterion for finding the DFT functionals as an alternative to obtain reliable results. The energy barrier between these two minima has been determined to be 0.58 eV at the CCSD(T) level of theory. Therefore, the water dimer cation dissociation is likely to occur via the oxonium channel.

The supplementary material contains energy calculations vs ROO bond distances for both the hemibonded and hydrogen bonded configurations of (H2O)2+ utilizing various functional methods (M06L, B3PW91-D3, PBE0-D3, BHandH, and M06HF) and the CCSD(T) level of theory. Additionally, we provide the optimized geometries, expressed in Å, for the neutral water dimer and water dimer radical cations.

Financial support from the Dr. Barabara-Mez-Starck Stiftung Foundation and CELEST (Center for Electrochemical Energy Storage Ulm-Karlsruhe) by the German Research Foundation (DFG) under Project ID 390874152 (POLiS Cluster of Excellence) is gratefully acknowledged. Computer time provided by the state of Baden-Wuerttemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 40/575-1 FUGG (JUSTUS 2 cluster) is gratefully acknowledged.

The authors have no conflicts to disclose.

Michael Busch: Data curation (equal); Formal analysis (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Mohsen Sotoudeh: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

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