From computational screening to the synthesis of a promising OER catalyst

The search for new materials can be laborious and expensive. Given the challenges that mankind faces today concerning the climate change crisis, the need to accelerate materials discovery for applications like water-splitting could be very relevant for a renewable economy. In this work, we introduce a computational framework to predict the activity of oxygen evolution reaction (OER) catalysts, in order to accelerate the discovery of materials that can facilitate water splitting. We use this framework to screen 6155 ternary-phase spinel oxides and have isolated 33 candidates which are predicted to have potentially high OER activity. We have also trained a machine learning model to predict the binding energies of the *O, *OH and *OOH intermediates calculated within this workflow to gain a deeper understanding of the relationship between electronic structure descriptors and OER activity. Out of the 33 candidates predicted to have high OER activity, we have synthesized three compounds and characterized them using linear sweep voltammetry to gauge their performance in OER. From these three catalyst materials, we have identified a new material, Co2.5Ga0.5O4, that is competitive with benchmark OER catalysts in the literature with a low overpotential of 220 mV at 10 mA cm−2 and a Tafel slope at 56.0 mV dec−1. Given the vast size of chemical space as well as the success of this technique to date, we believe that further application of this computational framework based on the high-throughput virtual screening of materials can lead to the discovery of additional novel, high-performing OER catalysts.


Introduction
2][3] The current state-of-the-art materials for OER catalysis oen contain materials like IrO 2 and RuO 2 . 2,3RuO 2 for example, typically has overpotentials between 250-350 mV at 10 mA cm −2 of current density and Tafel slopes between 50-70 mV dec −1 in basic conditions of 1 M KOH. 35,46,51owever, IrO 2 and RuO 2 are rare and expensive, necessitating the development of cheaper catalysts that are comparable to them on both metrics.][4][5] The crystal structure of spinel oxides is made up of oxygen anions arranged within a cubic close-packed sublattice while metallic cations are positioned within the tetrahedral and octahedral interstitial sites between the anions. 6,7The basic composition of a ternary-phase spinel oxide is A x B 3−x O 4 , where A and B are two different metals. 6,7In this structure, the distribution of both metals across both coordination geometries can be represented by the formula (A x−3 B 1−x+3 ) T d (A 3 B 2−3 ) O h where T d refers to tetrahedrally coordinated cations, O h refers to octahedrally coordinated cations and 3 refers to the inversion parameter (3 = 0x) 7 (Fig. 1B).If 3 equals 0, the structure is a normal spinel, and if 3 equals x, the structure is an inverse spinel. 7The presence of transition metals in both the AO 4 and BO 6 structural units likely allows both of them to contribute to OER activity.This could potentially increase the structural diversity and tunability of spinel oxide catalysts for OER catalysis, making them an interesting system to study.
In recent years, spinels have proven to be a rich source of earth abundant metal oxide catalysts capable of catalyzing OER efficiently. 5Most spinel oxides reported for OER tend to be Febased or Co-based, with some examples of highly active spinel oxides being NiFe 2 O 4 , CoFe 2 O 4 and CuCo 2 O 4 . 5For example, Liu and coworkers fabricated NiFe 2 O 4 nanospindles on an Fe 3 Ni foam substrate that exhibited an overpotential of 262 mV at a current density of 10 mA cm −2 , as well as a Tafel slope of 39.5 mV dec −1 in 1 M KOH. 8 Lu and coworkers synthesized CoFe 2 O 4 nanoparticles on carbon nanorods that had an overpotential of 240 mV at 10 mA cm −2 as well as a Tafel slope of 45 mV dec −1 in 1 M KOH. 9 Yadav and coworkers synthesized 3D CuCo 2 O 4 nanoowers on a carbon cloth substrate that had an overpotential of 288 mV at 10 mA cm −2 and a Tafel slope of 64.2 mV dec −1 in 1 M KOH. 10 The efficacy of these catalysts demonstrate that spinels could be a promising source of new materials for OER catalysis.However, while considerable work has been done so far to discover new spinel oxide catalysts for OER, the entire possible chemical space of spinels has yet to be explored.There could potentially be more efficient spinel oxide catalyst compositions that haven't been discovered yet.
It is not feasible to explore the entire chemical space experimentally, due to the immense number of permutations possible with each element in the periodic table.According to one estimate by Walsh and coworkers, the number of possible ternary phase inorganic materials that can be synthesized out of 103 elements in the periodic table is greater than 32 million, even aer constraining this space by imposing charge neutrality and electronegativity rules. 11However, it is possible to explore a greater fraction of this space computationally, in order to narrow down possible candidates before synthesizing them in a lab.This approach towards materials discovery has had a signicant impact in elds as diverse as catalysis, 12,13 Li-ion batteries, 14,15 thermoelectrics, 16,17 organic light-emitting diodes, 18 and transparent conducting oxides. 19Within the eld of catalysis, there have been high-throughput computational studies that have aimed to discover new OER catalysts.For example, Xu and coworkers utilized a bandcenter descriptor to screen 3d spinel oxides for OER activity. 13Ulissi and coworkers screened 2600 equimolar Ir-containing bimetallic oxides for acid-stable OER catalysts and identied 14 possible candidates predicted to be stable under acidic conditions. 12ørskov and coworkers screened 47 814 nonbinary metal oxides in the Materials Project and identied 68 possible acid-stable OER catalysts. 20In each of these studies, computational screening was used to predict the stability or activity of OER catalysts.However, to the best of our knowledge, no study has successfully utilized a high-throughput computational screen to discover novel, highly-active OER catalysts that have been experimentally demonstrated to compete with current benchmark catalysts reported in the literature.In this paper, we discuss a high-throughput computational workow that enabled us to achieve this goal.
In this study, we built a computational workow to screen for new highly active spinel oxides for OER in basic pH conditions.A computational database of ternary-phase spinel oxide materials comprised of 52 elements in the periodic table is rst constructed.These elements include alkaline, alkaline-earth, transition (with the exception of articial Tc), post-transition and some lanthanide metal elements (La, Ce, Nd, Gd and Lu) and are circled in Fig. S1.† Every possible combination of elements that could be constructed from this list was explored, for an overall spinel system A x B 3−x O 4 where A and B are two different metallic elements.For each system A x B 3−x O 4 , three different values of x (x = 0.5, x = 1, x = 1.5) are chosen in order to screen for three unique compositions A 0.5 B 2.5 O 4 , AB 2 O 4 and A 1.5 B 1.5 O 4 .Both normal spinel (3 = 0) and inverse spinel (3 = x) structures for each specic composition were explored.Once this database was created, the thermodynamic stability of these materials was assessed to isolate materials that are likely to exist in nature.Following that, the theoretical overpotential of OER on the surface of these materials was computed.We used the results of these computations to create a machine learning model to predict OER activity using electronic structure descriptors of bulk metal, bulk oxygen, surface metal, surface oxygen and adsorbate oxygen atoms as features.This allowed us to further probe the relationship between electronic structure descriptors and OER activity.The trends observed in the results of this high-throughput computational screen were applied to synthesize three catalysts, Co 2.5 Ga 0.5 O 4 , Co 1.5 Ga 1.5 O 4 and Co 1.5 Al 1.5 O 4 and characterize them for OER activity.We discovered that Co 2.5 Ga 0.5 O 4 had an overpotential of 220 mV at a current density of 10 mA cm −2 and a Tafel slope of 56.0 mV dec −1 , making it a new highly active OER catalyst successfully discovered using this novel computational workow.The success of this workow at predicting a real catalyst highly active for OER implies that it can be adapted to discover new highly active OER materials belonging to other classes of materials as well.

High throughput computational workflow
A computational database of 6155 ternary-phase spinel oxides was rst created from 52 transition, post-transition, alkaline, alkaline earth and lanthanide metal elements.For each combination of metal elements, three different ratios of metals were explored: A 0.5 B 2.5 O 4 , AB 2 O 4 and A 1.5 B 1.5 O 4 .Both normal and inverse spinels of each permutation of materials were also considered, leading to six different spinel structures: (A 0.5 -B 0.5 ) The spinel oxide Fe 3 O 4 belonging to the Fd3m spacegroup was used as a prototype structure for every single combination of elements except those utilizing Mn; for spinel oxides containing Mn in the octahedral positions, Mn 3 O 4 was used as a prototype structure instead in order to account for the Jahn-Teller distortion of octahedrally coordinated Mn 3+ ions. 21DFT calculations were subsequently performed to relax the structure of each compound and obtain its ground state energy.
In the next step of the high-throughput workow, the thermodynamic stability of each composition was assessed.First, the energies of the normal and the inverse spinel structures of each composition were compared to determine the most probable ground-state conguration of each composition, with the more stable one being chosen for the next step.Then, we constructed a convex hull based on the bulk energies of all known materials encompassing the phase space of the elements that constitute each spinel oxide composition.The Open Quantum Materials Database was used to construct the phase space of the elements that constitute each compound. 22,23The thermodynamic stability of each composition was subsequently determined by calculating the distance of the bulk energy of each composition from this convex hull.The total energy of each compound was compared to a linear combination of possible decomposition products that lie on the convex hull of the phase space in order to calculate this distance from the convex hull. 15,24If the difference between the energies was within 0.030 eV per atom (30 meV per atom), the compound would be considered thermodynamically stable. 25he heat maps in Fig. 2A-C show the formation energies of the different spinel oxide compositions.The higher density of formation energies below 0.030 eV per atom demonstrates that the AB 2 O 4 spinel composition is generally the most stable composition type, followed by A 0.5 B 2.5 O 4 and then A 1.5 B 1.5 O 4 .This is probably because in the latter two compositions, at least one of the cations should occupy both tetrahedral and octahedral sites in the spinel oxide structure.This cation should be able to assume both +2 and +3 oxidation states and should not have any strong preference for tetrahedral or octahedral coordination geometry.This constraint would limit the number of possible stable compositions for both A 0.5 B 2.5 O 4 and A 1.5 B 1.5 O 4 compositions relative to AB 2 O 4 .It can be observed that spinel oxides containing 3d metal cations are typically more stable than spinel oxides containing 4d and 5d metal cations.Mn, Co and Fe-containing spinels, in particular, are more likely than other metal cations to have formation energies that are within 0.030 eV of their convex hulls.This makes Mn, Co and Fecontaining spinels more likely to exist in nature. 25This observation can rst be partly rationalized by the fact that Mn, Co and Fe are stable in the +2 and +3 oxidation states. 26Furthermore, according to Kocevski and coworkers, Mn, Co and Fe ions do not have a strong preference for either tetrahedral or octahedral sites in spinels and are frequently seen in either of them. 27herefore, Co, Mn and Fe containing spinel structures have a greater amount of exibility in accommodating a wider variety of other cations as a result, making them more ubiquitous than other elements; an observation that was further corroborated by the calculations performed in this work. 27Spinels containing metals of higher oxidation states can be stabilized through the engineering of oxygen vacancies into the structure. 6This study, however, did not consider defects such as oxygen vacancies in order to keep computational costs reasonable.Aer this step had been performed, 290 materials were identied as being potentially thermodynamically stable.They were selected for the nal step of this high-throughput computational screen, where their OER activity was estimated.
In the nal step of this workow, the catalytic activity of each of these 290 spinel oxide materials is predicted by calculating the theoretical overpotential using the computational hydrogen electrode method proposed by Nørskov and coworkers. 28The computational hydrogen electrode states that the chemical potential of the proton and electron pairs can be described as half of the chemical potential of the H 2 molecule since these species are in chemical equilibrium at 0 V vs. the reversible hydrogen electrode (RHE). 28We can use this model to estimate the Gibbs free energies of the rst two reaction steps of the adsorbate evolution mechanism (AEM) by calculating the binding energies of the *O, *OH and *OOH intermediates of OER.To calculate these energies, we created slabs of each material, relaxed the surface of the slabs, and then added the relaxed *O, *OH and *OOH molecules on top of the metal sites of the slab surface.0][31] Higher index surfaces were not included in this step to keep the computational cost reasonable.The surface energies of all these slabs were calculated and the slab with the lowest surface energy was identied.Next, *O, *OH and *OOH molecules were xed onto the surface of the slab with the lowest surface energy.Each of these slabs has a 50% monolayer coverage of *O, *OH and *OOH intermediates respectively.This assumption was made in order to limit computational cost that could come with considering all the different permutations of intermediates possible.We note that surface coverages can change depending on the applied potential, and that the potential window of the stability of surface coverage should be taken into account in future studies where the predicted overpotential of catalysts is computationally assessed.The adsorption energy of each of these intermediates was calculated and used to determine the theoretical overpotential of OER on the surface of each material.
The distribution of the calculated theoretical overpotentials of the compounds is plotted in Fig. 2D.Out of the 290 materials for which these theoretical overpotential calculations were successfully performed, 33 compounds show potential to be useful OER catalysts with theoretical overpotentials below 0.5 eV.Spinel oxides like CuCo 2 O 4 and CoFe 2 O 4 that are known in the literature to be highly active for OER were also present amongst these 33 compounds. 5,9,10This further conrms the efficacy of this screen at discovering real materials that catalyze OER efficiently.To better understand the impact of elemental composition on catalytic activity, a frequency plot of the elements in all calculated compounds is plotted with respect to theoretical overpotential in Fig. 2E.It seems that Cobalt forms the largest number of catalytically active spinel oxides for OER, since they have the greatest number of low overpotential (<0.5 V) compounds compared to others.This insight corroborates existing trends in the literature on the use of spinel oxides for OER since most spinels reported for OER are either Co or Febased. 5e then use the results of these OER calculations to train a machine learning model to predict the binding energies of the *O, *OH and *OOH intermediates of OER.
Predicting *O, *OH and *OOH binding energies based on electronic structure descriptors Now that we have demonstrated that the success of this computational workow at discovering new catalysts, we have decided to use the database created by this workow to train a machine learning model to predict the binding energies of *O, *OH and *OOH intermediates.The objective of this task is to use this model to identify patterns, based on DFT-calculated electronic structure descriptors, that could aid with the swi screening of OER catalysts.7][68] However, the main weakness of solely using bulk descriptors is the likelihood that surface effects end up getting ignored.Stoerzinger and coworkers demonstrate this reality when they showed that OER activity on the surface of RuO 2 (100) and RuO 2 (101) was higher than their (110) and (111) counterparts in basic conditions due to the higher concentration of coordinatively unsaturated sites. 73Any comprehensive assessment of the link between OER activity and electronic structure descriptors would have to include electronic structure descriptors associated with both bulk and surface atoms as well.Inspired by a recent study by Lunger and coworkers which concluded that the Bader charge and the O2p bandcenter of surface oxygen was predictive of *O, *OH and *OOH binding energies on the surface of perovskite slabs, we decided to use the bandcenters and Bader charges of the atoms in each slab as features for our model. 71We rst calculated the O2p bandcenters of the surface oxygen, bulk oxygen and adsorbate oxygen atoms, the M3d bandcenters of the surface metal and bulk metal atoms and the Bader charges of the surface oxygen, bulk oxygen, adsorbate oxygen, surface metal and bulk metal atoms.These descriptors were then used as training data for a random-forest regression model.We used a random-forest model since they are known to be simple to use, accurate, and capable of dealing with small sample sizes and high-dimensional feature spaces. 75In addition, since random forests are also interpretable, we can also look at the relative importance of the input features to nd electronic structure descriptors that can correlate well with predicted OER activity. 75This model was able to predict the binding energies of *O, *OH and *OOH intermediates quite accurately, with an R 2 score of 0.83 for the training set and 0.78 for the test set (Fig. 3A).We also found that the outliers of our model, dened as a prediction which differs from the calculated value by more than 1 eV, were far more likely to contain novel elements like Nb, Pb, Cd, and Ir.The failure of our model to predict the theoretical overpotentials of these structures can be attributed to the fact that these materials were underrepresented in our dataset, since spinel oxides fabricated using these elements were far more likely to be unstable and were therefore ltered out by the convex hull stability screen.Such outliers, thankfully, are rare and can only be attributed to 4% of the calculations in our screen, further demonstrating the efficacy of our model at predicting OER binding energies.
The relative importance of each feature for the prediction of binding energies was then extracted from the random forest model.The Bader charge of adsorbate oxygen was the most important feature of all of them, with the O2p bandcenter of the adsorbate oxygen atom coming at a distant second (Fig. 3B).The signicance of electronic structure features associated with adsorbate oxygen for OER is corroborated by previous work by Nørskov and coworkers, which demonstrated a link between the O2p bandcenter of adsorbate oxygen and the calculated theoretical overpotential. 69We also decided to use LASSO to predict the binding energies of the OER intermediates, in order to determine how universal the features we identied were with respect to their relationship with the binding energies of OER (Fig. S14 †).Just like with the random-forest model, the Bader charge of adsorbed oxygen was determined to be the most important feature.The second most important feature in the LASSO model, however, was deemed to be the Bader charge of surface oxygen.The LASSO model performed poorly at predicting binding energies compared to the random-forest model; with an R 2 value of 0.65 for the training set, an R 2 value of 0.72 for the test set and a median absolute error of 0.56 eV for both the training and the test sets.As such, we decided to use the results of the random-forest model for further analysis instead.
We then further examine the relationship between the two most important features identied by our random-forest model, the Bader charge and the O2p bandcenter of adsorbate oxygen,  69 Since the optimal DE O -DE OH would be between 1.5 and 1.7 eV, the ideal catalyst is likely to have an oxygen Bader charge between −0.3e − and −0.6e − and an O-2p bandcenter between 2 eV and 4 eV below the Fermi level.Screening for catalysts based on both descriptors is likely to lead to the identication of more highly active catalysts for OER.

Synthesis and characterization of compounds
Out of the 33 possible candidates predicted to have theoretical overpotentials below 0.5 V, three Co-based spinel oxides, Co 1.5 Al 1.5 O 4 , Co 2.5 Ga 0.5 O 4 , Co 1.5 Ga 1.5 O 4 , were selected for synthesis and experimental characterization based on their novelty.They were synthesized via electrodeposition onto a carbon paper substrate.In order to obtain compounds that have a similar composition to these predicted materials, stoichiometric mixtures of the salts Co(NO 3 ) 2 , Ga(NO 3 ) 3 and Al(OH)(C 2 H 3 O 2 ) 2 were dissolved in the electrolyte solution used for electrodeposition in all three cases.Co 3 O 4 was also synthesized for benchmarking purposes since it is well-known to be active for OER. 49Aer synthesizing these materials, their morphologies and chemical compositions were characterized with SEM and energy dispersive X-ray spectroscopy (EDS) mapping.The SEM images of all four catalysts (Fig. 4) reveal the presence of an amorphous material evenly deposited on the surface of the carbon substrate.EDS maps (Fig. S2  heterogeneity in their distribution over the surface of the material.This seems to indicate the presence of another phase. To further investigate the structure of the synthesized catalysts, all four materials were characterized by XRD (Fig. 5A).The major peaks corresponding to the spinel phase are found at 31.3°, 36.9°,44.8°, 59.4°and 65.2°, which are the characteristic peaks of Co 3 O 4 . 32The peaks at 44.8°and 55.0°found in all four spectra can be attributed to the carbon paper substrate (Fig. S11 †).The presence of the spinel peaks in all four spectra demonstrates that spinel oxides have been successfully synthesized in all four cases.However, in addition to the peaks corresponding to the spinel oxide and the carbon substrate phases, new peaks corresponding to a different phase also appear in the Co-Ga and Co-Al oxide spectra.These can be attributed to b-Ga 2 O 3 in Co-Ga oxide spectra and q-Al 2 O 3 in the Co-Al oxide spectrum. 32,33Since neither Ga 2 O 3 nor Al 2 O 3 are redox-active, they are unlikely to make any signicant contribution to OER activity. 34,35he composition of the synthesized materials was characterized with XPS (Fig. S6-S9 †).The high-resolution Co2p spectra of each of these materials (Fig. 5B-E) show evidence of spinorbit splitting into 2p 3/2 and 2p 1/2 components at ∼780 eV and  ∼795 eV respectively, with shake-up satellites for each of these components present at ∼790 eV and ∼805 eV. 36Each of these spectra resembles the Co2p spectra of Co 3 O 4 , further corroborating the results of the XRD that have demonstrated the successful synthesis of cobalt spinel oxide. 36In order to determine the stoichiometry of each of these catalysts, the Co2p 2/3 spectra were deconvoluted based on the Co2p 2/3 peak tting parameters by Biesinger et al. and the ratio between the areas of the peaks associated with the Co 2+ and Co 3+ oxidation states were calculated. 36,37Since both Ga and Al are not stable in the +2 oxidation state, they are both likely to displace Co 3+ if successfully incorporated into the spinel oxide structure.Thus, comparing the Co 2+ and Co 3+ ratio will help determine the stoichiometry.The ratios of the peak areas of Co 2+ to Co 3+ are 1 : 2, 1 : 0.7, 1 : 1.

Catalytic characterization
Aer the synthesis and structural characterization of these catalysts, they were electrochemically characterized for OER activity by linear sweep voltammetry (LSV) in a 1 M KOH solution (Fig. 6B).The overpotential in this report is dened as the iR-corrected potential at 10 mA cm , it also has a much lower Tafel slope at 56.0 mV dec −1 than either catalyst (Fig. 6C).This allows Co 2.5 Ga 0.5 O 4 to start outcompeting the rest of the catalysts at around 1.55 V RHE .
The LSV plot of each catalyst was further corrected using the estimated electrochemically active surface area (ECSA) of each catalyst in order to determine the intrinsic activity of each  S5.† material for OER (Fig. S10 †).Co 3 O 4 has the highest ECSA out of all four catalysts, followed by Co 1.7 Al 1.3 O 4 , Co 2.5 Ga 0.5 O 4 and then GaCo 2 O 4 .Aer adjusting for the ECSA of all catalysts (Fig. 6A), Co 2.5 Ga 0.5 O 4 starts outcompeting Co 3 O 4 even earlier, at about 1.52 V RHE .

Benchmarking catalysts against stateof-the-art materials
The catalysts explored in this report were then compared to 38 other benchmark catalysts in the literature in order to assess their performance (Fig. 6D). [8][9][10]The catalysts explored in this report were compared to four families of catalysts: Nonspinel Co catalysts (marked as Others), spinel oxide catalysts, Ni-Fe catalysts and Ru/Ir-based catalysts.The ideal catalyst would be one with a lower Tafel slope and a lower overpotential than the state-of-the-art catalysts in the literature that currently exist.While Co 1.7 Al 1.3 O 4 has a low overpotential at a current density of 10 mA cm −2 , its high Tafel slope compared to the benchmark catalysts examined makes it uncompetitive.The performance of the GaCo 2 O 4 catalyst with respect to its overpotential and Tafel slope was also uncompetitive compared to the benchmark catalysts.On the other hand, Co 2.5 Ga 0.5 O 4 is as competitive as some of the better benchmark materials within the Ni-Fe family of catalysts on both the overpotential and Tafel slope metrics.It has also outperformed all the Co, spinel oxide and Ru catalysts examined in this report, with respect to both its overpotential and its Tafel slope.The high performance of Co 2.5 Ga 0.5 O 4 demonstrates the efficacy of this high-throughput computational workow at discovering promising OER catalysts.

Conclusions
In this work, a high-throughput computational framework was developed in order to screen for novel low-cost materials for OER in basic conditions.Out of the 6155 binary spinel oxides that were investigated for this work, 33 were predicted to have low theoretical overpotentials below 0.5 eV.This made them ideal candidates for further study through experiment.Based on an analysis of the results of the overpotential screen, Ga and Al-doped Co 3 O 4 were investigated using LSV to determine their OER activity.The data indicated that Ga 0.5 Co 2.5 O 4 was highly active for OER, surpassing many other state-of-the-art catalysts that have been reported in the literature.This catalyst is, to the best of our knowledge, novel and has not been investigated for OER before.These results further demonstrate the strength of this computational framework for facilitating the discovery of novel materials for OER.

DFT calculation details
DFT calculations were performed to optimize the structure of each spinel oxide in this database and calculate their energies.Every DFT calculation in this study was performed with the Vienna ab initio Simulation Package (VASP).The Projector-Augmented Wave (PAW) method was used to model the core electrons. 62The Perdew-Burke-Ernzerhof functional which utilizes the generalized gradient approximation approach was used to describe the exchange-correlation effects. 63For spinel oxides containing 3d metal elements with the exception of Zn, spin-polarized DFT calculations were performed; for all other materials non spin-polarized calculations were performed instead.The Hubbard U correction was employed for materials containing 3d transition metals; the exact values used in this study are described in ESI Table S4.† The energy cut-off used in these calculations was 520 eV.Each and every single one of these parameters was optimized to ensure that the right balance between computational cost and accuracy was achieved.
In the rst step of this screen, three structural optimizations were performed on the bulk structure of each spinel oxide within this constructed dataset before the energy of each structure was determined.A gamma-centered k-point mesh with a density of 8000 k-points/number of unit cell atoms was generated for each material using Pymatgen. 29Both the shape of the unit cell and the positions of the atoms were allowed to uctuate until the energy convergence criterion of 10 −5 eV between each self-consistent eld iteration step and the force convergence criterion of 0.01 eV Å −1 were met.Once the relaxation steps had been concluded, the shape of the unit cell and the positions of the atoms were xed in order to evaluate the energy of each structure.The Brillouin zone was integrated using the tetrahedron method with Blöchl corrections for this energy calculation. 64eating slabs and calculating the binding energies of the *O, *OH and *OOH intermediates When creating all slabs using Pymatgen, the width of each slab was set at 6 Å and the vacuum space was set at 20 Å.The energies of each slab were then calculated using DFT.Once the slab energies had been calculated, the surface energy of each slab was calculated using the following formula: 31 where E surf is the surface energy, E slab is the calculated energy of the slab, E bulk is the energy (per atom) of the bulk structure, n is the number of atoms in a unit cell and A is the surface area of the slab.The slab with the lowest surface energy was then isolated and utilized for the nal step of this screen.VASP was used in the DFT calculations of slab energies.The top two layers of atoms in each slab were allowed to relax for two steps while the bottom layers were frozen.A gamma-centered kpoint mesh with a density of 1000 k-points/number of unit cell atoms was generated for each material using Pymatgen. 29In each structural optimization step, the atoms in the top two layers were allowed to relax until the energy convergence criterion of 10 −5 eV and the force convergence criterion of 0.1 eV Å −1 were met.Once the relaxation steps had been concluded, the positions of all atoms in each slab was xed in order to evaluate the energy of the slab, where the tetrahedron method with Blöchl corrections was used to integrate the Brillouin zone. 64he energy cutoff used in these slab calculations was 400 eV.The slab surface with the lowest surface energy was used to calculate the binding energies of the *O, *OH and *OOH intermediates.
The binding energies of the *O, *OH and *OOH intermediates were then calculated by adsorbing each molecule on top of the metal sites on the relaxed surfaces of the slab with the lowest surface energy.Both A and B sites on the surface of the spinel catalysts were considered.The top two layers, along with the molecules, are allowed to relax for two steps while the bottom layers are kept xed.The energies of each slab were calculated based on the same parameters used to calculate the surface energies in the previous step.The adsorption energies of the OER intermediates on each slab are subsequently calculated using the following equations: 28 where DG ads is the adsorption energy of the adsorbate (*O, *OH or *OOH), DE ads is the electronic adsorption energy of the adsorbate, DE ZPE is the zero-point vibrational energy difference between adsorbed and gaseous species, TDS ads is the entropy difference between the gaseous and adsorbed species, E slab is the energy of the clean slab, E slab+OxHy is the energy of the slab with the adsorbate species on the surface, E O is the energy of an oxygen atom, E H is the energy of a hydrogen atom, E H 2 O is the energy of a water molecule and E H 2 is the energy of a hydrogen molecule.

Calculating the theoretical overpotential based on binding energies
Fig. S12 † shows the different steps of the adsorbate evolution mechanism (AEM) proposed for water oxidation and the equations necessary to calculate each step of the mechanism.According to the computational hydrogen electrode model proposed by Nørskov and coworkers, steps 1 and 2 are typically the potential limiting steps of the AEM mechanism of OER. 28,65herefore, in order to calculate the theoretical overpotential, DG 1 and DG 2 must be calculated based on the binding energies of the intermediates calculated in the previous section.1.23 eV shall then be subtracted from the larger of the two energies in order to calculate the theoretical overpotential. 65eating the machine learning model The random forest model was created using the Python package scikit-learn.Features such as the bandcenters of the oxygen and metal atoms were calculated, relative to the Fermi level, based on the formula below: where 3refers to the bandcenter, 3 refers to the energy level, r refers to the density of states at that energy level, 3 max refers to the maximum energy level within this range of integration and 3 min refers to the minimum energy level within this range of integration.The Bader charge was calculated using a Bader charge analysis code provided by the Henkelman group. 74The average bandcenters and Bader charges of all atoms were used to train the random forest model.To train and test the models, the dataset was split using train_test_split from scikit-learn, saving 20% as test set.The optimal model hyperparameters were identied using a Bayesian optimization approach (Bayes-SearchCV from scikit-optimize, on RandomForestRegressor from scikit-learn) on the training set with a 5-fold cross validation.The hyperparameters of the model that were optimized were the number of trees in the random forest model, the maximum depth of each tree, the minimum number of samples required to split an internal node, the minimum number of samples required to be at a leaf node and the number of features to consider when looking for the best split.For comparison, LASSO was also tested to predict the binding energies of the OER intermediates (Fig. S14 †).The same features used to train the random-forest model were used to train the LASSO model as well, and the hyperparameters were also optimized using Bayesian optimization.The dataset for the LASSO model was also split similarly to the random-forest model, with 20% of the dataset being saved for the test set.The hyperparameter tuned for the LASSO model is the alpha parameter that determines the magnitude of the L1 regularization term.
Synthesis of Co 3 O 4 catalysts and metal cation doped Co 3 O 4 catalysts on carbon paper substrates.These catalysts were synthesized using an electrodeposition method in a standard three electrode cell consisting of carbon paper as working electrode, carbon rod as a counter electrode and saturated calomel electrode(SCE) as reference electrode at the room temperature.Co(NO 3 ) 2 $6H 2 O (0.1 M) was dissolved in water as the electrolyte.Carbon paper with the size of 0.5 cm × 0.5 cm was then immersed into the electrolyte for the electrodeposition of Co 3 O 4 -species.Electrodeposition was performed at a potential range of −1.0 to +0.2 V/SCE for 50 scans using an Autolab PGSTAT302N workstation.During the deposition process, the stirring rate was kept constant at 1000 rpm.Aer that, the catalyst was calcined under vacuum at 250 °C for 1 hour and then annealed in air at 350 °C for 4 hours.The loading mass is about 1.6 mg.The Co 1.5 Ga 1.5 O 4 , Co 2.5 Ga 0.5 O 4 , and Co 1.5 Al 1.5 O 4 catalysts were synthesized following a process similar to that of the Co 3 O 4 catalyst, with addition of Co(NO 3 ) 2 $6H 2 O (0.05 M and 0.083 M), Ga(NO 3 ) 3 $xH 2 O (0.05 M and 0.017 M), and Al(OH)(C 2 H 3 O 2 ) 2 (0.05 M) precursors.The geometric surface area of the samples was ∼0.25 cm 2 .In order to conrm the identity of these compounds, they were further characterized with XPS, XRD and SEM.SEM was performed using a Hitachi SU7000.XRD was performed using a Bruker D8; all materials were exposed to Cu-Ka radiation (l = 0.15406 nm) and the data was collected with a point step of 0.02°.The XPS spectra in this study were analysed using the ThermoAvantage soware.All XPS spectra were calibrated based on the position of the C1s peak in each spectrum and the Co2p peak convolution was performed using the same soware as well.A Shirley background was used to determine the background of the XPS spectrum.
Electrochemical characterization.Electrochemical data were collected using a three-electrode system connected to an electrochemical workstation (Autolab PGSTAT302N) SCE and a carbon rod were used as reference and counter electrodes, respectively.1 M KOH without saturated O 2 was used as the electrolyte.Cyclic voltammetry (CV) measurements at 50 mV s −1 were performed for 3 cycles prior to recording linear scan voltammetry (LSV) at 5 mV s −1 for each sample.
Electrochemically active surface area (ECSA) determination.The ECSA of each catalyst was determined by measuring the electrochemical double-layer capacitance (C dl ) of each catalyst from the scan-rate dependence of the CV plot.Four different CV measurements were performed at a potential window between 0.72 to 0.82 V (vs. the reversible hydrogen electrode) on each catalyst at scan rates of 20,30, 40 and 50 mV s −1 respectively.The C dl was estimated at the average potential within this range by calculating the slope of the linear t at that point.A specic capacitance (C s ) of 40 mF cm −2 was used to calculate the ECSA of the catalyst using the following equation:

ECSA
The LSV of each plot was subsequently normalized by dividing the current density in each plot with the calculated ECSA of each catalyst.
Potential calibration and iR correction.The potentials versus SCE (E SCE ) were calibrated versus RHE (E RHE ) using the following equation: where R is the ideal gas constant (8.314J mol −1 K −1 ), T is the temperature (298 K), F is the Faraday constant (96 485 C mol −1 electrons), z is the number of electrons (1 mol) transferred, E 0 SCE is the standard potential of the SCE reference electrode that has been calibrated versus RHE.The pH is 14.
The potentials of the LSV were corrected by subtracting iR S , where i is the current measured at the corresponding potential and R S is the bulk and solution resistances.The resistance R S was calculated by tting electrochemical impedance data in the 0.1-1 Hz range with the Randles circuit model.The R S for each catalyst is given in Table S6 of the ESI.†

Fig. 1 (
Fig. 1 (A) Scheme describing the computational workflow to identify new catalysts for OER.In this workflow, Co 2.5 Ga 0.5 O 4 was discovered as a highly-performing spinel oxide catalyst.(B) Diagram describing the structure of a spinel oxide AB 2 O 4 , where A is tetrahedrally coordinated and B is octahedrally coordinated.

Fig. 2
Fig. 2 Heat map outlining the stability of compositions (A) AB 2 O 4 , (B) A 0.5 B 2.5 O 4 and (C) A 1.5 B 1.5 O 4 .Stable spinels in this study are considered to have formation energies below 0.03 eV per atom with respect to the convex hull.(D) Frequency plot of all thermodynamically stable materials based on the calculated theoretical overpotential.Materials deemed to have a low theoretical overpotential (less than 0.5 eV) marked in dark blue.(E) Frequency of each metallic element amongst thermodynamically stable spinels for which theoretical overpotentials for OER were also calculated.Elements making up low overpotential compounds (less than 0.5 eV) are also labelled in dark blue.
and the binding energies of *O, *OH and *OOH.We plot the relationship between the binding energies of *O and *OH and demonstrate that any increases in the Bader charge of adsorbate oxygen can be associated with an increase in the binding energies of the *O, *OH and *OOH intermediates.Furthermore, increases in the distance between the O2p bandcenter and the Fermi level are associated with decreases in the binding energies of *O (Fig. 3C) and *OH (Fig. S13 †).We then plot DE O -DE OH against both the oxygen bandcenter and the Bader charge of the adsorbate oxygen in the *O intermediate (Fig. 3D).We see that the relationship between DE O -DE OH and the Bader charge of oxygen resembles that of a volcano plot; DE O -DE OH increases when going from a Bader charge of −1.4e − to −0.5e − before decreasing against from −0.5e − to −0.1e − .DE O -DE OH decreases as the O-2p bandcenter moves further away from the Fermi level, a relationship that can be corroborated by Nørskov and coworkers.
-S5 †) of Co 3 O 4 and Co-Al oxide indicate that all the constituent elements of both materials are uniformly distributed.In contrast, the elements in the Co-Ga oxide samples show greater

Fig. 3 (
Fig. 3 (A) A graph created to illustrate the efficacy of the random forest model by comparing calculated binding energies of all the *O, *OH and *OOH intermediates with predicted binding energies on the test set.The model has an R 2 score of 0.83 for the training set and 0.78 for the test set.It also has a median absolute error of 0.16 eV for the training set and 0.21 eV for the test set.(B) A bar graph illustrating the relative importance of all the features used to predict the binding energies of *O, *OH and *OOH.The Bader charge of the adsorbed oxygen was plotted against the (C) binding energy of *O and (D) the difference in energy between O* and *OH, with the colorbar representing the bandcenter of the O2p band of adsorbate oxygen.

Fig. 4
Fig. 4 SEM images of all four catalysts with two different resolutions.The four images on the top show zoomed out images of catalysts presumed to be (a) Co 3 O 4 , (c) Co 1.5 Al 1.5 O 4 , (e) Co 2.5 Ga 0.5 O 4 and (g) Co 1.5 Ga 1.5 O 4 .The images at the bottom show zoomed in images of the same catalysts thought to be (b) Co 3 O 4 , (d) Co 1.5 Al 1.5 O 4 , (f) Co 2.5 Ga 0.5 O 4 and (h) Co 1.5 Ga 1.5 O 4 .

Fig. 5
Fig. 5 XRD and XPS spectra of all four compounds.(A) XRD spectra of all four catalysts.The spinel phase is present in all four catalysts, and in Co 2.5 Ga 0.5 O 4 , Co 1.5 Ga 1.5 O 4 and Co 1.5 Al 1.5 O 4 , the phases of Ga 2 O 3 and Al 2 O 3 are present as well.Co 2p XPS spectra of materials presume to be (B) Co 3 O 4 (C) Co 2.5 Ga 0.5 O 4 (D) Co 1.5 Ga 1.5 O 4 (E) Co 1.5 Al 1.5 O 4 .In order to determine the composition, the Co 2p 3/2 peak was deconvoluted.The peaks corresponding to Co 2+ and Co 3+ are labelled in each XPS diagram.
5 and 1 : 1 in the Co 3 O 4 , Co-Al, and the two Co-Ga catalysts respectively.This indicates that the exact structural formulae of the catalysts are Co 3 O 4 , Co 1.7 Al 1.3 O 4 , Co 2.5 Ga 0.5 O 4 and GaCo 2 O 4 .

Fig. 6 (
Fig. 6 (A) ECSA-adjusted iR-corrected LSV curves of all catalysts outlining the current density at different potentials (vs.RHE).Linear sweep voltammetry (LSV) measurements are performed in 1 M KOH solution at pH 14. (B) iR-corrected LSV curved of all catalysts.(C) Tafel Slopes of all catalysts shown on this semilog plot.Co 2.5 Ga 0.5 O 4 has the lowest Tafel slope of all the catalysts.(D) Benchmarking catalysts in this study with other state-of-the-art catalysts reported in the literature.Each experiment in our plot was done in 1 M KOH solution.A more detailed breakdown of each reference catalyst can be seen in TableS5.† Co 3 O 4 has the lowest overpotential of all the catalysts at 170 mV, followed by Co 1.7 Al 1.3 O 4 at 193 mV, Co 2.5 Ga 0.5 O 4 at 220 mV and nally GaCo 2 O 4 at 270 mV.However, while Co 2.5 Ga 0.5 O 4 has a higher overpotential than either Co 3 O 4 or Co 1.7 Al 1.3 O 4