Tuning metal single atoms embedded in N_xC_y moieties toward high-performance electrocatalysis

Abstract

Noble nanoparticle (NP)-sized electrocatalysts have been exploited for diverse electrochemical reactions, in particular, for an eco-friendly hydrogen economy such as water splitting. Recently, minimal amounts of single atoms (SAs) are exploited to maximize the active surface area and to tune the catalytic activity by coordinating the SAs in defect sites of N-doped graphene (G_N). For the hydrogen evolution reaction (HER) and oxygen evolution/reduction reactions (OER/ORR), we show high-performance 3d–5d transition metal (TM) SA catalysts using density functional theory (DFT) along with machine learning (ML)-based descriptors. We explore the stability and activity of TM–G_N from the view of structure/coordination, formation energy, structural/electrochemical stability, electronic properties, electrical conductivity, and reaction mechanism, which have not been seriously explored yet. Among various –N_nC_m moieties, the –N₂C₂ moieties tend to be more easily formed and show higher electrochemical catalytic performance and longer durability (without aggregation/dissolution) compared with the widely studied pure –C₄/C₃ and –N₄/N₃ moieties. We found that some TM(SA)s favor a new OER/ORR mechanism, completely different from any known mechanism. The ML-based descriptors showing super HER/OER/ORR performances better than those of bench-mark noble metal catalysts are assessed. In the N₂C₂ templates, Ni/Ru/Rh/Pt show low HER overpotentials. Here, the H adsorption sites are shared by both the metal and C (not N), which was undiscussed in most of the previous literature where the H is attached on top of a metal atom. Low OER overpotentials are found for Pt/Ni–N₂C₂, Ni/Pd–C₄, and Rh–N₄, while low ORR overpotentials are found for Ir/Rh-N₄, Pd–C₄, Ru–N₃C₁ and Ni/Pd/Pt–N₁C₃. The present findings should help in designing high-performance SA catalysts for other various electrocatalytic reactions such as the ammonia evolution reaction.

---

Broader context

A high catalytic performance in water splitting for hydrogen fuel is essential for an ever-clean hydrogen economy including hydrogen vehicles and hydrogen power generation. Transition metal d-orbitals can be modified to be highly active for specific electrochemical reactions by altering the coordination environment, for which we could utilize high-throughput computational screening. However, the complicated electrochemical reactions have shown wide-ranging different activities depending on the chemical environment and calculation methods, which has caused a lack of reliable clear-cut reference data. Thus, here we report the electrocatalytic properties of thirty 3d–5d transition-metal single-atoms [(TM)SAs] ligated by –N_nC_m, using density functional theory along with machine learning-based descriptors. High-performance (TM)SA catalysts in hydrogen evolution and oxygen evolution/reduction reactions are assessed. The –N₂C₂ moiety tends to show higher catalytic performance than those of the widely studied pure –C₄/C₃ and –N₄/N₃ moieties. A vast amount of information for electrochemical reactions including new reaction mechanisms will provide impacts on high-performance single-atom catalytic tuning. In addition to the data for electronic/geometric structures and electrocatalytic activities, their XANES spectra and conductivities, which are not available in the literature but are critically important for material characterization and faster reaction kinetics, will be a useful guide for developing new catalysts towards diverse electrocatalytic reactions.

Introduction

Recently, it has been of utmost importance to develop high-performance, energy-efficient and eco-friendly electrocatalysts. These can be utilized for various types of electrochemical reactions^1,2 such as the HER,³ OER/ORR,^4–7 water oxidation/splitting,⁸ fuel cell reactions,⁹ carbon-oxidation/CO-reduction–oxidation/CO₂-reduction reactions,^10,11 MeOH/EtOH oxidation reactions,¹² nitrogen-reduction/ammonia-oxidation reactions,¹³ peroxide/chlorine evolution reactions,¹⁴ acid oxidation reaction,¹⁵etc. Until recently, most of these catalysts were usually made of highly expensive and scarce novel metals (Pt/Ir/Ru).¹⁶ Hence, tremendous efforts have been paid to using a tiny amount of novel metal SAs/NPs^17,18 or replacing them with inexpensive alternatives.¹⁹ In this regard, minimal amounts of metal (M) SAs^20–26 are exploited to maximize the active surface area and to tune the catalytic activity by coordinating the SAs to glue atoms in defect sites of graphene (Gr) or other 2D materials.^17,27–36 However, while most metal NPs behave like conducting metals, metal SAs embedded in supports often lose some of the electric conductivity (σ), because once embedded metal SAs are ligated by C/N atoms in vacant sites, the π-conjugation through C/N atoms in the 2D sheet can be lost.^17,27,28,31 This would cause a reduction in the kinetic reaction rate of SA electrocatalysts.^31,37,38 To maximize the catalytic performance, (TM)SA-embedded electrocatalysts need to be optimized in view of coordination-structure, stability, conductivity, and catalytic activity. To this end, high-throughput computer screening is a powerful approach.^39–46

(TM)SA catalysts embedded in M–N₄/C₄ and M–N₃/C₃ moieties of G_N/Gr have been reported to be promising as HER/OER/ORR electrocatalysts.^31,34 The localized d orbital energy levels in (TM)SAs are modified by the coordination environments.⁴⁷ Thus, they dictate the stability and activity of SAs. Also, they dictate the σ when the (TM)SAs are embedded in π-conjugated G_N, which enables fast transport of charges (i.e., electrons) in electrochemical reactions.¹⁷ In this regard, we have considered 3d (Sc–Zn), 4d (Y–Cd), and 5d (La, and Hf–Hg) TMs in the mono/di-vacancy defect sites of G_N. G_N is particularly considered here because of its high σ compared with that of Gr and its high tunability towards better electrochemical activity by a variety of coordination environments around a (TM)SA with tetra-coordination (–N_n=0–4C_4−n) in di-vacancy sites or tri-coordination (–N_n=0–3C_3−n) in mono-vacancy sites. We thus studied 11 different pyridinic/graphitic-type defect configurations for each (TM)SA embedded in G_N. For high-throughput computational screening towards high-performance electrocatalysts of (TM)SAs embedded in G_N [M–N_nC_m–G_N], we have performed DFT calculations (see Method & Tables S1 and S2, ESI†).

Here we investigate the stability of each embedded (TM)SA, the X-ray absorption near-edge structure (XANES) (which helps in characterization of the coordination environment of SAs onto the G_N template), the electronic properties including σ, the reaction mechanisms and the HER/OER/ORR activities of M–N_nC_m–G_N. We find that the (TM)SAs, in particular, in the –N₂C₂ moiety, display superior durability and high σ in favour of fast HER reaction kinetics. These complicated electrochemical reactions have shown very different activities depending on different chemical environments and calculation methods, causing a lack of reliable clear-cut reference data. Thus, we report the lucid understanding of TM(SA) catalysts, in particular, for HER/OER/ORR towards water splitting and hydrogen fuel, which are vital for clean green energy environments such as hydrogen cars and electric power. We assess highly effective (TM)SAs ligated by various –N_nC_m, in particular, by –N₂C₂ for superb electrochemical performance. The machine learning-based descriptors describing high-performance HER/OER/ORR catalysts are also assessed.

Results and discussion

Formation and stability of M–N_nC_m–G_N

To analyse the degree of formation and stability of various –N_nC_m moieties in G_N/Gr defects, we study the representative pyridinic/graphitic-type defects and (TM)SA-embedded configurations and these formation energies including metal-embedding energies, based on the refined DFT method (Tables S1–S6, ESI†). Here we show that the (TM)SAs embedded in the –N₂C₂ vacant sites are consistently the most stable among possible conformations of M–N_nC_m–G_N.

SA-catalysts templated in G_N are often synthesized using Gr, G_N, N-containing carbon compounds (such as melamine) and other materials.^34,48–50 In the absence of metals, it is not energetically easy to form significant vacancy sites in G_N/Gr (based on Table S3, ESI†). M–N_nC_m–G_N is often formed in the presence of (TM)SAs and N sources (NH₃, N₂H₄, etc.) or N-containing carbon materials.^{17,27,34,48–50} According to our DFT-predicted formation energies of 11 different pyridinic/graphitic-type defect configurations (Fig. 1 and Fig. S1, ESI†) with optimal H atoms attached on edge C/N atoms (Fig. 2a and Table S3, ESI†), the –N_nC_m defect sites in G_N can be formed easily in the presence of H atoms (which is naturally present even in moderately high vacuum) with the formation energies of ∼2 eV. The formation of pyrrolic defect sites is less favoured due to their higher formation energies (∼5 eV) (Fig. S2 and Table S4, ESI†). We investigate the formation energies and metal embedding stabilities on the mono/di-vacancy pyridinic defect sites with tri/tetra-coordination (N_n=0–3C_3−n)/(N_n=0–4C_4−n) (Table S5, ESI†).

Fig. 1 Configurations of a (TM)SA embedded in M–N_nC_m–G_N [n + m = 4 (a di-vacancy site) or 3 (a mono-vacancy site)]. The superscripts “^b” and “^c” in N₂C₂^b and N₂C₂^c denote two different cis forms of two N atoms with respect to the central TM, which is Ni here [C: brown, N: blue, and Ni: dark green]. (cf. Fig. S1, ESI† where the TM is Zr located above the graphene sheet due to its slightly larger radius; Fig. S2 and Tables S8 and S9 (ESI†) which show the off-plane distances for various (TM)SA moieties before and after H-adsorption; Table S10, ESI† after O-adsorption.)

Fig. 2 Vacancy forming energies (E_v in eV) in –N_nC_mH_h moieties of G_N. (a) E_v in –N_nC_mH_h moieties of G_N in the presence of H (in atmosphere or acid medium) but in the absence of metal and their representative structures with the optimal H atoms attached on the edge C/N atoms at the RPBE(U) + TS level of theory. (E_v = E_total − mμ_C − nμ_N − hμ_H, where μ_C = energy of a C atom in graphene, μ_N = energy of a N atom in N₂, and μ_H = energy of a H atom in H₂), and [C: brown, N: blue, and H: cyan] (cf. Table S4, ESI† for E_v). (b) (TM)SA embedding stability (E_emb/coh ≡ E_emb − E_−coh; E_coh = −E_−coh in eV) color map of (TM)SAs in G_N over the TM cohesion. Bluish colors strongly favor SA-embedding over TM-cohesion, while reddish colors favour the opposite. “*” denotes unstable systems. The –N_nC_4−n (particularly –N₂C₂) divacancy sites tend to be more bluish, while the –N_nC_3−n (particularly M–N₃) mono-vacancy sites tend to be more reddish (cf. Table S6 for E_emb, Table S7 for E_emb/coh, ESI†). In –N_nC_3−n mono-vacancy sites, all (TM)SAs are off the G_N plane by >1.5 Å.

The stability criterion of “metal embedding (embedding a metal atom into/onto the support)” over “metal cohesion (adding a metal atom to the bulk)” is defined as the difference (E_emb/coh = E_emb − E_−coh) of the (TM)SA-embedding energy (E_emb) on G_N (Fig. 2b and Tables S6, S7, ESI†) from the bulk cohesive energy (E_coh= −E_-coh; Table S1, ESI†). Here, embedding energy is defined in a negative value. A more negative embedding energy means stronger embedding, while cohesion is generally defined by a positive value. If E_emb/coh < 0, the embedding into graphene is preferred over metal clustering. However, even for small positive E_emb/coh, the embedding into graphene can still be transiently stable unless the concentration of metal is significantly present to be able to form clusters.

Embedding energy shows that the N₂C₂ moiety is the most favorable for the M–N_nC_m–G_N formation with the trend of embedding strength (−E_emb) in the order of N₂C₂ > (N₃C₁, N₂C₂^b, and N₁C₃) > (N₂C₂^c) > (C₄, N₄, C₃) > (N₂C₁, N₁C₂) > N₃ (where N₂C₂^b and N₂C₂^c are the cis forms, while N₂C₂ is the trans form; Fig. 1 and Fig. S1, ESI†). The N₂C₂ case shows high embedding stability of (TM)SA over TM-cohesion, and thus their SAs can be easily formed (shown in blue-to-yellow color in Fig. 2b). In contrast, the cases of N₁C₂, N₂C₁ and N₃ are less stable than the (TM)-cohesion (shown in orange/red color). Thus, these SAs are not so stable that their clusters/NPs are formed. Most (TM)SAs are slightly large to fit into the di-vacant site of G_N, but several of the (TM)SAs (such as (Mn/Fe), Co/Ni/Cu/Zn/Pt) just fit into that site (Fig. S3 and Table S8, ESI†). The planar embedding structures are highly stable. However, in all –N_nC_3−n mono-vacancy moieties, all (TM)SAs are off the plane by more than 1.5 Å. Upon *H/*O adsorption, a slight/small structural change can be noted (Fig. S4 and Tables S9, S10, ESI†).

To understand the aggregation of SAs by clustering, we report the 2nd (TM)SA's clinging energy (E_2M) and the (TM)SA-stability over TM-cohesion (E_2M/coh) (Tables S11 and S12, ESI†). A few M–N_n=0–3C_3−n–G_N systems have smaller E_2M/coh (or even negative values) than M–N_n=0–4C_4−n–G_N due to their weak metal–ligand interaction. Then, in some cases 2nd (TM)SAs bind the C atoms in the cases of –N₁C₃–G_N and –C₄–Gr sites, showing small E_2M/coh. The Hg/Cd/Zn/Ag/Au/Pd TMs have rather small E_2M/coh, and Mn/Cu have moderate values, while other (TM)SAs show large stability (>3eV) over TM-cohesion. Since d-orbitals of (TM)SAs split due to the metal–ligand interaction by C/N atoms in G_N, the 1st (TM)SAs hardly interact with the 2nd (TM)SAs.

Electrochemical stability

To analyse the electrochemical stability of a (TM)SA embedded in a –N_nC_m moiety, we investigate its stability over metal-aggregation and metal-dissolution. We find that the –N₂C₂ moiety generally shows the highest durability against aggregation and dissolution.

The electrochemical stability is governed by dissolution potential ,⁵¹ where is the standard dissolution potential (pH = 0) of bulk metal in aqueous solution, e is the electron charge, and n_e is the number of electrons involved in dissolution. A more positive U_ds relative to the equilibrium potential indicates that the metal atoms bind the M–N_nC_m–G_N support more strongly and the dissolution of metal atoms can be better avoided under electrochemical reactions, which enhances the durability of (TM)SA in a harsh electrochemical environment. While there can be a different number of electron transfers in the dissolution process, we here use a broad range of U_ds (>0 V) based on the reported experimental studies on most of the synthesized (TM)SA-catalysts.⁵²

The stability against metal-aggregation and metal-dissolution in the M–N_nC_m–G_N system is plotted in the E_emb/cohvs. U_ds plot (Fig. 3). Negative E_emb/coh indicates the thermodynamical stability against clustering, while positive U_ds reflects the electrochemical stability against dissolution. Following these two stability criteria, the whole figure is divided into four quadrants, depicted by blue, green, yellow, and pink. We set the stability criteria as the unaggregatable and indissoluble system with E_emb/coh < 0 eV and U_ds > 0 V. It can be noted that most of the double vacancy systems are stable, and the –N₂C₂ moiety is particularly stable among –N_mC_n moieties (Tables S13 and S14, ESI†). In contrast, most TM elements (except for La and Y) in mono-vacancy site of –N_n=1–3C_3−n–G_N (excluding –C₃–Gr having weak stability) show poor stability over metal cohesion, and so these SAs tend to aggregate as clusters. This is a salient finding in that most studies for TM(SA) electrocatalysts in Gr/G_N have focused on pure C(C₄/C₃) or N(N₄/N₃) coordination except for very recent studies showing remarkable electrocatalytic performance in both activity and durability by using the –N₂C₂ moiety.^17,27–29 Among all the M–N_nC_m–G_N, we find that the –N₂C₂ moieties show favorable E_emb/coh and U_ds values for promising hosts, which provide high durability for metal atoms. This finding would help in accelerating various SA-driven electrochemical reactions.

Fig. 3 Stability of (TM)SAs in M–N_nC_m–G_N against metal aggregation and electrochemical dissolution in terms of (TM)SA-embedding over metal-cohesion (E_emb/coh) and dissolution potential (U_ds) of metal atoms for M–N_nC_m–G_N (cf. Tables S13 and S14, ESI†). The more negative E_emb/coh and the more positive U_ds indicate the better stability.

Generally, the early (TM)SAs (Sc/Ti/Y/Zr/La/Hf) are more likely to be stabilized on the defective sites of the substrates, as indicated by the more negative E_emb/coh values, while the late (TM)SAs (Co–Zn, Rh–Cd, and Ir–Hg) tend to be stabilized against electrolytic dissolution. Overall, (TM)SAs including Pt/Pd/Au/Ag/Cu/Ir/Rh/Ni in the –N₂C₂ moiety display superior durability.

Electronic structure and XANES

The band structure is reflected in the X-ray absorption near edge structure (XANES), conductivity, catalytic activity, etc. The coordination environments of defect sites in G_N affect the electronic properties of (TM)SAs and hence the thermodynamic adsorption/desorption strength of intermediates in electrochemical redox reactions. Here, we study the electronic properties of TM(SA) in M–N_nC_m–G_N and XANES showing the coordination environments of a (TM)SA to aid experimental characterization.

Fig. 4a–f show the projected density of states (PDOS) for (V/Ni/Zn)–N₂C₂–G_N and (V/Ni/Zn)–N₃–G_N defect sites. Square planar complexes are stable as 16-electron species.⁴⁷ The Ni on a Ni–N₂C₂–G_N (as Ni²⁺ ion) has eight d electrons after forming four bonds with four C/N atoms, resulting in a 16-electron d⁴ + dsp² electron configuration in the square-planar like structure (Fig. 4a). Early TM elements like V on di-vacancy sites cannot fill the 16-electron configuration, and so the orbitals are not occupied evenly, showing a magnetic property (Fig. 4b, Fig. S5, S6, Table S15, ESI†). Group 12 elements like Zn have fully occupied 10 d electrons in +2 oxidized states along with 8 electrons from –N₂C₂, resulting in 18 electrons in total (Fig. 4c). Meanwhile, (TM)SAs on mono-vacancy sites have tetrahedral like structures. As an example, in Ni–N₃–G_N, the Ni atom (3d⁸4s²) is tri-coordinated to N atoms, resulting in 16 electrons in total around the Ni atom. The calculation shows d⁹ electrons with one unoccupied minor spin state in the d-orbital (d_xz) due to the electron transfer to the N atoms, showing a magnetic moment (Fig. 4d). The magnetic properties can change upon *H/*O-adsorption (Fig. S7, S8 and Tables S16, S17, ESI†). The less oxidized metal atom can be found in early TM elements as in V–N₃–G_N (Fig. 4e) and group 12 elements as in Zn–N₃–G_N (Fig. 4f).

Fig. 4 Projected density of states (PDOS) of pristine M–N₂C₂–G_N (a–c) and M–N₃–G_N systems (d–f) and (g) XANES of Co–N₂C₂–G_N. (a and d) Case for M = V [early transition metal d³]. (b and e) Case for M = Ni [square planar d⁸; 16 electron rule]. (c and f) Case for M = Zn [group 12; fully occupied d orbital d¹⁰]. The major-spin (positive) and minor-spin (negative) PDOS of M (d orbitals) are plotted in each upper box, and those of M (s, p orbitals) and C, N and H are plotted in each lower box. Charge densities related to H-adsorption on each system are shown on top of PDOS plots. (Ni,Zn)–N₂C₂–G_N have planar structures, while V–N₂C₂–G_N has a non-planar structure. See the text for the planar d⁴ + dsp² electron configuration favoring 16-electron rule and the octahedral sp³d² electron configuration favoring 18-electron rule in the Ni case. (g) Co K-edge XANES of Co–N_nC_m–G_N showing edge shifts to a higher energy than that of Co bulk, while giving a rising pre-edge (an exemplary case to show the difference between various Co–N_nC_m–G_N coordination environments). See Fig. S9 (ESI†).

From these differences, XANES of M–N_nC_m–G_N systems show each characteristic peak. For example, the edge energy of XANES shifts according to the oxidization state of (TM)SAs. Also, the symmetry of the structure and p–d orbital hybridization increase the intensity of pre-edge peaks. The group of TM is another factor affecting XANES associated with the number of empty orbitals for electron excitation. The calculated K-edge XANES data of all 3d and 4d (TM)SAs are in Fig. S9 (ESI†), which are essential for analyzing the SA-coordination analysis but not available in the literature. In Fig. 4g, the Co K-edge XANES is provided as an exemplary case to show the difference between various Co–N_nC_m–G_N coordination environments. Co–N₂C₂–G_N and Co–N₂C₂^c–G_N show distinguishable pre-edge structures from each other (Supplementary Note, ESI†). The pre-edge peaks are usually broader and less intense in 4d (TM)SAs, and additional effects like symmetry and distortion of the structure affect XANES, which helps in characterizing the coordination environments for each (TM)SA. In-situ XANES analysis upon *H/*O–adsorption during the HER/OER can also be made, as seen from the PDOS (Fig. S8, ESI†).

Electrical conductivity

A good electrical conductivity is vital for high-performance electrocatalysts. This issue was not important in the past because the catalysts on electrodes were made with conductive metals or carbon materials, but nowadays the SA-embedded supports are not necessarily good conductors. Hence, σ is an important issue in determining the kinetic reaction rate, in particular, in fast reactions such as the HER. However, there is no serious study on this issue except for our recent work.^17,27,28

In the HER, the reaction rates in (TM)SAs are often very fast by the Tafel mechanism where H₂ gas is released by combining two neighboring adsorbed H atoms on catalytic sites; thus, fast charge transfer plays a critical role.³¹ This charge transfer effect is rather less important in the OER, which shows a rather sluggish reaction rate. Thus, we consider the fast H-adsorption in the HER. Here, the σ for M–N_nC_m–G_N before/after H-adsorption is calculated at 300 K and compared with the σ of 2.1% N-doped G_N [σ^o ≈ 7 × 10²σ(Gr)] (Fig. S10 and Tables S18, S19, ESI†). Overall, M–N₂C₂ tends to show higher σ than any other double vacancy M–N_nC_m–G_N system (Fig. 5). The change in electronic structure by coordination environment affects the σ significantly. The covalent bonding of Gr is disrupted by (TM)SA-doping and N-adsorption; hence, the σ of metal-doped G_N tends to be much reduced, showing poor conductivity. Among the (TM)SAs, early TM elements [group-3(Sc/Y/La), and group-4(Ti/Zr/Hf)] show moderately high σ (>0.2σ^o) regardless of N_nC_m-coordination/H-adsorption environments, while other (TM)SAs doped on –N_nC_m–G_N introduce band gaps, resulting in rather poor σ. As examples, we analysed the band structure of Ni/Pd/Pt–N₂C₂–G_N as a candidate for HER catalysts (Fig. S11, ESI†). Their pristine structures show poor conductivity with small band gaps near the Fermi energy level (E_F), while the electrons fill the valence band fully up to the Dirac-like point in the band structure. Upon adsorbing H*, the conduction band gets filled and the σ increases with one conducting channel at E_F = 0. Highly conductive systems have conducting channels at E_F and accelerate electrocatalytic reactions with fast charge transfer.

Fig. 5 Color code map for the electrical conductivity (σ) of M–N_nC_m–G_N with respect to that of G_N before and after H-adsorption. Each σ is calculated from the band structure. Assuming similar scattering rates (τ) for M–N_nC_m–G_N and G_N from the experimental value of ∼100 fs,⁶² the [σ(M–N_nC_m–G_N)/σ^o] values are reported in Tables S18 and S19 (ESI†) [σ^o/σ(Gr) = ∼7 × 10²]. “*” denotes an unstable state. Overall, among the divalent N_nC_4−n moities, the N₂C₂ moiety tends to show high σ for most of the TM(SA)s embedded in G_N. For example, Co, Ni, Pd and Pt show enhanced σ upon H-adsorption, which accelerates the HER rate towards the Tafel mechanism, which readily combines two neighboring adsorbed H atoms into H₂ gas. Since the OER/ORR rates are sluggish, the large overpotential is much more critical than the high σ. Thus, this high σ is particularly investigated for fast HER.

We also study the σ when O*, HO*, and HOO* are adsorbed on the active sites of M–N_nC_m–G_N, where the σ tends to decrease compared to the pristine case (Fig, S12 and Table S19b, ESI†). However, in the cases of OER/ORR, the overpotentials are large, and so lowering their high barriers, which exponentially forbid the charge transfer, is the critical issue, while the conductivity is not a significant issue unless it becomes very small.

Electrochemical activity: reaction mechanism

For the adsorption of intermediates (O*/HO*/HOO*-adsorption) in the OER/ORR, the well-known mechanism (path-I) for M–N/C prevails. However, we note a new mechanism (path-II) in the OER/ORR for certain (TM)SA cases of M–N_nC_m–G_N. Here, we discuss both mechanisms.

Based on the structural/electronic stability of M–N_nC_m–G_N, the candidates for HER/OER/ORR active catalysts are considered particularly among double vacancy systems. The schematics of HER/OER/ORR mechanisms are shown in Fig. 6. The H*/O*/HO*/HOO*-adsorption free energies (ΔG_{H*/HO*/O*/HOO*}) are investigated (Methods). For the HER, the reaction rate in Volmer–Tafel mechanism is fast by combining neighboring adsorbed H atoms on active sites into H₂, while that in the Volmer–Heyrovsky mechanism is slow because an additional H* should be near the H* adsorbed on the neighboring active site. In OER and ORR cases, path-I where each O*/HO*/HOO* is adsorbed on a metal SA site has generally been considered. However, we find a new path (path-II), which has OH⁻ on the metal–SA and beside O* (rather than on the O* attached to the metal-SA because of the difficulty in adsorbing OH⁻ on O*). Thus, the O-poisoned site (M–O*) works as a new active site, and each O*/HO*/HOO* is attached on the metal-SA site. The reverse steps of these two mechanisms are also applied to the ORR. Though path-I is generally preferred in most (TM)SAs, in some metals path-II is preferred. For example, V–N₄ follows path-II for the OER. The preferred mechanism depends on TM and the coordination environment.

Fig. 6 Schematic representation for the mechanisms and reaction pathways of electrocatalytic reaction cycles on a metal binding site. (a) HER. (b) Free energy profiles of V–N₄ for OER path-I and path-II. (c) OER: path-I. (d) OER: path-II. (e) ORR: path-I. (f) ORR: path-II. On most systems, the OER/ORR follows the conventional mechanism with the metal as the active site (path-I). On a few systems, surface poisoning by oxygen inhibits direct *OOH formation. They follow a different pathway with an O-poisoned metal site (O*–M) as the active site, while other intermediates are adsorbed on O*-M site [path-II: Re/Os for C₄; (Ta,Re)/(Mo,Tc,Ru,W,Os) for N₂C₂; (Nb,Mo,Tc,Ta,W)/(V,Ru,Re,Os) for N₄]. The direct O₂ OER mechanism (or dissociative ORR mechanism) is not preferable (see Tables S26 and S27, Supplementary Note, ESI†). See the equation for ΔG_1–8 in Methods.

HER activity

The HER of M–N_nC_m–G_N is achieved through their active sites such as metal atom M, and bridge sites between M and C, or C/N atoms. One should also consider if the Kubas active site⁵⁴ can help in accelerating H₂ desorption by combining both H atoms, leading to a Tafel-like reaction. Here we assess high-performance HER SA catalysts.

The catalytic activity is investigated in terms of H-adsorption binding/free energies (ΔE_H*/ΔG_H* in Tables S20 and S21, ESI†) for the intermediates during the HER. The solvent effect on the HER is rather insignificant within 0.03 eV for ΔE_H*/ΔG_H* (see Table S22, ESI†). The zero-point energy and entropy contributions to ΔG_H* when H is on a metal atom show an almost constant shift by 0.24 eV (for H₂),⁵³ but this shift changes significantly when H is on a C/N atom or a bridging bond between metal and C/N atoms, and thus such a correction is made by phonon calculations. Then, the refined HER overpotential (η^HER, Tables S30a and S31, ESI†) is calculated based on the ΔG_H*, which is used to find outstanding HER catalysts. Mo/Rh/La/Os/Pt in C₄–Gr, Ti/Zr/Ag/Cd in N₁C₃–G_N, Co/Ni/Ag in N₂C₂^b–G_N, Cr/Cu/Cd in N₃C₁–G_N, Rh/Cd in N₄–G_N, and Ni/Ru/Rh/Pt in –N₂C₂-G_N templates show −η^HER < 0.1 V (red colored in Fig. 7a and in filled red triangle in Fig. 7b) and are likely to be super-performing HER catalysts, given that the Pt(111) surface shows −η^HER = 0.2 V on the fcc-hollow site. Indeed, Pt–N₂C₂ (−η_expt/calc= 0.01–0.02/0.03 V),¹⁷ Ru–N₂C₂ (−η_expt/calc= 0.01–0.04/0.06 V)^27–29 and Rh–N₂C₂ (−η_expt/calc = 0.04/0.06 V)⁵⁵ are reported as remarkably active HER (TM)SA catalysts. In these cases, the H adsorption sites are shared by both the metal and C (not N) except for Ru (Fig. 7a).

Fig. 7 HER overpotentials (η^HER) of M–N_nC_m–G_N based on the refined values (Table S30 based on Tables S20, S21, ESI†). (a) Color code map for η^HER (upper triangle) and most stable H-adsorption sites (H*) (lower triangle), which are denoted by metal site (in white blank), metal-C bridge site (by a hat symbol “^”), C site (by “C”) or N site (by “N”). The solvent effect for ΔE_H*/ΔG_H* is generally very small (Table S22, ESI†). Also, it is worthwhile to consider the Volmer–Kubas reaction free energies (Table S23, ESI†). “*” denotes unstable systems. The reddish colors represent near-zero free energies (responsible for near-zero overpotential). Zr/Ti/Ag/Cd–N₁C₃, Pt/Ni/Ru/Rh–N₂C₂, Cd/Rh–N₄, and Rh/Os/Pt–C₄ would be high-performance HER catalysts (−0.1 V < η^HER < 0 in red color) and those within (−0.2 V < η^HER < −0.1 V in pink color) are likely to perform well in the HER, as compared with a Pt(111) surface (η^HER = −0.2 V). (b) Cases for −0.2 V < η^HER < 0, many of which belong to the –N₂C₂ type (denoted by filled red triangle) and the –N₂C₂^b/N₂C₂^c types (in open red triangles).

It is generally expected that the strongest H-adsorption site for M–N_nC_m–G_N is just metal-SA M, as noted in most of the previous literature, which did not consider other options. However, in certain cases it is a C or N site, while often sharing the metal site. The H adsorption sites can be shared by both metal and C (not N) (i.e., the bridge site denoted by a hat symbol in Fig. 7). The metal-adjacent C sites (but not N sites) are the most active in some cases of –N₃C₁, –N₂C₂, –N₂C₂^b, –N₂C₂^c, –N₁C₃, and –C₄, while the N sites can be the most active in certain cases of –N₃C₁ and –N₄. In certain cases where a C or N site is the most active, the metal site can be a coactive site. Among the active C/N atoms, the metals in Zr–C₄ and Pt–N₄ are nearly compatible in activity to the C and N atoms, respectively. In particular, when a H atom is adsorbed on 3d (Sc/Cr/Mn/Fe/Co/Ni/Cu/Zn), 4d (Y/Rh/Pd), and 5d (Pt) TM–G_N systems, the H tends to share a C site, while a H adsorbed on the (Ag/La/Pt)–G_N system tends to be on C/N sites instead of metal sites (Fig. 7a).

When a 2nd H atom is attached to the same TM atom embedded in G_N, the Kubas active site⁵⁴ can help in accelerating H₂ desorption by combining both H atoms, leading to a Tafel-like reaction, but its role is not dominating in the HER reaction. The Volmer–Kubas structures and free energy changes for selected (TM)SAs in (N₂C₂) systems are shown along with the PDOS of Ru upon H/2H-adsorption in Fig. S13, 14 and Table S23 (ESI†).

OER/ORR activity

We note that the (TM)SA–N_nC_m–G_N catalysts show good performance in the OER/ORR, based on their overpotentials. Here, their OER/ORR performance and their preferred pathways among two different paths (I and II) are investigated and the predicted high-performance OER/ORR materials are assessed.

We consider the catalytic activity in terms of the overpotential of the OER/ORR (η^OER/η^ORR). The HO*/O*/HOO*-adsorption free energies (ΔG_HO*/O*/HOO*) are calculated by adding zero-point energy and entropic contributions to binding energies of each intermediate for path-I and path-II (Tables S24–S27, Methods, ESI†). Then, η^OER/η^ORR values are obtained from the free energies of intermediates during OER/ORR (Methods). The lowest η^OER/η^ORR values between path-I and path-II are reported. The solvent effect is important in particular for d¹ or d² SAs, while it is much less in other cases (Table S28 and S29, ESI†). Here we discuss the activities by taking into account the solvent effect.

In the OER/ORR cases, only several TM(SA)s show good activity regardless of coordination environments (Fig. 8a). Among thermodynamically/electrochemically stable systems, the low η^OER/V values (compared with 0.43 ± 0.09 (upper limit: 0.52) for Pt(111)) are shown for Pt/Ni–N₂C₂ (0.27/0.40), Ni/Pd–C₄ (0.31/0.36), and Rh–N₄ (0.31) (Tables S30b and S31, ESI†). The low η^ORR/V values (compared with 0.47 ± 0.21 (upper limit 0.68) for Pt(111)) are shown for Ir/Rh–N₄ (0.31/0.36), Pd/Ni/Pt/Cd–C₄ (0.34/0.43/0.46/0.48), Ru–N₃C₁ (0.36), Ni/Pd/Pt/Hg/Zn–N₁C₃ (0.36/0.36/0.37/0.48/0.50), Zn–N₂C₂^b (0.48), Cd–C₄ (0.48), Rh–N₃C₁ (0.49), Cd–N₂C₂ (0.55), Mo/Tc/Co–N₄ (0.55/0.55/0.58), Co/Ni–N₂C₂^b (0.60/0.66), Ni/Ir/Ru–N₂C₂^c (0.64/0.64/0.68), and Zn/Pt/Ni/Co–N₂C₂ (0.65/0.65/0.66/0.67) (Tables S30c and S31, ESI†). Our predicted η^ORR[Pt–N₁C₃] = 0.37 V is in agreement with an experimental value (η^ORR = ∼0.3 V) of Pt in G_N as one of possible structures of Pt–N_nC_m.⁵⁶ Fe and Co atoms embedded in N-doped carbon templates were reported to be better in performance than those in pure carbon templates, respectively,^57–60 while Co behaves better than Fe in both N-doped and un-doped carbon templates, respectively,⁵⁷ the results of which are consistent with our predicted values. Co–N₄ is predicted to be the most active for the ORR among Co–N_nC_m–G_N, in agreement with experiments.⁶⁰

Fig. 8 OER/ORR overpotentials (η^OER/ORR) based on the refined values (Table S30 based on Tables S24–S29, ESI†). (a) Color code map for η^OER/ORR and their corresponding potential determining step (“1–8“ in each small box denotes ΔG_1–8, respectively; Fig. 6 & Methods). The lowest η^OER/ORR is reported among path-I and path-II. (b) Color code maps for η^OER < 0.6 V and η^ORR < 0.7 V (Table S30, ESI†). High-performance SA catalysts for the OER are predicted to be Pt–N₂C₂ (η^OER = 0.27 V), Ni–C₄ (0.31 V), Rh–N₄ (0.31 V), Pd–C₄ (0.36 V), and Ni–N₂C₂ (0.40 V), while those for ORR are Ir/Rh–N₄ (η^ORR =0.31/0.36 V), Pd–C₄ (0.34 V), Ru–N₃C₁ (0.36 V), Ni/Pd/Pt–N₁C₃ (0.36/0.36/0.37 V), Pt/Cd–C₄ (0.46/0.48 V), and Zn–N₂C₂^b (0.48 V) and Hg–N₁C₃ (0.48 V).

The HO* and HOO* are adsorbed on metal atoms much more strongly than on the surrounding C/N atoms. Although O* forms a bond with both M and C in certain cases, the active sites for the OER/ORR are metal sites in all investigated systems in contrast to those for the HER, which vary from metal to C/N or metal-C bridging site depending on the metal and ligands. As a result, the η^OER/η^ORR values of each (TM)SA are affected by the coordination environment.

In general, OH⁻ tends to be on top of O* to form HOO* (path I), known as the well-known OER/ORR mechanism. However, we find that in some systems the HO* adsorbs on the O-poisoned metal site (M–O*), which needs to be considered as the active site (path-II). Re/Os prefers path-II in all templates, while V/Nb follows path-II only in –N₄–G_N. Tc/Ru/Ta/W/Mo follows path-II only in –N₄–G_N and –N₂C₂^c–G_N. Also, the Cr/Fe follows path-II only in –N₂C₂^c–G_N. Tc/Ru/Ta/W/Mo follows path-II on most of the double vacancy –N_nC_m–G_N except –C₄–Gr. When O* tends to adsorb strongly on a metal atom, it tends to follow path-II. If the metal atom is more positively charged due to a higher coordination by N atoms, two O atoms can be further coordinated. For example, ΔG_O* on V–N₄ (−1.17 eV) is much stronger than those on V–N₂C₂ (−0.01 eV) and V–C₄ (0.26 eV). The adsorption of OH⁻ on top of O* makes the M–O* bonding weaker (i.e., ΔG_HOO* > ΔG_O*); then, the HOO* formation step along path-I tends to be less favoured. In contrast, the OH⁻ adsorption (beside O*) on the metal atom creates an additional M–O* bond [HO*–(M–O*)]; nevertheless, the repulsion between two O atoms around the metal atom arises. If this repulsion is not large enough, path-II would be favoured, and vice versa. The O* adsorption energy gets affected by the coordination environment around the metal atom. However, the optimal overpotential is governed by the balance among the four e⁻ reaction steps in the OER/ORR, and thus, the delicate preference between paths I and II is not easily describable. The stability of the M–O* site is confirmed by Pourbaix diagrams (Fig. S15–S17, ESI†).

Descriptors for HER/OER/ORR activity

As a supervised machine learning approach to the classification of di-vacancy configurations into good and bad catalysts, we searched for a descriptor taking the properties of atomic elements constituting the configuration, which works universally for all M–N_nC_m–G_N systems.

We considered the following machine learning algorithms: Gaussian naïve Bayes, logistic regression, k-nearest neighbors, radius neighbor classifier, support vector classifier, neural network, and decision trees as well as the following ensemble methods: random forest, extremely randomized trees, and gradient boosting. The best hyper-parameters were found using an exhaustive search implemented in GridSearchCV.⁶¹

We constructed the descriptors as Coulomb matrix elements (CM = Z₁Z₂/r₁₂ where Z = atomic number and r₁₂ = the distance from each other)⁶² between M and neighbouring C/N atoms, between M and adsorbent atoms (H/O), and between adsorbents and the closest non-metal surface atom (C/N), along with atomic numbers, atomic first ionization energy etc. (see Supplementary Note, ESI†). The classified activity data of M–N_nC_m–G_N systems at given criteria (|η^HER| < 0.2 V for HER and η^OER/ORR < 1 V for OER/ORR) are used for training and test. The Coulomb matrix for the distances from the TM to the coordinating C/N atoms and the first ionization potential of TM are found to describe high-performance HER/OER/ORR catalysts with ROC AUC⁶³ of 0.79–0.89 (Fig. S18–S21 and Tables S32–S34, ESI†). The importance of engineering coordination environments of (TM)SA–G_N catalysts is noticed based on the correlation between activity and coordination environment, which is described as a combination of the ionization energy of the metal atom and the Coulomb matrix elements with the ligating atoms.

Conclusion

TM(SA) electrocatalysts of M–N_nC_m–G_N are analyzed to characterize their stability/durability, activity, active sites, electronic properties, σ, and XANES depending on the coordination environment. This diverse coordination-driven information essential for the design of highly efficient SA-driven electrocatalysts is rather barely available in previous studies despite the development of diverse SA-driven electrochemical catalysts. The early (TM)SAs are more likely to be stabilized on the defective site of substrates, while the late (TM)SAs tend to be stabilized against electrolytic dissolution. The (TM)SAs in –N₂C₂ moiety display superior durability and high σ in favor of fast Tafel reaction kinetics. The present results show most suitable coordination environments for stability and activity towards the design of high-performance (TM)SA electrocatalysts. Overall, (TM)SAs including Pt/Pd/Au/Ag/Cu/Ir/Rh/Ni in the –N₂C₂ moiety display superior durability. Thus, instead of –N₄ or –C₄ moieties that have been well studied, the –N₂C₂ moiety is suggested to be highly useful. As for catalytic activity, the –N₂C₂ moiety performs very well in the HER and reasonably well in the OER, while –N₄, –C₄, –N₃C₁ and –N₁C₃ moieties perform well in the ORR. Ti/Zr/Ag/Cd–N₁C₃, Pt/Ni/Ru/Rh–N₂C₂, Cd/Rh–N₄, and Mo/Rh/La/Os/Pt–C₄ are likely to perform well in the HER (−η^HER/V < 0.1). The low η^OER/V values are shown for Pt/Ni–N₂C₂ (0.27/0.40), Ni/Pd-C₄ (0.31/0.36), and Rh–N₄ (0.31). The low η^ORR/V values are found for Ir/Rh–N₄ (0.31/0.36), Pd/Ni/Pt/Cd–C₄ (0.34/0.43/0.46/0.48), Ru–N₃C₁ (0.36), Ni/Pd/Pt/Hg/Zn–N₁C₃ (0.36/0.36/0.37/0.48/0.50), Zn–N₂C₂^b (0.48), Cd–C₄ (0.48), Rh–N₃C₁ (0.49), and Cd–N₂C₂ (0.50). Consequently, the diversity of –N_nC_m moieties can be exploited to achieve better performance for diverse types of electrocatlytic reactions. The calculated K-edge XANES of systems for (TM)SAs and conductivity would be beneficial in experimental characterization of their coordination environments and fast reaction kinetics of TM(SA). The present findings should be vital for the design of new (TM)SA catalysts for various electrocatalytic reactions.

Computational methods

Theoretical calculations

Spin-polarized DFT calculations using the Vienna ab initio simulation package (VASP)⁶⁴ are performed to optimize structures and to obtain electronic properties. Instead of the Perdew–Burke–Ernzerhof (PBE)⁶⁵ exchange–correlation functional (which is commonly used for bulk systems), the revised PBE (RPBE) type⁶⁶ (which was optimized for surface states) is used for these 2-dimensional (2D) systems with Tkatchenko–Scheffler (TS)⁶⁷ dispersion correction. For 3d metals, the Hubbard U_eff correction⁶⁸ (with the values of 2.11, 2.58, 2.72, 2.79, 3.06, 3.29, 3.42, 3.40, 3.87, and 4.12 eV in the order of Sc through Zn) is included based on the literature.³⁹ Given that our focus is more on the surface phenomena, the RBPE(U)+TS method is employed here. The RPBE+TS results without U are also provided for a guide even for 3d elements because the experimental cohesive energies and interatomic distances of TMs tend to be between RPBE(U)+TS and RPBE+TS values (Tables S1 and S2, ESI†). Thus, in the evaluation of η, we may consider both U=0 and U values as the upper/lower limits, or simply consider their mid value with the error bar (ε) as the half of their difference (see the legend of Table S30, ESI†).

Bulk transition metals are optimized with the conjugated gradient algorithm using a Monkhorst–Pack grid of 11 × 11 × 11 special k points and 700 eV cutoff for the kinetic energy and the Gaussian smearing method with an energy width of 0.02 eV. Crystal structures for all bulk transition metals are used for their primitive unit cells and their lattice vectors. Despite the fact that the dispersion correction (which is not present in PBE results) is essential, the TS dispersion correction (PBE+TS) tends to overestimate the experimental TM cohesive energies probably because the parameters of the PBE functional (without considering dispersion correction) were optimized directly for the experimental data. The Hubbard U correction is also necessary for self-energy correction in order to obtain more reliable electronic structure such as a band gap, but it tends to underestimate cohesive energies. Since the parameters of PBE were optimized for bulk materials, RPBE was introduced to give better surface energies.⁶⁹ In this regard, RPBE(U)+TS is expected to better describe the 2D systems including both 3d TMs and π systems (such as Gr and G_N). Thus, we have chosen the RPBE(U) + TS method to study free energies of (TM)SAs embedded in 2D G_N, while taking into account spin polarization. For the study of the catalytic phenomena, the projector augmented wave (PAW)⁷⁰ pseudopotentials are employed and the kinetic energy cutoff for plane-wave expansions is set to 500 eV with Gaussian broadening by 0.05 eV. Each energy of (TM)SA is optimized in the asymmetric supercell (19 Å × 20 Å × 21 Å) for a better description of their degeneracy states using the Γ point. Each isolated (TM)SA is embedded in a Gr supercell composed of 48 C atoms, and 7 × 7 × 1 k-point sampling is used. Geometric relaxations are done until the Hellman–Feynman forces converge to 0.01 eV Å.

Solvent effect

The solvation effect is considered using the pure implicit model implemented in VASPsol⁷¹ and the explicit model (which uses an implicit model in the presence of one water molecule added explicitly to the system). For the HER overpotentials, the solvent effect is found to be insignificant (no more than 0.03 eV for ΔG_H*) (Table S22, ESI†), so it was not necessary to include the solvent effect in HER analysis. In the case of the Mn system, the magnetic moment of the metal atom changes in the implicit solvation model, resulting in 0.09 eV in free energy change, but it does not change in the explicit model, which is more reliable, resulting in an insignificant solvent effect on the ΔG_H*. However, the solvent effect is somewhat significant in the OER/ORR, and so it is included in our results and discussion (Tables S27–S29, ESI†). In a bulk dielectric medium the implicit solvent model would be reliable, while in ~1 nm-sized confined medium the Debye screening is far from complete, and so the dielectric constant is much reduced due to the limited polarization arising from the rotationally forbidden surface water molecules at catalytic sites. Owing to the much reduced dielectric effect of interfacial water, the implicit solvent effect could overestimate the solvent effects on η of (TM)SA catalysts embedded on G_N/Gr. In this regard, we consider both vacuum and implicit solvent models as the upper/lower limits, namely the mid value with the error bar as the half of their difference (see the legend of Table 30, ESI†).

Reaction mechanism and Gibbs free energy calculation

While the HER is a fast reaction, the OER is a sluggish reaction, which generates oxygen gas through several electron/proton complex reaction pathways. The OER consists of four elementary steps, which strongly depend on the pH of the reaction medium. In acidic pH, the reaction operates via oxidation of two water molecules to form one oxygen molecule with four proton-coupled electron transfer steps, while in basic pH, the four hydroxyl groups are converted into O₂ and two water molecules with four electron transfer steps. ORR kinetics is also generally slow. In acidic conditions, O₂ reduces to H₂O with four proton-coupled electron transfer steps, while in basic conditions, O₂ reduces to OH⁻ with four electron transfer steps. The reactions for the HER, OER and ORR in acidic and basic pH are summarized as follows.^31,72,73 The detailed mechanisms we considered are described below.

HER mechanism

HER in acidic conditions with free energy equations:⁵³

* + H⁺ + e⁻ → H* → * + 0.5H₂ (1)

ΔG_H* = (E_H* − E_* − 0.5E_H2) − (ZPE_H* − ZPE_* − 0.5ZPE_H2) − T(S_H* − S_* − 0.5S_H2) (2)

HER in basic conditions is the same as below with the same free energy equation as eqn (2):

* + H₂O + e⁻ → H* + OH⁻ → * + 0.5H₂ (3)

OER and ORR mechanisms

We have considered various OER and ORR mechanisms, as shown below. The overpotential of each mechanism is calculated from the maximum difference between the Gibbs free energies in each step, as provided for path-I in (37) and (38). The lowest overpotential values for the best pathway are reported in Tables S28–S29 (ESI†).

(i) OER mechanism in basic conditions. a. Conventional mechanism on M (path-I):

* + OH⁻ ⇌ HO* + e⁻ (4)

HO* + OH⁻ ⇌ O* + e⁻ (6)

O* + OH⁻ ⇌ HOO* + e⁻ (8)

HOO* + OH⁻ ⇌ * + O₂ + e⁻ (10)

b. Conventional mechanism on O*-M (path-II):

O*–M + OH⁻ ⇌ O*–HO* + e⁻ (12)

O*–HO* + OH⁻ ⇌ O*–O* + e⁻ (13)

O*–HO* + OH⁻ ⇌ O*–HOO* + e⁻ (14)

O*–HOO* + OH⁻ ⇌ O*–M + O₂ + e⁻ (15)

c. Direct O₂ mechanism:

* + OH⁻ ⇌ HO* + e⁻ (16)

HO* + OH⁻ ⇌ O* + e⁻ (17)

O* + OH⁻ ⇌ O*–OH* + e⁻ (18)

O*–HO* + OH⁻ ⇌ O*–O* + e⁻ (19)

O*–O* + e⁻ ⇌ O₂ (20)

(ii) ORR mechanism in basic conditions. d. Conventional (associative) mechanism on M (path-I):

* + O₂ + H₂O + e⁻ ⇌ HOO* + OH⁻ (21)

HOO* + e⁻ ⇌ O* + OH⁻ (23)

O* + H₂O + e⁻ ⇌ HO* + OH⁻ (25)

HO* + e⁻ ⇌ OH⁻ + * (27)

e. Conventional (associative) mechanism on O*–M (path-II):

O*–M + O₂ + H₂O + e⁻ ⇌ O*–HOO* + OH⁻ (29)

O*–HOO* + e⁻ ⇌ O*–O* + OH⁻ (30)

O*–O* + H₂O + e⁻ ⇌ O*–HO* + OH⁻ (31)

O*–HO* + e⁻ ⇌ O*–M + OH⁻ (32)

f. Dissociative mechanism:

* + O₂ ⇌ O*–O* (33)

O*–O* + H₂O + e⁻ ⇌ O*–HO* + OH⁻ (34)

O*–HO* + H₂O + e⁻ ⇌ HO*–HO* + OH⁻ (35)

HO*–HO* + 2e⁻ ⇌ * + 2OH⁻ (36)

Then, the overpotential (η) of OER, and ORR for path-I are calculated as follows:^31,72,73

η^OER = max(ΔG₁, ΔG₂, ΔG₃, ΔG₄)/e − 1.23 V (37)

η^ORR = max(ΔG₅, ΔG₆, ΔG₇, ΔG₈)/e + 1.23 V (38)

The intermediates of the OER and ORR are the same in both acidic and basic conditions and the Gibbs free energies of adsorbates are calculated by adding zero-point energies (ZPEs) and entropic contributions (S) to the adsorption energies. ,, and equations are given for path-I and , , and can be calculated by changing G_* to G_O*–M in each equation, where the TS value of H₂ is 0.41 eV and that of H₂O is 0.67 eV.⁵³ For the active site and adsorbates, their TS values are 0 eV.

where the TS value of H₂ is 0.41 eV and that of H₂O is 0.67 eV.⁵² For active sites and adsorbates, their TS values are 0 eV.

Electrical conductivity and XANES calculation

Electrical conductivities (σ) are calculated using the Boltzmann transport theory⁷⁴ and the scattering rates τ are assumed to be the same as ≈100 fs for G_N.⁷⁵ X-ray absorption near edge structures for all optimized structures are calculated by using the FDMNES code.⁷⁶ Muffin-tin approximation and the Hedin–Lundqvist functional are used to conduct full multiple scattering calculations. The spectra are convoluted using the Lorentzian function.

Author contributions

M.H., D.Y.K. and C.W.M. initiated the project. M.H., D.Y.K., M.U., and C.W.M. performed DFT calculations and analyzed the results. V.G. and M.U. investigated machine learning-based descriptors. M.H., C.W.M. and K.S.K. wrote the manuscript. K.S.K. supervised the project.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by NRF (National Honor Scientist Program: 2010-0020414), A.I. Incubation Project Fund (1.210091.01) of UNIST, and KISTI (KSC-2018-CHA-0057, KSC-2019-CRE-0103, KSC-2019-CRE-0253, KSC-2020-CRE-0049, KSC-2020-CRE-0146 and KSC-2020-CRE-0185). C. W. M., M. H. and D. Y. K. acknowledge the support from KISTI (KSC-2019-CRE-0139, KSC-2020-CRE-0081, and KSC-2021-CRE-0009, respectively). D. Y. K. also acknowledges the NRF grant funded by the Korean Government (2018H1A2A1062086).

References

1. J. Zhu, L. Hu, P. Zhao, L. Y. S. Lee and K. Y. Wong, Chem. Rev., 2019, 120, 851–918.
2. L. Liu and A. Corma, Chem. Rev., 2018, 118, 4981–5079.
3. J. Mahmood, F. Li, S. M. Jung, M. S. Okyay, I. Ahmad, S. J. Kim, N. Park, H. Y. Jeong and J. B. Baek, Nat. Nanotechnol., 2017, 12, 441–446.
4. V. Vij, S. Sultan, A. M. Harzandi, A. Meena, J. N. Tiwari, W. G. Lee, T. Yoon and K. S. Kim, ACS Catal., 2017, 7, 7196–7225.
5. M. Lefèvre, E. Proietti, F. Jaouen and J. P. Dodelet, Science, 2009, 324, 71–74.
6. H. W. Kim, M. B. Ross, N. Kornienko, L. Zhang, J. Guo, P. Yang and B. D. McCloskey, Nat. Catal., 2018, 1, 282–290.
7. E. Jung, H. Shin, B. H. Lee, V. Efremov, S. Lee, H. S. Lee, J. Kim, W. Hooch Antink, S. Park, K. S. Lee, S. P. Cho, J. S. Yoo, Y. E. Sung and T. Hyeon, Nat. Mater., 2020, 19, 436–442.
8. S. Sultan, M. Ha, D. Y. Kim, J. N. Tiwari, C. W. Myung, A. Meena, T. J. Shin, K. H. Chae and K. S. Kim, Nat. Commun., 2019, 10, 5195.
9. M. K. Debe, Nature, 2012, 486, 43–51.
10. Y. Y. Birdja, E. Pérez-Gallent, M. C. Figueiredo, A. J. Göttle, F. Calle-Vallejo and M. T. M. Koper, Nat. Energy, 2019, 4, 732–745.
11. Y. Wu, Z. Jiang, X. Lu, Y. Liang and H. Wang, Nature, 2019, 575, 639–642.
12. J. N. Tiwari, W. G. Lee, S. Sultan, M. Yousuf, A. M. Harzandi, V. Vij and K. S. Kim, ACS Nano, 2017, 11, 7729–7735.
13. B. H. R. Suryanto, H. L. Du, D. Wang, J. Chen, A. N. Simonov and D. R. MacFarlane, Nat. Catal., 2019, 2, 290–296.
14. T. Lim, G. Y. Jung, J. H. Kim, S. O. Park, J. Park, Y. T. Kim, S. J. Kang, H. Y. Jeong, S. K. Kwak and S. H. Joo, Nat. Commun., 2020, 11, 412.
15. Z. Li, Y. Chen, S. Ji, Y. Tang, W. Chen, A. Li, J. Zhao, Y. Xiong, Y. Wu, Y. Gong, T. Yao, W. Liu, L. Zheng, J. Dong, Y. Wang, Z. Zhuang, W. Xing, C. T. He, C. Peng, W. C. Cheong, Q. Li, M. Zhang, Z. Chen, N. Fu, X. Gao, W. Zhu, J. Wan, J. Zhang, L. Gu, S. Wei, P. Hu, J. Luo, J. Li, C. Chen, Q. Peng, X. Duan, Y. Huang, X. M. Chen, D. Wang and Y. Li, Nat. Chem., 2020, 12, 764–772.
16. L. Zhang, K. Doyle-Davis and X. Sun, Energy Environ. Sci., 2018, 12, 492–517.
17. J. N. Tiwari, S. Sultan, C. W. Myung, T. Yoon, N. Li, M. Ha, A. M. Harzandi, H. J. Park, D. Y. Kim, S. S. Chandrasekaran, W. G. Lee, V. Vij, H. Kang, T. J. Shin, H. S. Shin, G. Lee, Z. Lee and K. S. Kim, Nat. Energy, 2018, 3, 773–782.
18. D. Liu, X. Li, S. Chen, H. Yan, C. Wang, C. Wu, Y. A. Haleem, S. Duan, J. Lu, B. Ge, P. M. Ajayan, Y. Luo, J. Jiang and L. Song, Nat. Energy, 2019, 4, 512–518.
19. I. Roger, M. A. Shipman and M. D. Symes, Nat. Rev. Chem., 2017, 1, 0003.
20. A. Wang, J. Li and T. Zhang, Nat. Rev. Chem., 2018, 2, 65–81.
21. J. Li, M. Chen, D. A. Cullen, S. Hwang, M. Wang, B. Li, K. Liu, S. Karakalos, M. Lucero, H. Zhang, C. Lei, H. Xu, G. E. Sterbinsky, Z. Feng, D. Su, K. L. More, G. Wang, Z. Wang and G. Wu, Nat. Catal., 2018, 1, 935–945.
22. X. Wan, X. Liu, Y. Li, R. Yu, L. Zheng, W. Yan, H. Wang, M. Xu and J. Shui, Nat. Catal., 2019, 2, 259–268.
23. H. Xu, D. Rebollar, H. He, L. Chong, Y. Liu, C. Liu, C. J. Sun, T. Li, J. V. Muntean, R. E. Winans, D. J. Liu and T. Xu, Nat. Energy, 2020, 5, 623–632.
24. X. Cui, W. Li, P. Ryabchuk, K. Junge and M. Beller, Nat. Catal., 2018, 1, 385–397.
25. E. D. Goodman, A. C. Johnston-Peck, E. M. Dietze, C. J. Wrasman, A. S. Hoffman, F. Abild-Pedersen, S. R. Bare, P. N. Plessow and M. Cargnello, Nat. Catal., 2019, 2, 748–755.
26. Y. Xiong, J. Dong, Z. Q. Huang, P. Xin, W. Chen, Y. Wang, Z. Li, Z. Jin, W. Xing, Z. Zhuang, J. Ye, X. Wei, R. Cao, L. Gu, S. Sun, L. Zhuang, X. Chen, H. Yang, C. Chen, Q. Peng, C. R. Chang, D. Wang and Y. Li, Nat. Nanotechnol., 2020, 15, 390–397.
27. J. N. Tiwari, A. M. Harzandi, M. Ha, S. Sultan, C. W. Myung, H. J. Park, D. Y. Kim, P. Thangavel, A. N. Singh, P. Sharma, S. S. Chandrasekaran, F. Salehnia, J. W. Jang, H. S. Shin, Z. Lee and K. S. Kim, Adv. Energy Mater., 2019, 9, 1900931.
28. A. M. Harzandi, S. Shadman, M. Ha, C. W. Myung, D. Y. Kim, H. J. Park, S. Sultan, W. S. Noh, W. Lee, P. Thangavel, W. J. Byun, S. Lee, J. N. Tiwari, T. J. Shin, J. H. Park, Z. Lee, J. S. Lee and K. S. Kim, Appl. Catal., B, 2020, 270, 118896.
29. J. N. Tiwari, N. K. Dang, S. Sultan, P. Thangavel, H. Y. Jeong and K. S. Kim, Nat. Sustainable, 2020, 3, 556–563.
30. K. Jiang, S. Siahrostami, T. Zheng, Y. Hu, S. Hwang, E. Stavitski, Y. Peng, J. Dynes, M. Gangisetty, D. Su, K. Attenkofer and H. Wang, Energy Environ. Sci., 2018, 11, 893–903.
31. S. Sultan, J. N. Tiwari, A. N. Singh, S. Zhumagali, M. Ha, C. W. Myung, P. Thangavel and K. S. Kim, Adv. Energy Mater., 2019, 9, 1900624.
32. J. N. Tiwari, A. N. Singh, S. Sultan and K. S. Kim, Adv. Energy Mater., 2020, 10, 2000280.
33. M. D. Hossain, Z. Liu, M. Zhuang, X. Yan, G. L. Xu, C. A. Gadre, A. Tyagi, I. H. Abidi, C. J. Sun, H. Wong, A. Guda, Y. Hao, X. Pan, K. Amine and Z. Luo, Adv. Energy Mater., 2019, 9, 1803689.
34. H. Fei, J. Dong, Y. Feng, C. S. Allen, C. Wan, B. Volosskiy, M. Li, Z. Zhao, Y. Wang, H. Sun, P. An, W. Chen, Z. Guo, C. Lee, D. Chen, I. Shakir, M. Liu, T. Hu, Y. Li, A. I. Kirkland, X. Duan and Y. Huang, Nat. Catal., 2017, 1, 63–72.
35. Y. Zhang, Nat. Rev. Chem., 2018, 2, 0151.
36. Y. Wang, J. Mao, X. Meng, L. Yu, D. Deng and X. Bao, Chem. Rev., 2018, 119, 1806–1854.
37. H. Wang, S. Min, Q. Wang, D. Li, G. Casillas, C. Ma, Y. Li, Z. Liu, L.-J. Li, J. Yuan, M. Antonietti and T. Wu, ACS Nano, 2017, 11, 4358.
38. S. Wang, D. Z. B. Li, C. Zhang, Z. Du, H. Yu, X. Bi and S. Yang, Adv. Energy Mater., 2018, 8, 1801345.
39. H. Xu, D. Cheng, D. Cao and X. C. Zeng, Nat. Catal., 2018, 1, 339–348.
40. Z. W. Seh, J. Kibsgaard, C. F. Dickens, I. Chorkendorff, J. K. Nørskov and T. F. Jaramillo, Science, 2017, 355, eaad4998.
41. J. Greeley, T. F. Jaramillo, J. Bonde, I. Chorkendorff and J. K. Nørskov, Nat. Mater., 2006, 5, 909–913.
42. J. K. Nørskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37–46.
43. H. Shin, H. Xiao and W. A. Goddard, J. Am. Chem. Soc., 2018, 140, 6745–6748.
44. C. Choi, S. Back, N. Y. Kim, J. Lim, Y. H. Kim and Y. Jung, ACS Catal., 2018, 8, 7517–7525.
45. M. Zafari, D. Kumar, M. Umer and K. S. Kim, J. Mater. Chem. A, 2020, 8, 5209–5216.
46. D. Y. Kim, M. Ha and K. S. Kim, J. Mater. Chem. A, 2021, 9, 3511–3519.
47. C. A. Tolman, Chem. Soc. Rev., 2004, 1, 337.
48. D. Deng, X. Pan, L. Yu, Y. Cui, Y. Jiang, J. Qi, W.-X. Li, Q. Fu, X. Ma, Q. Xue, G. Sun and X. Bao, Chem. Mater., 2011, 23, 1188.
49. D. Wei, Y. Liu, Y. Wang, H. Zhang, L. Huang and G. Yu, Nano Lett., 2009, 9, 1752.
50. J. N. Tiwari, Y. K. Seo, T. Yoon, W. G. Lee, W. J. Cho, M. Yousuf, A. M. Harzandi, D. S. Kang, K. Y. Kim, P. G. Suh and K. S. Kim, ACS Nano, 2016, 11, 742–751.
51. J. Greeley and J. K. Nørskov, Electrochim. Acta, 2007, 52, 5829–5836.
52. P. Vanýsek, in CRC Handbook of Chemistry and Physics 87th, ed. D. R. Lide, Taylor & Francis, Boca Raton, 2007, ch. 8, pp. 20–29.
53. J. K. Nørskov, T. Bligaard, A. Logadottir, J. R. Kitchin, J. G. Chen, S. Pandelov and U. Stimming, J. Electrochem. Soc., 2005, 152, J23.
54. G. J. Kubas, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 6901–6907.
55. H. Jin, S. Sultan, M. Ha, J. N. Tiwari, M. G. Kim and K. S. Kim, Adv. Funct. Mater., 2020, 30, 2000531.
56. J. Liu, M. Jiao, L. Lu, H. M. Barkholtz, Y. Li, Y. Wang, L. Jiang, Z. Wu, D.-J. Liu, L. Zhuang, C. Ma, J. Zeng, B. Zhang, D. Su, P. Song, W. Xing, W. Xu, Y. Wang, Z. Jiang and G. Sun, Nat. Commun., 2017, 8, 15938.
57. Y. Zhao, K. Watanabe and K. Hashimoto, J. Am. Chem. Soc., 2012, 134, 19528–19531.
58. W.-J. Jiang, L. Gu, L. Li, Y. Zhang, X. Zhang, L.-J. Zhang, J.-Q. Wang, J.-S. Hu, Z. Wei and L.-J. Wan, J. Am. Chem. Soc., 2016, 138, 3570–3578.
59. P. Yin, T. Yao, Y. Wu, L. Zheng, Y. Lin, W. Liu, H. Ju, J. Zhu, X. Hong, Z. Deng, G. Zhou, S. Wei and Y. Li, Angew. Chem., Int. Ed., 2016, 55, 10800–10805.
60. Q. Yang, Y. Jia, F. Wei, L. Zhuang, D. Yang, J. Liu, X. Wang, S. Lin, P. Yuan and X. Yao, Angew. Chem., Int. Ed., 2020, 59, 6122–6127.
61. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and É. Duchesnay, J. Mach. Learn. Res., 2011, 12, 2825–2830.
62. M. Rupp, A. Tkatchenko, K.-R. Müller and O. A. von Lilienfeld, Phys. Rev. Lett., 2012, 108, 058301.
63. D. K. McClish, Med. Decis. Making., 1989, 9, 190.
64. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 2002, 6, 15–50.
65. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 2002, 77, 3865–3868.
66. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 59, 7413–7421.
67. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
68. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 57, 1505–1509.
69. J. Wellendorff and T. Bligaard, Top. Catal., 2011, 54, 1143–1150.
70. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 50, 17953–17979.
71. K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. A. Arias and R. G. Hennig, J. Chem. Phys., 2014, 140, 084106.
72. J. Rossmeisl, A. Logadottir and J. K. Nørskov, Chem. Phys., 2005, 319, 178–184.
73. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
74. G. K. H. Madsen and D. J. Singh, Comput. Phys. Commun., 2006, 175, 67–71.
75. Y. W. Tan, Y. Zhang, K. Bolotin, Y. Zhao, S. Adam, E. H. Hwang, S. Das Sarma, H. L. Stormer and P. Kim, Phys. Rev. Lett., 2007, 99, 246803.
76. O. Bunău and Y. Joly, J. Phys.: Condens. Matter, 2009, 21, 345501.

---

Footnotes

† Electronic supplementary information (ESI) available: Supplementary note, Fig. S1–S21 and Tables S1–S34 as described in the main text. See DOI: 10.1039/d1ee00154j

‡ These authors contributed equally.

---