Efficient workflow for the investigation of the catalytic cycle of water oxidation catalysts: Combining GFN‐xTB and density functional theory

Abstract Photocatalytic water oxidation remains the bottleneck in many artificial photosynthesis devices. The efficiency of this challenging process is inherently linked to the thermodynamic and electronic properties of the chromophore and the water oxidation catalyst (WOC). Computational investigations can facilitate the search for favorable chromophore‐catalyst combinations. However, this remains a demanding task due to the requirements on the computational method that should be able to correctly describe different spin and oxidation states of the transition metal, the influence of solvation and the different rates of the charge transfer and water oxidation processes. To determine a suitable method with favorable cost/accuracy ratios, the full catalytic cycle of a molecular ruthenium based WOC is investigated using different computational methods, including density functional theory (DFT) with different functionals (GGA, Hybrid, Double Hybrid) as well as the semi‐empirical tight binding approach GFN‐xTB. A workflow with low computational cost is proposed that combines GFN‐xTB and DFT and provides reliable results. GFN‐xTB geometries and frequencies combined with single‐point DFT energies give free energy changes along the catalytic cycle that closely follow the full DFT results and show satisfactory agreement with experiment, while significantly decreasing the computational cost. This workflow allows for cost efficient determination of energetic, thermodynamic and dynamic properties of WOCs.


| INTRODUCTION
Chemical fuels produced by solar energy have shown potential as a clean energy alternative to carbon-based fossil fuels. Fuel production most commonly involves proton or CO 2 reduction. The electrons needed for this reduction can be obtained through the oxidation of water into protons, molecular oxygen, and electrons. Therefore, water splitting dye-sensitized photoelectrochemical cells and other photoelectrochemical devices for solar energy production have been intensely investigated in the recent decades. [1][2][3] Still, the efficiency of such devices remains quite low as the water oxidation usually requires a high overpotential leading to significant energy losses. Although high overpotential can in principle be lowered by molecular catalysts, their operation with a sufficiently high turnover number and turnover frequency is a challenge. The development of a combination of a photosensitizer and a (transition metal-based) stable and rapid catalyst in a photocatalytic complex is thus a crucial step forward in the realization of more efficient devices. Ru-based transition metal complexes have in this respect in the last decade emerged as promising water oxidation catalyst (WOC) candidates. [4][5][6][7][8] In addition to the catalyst itself, also the coupling of the WOC to a suitable photooxidative dye is challenging due to the different requirements for oxidation potentials and HOMO energies for the different catalytic steps that involve a number of oxidation states of the transition metal. Electronic configurations change through the catalytic cycle, which affects the nature of the HOMO and its energy. Computational investigations with for example, estimation of redox potentials, orbital energies, and excitation energies can help in finding suitable combinations of dye and WOC candidates. 9 Dynamics and solvent effects should also be taken into account to simulate both water oxidation as well as the initial photooxidation of a dye-WOC complex to obtain reliable insight in electron and hole transfer processes. 10 DFT-based ab initio molecular dynamics coupled with enhanced sampling techniques have been employed quite successfully in determining reaction barriers for water oxidation processes with oxidized WOCs, [11][12][13][14] redox mediators 15,16 or photo-oxidized dyes. [17][18][19][20] However, these simulations are computationally very demanding as they require a quantum mechanical description of the entire system, including explicit solvation, since the water solvent participates actively in the catalytic reaction.
Methods that well describe the organic molecular dyes, transition metal complexes and water are therefore needed, but this is complicated due to the open shell nature of the systems and the different spin states that need to be taken into account. Finding a computationally affordable and yet reliable method to determine thermodynamic requirements, oxidation potentials, and dynamic properties remains challenging.
A semiempirical method that has shown great potential in the structural characterization of transition metal complexes is the so called GFN-xTB (Geometries, Frequencies, Noncovalent interactions extended Tight Binding) method developed by Grimme et al. 21 This method is based on atomic and global parameters that are optimized on basis of DFT. GFN-xTB has been used rather successfully for the description of transition metal complexes [21][22][23] and lanthanoides. 24 This method has been also extended on several fronts, including the so called GFN-2xTB, that is founded on more physically relevant parameters and the GFN-FF, 25,26 a force field parametrized on GFN-xTB results. Both GFN-xTB and GFN-2xTB have already been used in the determination of redox potentials with low computational cost. [27][28][29] In this study, we determine whether GFN-xTB is a viable alternative to describe a Ru-based catalyst developed by Duan et al., 4 also with respect to potential molecular dynamics simulations. Therefore, we investigated the performance of GFN-xTB on this Ru-based WOC in comparison to DFT with different exchange correlation functionals and to available experimental data. We propose a computationally efficient workflow that combines GFN-xTB calculations for geometries and frequencies and B3LYP for energies, which leads to accurate relative Gibbs free energies along the catalytic cycle.  30 The other possible pathway involves two 2 2+ dimer through a radical coupling mechanism: this mechanism has been called interaction of two metal oxo species (I2M) or radical oxo coupling (ROC) in the literature. 32,33 This binuclear reaction pathway often shows higher turnover frequencies and lower overpotentials than the WNA pathway. It also circumvents the problem of scaling relations between the Ru-OH and Ru-OOH bond strength that makes optimizing catalysts far from trivial due to the interdependence of these two parameters. 34,35 However, the reaction rate in the I2M mechanism depends on the catalyst concentration, since the coupling involves two catalysts in the correct oxidation state, and can be accelerated by increasing the local concentration, for example, through accumulation in self-assembling nanospheres. 36 The investigated catalyst has been shown experimentally to perform via the I2M mechanism. 4 However, by exchanging both equatorial or axial ligands, the catalyst can be tailored and the I2M mechanism can be made either more dominant or less favorable, to the point where the WNA mechanism is enforced over the I2M route. [37][38][39][40] In general, complexes with excess spin density on the oxygen and attractive interactions between axial ligands of two complexes lead to I2M mechanisms, whereas more positive partial charge on the same oxygen (and thus a higher electrophilicity) and steric hindrance between two complexes leads to a preferred WNA. 32,38,39 Since both reaction mechanisms are possible for the investigated complex in different regimes, such as low concentration, immobilizing the catalyst on a surface and so forth, it is crucial that the computational methods used to investigate the catalyst describe both mechanisms reasonably well and predict the right mechanism. We therefore also evaluate the different methods with respect to predicting the correct catalytic pathway.   48,49 and OPBE. [48][49][50] All DFT-based simulations were performed with the Slater type TZP (triple-ζ polarized) basis set. 51 D3 dispersion corrections with BJdamping were used. 52 Scalar Relativistic effects were included via the Zero-Order Regular Approximation (ZORA). [53][54][55] The solvent environment was included through the COSMO implicit-water model. 56 Geometry optimizations were performed with unrestricted DFT, considering all possible spin states involving the d-orbitals. Vibrational analysis was done only for the geometries corresponding to the most stable spin state.

| GFN-xTB-based simulations
The first oxidation potential was determined by ΔSCF, 61,62 from the difference between the energy of the precatalyst 1 [Ru(II)-OH 2 ] and its oxidized form 2 [Ru(III)-OH 2 ] + at the same geometry. Since this oxidation is pH independent and does only involve an electron transfer, the relaxation of the environment can be neglected as it is slow in comparison to the fast electron transfer. We also note that intake of water for the Ru(II) oxidation state is quite challenging and could in our model only be achieved by breaking a coordination bond between one carboxylic acid group and the ruthenium, with formation of a hydrogen bond between water and this carboxylate. A representative geometry is shown in Figure S1. This finding supports the suggestion that the complex is not stable in the 7-coordinated state at this low oxidation state and water coordinates to the complex under breaking of one of the other coordination bonds to form a 6-coordinated complex. 6 Since the reaction involves the extraction of an electron without a coupled proton transfer, the energy difference is essentially determined against the absolute electrode potential (as in comparison to vacuum). Therefore, the energy needs to be converted to the nor- The other two reaction steps involve a proton coupled electron transfer. Here, the reorganization of the complex is necessary for the oxidation to take place, especially the coupled proton transfer. Therefore, Gibbs free energy differences between the reactants and products are used to obtain the oxidation potentials. As mentioned earlier, the proton and electron release is assumed to be equivalent to the release of half a hydrogen molecule. These potentials are therefore already versus NHE at pH 0 since this corresponds to the energy needed to reduce a proton at standard conditions. As the experimental values were determined at pH = 1, the computed oxidation potentials were adjusted to this pH value by shifting the H + potential using the Nernst equation as shown in Equation 5.
here, G 0 H þ ,pH¼0 is the standard Gibbs free energy at standard conditions with pH = 0, which is 0 eV by definition in NHE, k B the Boltzmann constant (~8.617 eV/K) and T is the temperature at standard conditions (T = 298,15 K). Following the Nernst equation for a single electron event, the computationally determined oxidation potentials involving a PCET step were shifted by À0.059 V to obtain the values versus NHE at pH = 1.

| GFN-xTB + DFT approach for free-energy calculations
For a good compromise between accuracy and computational cost, a combination of DFT and GFN-xTB calculations was performed to estimate the Gibbs free energy. Since the geometries obtained at the GFN-xTB and B3LYP level are remarkably similar (see Figure S2 and The procedure is also visualized in Scheme 3. This workflow was tested using B3LYP, OPBE and the double hybrid functionals rev-DOD-PBE, rev-DOD-PBEP86, and rev-DOD-BLYP. [63][64][65] The settings for the Double Hybrid calculations were the same as for the other DFT based simulations (thus including COSMO water, D3 [BJ] corrections, relativistic effects via ZORA), except for the higher basis set TZ2P. [51][52][53][54][55][56] We also note that in a very recent work by Spicher et al., the Single Point Hessian (SPH) method has been introduced, which combines GFN-xTB/GFN-FF methods with DFT to obtain reliable free energies in non-equilibrium geometries, for example, obtained with a different level of theory or from molecular dynamics snapshots. 66 3 | RESULTS

| Energetically preferred spin states
For all catalytic intermediates, the energetically favored state turned out to be the one with the lowest possible spin multiplicity. This is due to the 7-coordinated environment of the Ru, breaking the octahedral symmetry and thus removing the t2g orbitals degeneracy.
This result is consistently found using both GFN-xTB and DFT with all exchange correlation functionals considered in this work (see Tables S1-S5). The corresponding favored spin state is reported in Scheme 1 for all intermediates. It should also be noted that for all cases, the states of higher multiplicity lie significantly higher in energy than the ground-state (see Section S1) and are therefore neglected in further investigations. We note however, that an intersystem crossing (ISC) event is necessary for the release of oxygen from the dimer. Formation of the dimer is a radical coupling mechanism that necessitates two antiferromagnetically coupled 2 [Ru(V) = O] + species, since the formation of the new covalent bond is due to radical coupling of these two unpaired electrons. If the spins are parallel, thus ferromagnetically coupled, the bond formation has much less of a radical coupling character and thus requires a higher activation energy. This bond formation was studied computationally by Nyhlén et al., who also found a higher barrier for the ferromagnetically coupled 2 [Ru(V) = O] + pair. 44 For the formation of 3 O 2 , an ISC from singlet to triplet state is necessary. The authors showed that the dimer in the triplet state has a very low barrier toward oxygen release. 44 In the WNA mechanism, the last step also necessitates an ISC for 3 O 2 release.

| Relative Gibbs free energy along the catalytic cycle
The Gibbs free energy differences of all catalytic intermediates at pH = 0 with respect to the first intermediate The conclusion from these results is that GFN-xTB in the present form (not allowing for unrestricted calculations) appears to be unsuitable for correctly describing the energetic and thermodynamic properties of the different catalytic intermediates. We note that this is not surprising, since GFN-xTB has been developed with the goal of computing accurate geometries and frequencies, but not primarily for accurate energetics.

| Preferred reaction mechanism
Experimental investigations on the catalyst showed that it operates via the I2M radical coupling mechanism. 4 All investigated methods predict this correctly. As visible in Figure 1 predicting the I2M mechanism to be much more likely. All methods show the oxygen release from the dimer to be endergonic under standard conditions. However, we note that this endergonic behavior might be an artifact of the missing explicit water solvation, as the water uptake is highly sensitive to the hydrogen bonding network. 44 Furthermore, an ISC event is necessary before triplet oxygen can be released.
The triplet state of the dimer, as seen in the SI, is considerably higher in energy than the singlet for all investigated methods, making the oxygen release favorable and irreversible. While for the DFT based methods, the Gibbs free energy difference between the dimer and the final catalytic step to regenerate the 2 [Ru(III)-OH 2 ] + under the release of oxygen is relatively small, it is quite large for GFN-xTB. Although this does not change the preference for the I2M mechanism, it underlines again that GFN-xTB alone is not suited to describe the energetics in the catalytic cycle in a satisfactory manner.

S C H E M E 3
Workflow for the GFN-xTB + DFT combined method. The colors denote the different computational methods used: Blue represents GFN-xTB, red DFT. Purple denotes the combination of both. After a GFN-xTB based geometry optimization, vibrational analysis with GFN-xTB is performed, and the calculation is completed with a single point DFT at this geometry. The combination of these results gives an estimate of the Gibbs free energy (purple box)

| Relative Gibbs free energy using GFN-xTB + DFT
While GFN-xTB did not prove to be reliable for the energies, it provides geometries and frequencies that are quite similar to the B3LYP results (see root mean square displacements and the Internal Energy and entropic terms in Table S6). For this reason, a combination of GFN-xTB geometries and frequencies with single point energies from B3LYP were also tested (see Scheme 3). Using a geometry optimization and vibrational analysis from GFN-xTB in combination with a single point energy from B3LYP significantly reduced computational cost while still providing reliable results for the free energy. The relative Gibbs free energy of the catalytic intermediates for both the WNA and the I2M mechanism is given in Figure 2 and Table 1 for both B3LYP and GFN-xTB + B3LYP. As can be seen in Figure 2, the GFN-xTB + B3LYP method predicts the correct reaction mechanism and follows very closely the B3LYP results: the deviation from B3LYP is quite small, around 0.1 eV for all the steps except one deviating by 0.26 eV. While the energy differences between consecutive steps are in excellent agreement, the accumulation of small errors lead to larger deviations for the later intermediates, as they are taken relative to the first step. The relatively large structural difference for the 2 [Ru(III)- Table S6) between the two methods might explain the slightly larger deviation at this step.
While this combination of B3LYP and GFN-xTB gives almost quantitatively the same results as the B3LYP, the computational cost is drastically reduced as only one single point calculation with the hybrid functional is needed, thus avoiding the demanding geometry optimization, and more importantly, frequency calculation. As the computational cost for GFN-xTB is extremely small compared to B3LYP (see Table S6), this leads to a speed up of a minimum of 2 orders of magnitude for this system, as 2*3 N (>300) points need to F I G U R E 1 Gibbs free energies at pH = 0 relative to the first catalytic step along the catalytic pathways for (A) the WNA mechanism and (B) the I2M mechanism obtained with: B3LYP (red, solid line), BLYP (blue, dashed), PBE (green, dashed), OPBE (cyan, dashed), and GFN-xTB (purple, dotted) T A B L E 1 Gibbs free energy difference in eV at pH = 0 for all catalytic steps and tested methods relative to the first catalytic intermediate All in all, this combination of DFT and GFN-xTB shows great potential in obtaining results close to the DFT description while being much more efficient computationally.

| Comparison to experimental oxidation potentials
The experimentally determined oxidation potential of the 1 [Ru(II)-  Table 2 and compared to the calculated oxidation potentials versus NHE at pH = 1 for all investigated methods. The results are also shown graphically in Figure 3.
In addition to the results reported in Table 1  only one of the steps shows such a large deviation, the oxidation potential was also determined using B3LYP geometries and frequencies combined with rev-DOD-PBE single point energies.
Here, the deviation remained over 0.7 V as well (see Table S7).
The poor performance of the double hybrids, in this case, is most likely due to the fact that they are used on geometries not opti-