OER scaling relations for state-of-the-art molecular OER catalysts

Firstly, we set out to prove the scaling relation between the HO\({}^{\ast}\) and HOO\({}^{\ast}\) 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 OER22.

As shown in Fig. 3, the representation of the binding energies of the HOO\({}^{\ast}\) species (ΔG\({}_{{\mathrm{HOO}}^{\ast}}\)) against those obtained for the HO\({}^{\ast}\) species (ΔG\({}_{{\mathrm{HO}}^{\ast}}\)) 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 principle36—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 centers29,30,31,37,38,39,40,41,42,43,44,45,46,47. 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.

Fig. 3 Linear scaling relation between the HO∗ and HOO∗ intermediates for the molecular OER catalysts investigated in this work and depicted in Fig. 2. The shaded blue region represents a 99% confidence interval for the linear model Full size image

We 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\({}^{\ast}\) group and the metal ligands might offer that opportunity given the innate rigidity of the HO\({}^{\ast}\) group (Supplementary Fig. 1). This additional stabilization of the HOO\({}^{\ast}\) intermediate would shift down ΔG\({}_{{\mathrm{HOO}}^{\ast}}\) with respect to ΔG\({}_{{\mathrm{HO}}^{\ast}}\), 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\({}^{\ast}\) 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\({}^{\ast}\) intermediate to form intramolecular H-bonds decreases the value of ΔG\({}_{{\mathrm{HOO}}^{\ast}}\) and, consequently, reduces the energy difference ΔG\({}_{{\mathrm{HOO}}^{\ast}}\) − ΔG\({}_{{\mathrm{HO}}^{\ast}}\). This effect is conveyed in Supplementary Fig. 1, where two clear distinct scaling relations are obtained depending on whether the HOO\({}^{\ast}\) 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\({}^{\ast}\) intermediate show a value of ΔG\({}_{{\mathrm{HOO}}^{\ast}}\) − ΔG\({}_{{\mathrm{HO}}^{\ast}}\) 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\({}^{\ast}\) to HOO\({}^{\ast}\) step in the WNA mechanism, as recently shown by experiments31, while theoretical investigations have highlighted the possibility that the second coordination sphere can modulate the binding energy of the crucial O\({}^{\ast}\) intermediate48.

We 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 elsewhere49,50. 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 IrO 2 , showing a larger stabilization of the HOO\({}^{\ast}\) intermediate on a O-covered surface51. 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.

OER volcano plots and reaction descriptors

After demonstrating that molecular catalysts hold a scaling relationship between the HO\({}^{\ast}\) and HOO\({}^{\ast}\) intermediates, we tested whether the OER descriptor proposed for heterogeneous catalysts, i.e. ΔG\({}_{{\mathrm{O}}^{\ast}}\)−ΔG\({}_{{\mathrm{HO}}^{\ast}}\), 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.

The 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 formation52,53,54,55. 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)\({}^{\ast}\) and O(V)\({}^{\ast}\), 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

$$\begin{array}{l}{\mathrm{O}}({\mathrm{IV}})^ \ast \ \rightleftharpoons \ {\mathrm{O}}({\mathrm{V}})^ \ast + {\mathrm{e}}^ - \\ \Delta G_{3^{\prime}} = \Delta G_{O({\mathrm{V}})^ \ast } - \Delta G_{O\left( {{\mathrm{IV}}} \right)^ \ast } - {\mathrm{eU}}\end{array}$$ (8)

$$\begin{array}{l}{\mathrm{O}}({\mathrm{V}})^ \ast + {\mathrm{H}}_2{\mathrm{O}} \ \rightleftharpoons \ {\mathrm{HOO}}^ \ast + {\mathrm{H}}^ + \\ \Delta G_{4^{\prime}} = \Delta G_{{\mathrm{HOO}}^ \ast } - \Delta G_{{\mathrm{O}}\left( {\mathrm{V}} \right)^ \ast } + k_{\mathrm{B}}T\,{\mathrm{ln}}\, a_{{\mathrm{H}}^ + }\end{array}$$ (9)

Fig. 4 2D-volcano representations featuring different OER descriptors. a Volcano plot using the conventional OER descriptor. Full markers denote catalysts which are theoretically predicted to undergo an additional ET step before O–O formation (see Supplementary Table 3). The red shaded area indicates the optimal range for the conventional OER descriptor. b Volcano plot including only the molecular catalysts predicted to undergo an additional ET. In this volcano, the calculated Gibbs energy for the additional ET from the M(IV)-oxo to M(V)-oxo intermediate is represented in the x-axis as a new OER descriptor. Note, the same scale and the left leg of the volcano in a has been used for comparability. The dotted line in b is drawn to guide the eye towards an ideal OER catalyst. The potential limiting step on each side of the two volcano plots is indicated in red text Full size image

Interestingly, 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\({}^{\ast}\) to O\({}^{\ast}\), nor the O\({}^{\ast}\) to HOO\({}^{\ast}\) 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, \(\Delta {G}_{{\mathrm{O}}({\mathrm{V}}) \ast } - \Delta {G}_{{\mathrm{O}}({\mathrm{IV}}) \ast }\), 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 scrutiny56, 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.

Another 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)\({}^{\ast}\) to O(V)\({}^{\ast}\) 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\({}^{\ast}\) and HOO\({}^{\ast}\) intermediates imposed by the conventional scaling, i.e. 3.2 eV, between the three steps. Assuming that one step between the HO\({}^{\ast}\) and HOO\({}^{\ast}\) 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)\({}^{\ast}\) to O(V)\({}^{\ast}\) has been known to be of critical importance experimentally and theoretically for individual catalysts57, this is the first computational rationalization of exceptional catalyst activity, using this intermediate, across a set of active catalysts.

Fig. 5 Three-dimensional volcano representations including. a OER catalysts following the conventional 4-PET pathway and b catalysts theoretically predicted to undergo the additional ET step before O–O formation. The PLS on each region of the volcano plots is indicated Full size image

The 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 intermediates58,59,60.

Comparing WNA and I2M reaction mechanisms under OER conditions

Our 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)29,30. 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 WNA61. 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 intermediate62. 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.

Fig. 6 Gibbs energy diagram for Ru-1 calculated at 0 V (solid lines) and 1.98 V (dashed lines) vs. RHE. Blue lines represent the I2M mechanism while the pale orange lines represent the conventional WNA mechanism Full size image

As 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\({}^{\ast}\) to HOO\({}^{\ast}\) 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 previously15.

We 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 O 2 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.

Fundamental principles for the design of ideal molecular OER catalysts

With 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\({}^{\ast}\) to HOO\({}^{\ast}\) intermediates.

Future 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 electrolyte18. Therein, the optimal binding of the HO\({}^{\ast}\) 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.

In 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\({}^{\ast}\) and HOO\({}^{\ast}\) 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 recently31. 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 systems63,64,65,66 and theoretically proven to be feasible for model Ir catalysts with no apparent loss of activity67, although substantial work remains to be done in order to devise more efficient and robust hybrid OER materials based on abundant elements.

Finally, one has to realize that homogeneous catalysts are essentially infinitely tunable owing to the vastness of chemical space68. 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.

Fig. 7 Proposed approach for the accelerated discovery of ideal OER catalysts Full size image

Overall, 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\({}^{\ast}\) to HOO\({}^{\ast}\) 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.