Universal scaling relations for the rational design of molecular water oxidation catalysts with near-zero overpotential

1Reducing anthropogenic CO2 emissions is an urgent challenge facing civilization in the 21st century. Solar and wind energy are attractive carbon-free energy sources that could supply energy demand sustainably, but an important difficulty with reliance on these energy sources is their inherent intermittency and localized energy input. These two issues could be addressed by storing the excess energy in the form of chemical bonds, such as H2, formed via electrochemical water splitting in an electrolyzer. The hydrogen gas liberated at the cathode could then be stored, and eventually combined with oxygen in a fuel cell, providing an entirely renewable and clean energy supply with water as the only reaction product. Hydrogen produced via this method could also be utilized for the reduction of anthropogenic CO2 to produce chemical feedstocks and hydrocarbon fuels which are easier to transport than H2,.

2The bottleneck reaction in water splitting is the so-called oxygen evolution reaction (OER), occurring at the anode of the electrolyzer (Eq. (1)),. Catalysts which are cost-effective, highly active and robust for continued periods of time are yet to be found for this endergonic reaction, which demands large overpotentials (≥400 mV) to reach substantial current densities (≥10 mA/cm2) and oxygen conversion rates,.

3The OER process typically occurs via four elementary steps involving different reaction intermediates and the formation of an O–O bond, which is eventually released as molecular oxygen. The two primary pathways proposed for the O–O bond formation are the water nucleophilic attack (WNA) and the interaction of two metal-oxo entities (I2M), as depicted in Fig. 1. Both reaction paths begin with two proton electron transfer (PET) events, giving rise to a metal-oxo intermediate. In the case of WNA, the metal-oxo species subsequently undergoes the nucleophilic attack of a solvent water molecule and a further PET to generate the O–O bond, whereas in the I2M mechanism the O–O bond formation involves the coupling of two separate metal-oxo moieties.

4With this understanding of the OER, Nørskov et al. showed that the activity of metal and metal oxide surfaces via a WNA mechanism can be readily assessed by computing the binding energies of the different reaction intermediates,. Particularly, the Gibbs energy change associated with each elementary step can be calculated as in Eqs. (2)–(5), where denotes a metal active site, aH+ is the activity of the proton, e is the number of electrons involved, and U is the applied potential:

5By applying the computational hydrogen electrode model to the above expressions, the chemical potential of the H+/e– pair can be replaced by one-half of the chemical potential of molecular hydrogen, as these two species, by definition, are in equilibrium at 0 V versus the reversible hydrogen electrode (RHE). Furthermore, the introduction of the error arising from the modeling of the O2 molecule with density functional theory (DFT) methods can be evaded by fixing the overall Gibbs reaction energy,, to the experimental value of 4.92 eV.

6From the relative Gibbs energies in Eqs. (2)–(5), one can then determine the potential limiting step (PLS), defined as the most thermodynamically demanding step:

7And the theoretical overpotential by subtracting the thermodynamic redox potential of water:

8It is important to note that ηtheor does not lead to a direct prediction of the current density attained at that given potential, but rather an upper bound of the voltage at which OER activity is to be initiated. Hence, an ideal catalyst with ηtheor = 0 V would perfectly distribute the overall change in Gibbs energy throughout the four OER steps, i.e. ΔG1−4 = 1.23 eV.

9Further insight into the energetics of the OER intermediates was provided by Koper, showing that the HOO and HO∗ binding energies ( and ) on metal surfaces scale linearly with a constant energy difference of ca. 3.2 eV. From this OER scaling relation, the following important implications emerged. Firstly, the value of can be directly obtained from the calculated, and vice versa, which reduces the computational time needed to assess the OER activity of a given surface catalyst. Secondly, because the energy difference between the HOO and HO intermediates is larger than the expected value for an ideal catalyst (2 × 1.23 = 2.46 eV), the best possible material subject to this constraint would have a value exactly halfway between the and binding energies, resulting in a minimum theoretical overpotential of 370 mV ([3.2–2.46 eV]/2). Consequently, OER electrocatalysts are predicted to be limited by an “overpotential wall” of 370 mV, in excellent agreement with the experimental benchmark studies of state-of-the-art electrocatalysts. While the uncertainty associated with the calculation of can reduce that overpotential wall (see Supplementary Fig. 6 and Supplementary Note 3 for a detailed discussion), as has been experimentally observed for a few exceptional OER electrocatalysts–, hurdling or reducing this wall remains a major challenge.

10Another major implication of the above OER scaling is that since ΔG1 and ΔG4 are rarely PLS, the Gibbs energy difference between the O and HO intermediates (ΔG − ΔG) can be employed as a reaction descriptor to predict the OER activity of heterogeneous catalysts,. In fact, this descriptor has been successfully applied to rationalize the activity of a wide variety of metal oxide and single-atom electrocatalysts, wherein the most active ones exhibit a descriptor which is close to the optimal value predicted by the scaling relations, i.e. 3.2/2 = 1.6 eV,,. This descriptor also allows for the fast, high-throughput screening of materials from ab initio simulations, a procedure that has proven to be effective in the design of more efficient electrocatalysts for the other relevant reactions such as the hydrogen evolution reaction and CO2 reduction,.

11Recent theoretical studies have reported OER scaling relations for a few model molecular systems bearing a corrole, porphyrin ligands, and functionalized graphitic materials,,, but the generalization of these relations to catalysts featuring different metals and ligand scaffolds, as well as the applicability of the OER descriptor, is yet to be demonstrated. Herein we confirm the existence of universal scaling relations for a wide variety of well-established molecular OER systems. This investigation is prompted by the superior performance of certain homogeneous catalysts for which turnover frequency values (TOFs) of 2–3 orders of magnitude higher than the best heterogeneous systems have been reported–. To demonstrate the robustness of the established scaling relations, we have exhaustively interrogated 17 selected molecular OER catalysts featuring a wide variety of metals and ligands (Fig. 2)—all of them experimentally tested and reported in the literature—some of which are amongst the most active OER systems known to date. For catalysts with two potential active sites, only one of them was considered for simplicity. In selecting catalysts, we considered only those exhibiting a well-defined and well-established structure. For instance, Ir catalysts bearing pentamethylcyclopentadienyl (Cp) ligands, were discarded, despite some of them showing considerable OER activity, as the integrity of many of these complexes has been shown to be compromised under OER conditions,. More complex and unique systems with multimetallic centers were also discarded due to the multiple distinct valence states and active sites which would have to be considered prior to oxygen evolution.

12Interestingly, we find that the best molecular OER catalysts can circumvent the scaling between the HOO and HO intermediates by undergoing an additional one-electron oxidation prior to O–O bond formation. As a result, we predict molecular catalysts enabling this extra oxidation to potentially behave as ideal OER catalysts, exhibiting zero theoretical overpotential. Further, we demonstrate that the activity of such highly active systems is underestimated by the conventional OER descriptor and propose a new descriptor to screen and narrow down the search of promising molecular catalysts. Finally, we use all this knowledge to establish fundamental principles for the rational design of ideal OER catalysts to advance the development of commercial water electrolyzers.

13Firstly, we set out to prove the scaling relation between the HO and HOO intermediates observed for heterogeneous systems using the molecular OER catalysts depicted in Fig. 2. With this aim, we employed DFT calculations at the TPSSh level along with Grimme D3 dispersion corrections to obtain the Gibbs energies of the different OER intermediates, following the methodology described in the “Computational Methods” section and applying the computational hydrogen electrode model, previously used in studies of heterogeneous oxide catalysts for the OER.

14As shown in Fig. 3, the representation of the binding energies of the HOO species (ΔG) against those obtained for the HO species (ΔG) shows that the energetics of these two intermediates are strongly related, displaying a linear correlation with a slope close to unity—as expected based on the bond order conservation principle—and an intercept of 3.26 eV. This confirms that the same linear scaling relation observed in metal and metal oxide systems applies to any molecular OER catalyst. Also noteworthy is the robustness of this relation, which is derived from complexes showcasing many distinct ligands and metal centers–,–. In addition, the constraint imposed by the OER scaling implies that both homogeneous and heterogeneous catalysts are subject to a minimum overpotential of ca. 370 mV. The OER scaling relation can therefore now be labeled as “universal”, since it can be applied to any OER material regardless of their nature.

15We then examined in detail the data presented in Fig. 3 to identify catalyst features which could improve the OER scaling, thus reducing the predicted overpotential. Particularly, we hypothesized that intramolecular H-bonds between the HOO group and the metal ligands might offer that opportunity given the innate rigidity of the HO group (Supplementary Fig. 1). This additional stabilization of the HOO intermediate would shift down ΔG with respect to ΔG, closing the energy gap between the two intermediates, and therefore, improving the OER scaling. To prove our hypothesis, we classified our data points as catalysts with or without a H-bond by measuring the minimum distance from the H atom in the HOO group to the nearest N/O-ligand using a number of different cut-offs, chosen such that each split contained sufficiently different data points. As we anticipated, the ability of the HOO intermediate to form intramolecular H-bonds decreases the value of ΔG and, consequently, reduces the energy difference ΔG − ΔG. This effect is conveyed in Supplementary Fig. 1, where two clear distinct scaling relations are obtained depending on whether the HOO intermediate exhibited a H-bond or not. Importantly, we note that all linear trends display slopes close to unity; however, those catalysts with an intramolecular H-bond present in the HOO intermediate show a value of ΔG − ΔG that is ca. 0.1 eV lower than those without the H-bond. While this effect might not be decisive—being within our reported mean absolute error (MAE) of the chosen functional—we conclude that H-bonding can be leveraged to improve the scaling, and therefore, it should be considered as an important feature in the design and fine-tuning of molecular OER catalysts. H-bonding can also be leveraged to boost the kinetics of the O to HOO step in the WNA mechanism, as recently shown by experiments, while theoretical investigations have highlighted the possibility that the second coordination sphere can modulate the binding energy of the crucial O intermediate.

16We note that molecular dynamics simulations including explicit solvation would allow for a more complete characterization of the significance of intramolecular H-bonds on the molecular OER scaling. This methodology has been used to explain the mediation of O–O bond formation by solvent structures for the RuV–oxo intermediate of Ru-7 elsewhere,. Applying this level of analysis to all of our investigated catalysts (and OER intermediates) is, however, extremely computationally demanding. To circumvent this, DFT calculations in concert with minima-hopping have been employed to study solvent effects on the OER intermediates adsorbed on rutile IrO2, showing a larger stabilization of the HOO intermediate on a O-covered surface. Such investigations, however, are outside the scope of this work, which is attempting to establishing and generalizing OER-scaling relations for molecular catalysts and to set the guidelines for the accelerated discovery of high-performance OER catalysts using methods that are computationally efficient and less time consuming.

17After demonstrating that molecular catalysts hold a scaling relationship between the HO and HOO intermediates, we tested whether the OER descriptor proposed for heterogeneous catalysts, i.e. ΔG−ΔG, could also be applied to rationalize the experimental performance of the investigated complexes. By plotting this descriptor against the computed theoretical overpotential using in Eqs. (2)–(7), a volcano representation is obtained, where the most active catalysts should lie at the top. This representation also allows us to categorize molecular catalysts based on the PLS. More specifically, catalysts lying on the left-hand side of the volcano display a small OER descriptor value, and therefore, the O∗ to HOO∗ transition becomes the PLS. On the other hand, catalysts found on the right leg of the volcano exhibit a high descriptor value, which dictates the OER overpotential.

18The volcano plot obtained for the 17 molecular catalysts is shown in Fig. 4a, where, surprisingly, none of them are predicted to have an outstanding OER activity. This contrasts starkly with the reported experimental data, which have proven some of the complexes to be among the best molecular OER catalysts. In the following, however, we demonstrate that this unsuspected behavior is caused by an oversimplification of the reaction mechanism, typically assumed in most theoretical studies of heterogeneous OER catalysts. Notably, some metals, such as Ru, Mn, or Fe, have been experimentally proven to undergo an additional one-electron transfer (ET) before O–O bond formation–. Bearing this in mind, we computed (whenever possible) the Gibbs energies for the ET from the M(IV)-oxo to M(V)-oxo species with M = Ru, Mn, Fe, hereafter referred to as O(IV) and O(V), respectively (see Supplementary Tables 2 and 3 for details). The Gibbs energy change associated with this step, and the subsequent O–O bond formation via a PET event, are given by

19Interestingly, upon accounting for the additional ET in the M(IV)-oxo intermediates, we observed that catalysts which undergo this step exhibit an improved theoretical overpotential, deviating significantly from the lines defined by the volcano in Fig. 4a. This divergence from the scaling stems from the fact that those catalysts are neither limited by the HO to O, nor the O to HOO steps. Instead, the additional ET from the M(IV)-oxo species becomes the PLS, which demands less energy leading to a lower predicted overpotential. Based on this observation, we proceeded to apply the Gibbs energy change associated with the additional ET step,, as an alternative OER descriptor on the x-axis. The result of adopting this new descriptor can be appreciated in the volcano plot depicted in Fig. 4b, where the relevant catalysts now lie on the right leg of the volcano and where the best Ru-based catalysts cluster together around the top. Notably, DFT calculations predict Ru-1–3, Ru-7, and Ru-9 to exhibit the lowest theoretical overpotential amongst the 17 complexes investigated in this work, which agrees with the highest experimental TOFs reported for these complexes (see Supplementary Fig. 5 and Supplementary Note 2). It is important to note that the active species for Ru-9 has come under scrutiny, but our calculations predict that this catalyst would be a relatively active OER catalyst. The trend in the experimental TOFs (see Supplementary Table 4) is also captured by the new descriptor, highlighting its predictive (both qualitatively and semi-quantitatively) power and its potential for the rapid screening of molecular OER catalysts. We also note that when comparing energy differences between the PET and one-electron oxidation steps, it is important to consider whether this difference is within the inherent DFT error (Supplementary Table 1), in which case neither of the two pathways can be entirely discarded.

20Another important observation from Fig. 4b is that Ru-1 and Ru-2 are predicted to display a theoretical overpotential below 370 mV, suggesting that molecular catalysts undergoing the O(IV) to O(V) transition might be able to hurdle the overpotential wall. In fact, such a catalyst with the perfect distribution of the energy levels would present a Gibbs energy of only 1.07 eV for each of the steps in Eqs. (3), (8), and (9)—obtained by dividing the energy difference between the HO and HOO intermediates imposed by the conventional scaling, i.e. 3.2 eV, between the three steps. Assuming that one step between the HO and HOO intermediates equals the thermodynamic potential of water (i.e. 1.23 eV), the remaining two steps would sum to 1.97 eV (3.2–1.23). This implies that the design of ideal OER catalysts with 0 mV overpotential is theoretically within reach, provided that the 1.97 eV is properly distributed across those two steps. With this enlightening insight, we next examined the optimal range of binding energies which would lead to a near-zero or 0 mV overpotential using the 3D volcano representations in Fig. 5. We note that while the oxidation of O(IV) to O(V) has been known to be of critical importance experimentally and theoretically for individual catalysts, this is the first computational rationalization of exceptional catalyst activity, using this intermediate, across a set of active catalysts.

21The difference between the 3D volcanos derived for the catalysts following the conventional 4-PET mechanism (Fig. 5a) and that with the additional ET step (Fig. 5b) is striking, namely, the area of high activity is approximately two times larger for the latter volcano, while the activity itself is clearly much more favorable (i.e. ηtheor values of 370 and 0 mV, respectively). The shape of the volcano featuring the new OER descriptor can be rationalized by considering the equations that determine the boundaries separating the distinct regions (see Supplementary Fig. 2 and Supplementary Note 1), which denote the PLS at each region, as shown in Fig. 5b. This new volcano provides also a comprehensive illustration of the observed superior performance of Ru-based OER catalysts compared with other metal-based catalysts. Furthermore, it sets up a general guideline for the future design of complexes as cost-effective and high-performance OER catalysts, given the propensity of some earth-abundant metals (e.g. Mn and Fe) to stabilize high-valent oxo intermediates–.

22Our investigations have thus far assumed a WNA mechanism. However, it is important to note that even if the I2M pathway—or indeed any other mechanism—was found to be more favorable, it would only result in a better-predicted activity. Hence, the scaling relations and reaction descriptors established in this work are perfectly valid and relevant, since they provide an upper limit of the OER activity for a given molecular catalyst, requiring only a few DFT calculations. Besides, when comparing the WNA and I2M mechanisms, one should bear in mind that the preference for one or the other pathway will be strongly dependent on the reaction conditions, something that can be overlooked. In particular, the O–O bond formation via the WNA mechanism involves a PET step, and therefore, the Gibbs energy associated with this process will be reduced as potential increases according to Eq. (4). On the contrary, the formation of the M–O–O–M dimer through the I2M mechanism is a chemical step, implying that the Gibbs energy with respect to the monomer will remain invariant with the applied bias. To illustrate this, we have analyzed this mechanism for Ru-1, 2, and 3 for which the corresponding Ru(V)-oxo species has been kinetically confirmed to dimerize leading to the complex Ru(IV)–OO–Ru(IV),. Furthermore, anchoring Ru-1 on indium tin oxide has also been shown to increase the observed overpotential, presumably due to a change in mechanism from I2M to the WNA. It should be noted that the optimized geometries of the complexes as single molecules or dimers is octahedral, despite the fact that the favorable activity of some Ru catalysts have been explained by their ability to dynamically access a seven-coordinated intermediate. The effect of the applied voltage on the thermodynamics of both WNA and I2M mechanisms for Ru-1 is summarized in Fig. 6, while the data for Ru-2 and 3 dimers is presented in Supplementary Table 3.

23As can be observed from Fig. 6, DFT calculations predict that the I2M mechanism will dominate at applied voltages up to 1.98 V, in good agreement with the reported experiments using cerium ammonium nitrate as the sacrificial chemical oxidant (E0 = 1.61 vs. NHE). At more oxidizing potentials, however, simulations predict the WNA pathway to become accessible. Hence, the I2M mechanism is expected to have an important contribution to the observed OER activity at low to moderate potentials (or mild chemical oxidants), whereas the WNA is expected to compete at strong oxidizing conditions. As a result, the I2M path should be considered for catalysts with low predicted overpotentials, especially if the O to HOO transition is the PLS. The results presented in Supplementary Table 3 indeed suggest that the degree of exergonicity of this dimerization step is the deciding factor in the improved kinetics of Ru-1 with respect to Ru-2 and 3. The importance of circumventing the overpotential wall through O–O-based methods, such as the I2M mechanism has been discussed previously.

24We emphasize that a complete quantum chemical characterization of the underlying mechanism for catalysts is to remain vital due to the need for increased accuracy in predicting the OER activity and in forming an in-depth understanding of this complex reaction. We also note that catalysts which evolve O2 through more intricate pathways could exhibit lower overpotentials than predicted by our analysis, although the extrapolation of this knowledge to the design of novel molecular catalysts is not straightforward. Instead, our work sets the foundations for the fast and efficient screening of molecular OER catalysts and posits that catalysts, which are close to the top of volcano as seen in Fig. 4b have the potential to evolve oxygen through both the WNA mechanism as well as through more complex routes. This approach, yet to be exploited in homogeneous catalysis, allows for the exploration of a substantial portion of the vast potential chemical space, which would be otherwise intractable, and hence, it is expected to lead to the discovery of novel catalysts with an improved performance in a reasonable time.

25With all this knowledge, in the following we establish and discuss a series of catalyst design principles to accelerate the discovery of novel molecular systems based on earth-abundant elements exhibiting an ideal OER activity. Firstly, we propose that any high-throughput search should focus on transition metal complexes which can stabilize high oxidation states and thereby undergo the additional ET step to hurdle the OER overpotential wall. Although a complete characterization of the activity of these complexes would require the modeling of many more intermediates, the approach presented herein allows for the rapidly accelerated screening of molecular catalysts using thermodynamic OER descriptors. Special attention, however, should be paid when predicting the activity of catalysts with a low conventional OER descriptor—those lying on the left leg of the volcano in Fig. 4a—as their activity can be underestimated if they undergo an extra ET step. To prevent this, we propose to use the Gibbs energy change for the additional ET as a new OER descriptor. In addition, the I2M mechanism must also be taken into account for those complexes so as to avoid overestimating the overpotential by implicitly assuming the oxygen evolution proceeds via a WNA from the O to HOO intermediates.

26Future design of highly active molecular OER catalysts may also benefit from the knowledge available in the heterogeneous catalysis community, where OER materials with reduced overpotentials have been successfully predicted by means of DFT calculations and scaling relations. An instructive example is the recent computer-aided design of a ternary Co, Fe, and W oxide material which exhibited a record low overpotential in alkaline electrolyte. Therein, the optimal binding of the HO intermediate of the ternary system was tailored by combining pure metal oxides and oxyhydroxides with weak and strong binding energies. By analogy, in molecular systems, we envisage that a similar effect could be achieved by altering the transition metal center to produce a favorable binding energy. The blending of features of two homogeneous catalysts is, however, a more challenging process. Particularly, we anticipate that replacing the metal center in molecular systems may lead to a more drastic effect on the electronic structure compared to doping of heterogeneous systems, where the effect is more distributed throughout the system. This may completely change the OER descriptor value and lead to a different PLS, illustrating the innate difficulty in tuning catalysts involving many ETs and intermediates.

27In designing more efficient OER complexes, the choice of ligand is crucial. For example, H-bonding with the metal ligands can be leveraged to improve (or break) the OER scaling between the HO and HOO intermediates, while ligands with different electronic properties may play a role in stabilizing other high-valent intermediates or by acting as a proton shuttle to improve the overall reaction kinetics, as has been experimentally demonstrated recently. Therefore, the electronic and steric effects of the ligand environment need to be systematically investigated to illuminate potential strategies for merging the desirable features of distinct ligands to identify highly active OER catalysts. Importantly, advances in this area would enable the future implementation of ideal molecular OER catalysts into commercial electrolyzers via their immobilization onto conductive electrode supports leading to hybrid materials with greatly enhanced selectivity, stability, and electron transport. In fact, this has been experimentally achieved for a wide variety of molecular systems– and theoretically proven to be feasible for model Ir catalysts with no apparent loss of activity, although substantial work remains to be done in order to devise more efficient and robust hybrid OER materials based on abundant elements.

28Finally, one has to realize that homogeneous catalysts are essentially infinitely tunable owing to the vastness of chemical space. Hence, highlighting the fruitful regions within this space is a most demanding and exciting challenge which we can begin to address with these design principles. With ample amounts of data, we also envision machine learning (ML) algorithms to become an integral part of such processes to greatly enhance the speed of catalyst discovery in an automatized approach, as illustrated in Fig. 7.

29Overall, this work provides strong theoretical evidence that OER scaling relations can be generalized to any material regardless of their nature. Additionally, it demonstrates that the conventional OER descriptor can be applied to rapidly screen molecular catalysts and to predict their activity requiring only a few DFT calculations. This descriptor, however, underestimates the catalytic performance of certain metal complexes with the O to HOO transition as PLS, since it does not consider the additional one-electron oxidation that they undergo before O–O bond formation. To capture this behavior, we propose using a new OER descriptor based on the Gibbs energy of that extra oxidation, and demonstrate that it provides a good account for the catalytic performance of the most active complexes considered herein, according to their reported experimental TOFs. This proves the predictive power (both qualitative and semiquantitative) of the new OER descriptor, which is unprecedented. In addition, this new descriptor pinpoints the additional one-electron oxidation as the means by which certain catalysts can circumvent the “overpotential wall”. In fact, we predict that molecular systems undergoing that extra step can exhibit zero theoretical overpotential provided the energy levels of the relevant OER intermediates are properly aligned. Finally, all the knowledge acquired in this work has been used to establish the fundamental principles for the rational design of ideal OER catalysts. We expect that these guidelines can be applied in a fully automatized approach driven by ML algorithms and OER descriptors to accelerate the discovery of ideal OER materials based on earth-abundant elements. Altogether, this work represents a significant conceptual advance towards developing energy conversion and storage technologies based on water splitting by bringing together concepts from homogeneous and heterogeneous catalyses. Therefore, we foresee this work will inspire future research directions in the field and spur the development of commercial water electrolyzers.

30DFT calculations were performed at the TPSSh level using the Gaussian09 software. The choice of this functional was based on the satisfactory results obtained in a preliminary benchmark study of the redox potentials of four molecular OER catalysts investigated in this work using three different functionals with varying levels of Hartree–Fock (HF) exchange, namely B3LYP (20% HF), TPSSh (10% HF),, and BP86 (0% HF),. The results and details of these calculations are presented in Supplementary Table 1. Overall, the TPSSh functional provided the lowest average MAE (MAE = 303 mV) compared with the reported experimental values (see Supplementary Table 1), and therefore, it was chosen for subsequent theoretical investigations. We note that an error of ca. 300 mV is within the accuracy reported in previous theoretical works. This choice of functional, in concert with the chosen basis sets is further supported by the satisfactory results obtained in the computation of heats of formation of many different transition metal complexes.

31To describe the Ru, Mn, Fe, Cu, Co, and Ni metals, the Lanl2dz effective core potential was used along with f-polarization functions with exponents 1.235, 2.195, 2.462, 3.525, 2.78, and 3.13, respectively,. The more electronegative O, N, F, and Cl atoms were described using the 6−31+G(d) basis set, while the 6–31G(d,p) basis set was employed for C and H atoms. To account for non-covalent interactions, Grimme-D3 dispersion corrections were added to the Gibbs energy through single point calculations using the optimized geometries in solution (see below). We note that the primary conclusions of this work remain the same regardless of the addition of dispersion; data for calculations without dispersion corrections are included in the Supplementary Figs. 3 and 4. For catalysts with two potential active sites, only one was considered for simplicity. This approach has been successfully applied to the modeling of the OER activity of polyoxometalates,, providing results in good agreement with experiments.

32Molecular structures were optimized in solution using the implicit SMD solvation model with H2O (ε = 78.3553) as a solvent. Whenever high-spin and low-spin states were accessible, both were considered in geometry optimizations and the lowest energy structure was selected.

33Gibbs energies were calculated at the temperature of 298.15 K and pressure of 1 atm, except for the isolated H2O molecule that was computed at the temperature and pressure at which both the liquid and gas phases are in equilibrium, i.e. 300 K and 0.035 atm. Relative Gibbs energies are referenced to H2O and H2 in solution, to avoid introducing the error associated with the modeling of O2 with DFT methods, and the global reaction Gibbs energy was fixed to the experimental value of 4.92 eV. All the optimized structures were verified to be real energy minima by means of vibrational frequency analysis.