High-throughput screening and rational design to drive discovery in molecular water oxidation catalysis

The difficulty of the oxygen evolution reaction is a fundamental impediment to the sustainable production of hydrogen, wherein molecular catalysts show the most impressive activity in terms of turnover frequency for this reaction. In this work, we have interrogated 444 automatically generated molecular water oxidation catalysts composed of well-known ligand scaffolds and six different transition metals (Cr, Mn, Fe, Ru, Co, and Ni). This data confirms the method-independent universal scaling relationship for water oxidation catalysts, describes routes toward circumventing this relationship, and justifies the ascendency of Ru catalysts for this reaction. Leveraging this information while applying catalyst design principles, we summarize experimental results, giving credence to our prediction of 9 earth-abundant molecular catalysts with theoretical overpotentials ranging from 200 to 400 mV as promising leads for experimental investigation. We also establish insights into spin-dependent scaling relations for key OER intermediates. Altogether, this work outlines the first steps towards enabling inverse design for molecular OER catalysts.


Introduction
The transition toward sustainable energy would be greatly facilitated by developing feasible means to store intermittent energy. 1 One promising strategy is the production of hydrogen from renewably-driven water electrolysis, with the ability to use H2 as an energy-dense fuel to decarbonize sectors of industry and transportation. Aside from the merits of using it as a fuel, hydrogen is essential to fertilizer production, and finds industrial applications from oil refineries to electronics. 2 However, this hydrogen is predominantly made via steam reforming of methane, contributing significantly to CO2 emissions. 3 To date, the anodic oxygen evolution reaction (OER) is the principal challenge preventing electrochemical hydrogen production at scale. Cheap, efficient and performant catalysts for the OER would impel green H2 production and improve CO2 electroreduction, which also occurs at the cathode concurrently with the anodic reaction shown in Eq. 1. 4,5 To evolve oxygen, a variety of mechanisms are possible, including electron and proton transfer steps, which can be sequential or proton-coupled electron transfers (PCETs). The most commonly studied mechanisms are the water nucleophilic attack (WNA) mechanism and the interaction of two metal-oxo units (I2M), as depicted in Figure 1. Computational studies assuming the WNA mechanism typically describe Gibbs energy changes associated with each step in the WNA as functions of intermediate binding energies e.g. ∆ *+ * , as outlined in  ! + * → * + ' + - Where denotes the applied potential, * " the activity of the protons, and the Gibbs energy of the overall reaction has been replaced by 4.92 eV to address the well-documented issues with modelling the triplet ground state of the O2 molecule. 7 These binding energies allow for the calculation of the theoretical overpotential, 23456 , according to Eq. 6.
Linear scaling relationships have been established between the binding energies of the HO* and HOO* intermediates using density functional theory (DFT) for heterogeneous and homogeneous OER catalysts, [8][9][10][11][12][13][14] and can thereby be said to be universal. The essential consequence of these relations is that ∆ *+ * and ∆ *++ * cannot be varied independently, so when they are not thermodynamically ideal, i.e. ∆ *++ * = ∆ *+ * + 2.46, they place a fundamental restriction on catalytic activity, as any catalyst which follows the conventional WNA mechanism will have a lower limit placed on the possible overpotential, a constraint often referred to as the "overpotential wall".
Recently, we highlighted the importance of a one-electron oxidation occurring at the metaloxo unit of certain molecular complexes, prior to O-O bond formation, which must be considered to correctly describe the impressive activity of the best homogeneous OER catalysts. 10 Furthermore, we showed that this additional step allows for the circumvention of the restriction imposed by the above scaling relation, leading to the possibility for catalysts to exhibit near-zero overpotentials provided the energy levels of redox intermediates are appropriately distributed. The elementary steps involved in this process and their associated Gibbs energy change are shown in Eqs. 7 and 8. The oxidation state of the oxo intermediates is given in parenthesis, and each adsorbate binding energy, ∆ 7 * , is defined in the ∆ 08 = ∆ +(9) * − ∆ +(: With this insight, we proposed a series of catalyst design principles to allow for the discovery of catalysts based on earth-abundant metals which can stabilize high oxidation states such as Mn and Fe. Mononuclear first-row transition metal complexes based on Cr, 15 Fe, [16][17][18]19 Co, [20][21][22][23] Ni, 24 and Cu 25,26 have been shown to be active towards the OER. Despite these exciting developments, earth-abundant catalysts have yet to show comparable activity to those based on single-site Ru in terms of turnover frequency (TOF) and overpotential. [27][28][29] In this study, we have evaluated and refined our design principles by studying intermediates across metals, ligands and oxidation states. Traditional computational OER studies apply the conventional OER descriptor, ∆ + * − ∆ *+ * which is thought optimal when it is equal to half the value of the intercept found in the scaling relationship between ∆ *+ * and ∆ *++ * . While this is effectively a pH-independent descriptor, OER catalyst performance is known to be pHdependent, with activity varying significantly in alkaline or acidic conditions. This is especially relevant due to the exigency of developing catalysts which are active in acid. 30 Alterations in activity due to pH can be attributed to interwoven factors such as catalyst instability, ionic transport, switching mechanisms or a combination thereof. In this study, we examine distinct reaction mechanisms by modelling descriptors in differing oxidation states, which is typically not done in studies utilizing the OER descriptor. Since one can expect more hydrogenated intermediates in acidic conditions, one would expect the extra oxidation mechanism to be less favored than in base. A concrete example of this is provided in an experimental Pourbaix 7 diagram reproduced in Supplementary Figure 1. We believe such a simple consideration could allow for pH-dependence using simple descriptors.
To test our design principles, herein we report a high-throughput computational screening of mononuclear complexes with the aim of accelerating the discovery of earth-abundant OER catalysts by focusing on non-critical elements which can stabilize high-valent M(V) intermediates and deepening our fundamental understanding of this reaction. In particular, we present a fully-automated approach to create transition metal complexes with molSimplify 31 via permutation of a selected set of ligands which have been previously seen in experimental OER catalysts. 32 As shown in Figure 2, five different geometries are generated from these ligands depending on their denticity, each of them containing one of six metal centers (i.e. Cr, Mn, Fe, Ru, Co, and Ni), leading to a total of 444 molecular catalysts which are octahedral upon the addition of the adsorbate. For each catalyst, *+(:::) * and *++(:::) * are evaluated, with *+(:9) * , *++(:9) * , +(:9) * and +(9) * calculated for a subset of these complexes.

Insights from Molecular OER Scaling Relations
With the 444 catalysts formed via the automatic structure generation outlined in Figure 2  Notably, we find that this data can be used to begin to rationalize past decades of experimental research into monometallic OER catalysts. 32 Ru complexes exhibit a tendency to shift below the line of best fit, meaning these complexes stabilize the HOO* intermediate to a greater extent than other metals, on average. This phenomenon, which reduces the minimum theoretical overpotential via the conventional WNA mechanism, was also noted in a computational study of IrO2 and SrIrO3 surfaces, and formed part of the rationalization for these materials' impressive activity in acidic media. 38 We note that this effect is seen for each of the investigated DFT functionals to varying degrees, as displayed in the Supplementary Figure 2. Similar to Ru, Ni catalysts are found to deviate from the line of best fit, which we will rationalize later. mechanism, 27,28 and we recently showed that both facets must be considered to explain their favorable activity. 10 These complexes, therefore, circumvent the scaling relationship by sidestepping the HOO* intermediate.
To further confirm the robustness and method-independence of the OER scaling, we applied the error estimation framework utility from BEEF-vdW to represent that relation for GGA functionals which are within one standard deviation of the mean value of this ensemble ( Figure   3b). This was performed similarly in a study on the functional-dependence of scaling relations, 39 although to the best of our knowledge, this is the first time that the analysis in Figure 3b has been done for the OER. In the GGA space, we find that the distribution of the binding energies for each metal is the same, which reinforces the functional independence of the OER scaling relation.
Interestingly, Figure 3a also indicates that the *++(:::) * values for Ni(III) complexes deviate from the linear scaling relationship. We posit that this is because the HOO ( Ni, Fe and Zn oxides can significantly enhance the rate of the OER. 40 The authors proposed that this effect stems from the induced magnetic field which facilitates the parallel spin alignment along the OER to yield the triplet oxygen product. Furthermore, a recent perspective has proposed that the magnetic enhancement effect could be related to a deviation from scaling relationships. 41 This prompted us to investigate this hypothesis by inspecting the spin density in the HOO* intermediate and the deviation of the *++(:::) * values from the scaling shown in Figure 3a.
To do this, we describe the collinearity of oxygen atoms, , using the Mulliken spin densities, , of the individual oxygens in the HOO* intermediate as: The result of plotting along with the corresponding Mulliken spin density on the metal ( < ) in the HOO* intermediate is shown in Figure 4, where we identify a distinct region defined by > 0.1 and < > 1 in which Ni catalysts deviating favorably from the OER scaling lie. The theoretical basis for spin-dependent oxygen electrochemistry has been developed and examined by Gracia and co-workers. [43][44][45] In those works, this effect was claimed to be exchange-mediated, which is borne out by our analysis as we only observe this effect when using DFT functionals which include Hartree-Fock exchange, as shown in Supplementary

Rational Screening of Molecular Catalysts
We next sought to test the possibility of near-zero overpotential catalysts via an extra-oxidation mechanism. 10 With this aim, we computed the conventional OER descriptor, +(:9) * − *+(:::) * , for catalysts in Figure 3a which exhibit a value of *+(:::) * ≤ 1.5 eV, leading to the volcano plot presented in Figure 5a. Use of this volcano representation has found useful application across OER research. 46,47 This upper limit was chosen to place a lower limit on the overpotential due to *+(:::) * , where our aim is to find catalysts with overpotentials that circumvent the "overpotential wall". While catalysts which are close to the peak of the volcano in Figure 5a demand a considerable overpotential, they evolve oxygen via PCETs exclusively, thereby avoiding the building up of positive charge and fostering redox potential levelling. 48 The best leads for future experiments which exhibit a PCET-only mechanism are chosen so that < 0.55, (since this is close to the lower limit of 0.5 imposed by the scaling in Figure 3a) +(:9) * < *+(:9) * , and *++(:::) * < +(9) * , since we do not want to perform the extra oxidation in this case. This identifies a 32-Fe catalyst composed of the 3a tridentate and 2a bidentate ligands as seen in the labelled catalysts from Figure 2. The 32-Fe catalyst is particularly interesting as we expect that the bipyrimidine ligand (labelled 2a in Figure 2) could be used to immobilize the catalyst onto a solid support. 49 Also included in this plot as black-outlined datapoints are catalysts which are experimentally known to catalyze the OER.
We note that the volcano peak at (1.7 eV, 0.47 V) in Figure 5a is larger than the overpotential wall previously reported by our group, which we attribute to the fact that individual metals have slightly different scaling relationships and intercepts, as seen in Supplementary Figure 5.
This difference leads directly to a difference in overpotential walls. In a recent work, we The volcano plot for this higher oxidation OER descriptor is shown in Figure 5c. Interestingly, we find that Ru and Cr catalysts predominate in the left leg of the volcano, implying a relative ease in the formation of the M(V) oxidation state, while Fe candidates lie on either side, and Mn and Co exclusively lying to the right leg of the volcano. In calculating the overpotential for Creating lines of best fit through the antiferromagnetic complexes in Supplementary Figure 7, we observe a slope less than 1. In these cases, the slope can be predicted by using Eq. 11. The results of using Eq. 11 to predict the slopes of lines of best in Figure 5d are also outlined in Table 1 can also be confirmed from the calculated Wiberg bond orders in Supplementary Table 1. As such, one would expect them to exhibit a slope of 2 since this is the expected value for precisely this type of bond, based on the bond order conservation principle, which our equations do not predict well. However, we note that the cluster of Cr(IV)-O points (black circles with yellow outlines in Figure 5d) are the most closely-packed subset of data, which may obscure the true scaling relation due to the sensitivity of the lines of best fit to new points that lie outside the bounds of the considered points. Table 1. Predicted slopes using the spin-based descriptors from Eqs. 10 and 11. The slope m was found by fitting the line of best fit to the relevant subset. We note that the Co data in Figure 5d also agrees with our recent studies of Co-based polyoxometalates (POMs), which found that these catalysts cannot show overpotentials less than 0.8 V due to the restriction placed on oxygen evolution by the *+(K) * vs. With all the knowledge acquired above and using the data from Figure 5b, we can now propose catalysts for experimental realization on the basis of predicted overpotentials less than 400 mV and consisting of earth-abundant metals. For this, we filter candidates by enforcing that +(9) * < L;+(:::) * so the extra oxidation is possible, and that the Gibbs energy change associated with the atom-proton transfer (APT), L;+(:::) * − +(9) * , be less than 1eV, as this is the largest value exhibited by the set of transition metal complexes known experimentally to catalyze the OER (in particular, the red square in Figure 5b). This leads to 9, predominantly Cr-based, low overpotential catalysts, shown in Figure 6. While the low overpotentials in these catalysts are very encouraging, we also need to consider the conditions under which the extra oxidation mechanism is expected occur. In particular, to properly determine which mechanism is expected to occur, the p M for the HO(IV)*/O(IV)* pair needs to be calculated. Doing this using only the Gibbs energy differences of the transition metal species involved has proven challenging in protic environments such as water. 58 However, we  Each of the tailored catalysts shown in Figure 6 is based on a ligand which has been previously seen in molecular Ru OER catalysts. 32 Assuming these catalysts can be made, experimental investigation will lead either to a refinement of our understanding of OER by rebutting our model, in the worst case, or novel low-cost catalysts with high activity in the best case.
The primary issue foreseen for these catalysts would be the APT, and it must be noted that low In summary, this study presents the first automated approach to rational computational design of molecular OER catalysts based on earth-abundant elements. Through the application of key design principles and the fastidious choice of ligand scaffolds, we have solidified methodindependent scaling relations and provided fundamental insights into molecular water oxidation catalysts and the OER generally. We then proposed a link between deviations from our scaling relation in nickel catalysts to recent experimental OER observations. Further screening led us to develop metal-specific scaling relations between O* and HO* intermediates that did not exhibit typical behavior in expected slope, which we justified on the basis of a straightforward spin-based heuristic. As a result of the computational screening, we identified 9 low-overpotential, earth-abundant catalysts with theoretical overpotentials below 400 mV.
We further filtered these catalysts based on their ability to evolve oxygen by the assumed mechanism in acid, leading to 3 promising Cr-based catalysts. Our data also has the potential to be used as the springboard for the machine learning-driven discovery of OER catalysts, further boosted by the scaling relations presented in Figures 3 and 5d. If the full potential of such ML models were realized they would help to avoid explorations of the chemical landscape that are not fruitful, while allowing researchers to trawl this space more effectively. The set of ligands used for this study were chosen due to their presence in other oxygen evolving complexes, therefore, it is inherently restricted by what mononuclear OER catalysts already exist. This could be overcome when generative models, which have been used to faithfully reproduce organic molecules, 61 are applied to transition metal complexes. Then, machine learning approaches will become even more powerful, with the ability to amplify the effects of high-speed predictions while testing catalysts over a much broader region of chemical space.
Future work will focus on studying the stability of these complexes against ligand oxidation or replacement.
Going beyond trial-and-error catalyst discovery for renewable energy technologies will require different approaches. This work emphasizes that large-scale computational analyses are one such approach in homogeneous OER catalysis.

Computational Methods
DFT calculations were performed at the TPSSh, M06-L and B3LYP level using the Gaussian09  The set of ligands seen in Figure 2 were chosen as they had seen application in experimentally realized catalysts. The charges of these ligands were pre-assigned and used to define the total charge of the catalyst. The spin state of metals in each reaction intermediate was determined by investigating which multiplicity was the most stable for a representative set of molecules.
In particular, if for every metal in a certain oxidation state a specific multiplicity was always lower in energy by 0.3 eV, that spin state was assigned throughout for all calculations.
Following this criterion, all Mn oxidation states, Fe(II), Co(II), Cr(II) were assigned to be highspin, while all Ru oxidation states, Co(III), Ni(II), were assigned low-spin. For Fe(III) complexes, both high and low-spin structures were modelled as no prevalent spin-state was determined. The results of this investigation are shown in Supplementary Table S2. The spinstates exhibiting the lowest energy were subsequently modelled using the M06-L and B3LYP exchange-correlation functionals, while the low-spin state was modelled regardless of the lowest energy spin state from TPSSh as GGA functionals are known to favor low-spin states. 71 generous provision of computational facilities and support. Authors also thank Mr. Eric Mates-Torres for providing assistance with the graphical design.