Multifunctional Electrocatalysis on Single-Site Metal Catalysts: A Computational Perspective

Abstract

Multifunctional electrocatalysts are vastly sought for their applications in water splitting electrolyzers, metal-air batteries, and regenerative fuel cells because of their ability to catalyze multiple reactions such as hydrogen evolution, oxygen evolution, and oxygen reduction reactions. More specifically, the application of single-atom electrocatalyst in multifunctional catalysis is a promising approach to ensure good atomic efficiency, tunability and additionally benefits simple theoretical treatment. In this review, we provide insights into the variety of single-site metal catalysts and their identification. We also summarize the recent advancements in computational modeling of multifunctional electrocatalysis on single-site catalysts. Furthermore, we explain each modeling step with open-source-based working examples of a standard computational approach.

1. Introduction

From the standpoint of worldwide climate change, where almost every country still depends on fossil fuels, there is a clear need to provide clean and efficient energy production and conversion systems for the coming years. Water splitting electrolyzers, metal-air batteries, and regenerative fuel cells are the key devices of the sustainable and green hydrogen cycle, promising effective use of the growing amount of recurrent energy sources such as solar and wind energy, which are supplied to the electric grid. Efficient electrical energy transformation to fuels is becoming increasingly crucial to smooth green energy’s daily and seasonal oscillation. Thus, all kinds of the above-mentioned electrochemical energy conversion systems are essential for smart-grid management.

Electrochemical water splitting is the simplest, efficient, and cleanest process for high-purity hydrogen generation. Intermittent sources of renewable energy (e.g., solar and wind energy) can convert and store energy as chemicals. The overall water splitting reaction is a sum of electrocatalytic electrode reactions: anodic Oxygen Evolution Reaction (OER) and the cathodic Hydrogen Evolution Reaction (HER).

The electrochemical Oxygen Reduction Reaction (ORR) and the Hydrogen Oxidation Reaction (HOR) are two essential half-reactions in fuel cells. These electrochemical devices can collect, store, and gradually release electrical energy on demand. Regenerative fuel cells offer unique benefits, such as high power density, high specific energy density, high efficiency, long life, good cycling ability, and zero environmental impact, making them a prospective solution for the long-term energy storage and power source in smart-grid and space technology-related applications [1].

Eventually, from an energy storage perspective, the amount of energy stored in secondary batteries is directly related to the extent of active masses in the battery. Using oxygen as an active component at the positive electrode is a viable option since oxygen is a part of the atmospheric gases and does not increase the battery mass [2]. ORR is the key reaction in the secondary batteries in the discharge mode, while OER applies to the charging process of secondary batteries relying on oxygen electrodes.

There are still significant technological challenges along with fundamental studies of these principal electrode reactions such as ORR, OER, and HER. Currently, most of the chemical industry applies precious metal catalysts. So far, Pt metal was reported to be one of the most potent in terms of the intrinsic activity and turnover frequency for HER [3]. One strategy to decrease production costs is to replace the expensive Pt-group metals with a suitable cheaper alternative and/or to utilize the catalyst more effectively in such a way that every catalyst’s atom is an active player.

In this scope, an appropriate catalyst design is still open and needs to be explored in detail before any widespread commercialization appears to be feasible. The modern and economical application requires the use of reasonable cost and performance materials for electrocatalysis. Single-atom catalysts (SACs) are up-and-coming candidates for that purpose because of their tunability, selectivity, and suitability for modeling (Figure 1). SACs comprise isolated single metal atoms dispersed on various substrates, resulting in the most efficient utilization of the active catalytic sites. Most common examples of SACS include heteroatom-doped carbon materials (MNC) and metal-organic frameworks (MOFs) derived catalysts or even pure MOFs. Despite the excellent activity, the major drawback is that these single distributed atoms are not chemically stable in harsh, corrosive, and oxidizing electrocatalytic reaction conditions [4]. The right choice of the suitable substrate stabilizes the catalyst and prevents its aggregation [5].

Bifunctional catalysts are of high priority because of their capability of being used solely and directly for HER and OER as well as ORR and OER. The bifunctional oxygen/air electrodes consume oxygen during electricity generation (e.g., battery discharge, fuel cell operation) and produce oxygen during battery charge or water electrolysis. In water splitting, a bifunctional catalyst can simultaneously conduct both the oxidation and reduction reactions of water, i.e., OER and HER, therefore they can be combined into a single water splitting device, which significantly reduces the cost of the hydrogen energetics cycle [6].

Up to now, the only commercialized bifunctional electrocatalysts capable of both ORR and HER are Pt-based materials. However, due to their cost and limited supply, many efforts were made to reduce the quantity of Pt used in electrocatalysts, but at the same time increase their performance [4,7]. As a substitute to the noble metal catalysts, bifunctional SACs can be used, and they are currently intensively explored to enhance their intrinsic activity for water splitting [5,8,9,10], ORR and OER in metal-air batteries applications [11,12,13,14].

The first part of this work focuses on an overview of the SACs as bifunctional electrocatalysts and their experimental characterization. Second, a review of recent advances in modeling multifunctional electrocatalysis (HER, OER, ORR) on SACs is presented. Finally, an example of a readily accessible computational workflow with free and open-source software is provided, and the remaining challenges and prospects are addressed.

2. Experimental Insights

2.1. Different Bifunctional Single-Site Metal Catalysts

A wide array of catalyst classes is proposed to enhance the kinetics of ORR/OER/HER processes. However, the development of materials with bifunctional catalytic activity is still an active area of research. To achieve efficient bifunctional electrocatalytic performance, earth-abundant transition metals and nitrogen are typically incorporated into a graphitic framework to form MNC catalysts [15,16]. MNC catalyst usually comprises the isolated transition metal M (M = Fe, Co, Ni, Mn, etc.) atom coordinated by carbon (MC_x) and nitrogen (MN_x) [17,18]. These isolated sites optimize the adsorption energy of electrochemical reaction intermediates, which reduces the reaction energy barrier for ORR, HER, and OER [19,20]. Beyond the high activity of MNC sites, the performance of the materials also depends on the carbon material conductivity and porosity, which significantly affects the mass and electron transport. Therefore, the design of the transition metal-containing N-doped mesoporous carbons with a maximum distribution of exposed active centers is a promising route to produce efficient bifunctional electrocatalysts for ORR/OER and water splitting processes.

Single-atom catalysts have maximum metal dispersion and the highest atom-utilization efficiency, enabling the rational and economical utilization of metal resources [21,22]. To date, numerical chemical methods to prepare single-site metal catalysts toward HER, OER or ORR were proposed [23]. In sharp contrast, research on the exploration for particularly bifunctional electrocatalysis is still at an early stage [24]. This section focuses on general synthetic routes and on the analysis of the electrochemical performance of bifunctional single-atom MNC electrocatalysts.

2.1.1. MNC-Derived Bifunctional Single-Site Metal Catalysts

Based on synthesis precursors, single-atom MNC electrocatalysts can be roughly divided into two categories. One of them is MNC-derived SACs typically prepared by immobilizing transition metal and nitrogen atoms onto high surface area carbon supports (Figure 2) [25]. Different carbon materials were employed as supports for atomically dispersed transition metals [26]. One of the most suitable support for single-atom catalysts is defective graphene, where single atoms are trapped in defective sites of heteroatom-doped graphene [27]. In defective graphene, carbon atoms are localized near divacancies (carbon defects), and by strong interaction with single-atom MC_x sites are formed [28]. For example, Zhang et al. synthesized NiNC bifunctional catalyst by doping defective graphene with nitrogen and Ni salt precursors and obtained material with remarkable water splitting electroactivity (overpotential (η) of 0.07 V at 10 mA cm⁻² for HER and 0.27 V at mA cm⁻² for OER in 1 M KOH solution) [29]. It was demonstrated that among transition metals anchored on hollow N-doped graphene, the electrocatalytic ORR performance follows the sequence Fe > Co > Cu > Ni [30]. On the other hand, Fei et al. reported that the OER activity of various SA-doped graphene samples follows the order: Ni > Co > Fe [31].

Wang et al. prepared a highly active bifunctional ORR/HER single-site metal electrocatalyst by carbonization of graphdiyne with melamine and cobalt nitrate hexahydrate [32]. Qin et al. utilized Fe-containing ionic liquid for a self-sacrificing, template-assisted, and controlled pyrolysis to obtain active atomic Fe sites containing FeN₄ coordination and Fe³⁺ species anchored on hollow carbon nanospheres [33]. Very recently, a series of metal SACs anchored on nitrogen-doped carbon aerogel was reported by Cheng et al. [34]. Their study indicated that Co atoms facilitated the formation of pyrrolic-N sites, thus resulting in enhanced HER activity, while Ni atoms promote the production of pyridinic-N sites, which accelerate the OER.

Shang et al. synthesized manganese SAC by co-pyrolysis of manganous nitrate and chitosan, which delivered excellent bifunctional oxygen electroactivity ΔE of 0.67 V, where ΔE = E_j10 − E_1/2, E_j10 is the OER potential at a current density of 10 mA cm⁻² and E_1/2 is ORR half-wave potential [35]. By pyrolyzation of polyvinylpyrrolidone–Fe(NO₃)₃ nanocrystals with melamine, Du and co-workers prepared highly active and stable FeNC bifunctional SAC material with abundant accessible active centers and favorable mass transport channels (ΔE = 0.80 V) [36].

Although the monometallic single-atom catalysts have shown excellent activity for ORR, their OER and HER electrocatalytic performance are still unsatisfactory. Therefore, incorporating additional active sites into the N-doped carbon matrix could be a potential solution to endow the electroactivity of bifunctional catalyst materials [37,38,39,40,41]. Recently, Chen et al. reported incorporating Ni/Fe metal cations into defect-rich 3D carbon nanosheets via a modified salt-sealed strategy [42]. Dual atom Ni/Fe SAC displayed high catalytic activity for ORR and OER, with a ΔE value of 0.69 V in 1 M KOH solution. The introduction of a secondary metal atom alters the coordination environment of the active center. Moreover, an additional catalytic center located nearby the primary single atom site modulates an adsorption–desorption strength for different ORR/OER intermediates, such as OH*, O*, and OOH* [41].

One of the main challenges in MNC-derived SACs is the aggregation of metal into particles caused by a sharp increase in surface energy during atomic dispersion [43]. MNC-based SACs still suffer from aggregation caused by the lack of strong interaction between metal and substrate atoms. Another challenge that hinders the large-scale application of MNC-based SACs is low metal loadings (generally < 1 wt.%), which are caused by the insufficient surface area of carbon supports.

2.1.2. Metal-Organic Framework-Derived Bifunctional Single-Site Metal Catalysts

Porous solid materials, called metal-organic frameworks (MOFs), are another promising class of SAC precursors [44]. Mononuclear metal complexes are utilized as starting materials in the typical synthesis route of MOF-based single-site catalysts (Figure 3). The atomic dispersion and stabilization of metal sites are achieved by removing the organic ligands, typically by carbonization at high temperatures [45]. The most convenient precursors for mononuclear metal complexes are materials with uniform and regular pore structures, such as zeolites and MOFs [46,47,48,49].

These porous catalysts supporting materials also stabilize metals against migration and their aggregation into nanoparticles [50,51]. The recent progress and advances in the fabrication methods of MOF-derived SACs for electrochemical applications were recently reviewed by several research groups [52,53].

To date, numerous MOF-based MNC single-atom bifunctional catalysts were prepared and investigated for application in ORR/ORR/HER [54]. For example, Zhang et al. surveyed a series of CoNC catalysts synthesized by self-sacrificial templating of benzene-1,3,5-tricarboxylate linker-based MOFs [55]. According to their DFT investigation, the CoN_2–4 were primary active catalytic sites responsible for ORR activity. Zhang et al. employed direct carbonization of ZIF-67/graphene oxide precursors to prepare a highly active (ΔE of 0.67 V) CoNC bifunctional electrocatalyst with a high surface area and a hollow capsule-like nanostructure [56]. Chen et al. utilized the KCl-template method to synthesize a ZIF-67-derived CoN_x catalyst with intrinsic oxygen electroactivity: ΔE = 0.63 V (ORR E_1/2 = 0.91 V, and OER E_j10 = 1.54 V) [57]. ZIF-67 precursors are also frequently used for the fabrication of bifunctional water splitting CoNC catalysts [13,57], where a combination of CoN_x with enhanced porosity and the presence of metallic nanoparticles synergistically enhanced the HER/OER performance.

Lian et al. recently introduced a new method of preparing SACs without the carbonization step [58]. They used Co²⁺ and hexaiminotriphenylene (Co₃HITP₂) based conductive MOF directly as a catalytic material and revealed that remarkable ORR/OER activity can be ascribed to a presence of metal centers with unpaired 3d electrons.

Atomic Fe species with abundant Fe–N₄ sites were successfully synthesized by pyrolysis of amino-substituted Zr-MOF and phthalocyanine iron precursor [59]. Optimized catalyst demonstrated excellent oxygen electroactivity (ORR: E_1/2 = 0.89 V and OER: E_j10 = 1.58 V). Pan et al. reported a novel polymerization–pyrolysis–evaporation method to prepare atomically dispersed Fe sites on N-doped porous carbon [60]. During the carbonization of Zn/Fe polyphthalocyanine precursor, the evaporated Zn²⁺ centers resulted in the formation of FeN₄ sites, which granted exceptional bifunctional ORR/OER activity (ORR E_1/2 = 0.89 V, OER E_j10 = 1.66 V). Luo et al. obtained FeN₄ and CoN₄ single-atom sites by pyrolysis of ZIF-8/67 precursor supported by SiO₂ nanospheres [61]. Promising bifunctional electroactivity was assigned to the electronic interaction between metal atoms, whereas Co modified FeN₄ and Fe modified CoN₄ were assigned as active ORR and OER sites, respectively. It was experimentally revealed that Fe single-atom sites, which are located close to Co sites, facilitate reactant adsorption and charge transfer on Co-active sites [39]. Zhu et al. synthesized atomically dispersed Fe–Ni pairs embedded in N-doped carbon hollow spheres [62]. According to their studies, the charge redistribution between Fe and Ni facilitated the adsorption of reaction intermediates, and the electronically tuned Fe and Ni atoms were responsible for ORR and OER, respectively.

Although the research of MNC and MOF-derived SACs is at a very early stage, it has already indicated very promising results. A better understanding of the synthesis processes will enable the fabrication of advanced SACs with unique structures and excellent multifunctional electroactivity.

2.2. Experimental Way of Identifying Single-Sites

The presence of single-atom metal sites on the catalyst surface can be characterized by electron microscopy techniques, such as high angle annular dark-field scanning transmission electron microscopy and scanning tunneling microscopy [17,19,23].

The most frequently used tool to confirm the presence of MN_x sites is X-ray photoelectron spectroscopy. However, for the rational design of efficient bifunctional SACs, it is essential to get direct insights into the chemical states of active sites and the evolution of atomic structure in realistic operation conditions.

Operando X-ray absorption spectroscopy is a powerful technique that provides valuable information for fundamental electrocatalysis [63]. Extended X-ray absorption fine structure technique enables monitoring of the local coordination environment of the active centers, and X-ray absorption near-edge structure is used to monitor the dynamic electronic structure of catalysts [64,65]. Shang et al. performed operando X-ray absorption spectroscopy measurements to monitor the local atomic coordination evolution of MnN₂ characteristics during the ORR and OER processes [35]. They observed that the average Mn oxidation state decreased from 2.9 to 2.2 during ORR and increased from 3.1 to 3.8 during the OER process. The real-time extended X-ray absorption fine structure measurements revealed that the real catalytic forms for ORR and OER were bond-length-extended Mn²⁺−N₂ and bond-length-shortened Mn⁴⁺−N₂, respectively.

Although physical characterization techniques are rapidly developing, it is still a big challenge to directly observe the atomic structures of SACs using just one method [66]. The combination of several imaging techniques with density functional theory (DFT) is necessary to verify the structure of SACs and to reveal the exact mechanism of the electrocatalytic process.

3. Computational Background

Computer modeling is a means of mirroring the real world. On the one hand, such reflection aims at providing a fundamental understanding of electrochemical reactions. For a clearer image, the experiment should be performed at precisely known conditions, e.g., using well-defined electrodes and extra pure electrolytes and reagents. At the same time, the simulation should be run using a realistically large interfacial model and for a sufficiently long time to observe the flow of electrochemical reactions with atomistic resolution. At the time being, such experiments and computations require unprecedented effort to be performed, so there are very few reactions that can be precisely characterized and compared at empirical and theoretical levels [67]. One of such reactions is HER, which stands still as a state-of-the-art research object [3,68,69,70,71,72,73]. On the other hand, modeling can outline the real world recognizable enough to provide guidelines for designing novel catalysts. The difference between the two approaches is in the detailing. Let us highlight three directions for development as shown in Figure 4: System complexity (type, size, time, and theory level), reaction mechanism (species and pathways), and free energy landscape sampling (thermodynamics and kinetics).

3.1. Free Energy Landscape

Electrochemical reactions proceed on the free energy landscape, described with minima and maxima corresponding to chemical species and transition states, respectively. An example of the landscape for the ideal and realistic catalyst is shown in Figure 5. Energy differences between minima are thermodynamic quantities, while differences between minima and transition states play a very crucial role in kinetics.

A minimum can be defined by the adsorption free energy of the corresponding chemical species, which is related to the reference-free energy of reagents and products:

ΔG_A = G_A* − G_* − G_A,

(1)

where the free energy of each chemisorbed species is estimated as

G_A* = E_A* + ZPE − TS – NeU,

(2)

where E is the potential energy obtained in calculations, ZPE is zero-point energy correction, TS term is calculated using vibrational analysis and statistical thermodynamics, N is the number of H⁺ + e⁻ pairs involved in the reaction (see below), e is the elementary charge, and U is the potential relative to the reversible hydrogen electrode (RHE) [74].

For instance, the adsorption free energies of OH, O, and OOH are defined relative to the reagents and products as:

ΔG_OH* = G_OH* − G_* − (G_H₂O − ½ΔG_H₂),

(3)

ΔG_O* = G_O* − G_* − (G_H₂O − ΔG_H₂),

(4)

ΔG_OOH* = G_OOH* − G_* − (2G_H₂O − 3/2ΔG_H₂),

(5)

At zeroth-order approximation, the adsorption energies are correlated with catalytic activity. Here, the Sabatier Principle is commonly applied to distinguish catalysts that bind a specific chemical species neither too strong nor too weak [75]. For example, Trassati’s volcano plot relates the HER exchange current density with the hydrogen adsorption energy on the metal surface. Similar volcano plots provide a helpful guideline with a single descriptor [76] as illustrated in Figure 6.

At first-order approximation, the free energy landscape is sampled by accounting for several selected minima. The energy differences are equalized to the reaction overpotentials. The classic example is OER and ORR presented as a (reversed) sequence of four reactions forming a pathway from H₂O to O₂ and back:

H₂O + * ⇌ OH* + (H⁺ + e⁻), ΔG₁ = G_OH*,

(6)

OH* ⇌ O* + (H⁺ + e⁻), ΔG₂ = G_O* − G_OH*,

(7)

H₂O + O* ⇌ OOH* + (H⁺ + e⁻), ΔG₃ = G_OOH* − G_O*,

(8)

OOH* ⇌ * + O₂ + (H⁺ + e⁻), ΔG₄ = G_O₂ − G_OOH*,

(9)

so that the potentials and overpotentials are defined in regard to the minimum and maximum free energy differences among ΔG₁, ΔG₂, ΔG₃, and ΔG₄:

U_ORR = min(ΔG_1–4)/e and U_OER = max(ΔG_1–4)/e,

(10)

η_ORR = 1.23 V − U_ORR and η_OER = U_OER − 1.23 V,

(11)

If only the spacing on the energy diagram were exactly 1.23 V, the overpotential would be zero. However, the adsorption free energy is interrelated by linear scaling relations:

G_A = γ_ABG_B + ξ_AB,

(12)

where slope γ_AB is proportional to the valencies of the species A and B [77,78], and intercept ξ_AB is in addition dependent on the adsorbate coordination [79]. Many of the studied classes of materials have the same slope, but different intercepts, thus showing distinct catalytic activity.

Under these assumptions, scaling relations set thermodynamic limits for overpotentials. For metals, oxides, and MNC materials, linear scaling relation ΔG_OOH* − ΔG_OH* = 3.2 ± 0.2 V [80,81] holds. Thus, η_ORR/ORR = 3.2/2 − 1.23 = 0.37 V, and the theoretical window for the bifunctional catalysts is 0.37 × 2 ≈ 0.74 V due to the scaling relations. Strategies for both utilizing and breaking scaling relations are commonly used to design electrocatalysts.

At second-order approximation, a more advanced free energy landscape sampling is utilized to formulate microkinetic models. Here, the focus is on the energy barriers that define the reactions’ rates, e.g., via transition state theory [82]. Although more advanced approaches can be used [67]. The Brønsted–Evans–Polanyi (BEP) relationship is commonly utilized for relating the energy barriers to the adsorption energies as:

ΔG_A^# = αΔG_A + β,

(13)

where α is the proximity factor ranging from 0 to 1 and β is an intrinsic reaction barrier [83].

3.2. Reaction Pathways

An optimal trajectory over the free energy landscape towards the desired products implies avoiding too low minima and too high maxima with correction for the applied potential. The variety of possibilities includes electron transfer (ET), proton transfer (PT), and concerted proton–electron transfer (CPET), among others, with possible pathways for OER and ORR shown in Figure 7. Within most computational approaches, it is most convenient to study CPET routes. That excludes ionic species from consideration, allows focus on the covalent-type bond energies, and is compatible with computational hydrogen electrodes.

3.3. Electrode–Electrolyte Interface Model

Besides the free energy landscape, the catalytic reactions are naturally happening in real space, namely, at interfaces. In electrocatalysis, the primary interface is between a solid electrode and an aqueous electrolyte. It is commonly referred to as the electrical double layer (EDL). The EDL sets up a field that affects the adsorbate. Thus, the more realistic the EDL model is, the more reliable description of the reaction can be obtained. The most rigorous treatment of that effect would involve an atomistic simulation of the two electrode models. Such simulations are standard for classical molecular dynamics at ideally polarizable electrodes with the recently developed contact potential method [84,85,86].

While constant potential simulations of reactions are not feasible, a Born–Haber cycle-based approach is highly effective in understanding trends in electrocatalysis. That allows for defying the absolute potential for reversible reactions taken as a reference. Figure 8 illustrates the cycle used for the most common reference—a generalized computational hydrogen electrode.

3.4. Computational Hydrogen Electrode

At the heart of most of the computational investigations of electrocatalytic activity lies the reference electrode, which was formulated by Norskov et al. in 2004 [87] and generalized in subsequent works by Rossmeisl et al. [71,88,89]. We briefly describe the computational hydrogen electrode (CHE) and its implications on calculations.

The CHE model assumes the following: (1) The interface is in equilibrium with the electrode and electrolyte; (2) The energy of the model is independent of the electrostatic field and boundary conditions. More generally, the interface must be charge neutral and large enough to screen all charges. Next, the CHE model relates the free energy of (H⁺ + e⁻) reaction to the 1/2H₂ energy. This assumption works only in cases when we consider the concerted proton–electron transfer (CPET) mechanism in each individual step [90], which was shown to be the most prevalent mechanism [91]. Notwithstanding, it is vital to keep in mind corner cases, and as was shown by Koper [90], in case of weak oxygen adsorption, the decoupling proton–electron pathway is more significant, indicating that using CHE to model such reactions is not appropriate. Experimentally, this pathway was confirmed for some metal oxides [92], observable as pH-dependent OER activity on the RHE scale [90]. However, as we will note later, the weak oxygen adsorption is undesirable for low overpotentials; CPET is a generally appropriate model.

3.5. Breaking Scaling Relations

Historically, a trial and error approach was used in searching for catalysts. However, recently efficient screening methods were developed utilizing the reactivity descriptors and scaling relations.

A variety of adsorbed intermediates at the catalyst surface is a complex system in heterogeneous catalysis and can be challenging to model in a complete ensemble using theoretical approaches. Although theoretical work provides deep insight into the many surface reactions and the adsorption energies of reaction intermediates and their assignee species, it can predict significant catalytic performance trends [93].

The scaling relation is a popular method, where linear relationships between adsorption energies are related to chemical species. Scaling relations allow predictions of adsorption energies of several species’ adsorption energies on a given surface based on the adsorption energy of a single species on that surface.

Craig et al. examined 17 of the most active molecular OER catalysts, based on different transition metals (Ru, Mn, Fe, Co, Ni, and Cu) containing MOFs, and demonstrated that they obey similar scaling relations to those established for heterogeneous systems. The authors also found that the conventional OER descriptor underestimates the activity of some very active OER complexes because the standard approach neglects the crucial one-electron oxidation step that many molecular catalysts undergo before O–O bond formation. This additional oxidative step descriptor allowed specific molecular catalysts to circumvent the “overpotential wall”, leading to enhanced performance in OER [94].

On the other hand, Huang et al. state that the adsorption free energy of intermediates, including OOH*, O*, and OH*, exhibit a strong linear correlation and cannot be optimized independently on a single type of catalytic site. Authors provide certain guidelines for the rational design of efficient catalysts for ORR and OER beyond the limitation presented by scaling relations, which includes: (a) introducing p states; (b) introducing a second adsorption site; (c) introducing a proton acceptor group; (d) the strain effect; (e) changing the solvent composition; (f) nanoscopic confinement; (g) O–O direct coupling in the absence of OOH* species [95].

In the recent studies, it was shown that the overpotential for ORR of 0.3–0.4 V (explained above) can be reduced below that limit V (for catalysts with protonable group such as hangman porphyrins) [96] or even below 0.2 V (for diporphyrins also known as pacman) [97,98]. In the computational models, the scaling relation is broken by (de)stabilizing one of the reaction intermediates.

Even though scaling relations are important for the development of efficient catalysts, blind breaking of the scaling relationships does not guarantee better catalytic activity, as was noted by Gogvindarajan et al. [99], also highlighting the importance of optimizing catalysts first of all for optimal performance with existing scaling relations, and only then performing more manageable breaking of the scaling relations.

3.6. Computational Models

In the following section, we provide an overview of computationally studied structures for single-atom catalysts. This includes usual MNC type structures, which attempt to model metal and nitrogen doping in graphene, additionally, we also discuss further modifications made to the MNC type structures and modeling of MOF-based catalysts.

3.6.1. MNC Type Structures

First studies of MNC type materials for OER and ORR date back to 2011, when Calle-Vallejo et al. [100] investigated a range of graphitic materials with active sites composed of four nitrogen atoms and various transition metal atoms corresponding to structures shown in Figure 9a,b. The study revealed multiple significant results; first of all, concerning scaling relations, it was found that trend-wise MNC catalysts have the same behavior as metal oxides, so the same scaling relations bind their catalytic activity. Nevertheless, it was noted that since active sites are formed in the interstices of the graphitic layers, the porous structure could facilitate additional 3D interactions, which could improve ORR and OER activities. Since graphitic materials could be more stable in an acidic environment and are generally more conductive than their oxide counterparts, further investigation of MNC catalysts is warranted. Second, regarding catalytic properties, it was observed that with an increase in the number of d-electrons of the transition metals, the interaction strength with adsorbates in the active sites decreases. The optimum adsorption strength for catalysis was noted and found out that for ORR, efficient catalysis can happen with Fe, Ir, Mn, Ru, and Rh, and for OER on Co, Rh, and Ir, indicating those transition metal atoms belonging to groups 7 to 9 are promising catalysts towards both OER and ORR, hence allowing bifunctional catalysis. Compared to metal oxide catalysts, multiple non-noble options were observed, showing that MNC materials are promising candidates for cheaper catalysts. Furthermore, spin analysis was performed to show that in general, the metal oxidation state in most active sites is +2.

More recent studies investigated more detailed aspects of MNC catalysts, such as the effect of metal oxidation states and the environment of active sites. Amiinu et al. [101] synthesized a Co-based MNC catalyst with the remarkable bifunctional activity of ΔE ≈ 0.65 V and rationalized the observations with DFT calculations. In this study, Co center was considered in both pyridinic and pyrrolic nitrogen environments. Besides the free energy calculations based on the CHE model we described earlier, activity was explained with the help of the analysis of the frontier orbitals and projected density of states. It was found that the Co atom in the pyridinic environment induced higher HOMO energy (−4.647 eV, as compared to −4.765 eV in the pyrrolic environment), which suggested more favorable activity towards ORR.

Additionally, the projected density of states analysis revealed that Co in a pyrrolic environment-induced localized electronic states in both up-spin and down-spin channels with low density of states at the Fermi level, which suggested lower frontier electron density and activity when compared to the pyridinic environment. Effect of the ease of Co²⁺ oxidation on the metal surface towards bifunctional activity was also noted, indicating that Co²⁺ plays a significant role for ORR, while Co³⁺ is active for OER (see also [102]). As a conjecture, it was suggested that the electronegativity of oxygen coupled with the excellent electron affinity of N could also introduce positively charged C-atoms as active sites by OER. Such a mechanism with adjacent C participation was observed in computational studies by Fei et al. [31] for OER reaction on Ni MNC catalyst and labeled as a dual-site mechanism. The comparison between single-site mechanism (adsorption only on metal center) and dual-site (adsorption on metal center and adjacent C atom) for Ni MNC showed the smallest limiting barrier difference of 1.24 eV for single-site mechanism and 0.42 eV for dual-site mechanism, displaying that depending on the considered mechanism, Ni MNC is better OER catalyst than Co MNC, for which single-site mechanism is more favorable. As such, it was found that the catalytic activity and reaction pathways depend strongly on the metal center in MNC, which further allows decoupling of ORR and OER reactions for bifunctional catalysts, as ORR always proceeds via a single-site mechanism, but OER can take place via double-site mechanism under suitable conditions.

Investigation on the effect of catalytic performance, depending on N-dopant concentration and configuration, was performed by Zhang et al. [103], who calculated the activity of CoN_1–4 catalysts. It was found that the overpotential for both OER and ORR decreases with more nitrogen atoms in the catalytic center, therefore CoN₄ exhibits the lowest overpotential. This was further rationalized with a density of states calculations, showing that hybridization interaction between Co-d and O-p states is stronger on CoN₄ when compared to other catalysts with lower nitrogen concentrations. As a further extreme of nitrogen doping concentration, Cai et al. [104] investigated 3d transition metal doping of pure graphene materials (that is MeN₀ type). In this case, Ti and Cu doped graphene were predicted to have excellent ORR catalytic performance, while Ni and Cu doped graphene catalysts were predicted to be good OER candidates, overall indicating that in case of no nitrogen doping, Cu doped graphene is identified as bifunctional electrocatalyst. However, we note that the observed overpotentials are significantly higher than those calculated for MeN₄ systems.

Having confirmed that most active sites in MNC type materials are found in MeN₄ type sites, more recent studies have investigated the effect of the surrounding environment of the active site and its role on the catalytic activity. Therefore, Li et al. [12] investigated the effect of the location of the ZnN₄ site including those shown in Figure 9a,b,e,f. In this study, they found out that the ZnN₄-pyridine configuration has the lowest overpotential (0.61 V) for ORR. Meanwhile, ZnN₄-pyrrole and ZnN₄-edge show lower overpotentials (0.73 V and 0.63 V, respectively) for OER. As such, bifunctional activity can be controlled with the closest environment of the MeN_x site. Furthermore, an AIMD simulation was performed to investigate the stability of various configurations, which found that the formation energies decreased with the increase in ZnN coordination numbers for all calculated structures, indicating that they are more active due to their stability and also more likely to be observed in experiments. Instead of considering different MNC type site locations, it is also possible to consider adjacent defects in graphene, as discussed by Zhang et al. [59] and shown in Figure 9d, who considered the effect of structural defects that are formed by the mesoporous framework. In this case, they constructed structures with “small pores” and “big pores” to observe the effect of the adjacent graphene defects on OER and ORR activity. It was found that due to the high electronegativity of structures with the largest pores, the defects form highly active ORR sites and that the interfacial charge redistribution at the interface structure also simultaneously facilitates OER in addition to the highly active ORR site. Another type of defect that was studied by Ji et al. [11] for MeN₄ structures is the formation of small Co metal clusters on top of the active CoN₄ site. In this study, OER and ORR activity was compared for freestanding CoN₄ site, and CoN₄ site with Co cluster above it, and the DFT calculations revealed that the single-atom sites have lower adsorbed OH hydrogenation barrier, compared to the state with Co cluster, and therefore exhibiting significantly better OER and ORR performance.

Moreover, as a last structural consideration, we mention the engineering of MeN₄ sites on carbon nanotubes with an example shown in Figure 9c, which can coexist with more conventional MNC sites on 2D graphene. In a study by Ban et al. [40], the activity of CoN₄-tube and CoN₄-sheet was compared and they found that CoN₄-tube reduces the electron density of CoN₄ atoms during the adsorption of OH, and therefore leading to accelerated reaction to form O*, hence facilitating ORR. Meanwhile, the transformation from O* to OOH* occurs more easily on CoN₄-sheet, facilitating OER. Therefore, it was found that the synergistic effects established for CoN₄-tube and CoN₄-sheet catalysts can increase the bifunctional catalytic activity.

3.6.2. Modified MNC Structures

Due to the versatile nature and the ease of MNC type system modification, further extensions of the usual systems were studied. In this case, we will consider three variations of modifications that have been investigated for MNC materials.

Starting with the simplest modification, additional co-doping using S or P atoms, catalytic activity was studied by Chen et al. [15]. It was discovered that co-doping of FeN_x/C with S atoms yields an extremely low ORR overpotential of 0.24 V, which based on computational studies can be attributed to the uneven charge distribution due introduction of S and N atoms. Further studies on co-doping were undertaken by Sun et al. [105], who considered doping with P atoms. Corroborating previous results on doping, the DFT calculations revealed that the doping alters the charge density distribution on the active Fe sites and leads to improved ORR and OER activity. Furthermore, the study notes that since C vacancies are more favorable near the FeN₄ centers, it provides easier doping close to the active site.

Further modifications of MNC systems include the addition of neighboring structures that induce additional electronic effects. Studies of this type of modification were reported by Wan et al. [63], who considered metal porphyrin deposition on graphitic-oxide nanosheets, and by Liu et al. [38] who considered the effect of MoC cluster near the active CoN₄ site (CoNC-MoC). With porphyrin on graphitic-oxide nanosheets, increased ORR and OER activity was explained in terms of optimized electronic structure and local coordination environment due to the interaction between metal porphyrin molecules and the graphitic-oxide nanosheet. However, on the CoNC-MoC catalyst, an additional "single-site double adsorption" mechanism becomes possible, which accounts for facilitating both ORR and OER performances simultaneously. The effect of MoC also includes the modulation of the d-band center for the active Co site.

As a significantly promising approach to increase the bifunctional activity has emerged bimetallic catalysts. In this case, two MeN₄ motifs are modified and put together to form one single-site motif. In particular, Li et al. [106] compared the activity of isolated FeN₄ and CoN₄ with combined FeCoN₆ (called FeN₄/CoN₄). First, the formation of the bimetallic FeCo system was calculated to be more favorable than the formation of individual Fe and Co sites. Therefore, it is expected to form more easily. Second, increased activity towards both ORR and OER was observed for the FeN₄ active site when adjacent to CoN₄ compared to the individual FeN₄ site. This increase in activity is attributed to the slightly lower adsorption energies for all intermediates due to the synergistic contribution. Further study of the bimetallic site synergy was carried out by Lin et al. [107] for the same FeCo bimetallic system, this time on nanotubes. For the developed material, the synergy was observed employing electron distribution calculations, which revealed stronger electron redistribution in the FeCo system when compared to single-doped CoN₄ or FeN₄ systems. Furthermore, when considering the density of states and d-band centers of the bimetallic system, two new peaks emerged around −2.7 eV, corresponding to the 3d orbitals of Fe and Co atoms due to their bonding. As such, the increased bifunctional activity of bimetallic MNC systems can be attributed to more favorable electronic interactions.

Nevertheless, in the case of bimetallic systems involving the Zn center, Wang et al. [108] reported different observations. First of all, the formation energies of bimetallic catalysts were predicted to be higher than that for their SAC counterparts, and the activity of both single-metal and bimetallic models involving Zn was poor, indicating their unsuitability for bifunctional catalysis. However, it was shown that involving electron-withdrawing groups (e.g., hydroxyl groups) in the ZnCoN₆ system improves the bifunctional activity due to adjusting the position of ΔG_O* relative to ΔG_OOH* and ΔG_OH*.

3.6.3. MOFs

Although MNC type materials are the dominant single-site metal atom catalysts, a more convenient approach involves using pure MOFs as a catalyst. Unfortunately, very few studies have reported MOFs with multifunctional catalytic characteristics comparable to that of MNC materials.

As one example prediction of a bifunctional active catalyst, we mention work by Mao et al. [109], who considered a series of TM-DBQ-CP MOFs with TM = Ti, V, Cr, Mn, Fe, and Co centers. Out of all metals considered, only MOF with Ti center displayed promising catalytic activity with OER overpotential of 0.41 V and ORR overpotential of 0.46 V. In this study, the predicted activity was highly correlated with d-band center (ε_d), showing that for optimal OER, the ε_d should be around 0.2 eV, meanwhile for optimal ORR it should be around 0.84 eV. Further investigation of multifunctional catalytic MOFs was carried out by Wang et al. [110], discussing Fe-PTC as a promising trifunctional catalyst with overpotentials of −0.23, 0.89, and 0.44 V for HER, OER, and ORR, respectively. Furthermore, it was also calculated that the H₂O₂ production process competing with ORR has a lower overpotential of 0.22 V, thus further hindering the ORR. Nevertheless, it provided a conclusion that with respect to HER, the active site is the S atom adjacent to the metal center, instead of the TM center for all considered transition metals in TM-PTC.

4. Computational Modeling Example

As a result of the development of freely accessible software and improvement in modern computer processing power, computational modeling of single-site catalysts has also evolved and become easier to implement. Here, we provide an example for the use of ASE [111] and GPAW [112] modeling of MeN₄ type catalysts based on python scripts described in Supplementary File S1. We envision that a modified application of the provided example is attainable for researchers with a brief background in computational modeling, and can assist in further development and comparison of single-atom catalysts.

In the example we employ the following DFT setup: PBE functional [113], TS09 dispersion correction [114], 2 × 2 × 1 k-point spacing, 0.18 Å grid spacing, force minimization below 0.01 eV/Å, implicit solvation using continuum solvation model using water solvation parameters [115] as implemented in GPAW with transition metal vdW radii [116] specified manually, vibrational analysis performed as implemented in ASE.

The first stage in the modeling is to obtain a reasonable model of the catalyst surface. For this, we start with a simple sheet of graphene and then substitute two C atoms with the metal center, and furthermore change the four surrounding C atoms with nitrogens. For the constructed structure, we perform optimization of the geometry to yield our model of the catalyst surface. The structures considered are shown in Figure 10.

As the next stage, after having obtained the model catalyst surface, we add adsorbates on the top of the metal atoms. Thus, considering ORR and OER processes, we consider intermediates such as OOH*, O* and OH* directly on top of the metal. Further, for the dual-site mechanism occurring on NiN₄, we also consider the adsorption of O* and OH* intermediates on the pyridinic C atom. Constructed adsorbates are shown in Figure 11.

With the absorbates discussed, we optimize the constructed structures to obtain the DFT energies of our reaction intermediates. Further, we perform vibrational analysis for each of the adsorbates to obtain ZPE, entropy and Cv corrections. The vibrational analysis also allows us to determine surface structures that are not yet in a stable optimized state, by observing imaginary frequencies. In the case of imaginary frequency, we take a new structure from the vibrational mode with the highest imaginary frequency and repeat the optimization and vibrational analysis until no imaginary frequency is present.

As the last step, we consider the small molecules involved in the elementary steps, H₂ and H₂O, by placing them in the same-sized cell as in the case of the surface model developed and optimized, and gather vibrational frequencies.

Having optimized all structures and obtained their DFT energies and applying corresponding corrections, we can construct the free-energy diagram and determine the overpotentials for each of the catalytic reactions at the corresponding catalyst we modeled. The raw calculated DFT energies and corrections are shown in Table S1 and all final optimized structures are compiled in an ASE Database S1. The constructed diagram is shown in Figure 12, and from it, we can extract the corresponding results from Section 3.6.1.

First of all, observing states for potential U = 0 V, we observe a general decrease in adsorption strength as the number of d electrons on the metal center increases. A more interesting result is observable from comparing single-site and dual-site mechanisms on NiN₄ catalyst surface. From first sight, it seems that the dual-site mechanism provides lower overpotential for the OER process, as the predicted overpotential for the dual-site mechanism is 0.57 V, while for the single-site mechanism it is 0.83 V. Nevertheless, an important observation, is that the predicted adsorption strength of the OH* intermediate is stronger on metal-center, hence it is more likely to observe the single-site mechanism in practice. With respect to efficient bifunctional catalysis, we can generalize this observation further, and link it with scaling relations. If we consider single-site and double-site mechanisms, which both share the same OOH* adsorbate state, we expect the scaling relation to hold for the single-site mechanism, and so for that pathway, the overpotential for the OER process can be improved to approximately 0.37 V. It would seem from first sight, that if the double-site mechanism breaks the scaling relation, it could provide lower overpotential, and hence increased bifunctional activity. However, since *OOH state is shared, the breaking of scaling relation would necessarily imply, that the *OH state for double-site adsorption is higher in energy than the corresponding state for single-site adsorption, and therefore in practice the limiting reaction will still be the single-site adsorption.

5. Perspectives

Single-site metal catalysts have proven to be a breakthrough in providing multifunctional electrocatalysis to substitute more expensive Pt-group metal-based catalysts. Their main strengths arise from their intrinsic ability to be easily modified and tuned according to the reaction needed. Nevertheless, many significant problems remain to be addressed both in the experimental and theoretical studies of such catalysts. The most significant experimental complications arise from the lack of understanding of the synthesis effect on the final structure and the final composition of the catalyst structure. As a result, even if theoretical studies suggest a particularly promising catalyst structure, there is no clear route for the synthesis of these exact surface structures that guarantee the formation of the desired catalyst. Therefore, significant underexplored research lies in identifying the details of processes happening during real synthesis. Main techniques to address the formation of single-site catalysts include advancing in situ imaging and coupling it with molecular dynamics simulations.

From the currently examined structures, it becomes apparently clear that simple catalysts featuring only single-sites have intrinsic limitations for bifunctional ORR and OER catalysis, which from a theoretical perspective, arises due to the nearly similar adsorption of OOH* and OH* intermediates, leading to the scaling relationships. In this regard, the main advantage of single-site catalysts—their ease of tuning, should be exploited more, to decouple the adsorption of the species concerned. Even so, the more complicated structures to consider lead to increased complexity in understanding experimental outcomes. To address the upcoming challenges, the comparative analysis of experimental and theoretical results is increasingly more and more important.

Notable steps in consolidating experimental and theoretical results were carried out by combining DFT predicted EXAFS data with experimental measurements to confirm the existence of distinct MeN₄ type sites in some catalysts. However, other combinations of the agreement are still lacking. As an important direction, we highlight microkinetic modeling, which would allow for a more systematic understanding of the kinetic process mechanism occurring during catalysis and furthermore would allow us to theoretically predict more complicated experimental results, including rotating-disk electrode experiments. Starting from an earlier stage, we identify that the formation of catalysts during synthesis can be investigated with larger-scale molecular dynamics simulations, e.g., applying simulated annealing.

In terms of the theoretical approach, most predictions depend on the accuracy and setup of DFT calculations. The most commonly employed CHE model provides a good starting point, but for finer details, it lacks the ability to deal with pH dependencies and constant potential simulations. To amend the pH dependency the application of GCHE extension appears to be a significant breakthrough that requires further research and application. Nevertheless, the constant potential simulations are at large an unsolved problem, partially due to being computationally expensive and partially due to being in general a less explored method.

Despite the multiple avenues of possible cutting edge research, the conducted studies and developed open-source packages, including catalysts-oriented databases (e.g., CatApp [117], CMR [118] katlaDB [98], OC20 [119],) provide a great background and base for further work in exploring the single-site atom catalysis as a substitute to more expensive materials and yields a promising path towards better sustainable energetics.