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

A major roadblock in realizing large-scale production of hydrogen via electrochemical water splitting is the cost and inefficiency of current catalysts for the oxygen evolution reaction (OER). Computational research has driven important developments in understanding and designing heterogeneous OER catalysts using linear scaling relationships derived from computed binding energies. Herein, we interrogate 17 of the most active molecular OER catalysts, based on different transition metals (Ru, Mn, Fe, Co, Ni, and Cu), and show they obey similar scaling relations to those established for heterogeneous systems. However, we find that the conventional OER descriptor underestimates the activity for very active OER complexes as the standard approach neglects a crucial one-electron oxidation that many molecular catalysts undergo prior to O–O bond formation. Importantly, this additional step allows certain molecular catalysts to circumvent the “overpotential wall”, leading to enhanced performance. With this knowledge, we establish fundamental principles for the design of ideal molecular OER catalysts.

where a catalytic Tafel plot is shown for a few complexes.
It is interesting to compare here the kinetic performance of complexes [2 IV (OH)] + and [10 IV (OH ax )] + , both containing the tda 2− ligand but with the Ru-hydroxido group located in the equatorial or axial positions, respectively. At pH = 7 the Ru(V/IV) redox couples appear at 1.43 V for 2 and 1.49 V for 10 and thus they represent overpotentials of 610 and 680 mV, respectively. However, even though 10 has a 60 mV higher thermodynamic driving force than 2 , the latter is more than 4 orders of magnitude faster (TOF max : 0.4 s −1 for 10 and 7,700 s −1 for 2 ). This enormous increase in the water oxidation kinetics is associated with the capacity of the dangling carboxylate in 2 that acts as an intramolecular base thus Figure 14. Catalytic Tafel plots for complexes 2 (dark blue), 8 (brown), 9 (light blue), 10 (black), and 11 (gray). The kinetic data of complexes 2 , 8 , 9 , and 10 correspond to pH 7 whereas complex 11 is at pH 1. Complexes 2 , 9 , and 10 oxidize water following a first order mechanism in respect to the concentration, and therefore their plots are concentration independent. In the case of complex 8 , its mechanism is second order in respect to the concentration, and therefore its Catalytic Tafel Plot is concentration dependent. A concentration of 10 mM is used for the Catalytic Tafel Plots of 8 .  Ru-7 50,000 10 1.7 × 10 -3 M CAN (n.a.) 31 S13

Supplementary Notes
Supplementary Note 1 The shape of the 3D volcano plot using the novel OER descriptor, ∆ #(%) * − ∆ #()%) * , depicted in Figure 5 and Supplementary Figure 2 can be explained by considering the equations that determine the boundaries separating distinct regions which denote the potential limiting step (PLS).
To begin the analysis, we define the Gibbs binding energies of the different OER intermediates within the computational hydrogen electrode model as follows: Where, DEF * is the value assigned to the normal hydrogen electrode potential in water, for which we used a previously calculated value of -4.28 V. The method applied to arrive at this value uses the absolute solvation free energy of the proton, which is calculated using the cluster pair approximation method, enthalpy and entropy corrections are then added by numerically solving equations from Fermi-Dirac statistics. This method was shown to give consistent and sufficient agreement with experimental data, as reported in Ref.
[55] of the main text.
We then use the above binding energies to determine the PLS using the established linear scaling relationship for molecular catalysts, shown in Figure  Step 5' is not seen on the volcano plot in Supplementary Figure 2 because it is not limiting until ∆ *# * is sufficiently negative. It should also be noted that there is an important distinction between the two volcano plots using the conventional and new OER descriptors shown in Figure 5. In particular, the overpotential can be directly obtained from the conventional volcano plot by computing the ∆ *# * and ∆ # * values -assuming that the first and last OER steps are not PLS, which is generally the case. On the contrary, the volcano plot using the novel OER descriptor requires three values, ∆ *# * , ∆ #(%) * and ∆ #()%) * (as can be seen above for ∆ GUL2 ).

S15
Supplementary Note 2 Notably, the catalyst we predict to exhibit the lowest theoretical overpotential is Ru-2, followed by Ru-1 and Ru-3. In experiments, this trend is respected also, with the Ru-1 and Ru-2 catalysts showing similar yet undoubtedly superior performance than Ru-3. The reason for the apparent difference between the experimental TOFs reported for Ru-3 and Ru-1-2 and our theoretical overpotentials may be ascribed to the fact that these Ru catalysts have been experimentally proposed to follow a I2M mechanism rather than the WNA assumed in our theoretical overpotentials. On the other hand, Ru-7 shows remarkable activity at high pH, as seen in Supplementary Table 4. However, to the best of our knowledge, this catalyst shows no activity in acidic media under which the other catalysts were tested. Furthermore, when Ru-3 is tested at pH = 7 (Supplementary Figure 5), it displays a better activity at low potentials than Ru-7, which is concurrent with our calculations predicting a lower overpotential for Ru-3. Hence, the predictive (both qualitatively and semi-quantitatively) strength of our approach is further demonstrated.

S16
Supplementary Note 3 In the review process of this paper, one referee suggested a more detailed explanation about the uncertainty in the overpotential wall. The source of this uncertainty is in the deviation from the line of best fit in Figure 2, and the uncertainty can be calculated as the standard deviation from the mean value of (∆ *## * − ∆ *# * ). In the following we provide a generalized explanation for the uncertainty of the overpotential wall to illustrate this is more detail.
Given an uncertainty , and a mean ∆ *## * − ∆ *# * value of 3.2 eV, it may be tempting to assume that this uncertainty is the same as the uncertainty associated with the overpotential wall.
However, we will now set out our interpretation to show that the interval of possible overpotential walls is in fact given by: W 3.2 − 2 − 1.23, 3.2 + 2 − 1.23Y As an instructive demonstration, in Supplementary Figure 6 we have constructed three putative volcano plots assuming different constraints imposed by the three putative values of (∆ *## * − ∆ *# * ), given by the OER scaling relations. The y-axis is the theoretical overpotential, which is indirectly calculated using DFT according to Eq. 7 in the main text, while the x-axis is the socalled OER descriptor, (∆ # * − ∆ *# * ). Because this latter energy is calculated by DFT, it has an associated uncertainty. Importantly, this uncertainty is the same for all the points; however, due to the constraint imposed by the (∆ *## * − ∆ *# * ) value obtained from the scaling, any shift away from the peak of a given volcano (located at [∆ *## * − ∆ *# * ]/2 on the x-axis) can only increase the theoretical overpotential. These uncertainties, represented by the red error bars in Supplementary Figure 6, are distinct from , which is the scale of the difference between the maximum and minimum heights of the overpotential wall (given by the expressions above), also S17 represented in Supplementary Figure 6. Hence, the uncertainty associated with the overpotential wall is simply given by [ /2].