Journal of Catalysis 412 (2022) 78–86 Contents lists available at ScienceDirect Journal of Catalysis journal homepage: www.elsevier.com/locate/jcat Metal coordination determines the catalytic activity of IrO2 nanoparticles for the oxygen evolution reaction Danilo González, Mariona Sodupe, Luis Rodríguez-Santiago, Xavier Solans-Monfort ⇑ Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain a r t i c l e i n f o a b s t r a c t Article history: H2 production through water electrolysis is a promising strategy for storing sunlight energy. For the oxy- Received 14 February 2022 gen evolution reaction, iridium oxide containing materials are state-of-the-art due to their stability in Revised 10 May 2022 acidic conditions. Moreover, precious metal content can be reduced by using small nanoparticles that Accepted 31 May 2022 show high catalytic activities. We performed DFT calculations on a 1.2 nm large IrO2 Wulff-like stoichio- Available online 4 June 2022 metric nanoparticle model (IrO2) with the aim of determining the factors controlling the catalytic activity of IrO2 nanoparticles. Results show that at reaction conditions tetra- and tricoordinated iridium centers Keywords: are not fully oxidized, the major species being IrO(OH) and IrO(OH)2, respectively. Although the com- Water electrolysis Oxygen Evolution Reaction puted overpotential show that all centers present relatively similar reactivities, low coordinated iridium Iridium oxide centers tend to be more active than the pentacoordinates sites of the well-defined facets. These low coor- Nanoparticles dination sites are likely more abundant on amorphous nanoparticles, which could be one of the factors DFT explaining the higher catalytic activity observed for non-crystalline materials. Ó 2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction lyst amount[21–33] Indeed, nanoparticles of about 1.5 – 2.0 nm have been shown to be within the most active species, The use of H2 as sunlight energy storing source is seen as one of [21,26,29,32,33] the origin of their high activity being associated the most promising strategies for obtaining clean energy.[1,2] This to three factors: i) the Ir(III)/Ir(IV) ratio; ii) the amount of surface can be achieved through the (photo)electrochemical splitting of hydroxylated species and iii) the degree of amorphization. water that converts water in oxygen at the anode through the oxy- Three main reaction mechanisms have been proposed for the gen evolution reaction (OER) and H2 at the cathode by means of the OER (Scheme 1):[10,11,34–37] i) the water nucleophilic attack hydrogen evolution reaction (HER).[3–9] The OER is a four-electron (WNA) in which the O-O bond of O2 is formed by the reaction of process that presents high overpotentials and slow kinetics and is an Ir=O species of the surface and a water molecule of the solution; considered the limiting half reaction for the application of water ii) the oxo-coupling mechanism (I2M) in which the O-O bond is electrolysis as source for sustainable energy conversion and stor- formed by the coupling of two Ir=O species on the surface and age.[10–12] Several OER promising catalysts have been reported iii) the lattice oxygen evolution reaction mechanism (LOER) that in the literature.[13,14] Within the large list of catalytically active involve oxygen atoms from the material in the O2 formation. Iso- species, IrO2 is considered to be the most efficient one in acidic tope labeling experiments showed that the LOER mechanism is media due to both its catalytic activity and stability in harsh reac- important in amorphous materials but does not contribute in crys- tion conditions.[11,12,15]. talline rutile-like materials.[11,38] Based on that, we decided not Iridium is an expensive precious metal and consequently lower- to consider the LOER pathway in the present contribution, since ing the iridium content on the catalyst is an essential issue for the our models (see below) are constructed from crystalline IrO2. practical implementation of the water splitting.[15–20] Several Indeed, for crystalline IrO2, the WNA pathway is the most accepted research groups have synthetized small IrO2(IrOX) nanoparticles mechanism, except in very flexible environments, where the I2M of different sizes and shapes showing intrinsic activities that are mechanism becomes competitive.[39,40] Remarkably, the local slightly higher than with crystalline films. Such improved intrinsic morphology of the material seems to tune the catalytic activity activity becomes remarkably high when normalizing by the cata- of IrO2 for OER.[41–44]. From a computational point of view, the catalytic activity of ⇑ Corresponding author at: Departament de Química, Universitat Autònoma de iridium-based materials for the OER has been addressed either Barcelona, Carrer dels Tillers s/n, 08193 Bellaterra, Spain. with slab models of the (110) surface[39,45–56] or small clusters E-mail address: xavier.solans@uab.cat (X. Solans-Monfort). up to 13 Ir centers.[57,58] Regarding the studies with slab models, https://doi.org/10.1016/j.jcat.2022.05.023 0021-9517/Ó 2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 Scheme 1. Reaction mechanisms for the oxygen evolution reaction proposed in the literature. most contributions limit the study to the thermodynamic cost of each PCET step. Only in a few cases, the energy barriers associated with some of the elementary steps have been computed (mostly the chemical ones) and the associated energy barriers tend to be low,[47,48,52,53,56] thus suggesting that the thermodynamics of the PCET is sufficient to get overpotential trends. Moreover, while most of these contributions agree that the most favorable mecha- nism on the (110) surface is the WNA, some discrepancies exist on the nature of the overpotential defining step. Indeed, three differ- ent steps have been proposed to be rate limiting depending on the level of theory used and the definition of the elementary steps of the reaction: i) the Ir-OH to Ir-O oxidation; ii) the electrochem- ical step associated with the Ir-OOH intermediate formation and iii) the O2 release. On the other hand, regarding the OER on cluster models, Auer and co-workers analyzed the effect of the applied potential and pH on the structure of two models containing either Fig. 1. (IrO2)33 nanoparticle model with the undercoordinated iridium centers 3 or 13 iridium centers.[58] Results for the larger cluster show that (coloured atoms) where OER activity has been studied. at acidic reaction conditions the surface presents mainly Ir-O and Ir-OH groups, the number of Ir-OH decreasing when increasing the applied potential and the pH. Moreover, the O-O bond forma- the tip site. This reconstruction implies the breaking of Ir-O bonds tion through the WNA mechanism is predicted to present a non- and the formation of Ir-Ir interactions as described previously.[61] negligible energy barrier that decreases while increasing the However, since IrOX nanoparticles are experimentally synthetized applied potential. in water solution, it is very likely that the reconstruction observed Remarkably, nanoparticles, even being crystalline, present par- in vacuum does not represent the experimental structure. With the ticular sites such as tips, edges and corners whose representation aim of better describing the aqueous environment in which IrO2 with extended models may not be appropriate, even if other sur- nanoparticles are usually synthetized, we constructed a model faces other than the (1 1 0) are taken into account.[48] In addition, including a water molecule in each surface iridium vacant site. clusters may not properly represent the electronic structure of the The full optimization of this model lead to some surface relaxation commonly used nanoparticle sizes. In this contribution, we address but the reconstruction around the tip was prevented. After full the OER reaction on all potentially active sites present in the stoi- optimization, the added water molecules were removed, and all chiometric Wulff-like (IrO2)33 nanoparticle model. Results show nanoparticle atoms, except the singly coordinated oxygen atoms, that although the reactivity of all considered sites is relatively sim- were kept fix in all subsequent calculations. ilar, the oxygen evolution reaction on tri- and tetracoordinated At reaction conditions, it is expected that the nanoparticle iridium centers tends to present lower overpotentials than on pen- would be highly oxidized, thus presenting several Ir-OH and Ir=O tacoordinated sites. This is associated with the presence of Ir(OH)n groups at the vacant sites of the surface iridium atoms. Inclusion species at reaction conditions. These low coordinated sites are of these groups in each iridium center in combination with the spin likely more abundant in amorphous materials, thereby providing polarized formalism makes SCF convergence more time consum- a plausible explanation on the high catalytic activity of non- ing. Moreover, determining the surface termination at 1.5 V poten- crystalline IrO2 and its dependence with the presence of Ir(OH)n tials including explicitly the 38 water molecules and analyzing all species on the surface. potential PCET associated to these water molecules is prohibitive. In this context, we decided to use a model where we only consid- ered the oxidation of water molecules located at adjacent iridium 2. Computational details centers, with the final goal of having two Ir-O centers to carry on the reaction. One of the two Ir-O centers acts as the active site A Wulff-like stoichiometric nanoparticle model of 1.2 nm size and the neighbor one either actively participate in the oxo- was constructed with the BCN-M computational tool (Fig. 1).[59] coupling mechanism or it assists the chemical water nucleophilic The surface energies provided to construct the model are 0.150, attack step. In addition, we constructed a second model in which 0.158, 0.185 and 0.206 eV A-2 for the (1 1 0), (0 1 1), (1 0 0) and all surface iridium centers are oxidized to the most stable structure (0 0 1) crystallographic surfaces,[60] respectively. The final model at 1.5 V (according to the results obtained with the first model). contains 33 IrO2 units, its size is close to some of the experimen- This model was used to analyze how the imposed constraints asso- tally reported active nanoparticles,[21,26,29,32,33] and it is com- ciated to the partial optimization and the degree of oxidation of the putationally affordable. Remarkably, the full optimization of nanoparticle influences the reactivity of the A5ax-A5ax site. That is, Ir33O66 in vacuum implied an important reconstruction around with the fully oxidized model we computed the energetics of the 79 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 WNA mechanism assuming two approaches. In the first one the Therefore, the minimum required overpotential (gOER) is com- nanoparticle was kept fixed as in the partial oxidized model and puted as the DG° difference between the highest in Gibbs energy in the second one the nanoparticle was fully optimized. That is, PCET process and 1.23 V and it is used to identify the most active comparison between the three models, partial oxidation with con- sites. Note that the reaction energy is imposed to be 4.92 eV due strained optimization, full oxidation with constrained optimization to the way the G° of O2 is computed (see below). and full optimization of the oxidized model allows distinguishing the effects arising from the degree of oxidation of the nanoparticle and the effects of the constrained optimization. Results are shown gOER = max[DG (PCET)]/e —1.23 ð1Þ in Table S1 of the supplementary materials, where we also added Values reported along the text are based on Gibbs energies at the data for the fully optimized (1 1 0) surface considering two 1 atm and T = 273.15 K assuming the computational standard hydro- supercell models, one fully oxidized (1 1 0)2x1 and one half oxidized gen electrode for the proton coupled electron transfer steps (PCET). (1 1 0)2x2. Results show that the degree of oxidation of the neighbor [45,72] Thermal contributions were computed through different iridium centers has very little influence on the overall overpoten- approximations: i) For the OER intermediates, the thermal correc- tial. However, the effect of relaxing the geometry constraints is tions were obtained by considering the contributions to the vibra- more pronounced, the overpotential decreasing by about 0.20– tional modes of the adsorbed species and the Ir-Oads stretching 0.25 V. All PCET steps become easier and close to the results with mode. This is 6, 3, 1, 5 and 3 vibrational modes per adsorbed *H2O, the surface models. Indeed, according to the results provided *OH, *O, *OOH and *OO species, respectively; ii) The entropy term below, it is very likely that slab models of the sites existing on of H2O and H2 were taken from tabulated values for liquid H2O the surface may lead to similar results to those obtained for the and H2 gas, respectively; and iii) In analogy with previous contribu- realistic nanoparticle model. The decrease of the PCET reaction tions,[39,45,53] the O2 Gibbs energy (in eV) was computed as energies is associated to a relative stabilization of the most oxi- 4.92 + 2 G0H2O – 2G0H2 where 4.92 is the experimental reaction energy dized species with respect to the less oxidized ones that may be and G0H2O and 2G0H2 are the computed Gibbs energies for water and related with the fact that the restricted structure corresponds to hydrogen, respectively. that in which water is adsorbed on each vacant site, i.e. the struc- ture of the less oxidized species. A similar effect should be expected with the other sites, but due to the limitations of the 3. Results and discusion model, small difference between sites should be taken with care. All calculations were performed with the VASP code,[62,63] We have studied the catalytic activity for the oxygen evolution using the spin polarized formalism, the PBE density functional reaction of the six different undercoordinated iridium centers pre- [64] and the Grimme’s D2 empirical correction to account for dis- sent in the (IrO2)33 nanoparticle (Fig. 1): A5ax (brown), B5ax (dark persion forces.[65] PBE functional has been widely used in model- blue), C5eq (orange), D4ax/eq (green), E4eq/eq (purple) and F3ax/eq/eq ing the OER with iridium based materials including finite clusters (light blue). We have only considered the I2M and WNA proposed [39,40,48,52,53,55,58] and according to Goddard and co-workers mechanisms since our nanoparticle model is based on the crys- represents the metallic nonmagnetic electronic structure of the talline material, for which the LOER reaction mechanism bulk in better agreement with experiments than hybrid function- (Scheme 1c) has been reported to be non-operative. The six centers als or the PBE-U approach.[66] Unfortunately, it is also well differ on the coordination number and environment. Taking into known that GGA functionals tends to delocalize the electron den- account that iridium in rutile shows two short Ir-O axial distances sity of unpaired electrons[67] as it is the case of the Ir-O oxo spe- and four long equatorial ones,[60] the six centers can be divided cies. To ensure that the electron delocalization has no effect on in five different types: i) pentacoordinated atoms with an axial the O-O bond formation step, we performed selected calculations vacant site (A5ax and B5ax); ii) pentacoordinated sites with an equa- with the PBE + U approach. As shown in Table S2, inclusion of the torial vacancy (C5eq); iii) tetracoordinated centers with one axial Hubbard term marginally increases the spin density on the oxyl and one equatorial vacant site (D4ax/eq); iv) tetracoordinated sites species of the IrO2 (1 1 0) main surface and this makes the I2M with two equatorial vacant sites (E4eq/eq); and v) iridium centers mechanism slightly less unfavorable. However, the WNA mecha- presenting only three bonds with lattice oxygens (F3ax/eq/eq). While nism is always preferred by at least 1.0 eV. PAW pseudopotentials other sites exist in larger Wulff-like nanoparticle models, they can were used to describe the ionic cores.[68,69] The valence elec- all be classified in one of these five types of centers.[59,61] Indeed, trons were represented with a plane wave basis with a kinetic we found that the water adsorption in the existing sites of (IrO2)33 energy cutoff of 500 eV. Moreover, single point calculations at and (IrO2)115 nanoparticle models is mainly determined by the nat- the final optimized geometries were performed to include solvent ure of the vacant site and nanoparticle morphology and not by the effects (water, e = 80.0) with VASPsol package.[70] It has been nanoparticle size.[61] Therefore, the results reported here should recently shown[71] that the inclusion of solvent through implicit be indicative of how the different sites of the nanoparticle react. solvent models does not necessarily leads to a better description The results and discussion section is organized in two parts. We of adsorbed species related to electrocatalysis such as OH or first focus on the oxidation of the surface iridium centers in the OOH. Therefore, we decided to compare the results with and presence of water as a function of the applied potential[61] and without implicit solvation (Table S3 of the supplementary mate- afterwards, we explore the OER reaction starting from the most rials). Inclusion of solvent through an implicit model has a non- stable species at 1.5 V, through either the I2M and WNA mecha- negligible effect on processes involving the adsorption of water, nisms. The choice for 1.5 V is based on the fact that this potential but it has small effect on the energetics of the proton coupled is close to the optimal value for performing the OER reaction with electron transfer steps (PCET). Consequently, the main trends different Ir-based materials.[26,53,73]. are preserved both with and without the inclusion of implicit sol- vent. The nanoparticle model was placed in a cubic box of a 30 Å 3.1. Oxidation of surface Ir-(H2O)X edge and the calculations were performed at gamma point due to the model dimensionality. The oxidation of the six iridium centers of (IrO2)33 model (Fig. 1) Since the global process has been experimentally determined to in the presence of water as a function of the applied potential has be endothermic by 4.92 eV and it implies four PCET, an ideal elec- been analyzed with the aim of establishing the most stable struc- trocatalyst should catalyze these reactions at potentials of 1.23 V. ture at potentials around 1.5 V. Since we considered both the 80 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 I2M and WNA mechanisms and the former path requires the par- iii) the blue series accounts for structures with at least one OH ticipation of two Ir-O groups, the active site oxidation was modeled group as the most oxidized species at each iridium site; iv) orange by considering either one single center that could lead to the for- color describes species with one Ir-O group in one of the metal cen- mation of two oxo groups (D4ax/eq, E4eq/eq or F3ax/eq/eq) or two vici- ters and; v) the pink colors indicate the presence of two Ir-O spe- nal iridium centers (A5ax-A5ax, D4ax/eq-B5ax, D4ax/eq-Ceq, F3ax/eq/eq- cies (Scheme 2). According to our results on surface models,[48] A5axor F3ax/eq/eq-D4ax/eq). The oxidation process takes place through the oxygen evolution catalysis requires the formation of species several proton coupled electron transfer steps (PCET) and trans- represented in pink in Fig. 2, at around 1.5 V potentials. forms the adsorbed water molecules to either Ir-OH or Ir-O species. Results show that the required potential for the H2O to OH oxi- The initial structure arises from adding one water molecule per dation is sensitive to the neighbor metal that is also being oxidized. vacant site in its preferred form (molecular or dissociated) as dis- The computed potentials for the H2O to OH oxidation for one speci- cussed in our previous contribution.[61] Then, we considered the fic center can vary up to 0.2 eV and this appears to be related with PCET processes until all metal centers present Ir-O species only. the formation of hydrogen bonds. Centers able to establish strong We explored all possibilities for the individual tri- and tetracoordi- hydrogen bonds with neighbor adsorbed species tend to present nated sites and for systems involving two centers the combinations higher potentials for the H2O to OH oxidation. However, this differ- of the most stable structures of each site only. The optimized struc- ence will likely be smaller if explicit solvent water molecules are tures of the most stable isomers involved in the surface oxidation added in the simulation due to the formation of hydrogen bonds (Figures S1 to S8), the relative Gibbs energies between the different between the adsorbed species and the solvent. Moreover, the isomers of the intermediates involved in the oxidation process and H2O to OH transition always takes place at low potentials and thus, the reaction Gibbs energies for all steps (Tables S4 to S18) are this oxidation step does not appear to play a key role in the OER reported in the supplementary material. with iridium-based metals. Relative stabilities of the different species as a function of the The Ir-OH to Ir-O oxidation is much less sensitive to the neigh- applied potential were computed following equation (2). bor iridium center that it is being oxidized and this is related with the fact that Ir-OH species do not form strong hydrogen bonding X(H2 O)-Y(H2 O) ! X(L1 )-Y(L2 ) + ne + nHþ with species adsorbed in the vicinity. Consequently, the values DG0 = G0 L1=L2 —G0 H2O=H2O —nU. ð2Þ for each iridium center are essentially not affected by the nature of the other metal that is being oxidized (less than 0.08 eV). The where X(L1)-Y(L2) stands for the oxidized species arising from X formation of the first Ir-O species at tri- and tetracoordinated irid- (H2O)-Y(H2O), n is the number of PCET processes required to go ium centers occurs at potentials between 1.18 and 1.50 V and the from X(H2O)-Y(H2O) to X(L1)-Y(L2) and U is the applied potential. computed values tend to be lower than the required potentials for Fig. 2 summarizes the results. The color labeling indicates the most oxidizing pentacoordinated centers (between 1.45 and 1.54 V). oxidized species in each center: i) grey indicates that all vacant Despite these small differences may be sensitive to the modeling sites are occupied by molecular or dissociated water molecules; approach, the general trend of the required potential for the forma- ii) the green color indicates that at least one OH group as most tion of the first Ir-O species follows the order E4eq/eq < F3ax/eq/eq <- oxidized species is present in one of the two considered centers; D4ax/eq < A5ax  B5ax < C5eq and this agrees with the observed trends for extended surface models.[48] Remarkably, the Ir-OH to Ir-O oxidation of the equatorial and axial sites of D4ax/eq with the other vacancy containing a Ir-OH group occurs at very similar potentials (difference of less than 0.03 V) and thus, both D4eq/ax(O/OH) and D4eq/ax(OH/O) structures should coexist. Full oxidation of the tri- and tetracoordinated iridium centers (D4eq/ax, F3ax/eq/eq) requires high potentials (over 1.6 V). This sug- gests that at reaction conditions these centers will not be fully oxi- dized and they will remain as Ir(O)(OH)x (x = 1 or 2 for D4eq/ax and F3ax/eq/eq sites, respectively). E4eq/eq located at the tip of the nanoparticle is an exception. According to calculations, the full oxi- dation of this center occurs at 1.29 V and thus, it will present two Ir-O groups in one single iridium atom at reaction conditions as already found for the analogous centers of the (0 0 1) surface and when using a slab model of the tip site (E4eq/eq-slab, Figure S9). [48] Overall, at 1.5 V the most stable structures for all considered systems are: D4ax/eq(OH/O), E4eq/eq(O/O), F3ax/eq/eq(O/OH/OH), A5ax(O)-A5ax(O), D4ax/eq(OH/O)-B5ax(O), D4ax/eq(OH/O)-C5eq(O), F3ax/eq/eq(O/OH/OH)-A5ax(O) or F3ax/eq/eq(O/OH/OH)-D4ax/eq(OH/O). 3.2. Oxygen evolution reaction catalytic activity In a second step, we studied the OER reaction starting from the most stable structure at potentials around 1.5 V. We performed the study in the pairs of sites having two Ir-O groups on the surface. We also considered E4eq/eq(O/O) that already has two oxyl species at 1.5 V. Despite being slightly less stable (less than 0.06 eV), we considered D4ax/eq(O/OH)-B5ax(O) and D4ax/eq(O/OH)-C5eq(O) instead of D4ax/eq(OH/O)-B5ax(O) and D4ax/eq(OH/O)-C5eq(O), since Fig. 2. Phase diagram of the different sites of (IrO2)33 at pH = 0, T = 298 K and the OER reaction can only take place when the oxyl group is in P = 1 atm as a function of the applied potential. Color labeling according to the axial site (see Supplementary material). We considered the Scheme 2. I2M and the WNA mechanisms shown in Scheme 3. 81 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 Scheme 2. Species considered for the H2O oxidation at the different nanoparticle sites. Color labeling indicate the most oxidized species at each Ir center (See text for details). Scheme 3. Reaction mechanisms considered in this work. Noteworthily, and similarly to the mechanistic proposal made regardless how the WNA process takes place (concerted or step- by other authors[47,53,58] for the OER on the IrO2 (1 1 0) surface, wise), the main intermediate involved in the O2 formation is Ir- we decoupled the water attack in two steps: i) the splitting of OO/Ir-OH rather than Ir-OOH/Ir-O. Comparison of the two approx- water over two Ir-O groups (1 to 2WNA process in Scheme 3) and imations, assuming that the two routes involve the Ir-OO/Ir-OH ii) a PCET process from the Ir-OOH/Ir-OH species (2WNA to 3WNA). species, suggests that the description of the WNA mechanism does Nevertheless, and for comparison with the results reported in the not modify the general trends. Finally, since the two iridium cen- literature,[56] the global reaction in which the water attack is cou- ters are not equivalent in most of the pair of sites considered, there pled with the PCET has also been considered (see Table S18 of the are two different WNA attack mechanisms depending on the irid- supplementary material). Results show that the Ir-OOH/Ir-O spe- ium presenting the Ir-OOH species. We considered the two possi- cies usually considered in the global reaction (30 WNA in bilities (WNAX and WNAY in Scheme S2 of the supplementary Scheme S1 of the supplementary material) is less stable than the material), the associated Gibbs energies for each pathway are Ir-OO/Ir-OH intermediate involved in the decopled mechanism reported in Table S18 and the optimized structures of all species reported here, in agreement with previous results on surfaces. considered in the OER processes are shown in Figures S10-S25 of [48,53] Moreover, similar results are obtained for RuO2 surfaces: the supplementary material. The values reported in the text corre- the Ru-OO/Ru-OH intermediate is more stable than the Ru-OOH/ spond to the most favorable pathway for each pair of sites. The DG Ru-O one and the geometrical features of the former matches with values for the individual steps are reported in Table 1 and Fig. 3 the experimental observations.[74] Therefore, it is very likely that shows the global energy profile for the WNA mechanism. 82 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 The I2M mechanism starts with the coupling of two Ir-O species ing that the oxo-coupling process is not only unfavorable, but also on the surface (1 to 2I2M in Scheme 3 and Table 1). The process is in slow. In any case, after O-O bond formation, O2 release by adsorp- all cases endergonic, and the reaction Gibbs energies varies from tion of two water molecules is thermodynamically favorable (DG° 0.63 to 1.26 eV. These values tend to be slightly higher than those between 0.92 and 1.33 eV) in all sites except at E4eq/eq, where obtained with extended models,[48] although this may originate the reaction occurs in one single metal center. The catalytic cycle from the restrictions imposed for the optimizations with the is closed through the oxidation of the adsorbed water species nanoparticle models (See computational details). The energy cost through processes that are equivalent to those described before for the O-O bond formation correlates reasonably well with the and require potentials between 1.29 V and 1.58 V to recover the OO distance in the initial structure (Table 1), since the more sep- two Ir-O species. Overall, the I2M mechanism presents relatively arated the oxygens initially are, the more unfavorable the process low overpotentials (between 0.06 and 0.35 V). Indeed, these values is. Moreover, we located the transition states for the oxo-coupling are similar or even lower than those obtained for the WNA path- step involving A5ax-A5ax, E4eq/eq and F3ax/eq/eq-D4eq/ax sites as way (see below). However, the feasibility of the I2M mechanism selected cases presenting very different environments. The com- is highly controlled by the oxo-coupling process, an endergonic puted energy barriers are in all cases large (>1.16 eV), thus suggest- chemical step with high energy barriers. Table 1 Reaction Gibbs Energies (in eV) for the chemical and electrochemical steps of the I2M and WNA reaction mechanisms. Values in parenthesis correspond to the Gibbs energy barriers. For those sites including Ir centers of different nature, only the most favorable pathway is reported. Reaction conditions are 0 V potential versus SHE, pH = 0, T = 298 K and P = 1 atm. Reaction A5ax-A5ax E4eq/eq D4ax/eq-B5ax D4ax/eq-C5eq F3ax/eq/eq-A5ax F3ax/eq/eq-D4ax/eq I2M 1? 2I2M 0.81 (1.16) 0.80 (1.78) 0.63 1.02 1.17 1.26 (1.41) 2I2M + 2 H2O ? 3I2M + O2 1.26 0.33 0.92 0.89 0.36 1.33 3I2M ? 4I2M + H+ + e- 0.96 0.57 0.93 0.55 0.55 1.38 4I2M ? 5I2M + H+ + e- 1.44 0.75 1.35 1.22 0.74 0.82 5I2M ? 6 + H+ + e- 1.45 1.18 1.44 1.49 1.28 1.31 6? 1 + H+ + e- 1.52 1.29 1.49 1.58 1.53 1.48 WNA 1 + H2O ? 2WNA 0.53 (0.13) 0.47 (0.84) 0.17 0.27 0.09 0.03 (0.43) 2WNA ? 3WNA + H+ + e- 1.31 0.88 1.19 1.01 0.93 1.16 3WNA ? 4WNA + H+ + e- 1.62 1.25 1.54 1.53 1.47 1.41 4WNA + H2O ? 5WNA + O2 0.40 0.47 0.11 0.44 0.26 0.09 5WNA ? 6 + H+ + e- 1.40 0.56 0.99 0.68 0.83 0.99 6 ? 1 + H+ + e- 1.52 1.29 1.48 1.52 1.34 1.48 Fig. 3. Calculated Gibbs energy profile for the most favorable OER catalytic cycle through the WNA mechanism at: a) the reference (1 1 0) sites of the nanoparticle; b) the tip site of the nanoparticle; c) sites including one penta and one tetracoordinated metal center and d) sites including one tricoordinated site. The blue lines correspond to chemical steps, the black lines correspond to the energetics of the PCET at an applied potential of 0 V and the red values are for an applied potential of 1.5 V. Labelling taken from Fig. 1 and Scheme 3. 83 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 The WNA mechanism starts with the water attack to the two Ir- reports the required overpotentials (gOER) and outlines the poten- O of the surface, which leads to the formation of Ir-OOH and Ir-OH tial determining step. The Figure is divided in four groups. A5ax- species (1 to 2WNA in Scheme 3 and Table1). The reaction Gibbs A5ax site, which corresponds to two iridium centers on the most energies for this step varies from 0.53 to 0.47 eV, but in most stable (1 1 0) surface and thus, its reactivity is taken as the refer- cases the process is essentially isergonic. In general, the reaction ence value (Fig. 3a). The energetics associated with E4eq/eq, which Gibbs energies for the WNA step at the different sites of the is the sole site involving only one iridium center, is shown in nanoparticle and the values obtained with extended surface mod- Fig. 3b and those of all other sites, organized according to the coor- els are very similar and they all suggest that the process is easy. dination number of the least saturated site, are given in Fig. 3c for The sole exception is the reaction at E4eq/eq. In this case, the process tetracoordinated metals and in Fig. 3d for tricoordinated species. is unfavorable and the computed value (DGo = 0.47 eV) is larger The computed overpotentials vary from 0.06 V to 0.39 V and they than that of the closely related (0 0 1) surface (DGo = 0.25 eV). are mostly between 0.24 and 0.39 V. The differences are less than [48] Moreover, the computed Gibbs energy barriers for A5ax-A5ax, 0.15 eV, and fall within the DFT accuracy. However, since errors E4eq/eq and F3ax/eq/eq-D4ax/eq indicate that the process should be fast due to the level of theory are likely to be similar for all sites, trends (DGà less than 0.5 eV) with the exception for the E4eq/eq tip site for are expected to be properly described. Indeed, the computed over- which the transition state is 0.84 eV above the bisoxyl species + H2- potential for the reference A5ax-A5ax site corresponds to the highest O. Remarkably, even for this E4eq/eq site, the comparison between value (0.39 V), suggesting that most of the sites of the nanoparticle the WNA reaction Gibbs energy (0.47 eV) and the associated Gibbs would be more reactive than the sites on the most stable (1 1 0) energy barrier (0.84 eV) with those computed for the O-O bond for- surface in agreement with the higher activity found for nanoparti- mation in the I2M mechanism (0.80 eV and 1.78 eV, respectively) cles.[22,26,27] The case of E4eq/eq is special: it presents the lowest suggests that the WNA pathway is the preferred chemical step. overpotential (only 0.06 V) but the corresponding chemical pro- Analysis of the spin densities over the oxygen atoms of Ir-O species cesses are unfavorable and present high energy barriers. This (Table 2) shows that they have a marked oxyl radical character and seems to be detrimental in the final catalytic activity. Indeed, the thus, according to electronic density distribution, the water attack experimentally reported Tafel plots show almost not activity at should be viewed as homolytic instead of nucleophilic. This is con- potentials below 1.4 V,[26,75,76] which would be consistent with sistent with the fact that E4eq/eq is the center presenting the lowest the fact that this site is not active. Noteworthily, the overpotentials spin densities over the oxygen atoms. computed for sites including one tetracoordinated metal as the After the O-O bond formation through the WNA mechanism, least saturated center (D4ax/eq-Bax and D4ax/eq-Ceq, Fig. 3c) are larger two PCET are required before O2 is released (2WNA ? 3WNA + 1H+ than those computed for sites including initially tricoordinated + 1e- and 3WNA ? 4WNA + 1H+ + 1e- in Scheme 3 and Table 1). iridium atoms (F3ax/eq/eq-A5ax and F3ax/eq/eq-D4ax/eq, Fig. 3d) by For all centers, the 3WNA to 4WNA PCET is more challenging than 0.05 – 0.07 eV. The two groups show overpotentials that are lower the 2WNA to 3WNA one, the computed values being similar to those to those of the A5ax-A5ax site of the (1 1 0) facet. Despite the differ- of the Ir-OH to Ir-O oxidation in the vicinity of a Ir-O group ences are small and fall within the accuracy of the computational (6 ? 1 + H+ + 1 e- step). Thus, the potential determining step approach, this is indicative that initially low coordinated sites involves in all cases a Ir-OH to Ir-O oxidation and depending on would be more active. Unfortunately, according to Alexandrov the site, this step is either the 3WNA to 4WNA oxidation or the 6 to and co-workers,[46] these Ir-OH containing species are likely less 1 process (Scheme 3 and Tables 1 and S17). The resulting overpo- stable and would readily dissolve in acidic conditions. This would tentials are similar to those of the I2M pathway and range from correlate with the number of surface OH groups at reaction condi- 0.06 to 0.39 V. tions. Indeed, at working conditions it is expected that all vacant Comparison between the O-O bond formation steps of the I2M sites will be saturated with adsorbed species and according to and WNA mechanisms (Table 1) shows that the water attack is the required overpotentials for the Iridium oxidation, these will always more favorable than the oxo-coupling process. Moreover, mostly be OH groups. Therefore, the tricoordinated metals would the Gibbs energy barriers computed for some selected sites present two OH groups on the surface and the tetracoordinated (A5ax-A5ax, E4eq/eq and F3ax/eq/eq-D4ax/eq) show that the barriers are species one single OH. Overall, our results seem to suggest that low for the WNA mechanism and significantly higher for the I2M the highest activity of IrO2 nanoparticles when compared with pathway. Therefore, assuming a similar value for the other centers crystalline materials and particularly the (1 1 0) facet arise both of the nanoparticle, the WNA seems to be the preferred route in all by the increase of active sites per mass of catalyst and the higher sites. Consequently, the analysis of the overpotentials required to activity of low coordinated sites. Remarkably, these low coordi- perform the reaction in each site is made by comparing the values nated centers are likely more abundant on non-crystalline materi- of the WNA mechanism. als, thus suggesting that this could be one of the key factors for the Fig. 3 summarizes the Gibbs energies of the most favorable highest activity of amorphous materials and, particularly, nanopar- WNA pathway for the 6 different sites of the nanoparticle. It also ticles.[24,26,30]. Table 2 Ir-OL distances (in Å) and spin moment for the most stable structures at working potentials. Species Ir1-O11 Ir1-O12 Ir1-O13 Ir2-O21 Ir2-O22 O1-O2a SO1b SO2c SIr1 SIr2 AO-AO 1.805 1.808 3.067 0.626 0.640 0.533 0.527 EO/O 1.772 1.770 – 2.900 0.233/0.255 0.406 DOH/O 1.941 1.794 N/A 0.582 0.565 DO/OH-BO 1.772 1.934 1.796 2. 889 0.441 0.471 0.240 0.281 DO/OH-CO 1.939 1.794 1.950 1.804 4.434 0.587 0.698 0.582 0.532 FO/OH/OH 1.778 1.917 1.930 1.804 0.314 0.314 FO/OH/OH-AO- 1.776 1.916 1.931 1.787 3.249 0.412 0.646 0.389 0.538 FO/OH/OH-DOH/O 1.768 1.927 1.921 1.804 3.234 0.692 0.242 0.700 0.239 a Distance between oxygens involved in the oxo-coupling step. b Oxygen atom of the Ir-O species at Ir1. c Oxygen atom of the Ir-O species at Ir2. 84 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 4. Conclusions [16] J. Cheng, J. Yang, S. Kitano, G. Juhasz, M. Higashi, M. Sadakiyo, K. Kato, S. Yoshioka, T. Sugiyama, M. Yamauchi, N. Nakashima, ACS Catal. 9 (2019) 6974– 6986. Calculations on the oxygen evolution reaction catalytic activity [17] H.N. Nong, T. Reier, H.-S. Oh, M. Gliech, P. Paciok, T.H.T. Vu, D. Teschner, M. of several iridium sites present on the (IrO2)33 Wulff-like stoichio- Heggen, V. Petkov, R. Schlögl, T. Jones, P. Strasser, Nat. Catal. 1 (2018) 841–851. [18] D. Böhm, M. Beetz, M. Schuster, K. Peters, A.G. Hufnagel, M. Döblinger, B. metric nanoparticle model show that pentacoordinated sites Böller, T. Bein, D. Fattakhova-Rohlfing, Adv. Funct. Mater. 30 (2020) 1906670. located at the (1 1 0) surface are slightly less active when compared [19] Q. Wang, Z. Zhang, C. Cai, M. Wang, Z.L. Zhao, M. Li, X. Huang, S. Han, H. Zhou, to other centers of the nanoparticle. In contrast, tri- and tetracood- Z. Feng, L. Li, J. Li, H. Xu, J.S. Francisco, M. Gu, J. Am. Chem. Soc. 143 (2021) 13605–13615. inated sites mainly located at tip, corner or edges (in larger mod- [20] M. Ledendecker, S. Geiger, K. Hengge, J. Lim, S. Cherevko, A.M. Mingers, D. els) of the nanoparticle present lower overpotentials. Despite the Göhl, G.V. Fortunato, D. Jalalpoor, F. Schüth, C. Scheu, K.J.J. Mayrhofer, Nano computed differences are small and may fall within the modeling Res. 12 (2019) 2275–2280. approach accuracy, present work suggests that using nanoparticles [21] Y. Zhao, E.A. Hernandez-Pagan, N.M. Vargas-Barbosa, J.L. Dysart, T.E. Mallouk, J. Phys. Chem. Lett. 2 (2011) 402–406. instead of films has two main advantages. On one side, the use of [22] T. Reier, M. Oezaslan, P. Strasser, ACS Catal. 2 (2012) 1765–1772. nanoparticles increases the number of active sites per mass of cat- [23] H. Jiang, S. Ge, Y. Zhang, M. Dong, S. Wu, M. Wu, J. Zhang, R. Ge, Z. Guo, J. Phys. alyst and, on the other side, nanoparticles present low coordinated Chem. C 124 (2020) 2–8. [24] P. Jovanovič, N. Hodnik, F. Ruiz-Zepeda, I. Arčon, B. Jozinović, M. Zorko, M. Bele, sites at corners and edges that would be more active than sites M. Šala, V.S. Šelih, S. Hočevar, M. Gaberšček, J. Am. Chem. Soc. 139 (2017) found in crystalline facets. These low coordinated sites are not fully 12837–12846. oxidized at potentials around 1.5 V and the most stable structure at [25] R.M. Kluge, R.W. Haid, A.S. Bandarenka, J. Catal. 396 (2021) 14–22. [26] D.F. Abbott, D. Lebedev, K. Waltar, M. Povia, M. Nachtegaal, E. Fabbri, C. reaction conditions is Ir(O)(OH)X (X = 1 or 2 for tetra- or tricoordi- Copéret, T.J. Schmidt, Chem. Mater. 28 (2016) 6591–6604. nated centers, respectively). These low coordinated sites are likely [27] Y. Lee, J. Suntivich, K.J. May, E.E. Perry, Y. Shao-Horn, J. Phys. Chem. Lett. 3 more abundant in amorphous materials, which could be one of the (2012) 399–404. [28] P. Lettenmeier, L. Wang, U. Golla-Schindler, P. Gazdzicki, N.A. Cañas, M. Handl, factors contributing to the higher catalytic activity of amorphous R. Hiesgen, S.S. Hosseiny, A.S. Gago, K.A. Friedrich, Angew. Chemie Int. Ed. 55 materials. (2016) 742–746. [29] L. Fu, X. Zeng, C. Huang, P. Cai, G. Cheng, W. Luo, Inorg. Chem. Front. 5 (2018) 1121–1125. Declaration of Competing Interest [30] E. Oakton, D. Lebedev, M. Povia, D.F. Abbott, E. Fabbri, A. Fedorov, M. Nachtegaal, C. Copéret, T.J. Schmidt, ACS Catal. 7 (2017) 2346–2352. The authors declare that they have no known competing finan- [31] R. Badam, M. Hara, H.-H. Huang, M. Yoshimura, Int. J. Hydrogen Energy 43 (2018) 18095–18104. cial interests or personal relationships that could have appeared [32] F. Bizzotto, J. Quinson, A. Zana, J.J.K. Kirkensgaard, A. Dworzak, M. Oezaslan, M. to influence the work reported in this paper. Arenz, Catal, Sci. Technol. 9 (2019) 6345–6356. [33] D. Lebedev, C. Copéret, A.C.S. Appl, Energy Mater. 2 (2019) 196–200. [34] H. Dau, C. Limberg, T. Reier, M. Risch, S. Roggan, P. Strasser, ChemCatChem 2 Acknowledgement (2010) 724–761. [35] E. Fabbri, T.J. Schmidt, ACS Catal. 8 (2018) 9765–9774. [36] S. Divanis, T. Kutlusoy, I.M. Ingmer Boye, I.C. Man, J. Rossmeisl, Chem. Sci. 11 This work was supported by the Ministerio de Ciencia e Inno- (2020) 2943–2950. vación, SpainGeneralitat de Catalunya, Spain, MICINN [CTQ2017- [37] O. Piqué, F. Illas, F. Calle-Vallejo, Phys. Chem. Chem. Phys. 22 (2020) 6797– 89132-P and PID2020-112715GB-I00], the Generalitat de Catalu- 6803. nya [2017SGR1323] and Red Española de Supercomputación [38] K.A. Stoerzinger, O. Diaz-Morales, M. Kolb, R.R. Rao, R. Frydendal, L. Qiao, X.R. Wang, N.B. Halck, J. Rossmeisl, H.A. Hansen, T. Vegge, I.E.L. Stephens, M.T.M. [QCM-2019-1-0043, QS-2019-2-0031]. Koper, Y. Shao-Horn, ACS Energy Lett. 2 (2017) 876–881. [39] M. García-Melchor, L. Vilella, N. López, A. Vojvodic, ChemCatChem 8 (2016) 1792–1798. Appendix A. Supplementary data [40] K. Klyukin, A. Zagalskaya, V. Alexandrov, J. Phys. Chem. C 122 (2018) 29350– 29358. Supplementary data to this article can be found online at [41] L.C. Seitz, C.F. Dickens, K. Nishio, Y. Hikita, J. Montoya, A. Doyle, C. Kirk, A. Vojvodic, H.Y. Hwang, J.K. Norskov, T.F. Jaramillo, Science 353 (2016) 1011– https://doi.org/10.1016/j.jcat.2022.05.023. 1014. [42] G. Buvat, M.J. Eslamibidgoli, A.H. Youssef, S. Garbarino, A. Ruediger, M. Eikerling, D. Guay, ACS Catal. 10 (2020) 806–817. References [43] J. Gao, C.-Q. Xu, S.-F. Hung, W. Liu, W. Cai, Z. Zeng, C. Jia, H.M. Chen, H. Xiao, J. Li, Y. Huang, B. Liu, J. Am. Chem. Soc. 141 (2019) 3014–3023. [1] L. Hammarström, Acc. Chem. Res. 42 (2009) 1859–1860. [44] W. Sun, Z. Wang, W.Q. Zaman, Z. Zhou, L. Cao, X.-Q. Gong, J. Yang, Chem. [2] S. Anantharaj, V. Aravindan, Adv. Energy Mater. 10 (2020) 1902666. Commun. 54 (2018) 996–999. [3] M.G. Walter, E.L. Warren, J.R. McKone, S.W. Boettcher, Q. Mi, E.A. Santori, N.S. [45] J. Rossmeisl, Z.-W. Qu, H. Zhu, G.-J. Kroes, J.K. Nørskov, J. Electroanal. Chem. Lewis, Chem. Rev. 110 (2010) 6446–6473. 607 (2007) 83–89. [4] X. Sala, S. Maji, R. Bofill, J. García-Antón, L. Escriche, A. Llobet, Acc. Chem. Res. [46] A. Zagalskaya, V. Alexandrov, J. Phys. Chem. Lett. 11 (2020) 2695–2700. 47 (2014) 504–516. [47] H.N. Nong, L.J. Falling, A. Bergmann, M. Klingenhof, H.P. Tran, C. Spöri, R. Mom, [5] J.D. Blakemore, R.H. Crabtree, G.W. Brudvig, Chem. Rev. 115 (2015) 12974– J. Timoshenko, G. Zichittella, A. Knop-Gericke, S. Piccinin, J. Pérez-Ramírez, B.R. 13005. Cuenya, R. Schlögl, P. Strasser, D. Teschner, T.E. Jones, Nature 587 (2020) 408– [6] B.M. Hunter, H.B. Gray, A.M. Müller, Chem. Rev. 116 (2016) 14120–14136. 413. [7] J.H. Montoya, L.C. Seitz, P. Chakthranont, A. Vojvodic, T.F. Jaramillo, J.K. [48] D. González, J. Heras-Domingo, M. Sodupe, L. Rodríguez-Santiago, X. Solans- Nørskov, Nat. Mater. 16 (2017) 70–81. Monfort, J. Catal. 396 (2021) 192–201. [8] C. Spöri, J.T.H. Kwan, A. Bonakdarpour, D.P. Wilkinson, P. Strasser, Angew. [49] Z. Xu, J. Rossmeisl, J.R. Kitchin, J. Phys. Chem. C 119 (2015) 4827–4833. Chemie Int. Ed. 56 (2017) 5994–6021. [50] S. Siahrostami, A. Vojvodic, J. Phys. Chem. C 119 (2015) 1032–1037. [9] H. Ding, H. Liu, W. Chu, C. Wu, Y. Xie, Chem. Rev. 121 (2021) 13174–13212. [51] I.C. Man, H. Su, F. Calle-Vallejo, H.A. Hansen, J.I. Martínez, N.G. Inoglu, J. [10] E. Fabbri, A. Habereder, K. Waltar, R. Kötz, T.J. Schmidt, Catal, Sci. Technol. 4 Kitchin, T.F. Jaramillo, J.K. Nørskov, J. Rossmeisl, ChemCatChem 3 (2011) 1159– (2014) 3800–3821. 1165. [11] T. Reier, H.N. Nong, D. Teschner, R. Schlögl, P. Strasser, Adv. Energy Mater. 7 [52] Y.-H. Fang, Z.-P. Liu, J. Am. Chem. Soc. 132 (2010) 18214–18222. (2017) 1601275. [53] Y. Ping, R.J. Nielsen, W.A. Goddard, J. Am. Chem. Soc. 139 (2017) 149–155. [12] L. An, C. Wei, M. Lu, H. Liu, Y. Chen, G.G. Scherer, A.C. Fisher, P. Xi, Z.J. Xu, C. [54] J.A. Gauthier, C.F. Dickens, L.D. Chen, A.D. Doyle, J.K. Nørskov, J. Phys. Chem. C Yan, Adv. Mater. 33 (2021) 2006328. 121 (2017) 11455–11463. [13] C.C.L. McCrory, S. Jung, I.M. Ferrer, S.M. Chatman, J.C. Peters, T.F. Jaramillo, J. [55] L.G.V. Briquet, M. Sarwar, J. Mugo, G. Jones, F. Calle-Vallejo, ChemCatChem 9 Am. Chem. Soc. 137 (2015) 4347–4357. (2017) 1261–1268. [14] J. Yu, Q. He, G. Yang, W. Zhou, Z. Shao, M. Ni, ACS Catal. 9 (2019) 9973–10011. [56] C.F. Dickens, C. Kirk, J.K. Nørskov, J. Phys. Chem. C 123 (2019) 18960–18977. [15] D. Lebedev, R. Ezhov, J. Heras-Domingo, A. Comas-Vives, N. Kaeffer, M. [57] A. Halder, C. Liu, Z. Liu, J.D. Emery, M.J. Pellin, L.A. Curtiss, P. Zapol, S. Vajda, A. Willinger, X. Solans-Monfort, X. Huang, Y. Pushkar, C. Copéret, A.C.S. Cent, Sci. B.F. Martinson, J. Phys. Chem. C 122 (2018) 9965–9972. 6 (2020) 1189–1198. 85 D. González, M. Sodupe, L. Rodríguez-Santiago et al. Journal of Catalysis 412 (2022) 78–86 [58] K. Bhattacharyya, C. Poidevin, A.A. Auer, J. Phys. Chem. C 125 (2021) 4379– [69] P.E. Blöchl, Phys. Rev. B 50 (1994) 17953–17979. 4390. [70] K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T.A. Arias, R.G. Hennig, J. [59] D. González, B. Camino, J. Heras-Domingo, A. Rimola, L. Rodríguez-Santiago, X. Chem. Phys. 140 (2014) 084106. Solans-Monfort, M. Sodupe, J. Phys. Chem. C 124 (2020) 1227–1237. [71] H.H. Heenen, J.A. Gauthier, H.H. Kristoffersen, T. Ludwig, K. Chan, J. Chem. [60] D. González, J. Heras-Domingo, S. Pantaleone, A. Rimola, L. Rodríguez-Santiago, Phys. 152 (2020) 144703. X. Solans-Monfort, M. Sodupe, ACS Omega 4 (2019) 2989–2999. [72] J.K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J.R. Kitchin, T. Bligaard, H. [61] D. González, M. Sodupe, L. Rodríguez-Santiago, X. Solans-Monfort, Nanoscale Jónsson, J. Phys. Chem. B 108 (2004) 17886–17892. 13 (2021) 14480–14489. [73] H.N. Nong, L. Gan, E. Willinger, D. Teschner, P. Strasser, Chem. Sci. 5 (2014) [62] G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169–11186. 2955–2963. [63] G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 558–561. [74] R.R. Rao, M.J. Kolb, L. Giordano, A.F. Pedersen, Y. Katayama, J. Hwang, A. Mehta, [64] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. H. You, J.R. Lunger, H. Zhou, N.B. Halck, T. Vegge, I. Chorkendorff, I.E.L. [65] S. Grimme, J. Comput. Chem. 25 (2004) 1463–1473. Stephens, Y. Shao-Horn, Nat. Catal. 3 (2020) 516–525. [66] Y. Ping, G. Galli, W.A. Goddard, J. Phys. Chem. C 119 (2015) 11570–11577. [75] K.A. Stoerzinger, L. Qiao, M.D. Biegalski, Y. Shao-Horn, J. Phys. Chem. Lett. 5 [67] X. Solans-Monfort, V. Branchadell, M. Sodupe, M. Sierka, J. Sauer, J. Chem. Phys. (2014) 1636–1641. 121 (2004) 6034–6041. [76] M. Kan, X. Qian, T. Zhang, D. Yue, Y. Zhao, A.C.S. Sustain, Chem. Eng. 5 (2017) [68] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758–1775. 10940–10946. 86