Critical assessment of theoretical modelling of single-atom catalysts

Abstract

The rational design of catalysts from first principles remains a central but still elusive goal of modern quantum chemistry. Although advances in computing power and electronic-structure methods have made predictive modelling increasingly feasible, accurately forecasting catalytic activity from theory alone remains challenging. Single-atom catalysts (SACs), with their more defined active sites, offer in principle simplified models compared with conventional heterogeneous systems, yet in many cases discrepancies between theoretical predictions and experimental results persist. Using the hydrogen evolution reaction (HER) as a prototypical case, this study analyses the limitations of current computational approaches, particularly those based on the computational hydrogen electrode (CHE), and the problems arising when theoretical predictions are compared with experimental evidence. Factors contributing to these discrepancies include the sensitivity of reaction thermodynamics to the local atomic environment, the often-unknown experimental structure of SACs, and the neglect of reaction intermediates that form on SACs and not on extended metal surfaces. Additional challenges arise from solvent effects, catalyst evolution under operating conditions, and the potential instability of theoretically designed materials under realistic electrochemical environments. Furthermore, intrinsic approximations in density functional theory introduce uncertainties that hinder quantitative accuracy. Overall, the CHE model, while valuable for identifying general trends, does not include many critical terms that contribute to the activity of untested catalysts. Progress toward the true rational design of catalysts will require integrating these chemical complexities and uncertainties, potentially through artificial intelligence and data-driven methods, to develop more robust descriptors and predictive frameworks.

---

1 Introduction

One of the most challenging goals of computational chemistry is the design in silico of new materials with unprecedented properties. In the field of catalysis, this implies the ability to predict the activity (and the stability!) of a novel catalyst in a given chemical reaction. Over the past 15 years, there has been an explosion of computational studies devoted to single-atom catalysts (SACs), systems where the active sites correspond to highly dispersed transition-metal (TM) atoms anchored on solid supports.^1–3 SACs are particularly attractive for several reasons: (a) every metal atom is isolated and exposed as an active site, meaning 100% utilization of precious metals (e.g., Rh, Pd, Pt); this drastically reduces cost and waste; (b) single atoms exhibit different electronic configurations compared to bulk metals or nanoparticles; this can create unusual reactivity and selectivity not accessible in conventional catalysts; (c) since the active sites are better defined and more uniform, SACs can achieve higher selectivity in chemical transformations, avoiding side reactions common with nanoparticles; (d) the coordination environment (support interactions, neighbouring atoms, defects) can be precisely engineered to control catalytic performance; (e) properly anchored single atoms (on supports like graphene, metal oxides, MOFs) can resist sintering or leaching better than nanoparticles, improving catalyst lifetime; (f) finally, SACs combine advantages of homogeneous catalysts (single sites, high selectivity) with heterogeneous catalysts (stability, reusability, easy separation).⁴

Another key reason for the strong interest of theoretical chemists in SACs is that they provide ideal model systems for studying catalytic mechanisms at the atomic scale, enabling direct correlations between structure and reactivity. This is especially true for SACs in which TM atoms are stabilized on two-dimensional supports such as graphene, MoS₂, etc. In such cases, modelling is greatly simplified: the bulk properties of the material need not be reproduced, and a single atomic layer is often sufficient to capture the relevant chemistry. This structural simplicity, compared to the complexity of supported metal nanoparticles, has spurred a surge of first-principles studies aimed at screening and designing new catalysts. By reducing computational complexity while retaining essential chemical accuracy, SACs serve as powerful platforms for theory-driven catalyst discovery.

Another factor driving the growth of computational studies in this field is the (apparent) relative simplicity of key reactions central to the energy transition: CO₂ and N₂ reduction, and above all, water splitting, where the hydrogen evolution (HER), oxygen evolution (OER), and oxygen reduction (ORR) reactions play pivotal roles. The popularity of these topics is closely linked to the computational hydrogen electrode (CHE) model, an elegant and powerful framework introduced by Nørskov and co-workers⁵ that greatly simplifies the theoretical treatment of electrochemical processes. Thanks to this approach, electrochemical reaction energetics can be evaluated without explicitly modelling the electrode potential or the solvated proton, which dramatically reduces computational complexity while preserving predictive power. The impact of the CHE model is evident: more than 100 000 DFT studies on the HER alone have been published to date, of which around 5000 focus on SACs.⁶

The large volume of computational studies has naturally produced a high number of predicted catalytic systems. The key question, however, is how many of these predictions are actually confirmed by experimental validation. To explore this point, we focus on a specific case of SACs, TM atoms stabilized in nitrogen-doped graphene (TM@N-Gr).⁷ Graphene itself is chemically inert but possesses excellent electrical conductivity. When doped with nitrogen, it provides strong anchoring sites capable of stabilizing isolated metal atoms. Moreover, different nitrogen configurations, pyridinic or pyrrolic, create distinct electronic environments, making nitrogen-doped graphene (N-Gr) a versatile support. The literature already reports hundreds of DFT studies on TM@N-Gr applied to the HER.

The number of experimental studies on TM@N-Gr is considerably smaller. Nevertheless, we identified at least a dozen cases where the activity of TM@N-Gr SACs has been reported. We discuss the applied overpotential (Table 1) at 10 mA cm⁻² as a proxy to quantify the catalytic activity. This is a common practise, although one should keep in mind that a proper treatment would require the study of the dependence of the measured current density on the applied overpotential. Nevertheless, the analysis of the overpotential at a fixed current density provides, in principle, a basis for assessing the reliability of theoretical predictions. For the sake of comparison, Table 1 lists the computed overpotentials corresponding to the experimentally known systems. These values were obtained at the DFT level using the PBE functional (see below for further details), which remains the most widely adopted in this field. For the structural model we worked under the assumption that the TM is incorporated into a pyridinic cavity. Both the choice of functional and the assumed coordination environment introduce uncertainties in the reported values, as we will discuss in detail in the rest of the paper. Nevertheless, pyridinic coordination is by far the most common structural model employed in theoretical studies of SACs based on N-Gr.

Table 1 Measured versus calculated overpotentials, η, for the HER on various TM@4N-Gr systems. The calculated overpotentials correspond to Gibbs free energies of H adsorption on the TM site by considering the formation of the MH intermediate (η^MH) (see next section for further details). Measured overpotentials are taken at current density equal to 10 mA cm⁻²

	Exp η^exp/mV	Ref.	DFT(PBE) η^MH/mV
V@N-Gr	Unreactive	8	160
Co@N-Gr	230	9	130
Co@N-Gr	30	10	130
Ni@N-Gr	530	9	1650
Rh@N-Gr	29 (acid)–45 (alkaline)	11	280
Ru@N-Gr	32	12	460
Ru@N-Gr	34	13	460
Ru@N-Gr	11 (acid)–54 (alkaline)	14	460
W@N-Gr	590	9	1050
Pt@N-Gr	30	15	1630
Pt@N-Gr	47	16	1630
Ce@N-Gr	≈40	17	270
N-Gr	508	16	1630

---

Experimental overpotentials for the HER have been reported for V, Co, Ni, Rh, Ru, W, Pt, and Ce @N-Gr.^8–17 In one case, even the activity of the N-Gr support itself has been investigated.¹⁶ For Co@N-Gr, different overpotentials have been reported, 230 mV⁹ and 30 mV,¹⁰ probably due to voltammograms being obtained in H₂SO₄ 0.5 M at different scan rates, 10 versus 2 mV s⁻¹, respectively. Overpotentials close to 0 mV correspond to active catalysts, and vice versa. The comparison between theory and experiment, however, proves disappointing. With the exception of Co, and, to some extent, Ce, for which DFT(PBE) calculations predict relatively high activity, the results diverge significantly. The most striking discrepancies concern Pt and Ru: both are experimentally highly active, yet DFT calculations classify them as completely inert. Conversely, V is predicted to be an efficient catalyst but has been found experimentally to be inactive. For Ni and W, large quantitative differences emerge between computed and measured overpotentials (Table 1).

The focus of this article is to understand the origin of these significant deviations. We will argue that applying the CHE approach to SACs without accounting for additional key factors can lead to misleading predictions of catalytic activity. These limitations come on top of well-known challenges in modelling solid systems, including uncertainties in identifying the true active site, possible structural and chemical evolution of the catalyst under operating conditions, stability issues, and the inherent variability introduced by different exchange–correlation functionals. Only by carefully addressing all these effects together can reliable predictions be made regarding the structure and activity of a given catalyst.

2 Modelling the hydrogen evolution reaction

2.1 The computational hydrogen electrode (CHE) model

In electrochemistry, the overpotential (η) of a reaction is defined as the extra potential (voltage) that must be applied to an electrode, beyond the thermodynamic equilibrium potential, in order to drive the electrochemical reaction at a desired rate:

η = U_onset − U_eq (1)

where U_onset is the actual electrode potential under operating conditions and U_eq is the equilibrium (or Nernst) potential, which is the thermodynamically required potential for the reaction to occur without kinetic hindrance. The CHE approach⁵ sets the chemical potential of a proton–electron pair (H⁺ + e⁻) equal to half the chemical potential of an H₂ molecule under standard conditions (T = 298 K, p(H₂) = 1 bar):

µ(H⁺ + e⁻) = 1/2 µ(H₂) (2)

This provides a convenient computational reference for electrochemical reactions involving protons and electrons, such as the HER. The hydrogen evolution reaction (in acidic medium) is:

2 H_(aq)⁺ + 2e_(s)⁻ → H_2(g) (3)

Using the CHE, instead of modelling solvated protons and the electrode explicitly, one references the free-energy change of adding one (H⁺ + e⁻) couple to half the energy of H₂. Thus, the free-energy change for hydrogen adsorption, H*, on a catalyst surface is written as:

ΔG_H = ΔE_H + ΔE_ZPE − TΔS_H (4)

where ΔE_H is the DFT-calculated adsorption energy of H, and the reference for the proton–electron pair is 1/2 E(H₂). In this way, the H adsorption energy becomes the key descriptor of the activity of a catalyst in the HER. In particular, if ΔG_H = 0, there is no thermodynamic barrier for the HER, and the catalyst is expected to be optimal, with its ΔG_H at the apex of the classical volcano curve. Catalysts on the left-hand side (strong-binding region) have limited activity because intermediates bind too strongly, poisoning the surface; catalysts on the right-hand side (weak-binding region) have low activity because the intermediates bind too weakly to stabilize transition states or form reactive species. In this framework, only species at the apex of the volcano correspond to an optimal binding energy where the trade-off is balanced, leading to the highest catalytic activity, following the Sabatier principle.

To simulate the electrode potential U, one shifts the free energy of each transferred (H⁺ + e⁻) couple by −eU. This makes it possible to compute reaction free-energy diagrams for the HER and other electrocatalytic reactions from DFT data without directly simulating an electrochemical interface. For a reaction involving n proton–electron coupled transfer steps, the onset potential is given by

U_onset = max{ΔG_i,…, ΔG_n}/e (5)

This is referred to as the potential determining step (PDS). Determining η and the PDS allows one to predict the performance of a given electrochemical reaction, such as the HER. In the following we discuss more in detail the factors that affect the calculation of ΔG_H, the key quantity in the CHE model of HER.

2.2 Details of the calculations

For the results reported in this paper, we performed spin-polarized density functional theory (DFT) calculations as implemented in the Vienna Ab-initio Simulation package (VASP).^18,19 If not differently mentioned, the Perdew–Burke–Ernzerhof (PBE) exchange and correlation functional was adopted.²⁰ The effect of core electrons on the valence electron density was accounted for by means of the projector augmented wave (PAW) method.²¹ The valence electron density was expanded on a plane-wave basis set with a kinetic cutoff equal to 400 eV. The threshold criteria for electronic and geometry loops were set to 10⁻⁵ eV and 10⁻² eV Å⁻¹, respectively. A 2 × 2 × 1 Monkhorst–Pack k-point grid²² was used to sample the reciprocal space. Grimme’s D3 parametrization scheme was adopted to account for dispersion interactions.²³

The N-Gr support was modelled starting from a graphene nanosheet, with a (4 × 4) supercell, where a carbon divacancy was created and four carbon atoms were replaced by four nitrogen ones. The metal atoms have been incorporated in the resulting cavity. The atomic coordinates were fully relaxed. Artificial interactions within the slab plane, along with their effects on the electronic structure and atomic relaxations, are minimized, though not entirely eliminated, by the 4 × 4 unit cell. A vacuum layer of 15 Å was added to avoid spurious effects due to interaction between periodic replicas of the system.

3 Results and discussion

3.1 Role of the local coordination

In the Introduction, we have shown that the use of ΔG_H as a proxy of the activity of a catalyst in the HER does not produce reliable results when one compares the theoretical values with the limited existing experimental data (Table 1). Several reasons can contribute to this discrepancy, but we will show in the following that the application of the CHE approach for SACs has intrinsic limitations that can be surmounted if one carefully considers other contributions and possible situations arising during the reaction.

A major challenge in comparing theory with experiment in heterogeneous catalysis lies in the fact that the true nature of the active site is rarely known with certainty. This has nothing to do with the CHE model, yet it remains a critical aspect. In recent years, several advanced techniques have been employed to achieve atomistic characterization of the local environment around active species. These include microscopy methods such as high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) and aberration-corrected STEM, as well as spectroscopic approaches such as X-ray absorption spectroscopy (XAS), X-ray photoelectron spectroscopy (XPS), electron paramagnetic resonance (EPR), and infrared (IR) spectroscopy with probe molecules.^24,25 Despite these advances, the precise coordination of a single atom on a solid support often remains unresolved, complicating direct comparisons with theoretical models. This is even more challenging if one considers that the observed catalytic activity often arises from the combination of different factors due to a non-homogeneous dispersion of the active phase.²⁶

Focusing on the case of N-Gr, the TM is coordinated to four nitrogen atoms (4N-Gr), typically forming a square-planar or slightly distorted geometry. Two distinct arrangements are possible in this configuration: pyridinic and pyrrolic. In both cases the TM is coordinated to four N atoms, but in the pyridinic structure the N atoms are at the apex of a six-membered ring of C atoms, while in the pyrrolic form they are part of a five-membered pentagonal ring. Other coordination motifs are also plausible, in which the metal atom is bound to only three or two nitrogen atoms. Experimental measurements may in fact probe a mixture of such sites, adding further complexity to the interpretation of results.

Even when restricted to the square-planar cases (pyridinic vs. pyrrolic), the structural differences are sufficient to influence catalytic performance. Ye et al.²⁷ systematically compared these two arrangements, demonstrating that they can yield distinct values of ΔG_H (Table 2). For some TMs (e.g., Co, Fe, Ni), ΔG_H remains nearly identical across both structures. However, for others (e.g., V, Ti, Mn, Cu, Rh), the difference can exceed 0.3 eV. Within the CHE framework, this translates into dramatic differences in catalytic activity. For example, V@4N-Gr in the pyrrolic configuration exhibits an overpotential as low as 16 mV, qualifying as an excellent catalyst, whereas the pyridinic counterpart performs poorly with an overpotential of 387 mV (Table 2).²⁷

Table 2 DFT Gibbs free energies of H adsorption (ΔG_H, eV) on SACs with pyridine-4N-Gr and pyrrole-4N-Gr coordination. From ref. 27

SAC	Pyridine-4N-Gr	Pyrrole-4N-Gr
Co	0.147	0.185
Cr	0.290	0.462
Cu	1.533	1.913
Fe	0.440	0.420
Mn	0.335	0.611
Ni	1.394	1.353
Sc	−0.394	−0.314
Ti	−0.662	−0.313
V	−0.387	−0.016
Rh	−0.307	−0.547
Ru	−0.606	−0.524

---

Even more important is the role of the coordination sphere if one replaces the N ligands with other light atoms, such as C or O. In a computational experiment,²⁸ 24 different TMs were incorporated into 4N-Gr and their catalytic activity toward the HER was evaluated. Then, a Pt atom was chosen as a representative metal centre and embedded in the same 4N-Gr framework but the N atoms surrounding Pt were then systematically substituted with C or O atoms, generating more than 20 distinct structures. Thus, in the first case, the HER activity was assessed as a function of the TM by keeping the support fixed, while in the second, the metal centre was held constant while varying its local coordination environment (Fig. 1).

Fig. 1 Free energy (ΔG_H) of H adsorbed on: (top) transition metal (TM) atoms stabilized on 4N-Gr, TM@4N-Gr; (bottom) Pt@C_xN_yO_z-Gr (Pt stabilized in 4N-Gr where C and O atoms are replacing some of the N atoms). The results are plotted in order of decreasing ΔG_H. Reproduced from ref. 28 (distributed under the terms of the Creative Commons CC BY license).

These two strategies reveal an important insight: comparable variability in catalytic performance can arise either from changing the active metal centre or from altering the chemical nature of its surrounding environment. In particular, ΔG_H was observed to shift from +2.1 eV (Au@4N-Gr) to −1.0 eV (Hf@4N-Gr) as the TM was varied, with Co, V, and Rh exhibiting especially low overpotentials, an outcome consistent with the wide variability of TM properties. Similarly, when changing the nearest neighbours, Pt@C_xN_yO_z-Gr displays ΔG_H values ranging from +1.5 eV to −1.5 eV, covering a span of 3 eV, comparable to the range obtained by changing the TM (Fig. 1). These results highlight that similar reactivity trends can be achieved either by varying the TM while keeping the support fixed, or by fixing the metal centre and tuning the support through local doping.²⁸

These results clearly demonstrate that the specific coordination environment of the TM strongly influences its chemical properties.²⁹ Even subtle variations in the surrounding atoms can lead to significant shifts in the calculated overpotentials.

The role of the support in determining the activity of supported TM atom is shown also by the analysis of bonding and activation of CO₂ to nine TMs (Fe, Co, Ni, Ru, Rh, Pd, Os, Ir, Pt) embedded on three different supporting materials: nitrogen-doped graphene (4N-Gr), the Au(111) surface, and titanium nitride (TiN).³⁰ The O–C–O bond angle provides a particularly robust descriptor of CO₂ activation: the greater the deviation from 180°, the stronger the activation. In general, CO₂ activation on SACs is challenging, as only a limited number of metal atoms can effectively bind and activate the molecule. SACs supported on 4N-Gr and Au(111) displayed poor binding ability, whereas TiN emerged as a highly promising support, markedly enhancing CO₂ activation (Fig. 2). Notable trends include the following: Ru and Os can activate CO₂ on all supports; Rh and Ir are inactive on 4N-Gr but become active when stabilized on Au(111) or TiN (in either Ti or N vacancies); Fe, Co, Ni, Pd, and Pt show activity exclusively when supported on TiN.³⁰ These findings highlight the strong synergistic effect between the TM atom and its surrounding environment in governing the overall catalytic activity.

Fig. 2 OCO bond angle of CO₂ adsorbed on TM@N-Gr, TM@Au(111) and TM@TiN (TM = Fe, Co, Ni, Ru, Rh, Pd, Os, Ir, Pt). The deviation from 180° is a proxy of the extent of CO₂ activation. Reproduced from ref. 30 (distributed under the terms of the Creative Commons CC BY license).

To summarize, in this section we have provided strong evidence that the catalytic activity of a SAC depends delicately on the local structural arrangement, the atoms directly bound to the SAC, and the electronic nature of the support. In this sense, there is no truly “isolated” single-atom catalysis, since chemical activity arises from the interaction of the active site with its surroundings. A meaningful comparison between theory and experiment is only possible when the theoretical models used in simulations closely resemble the actual experimental structures. At the same time, even the prediction of highly active catalysts can be difficult to translate into practice, as it requires controlling the outcome of synthetic procedures with atomistic precision.

3.2 What is the active site?

SACs can be supported on a wide range of materials, including oxides, sulfides, MXenes, nitrides, and metals.³¹ In certain cases, the term single-atom alloy is used to describe systems where a TM atom substitutes host metal atoms within a metallic support.³² In such situations, it is expected that the TM impurity is not the sole active site for stabilizing reaction intermediates. For example, we have shown that hydrogen atoms can bind with comparable strength to neighboring atoms surrounding the TM site in both metals and TiN surfaces.^33,34 Similarly, on oxide supports, hydrogen atoms can interact with surface oxygen atoms to form hydroxyl (OH) groups. This leads to the occurrence of direct and reverse hydrogen spillover phenomena,³⁵ resulting in a complex and dynamic behavior of surface species.

An important question is whether the same effect occurs in SACs stabilized on carbon-based matrices, such as graphene or doped graphene. Some studies have suggested that this may indeed be the case.²⁷ To investigate this, we computed the free energy of H adsorption on the bare support (4N-Gr) and on two SACs, Ti@4N-Gr and Rh@4N-Gr (see SI, Tables S1–S5 and Fig. S1–S20). We then compared the binding strength of H at different surface sites (Fig. 3). Specifically, in addition to the TM site, we examined adsorption on the nearest-neighbor N atom and on three C atoms located in the second and third coordination shells (Fig. 3).

Fig. 3 Various adsorption sites for H atoms near the TM in TM@4N-Gr.

Hydrogen binds strongly to the N site of the bare 4N-Gr support, with a free energy of adsorption ΔG_H = −1.63 eV (Table 3). This strong binding is not surprising due to the undercoordinated nature of the N atom (see also Fig. S1). In contrast, all C sites exhibit large positive ΔG_H values. Overall, DFT calculations predict that 4N-Gr would require a prohibitive overpotential of 1630 mV for the HER, according to the CHE model. Experimentally, an overpotential of 508 mV has been reported for N-doped graphene,¹⁶ which corresponds to an inactive catalyst. Direct comparison between theory and experiment is challenging, however, because the precise atomic structure of the material is not well defined. Moreover, processes such as electrochemical oxidation or defect formation can generate new active sites during operation. In such cases, the observed exchange current may reflect not only the intrinsic activity of the pristine carbon matrix but also that associated with potential-driven structural modifications and degradation.

Table 3 DFT Gibbs free energies of H adsorption (ΔG_H, eV) on different sites of SACs with pyridine-4N-Gr structure

	4N-Gr	Ti@4N-Gr	Rh@4N-Gr
TM-site	—	−0.85	−0.29
N-site	−1.63	0.82	Unbound
C-site	1.55	0.38	1.31
C1-site	1.31	0.36	0.82
C2-site	0.94	0.89	0.74

---

Next, we analyzed the case of Ti@4N-Gr (Table 3). In this system, H binds preferentially to the Ti site, with ΔG_H = −0.85 eV. There is no competition with the N or C sites of the surrounding matrix, as ΔG_H at those sites remains positive. Interestingly, the C and C1 sites near the Ti active centre exhibit significantly reduced ΔG_H values compared to the pristine 4N-Gr matrix (decreasing from about 1.3–1.5 eV in 4N-Gr to around 0.4 eV in Ti@4N-Gr), highlighting the strong electronic influence of the Ti dopant on the surrounding atoms.

A particularly interesting case is Rh@4N-Gr. It has been proposed that, on this SAC, H binding at the Rh site competes with binding at the adjacent N site, with both exhibiting nearly identical theoretical overpotentials.²⁷ However, our new results do not support this conclusion. We investigated the binding behavior of Rh@4N-Gr toward H and found that H indeed binds to Rh with ΔG_H = −0.29 eV (in excellent agreement with the previously reported value of −0.31 eV, ref. 27) (Table 3). In contrast, no local minimum is observed for H adsorption at any of the four N sites. When H is initially placed near one of these sites, H spontaneously relaxes toward the Rh centre. It is possible that in ref. 27, the calculation got stuck in an unphysical local minimum.

This finding is significant, as it indicates that in carbon-based matrices functionalized with TM atoms, the TM itself serves as the primary active centre. While it cannot be ruled out that in some cases the matrix might exhibit higher activity than the TM, such scenarios for carbon-based supports are likely exceptions rather than the norm.

3.3 HER mechanism on SACs: binding multiple H atoms

When a proton is reduced on a metal electrode, the reaction is:

H⁺ + e⁻ + * ⇄ H* (6)

where * is the adsorption site, and H* denotes an adsorbed H atom on the surface (Volmer step). In order to form H₂, the reaction can follow two different paths. According to the Heyrovsky mechanism, a second proton is reduced on the same metal site:

H* + H⁺ + e⁻ → * + H₂ (7)

According to the Tafel mechanism, two H atoms bound to two different * sites diffuse and recombine:

H* + H* → 2* + H₂ (8)

Once H₂ is formed, it interacts only weakly with the surface through dispersion forces. In this physisorbed state, its energy is nearly identical to that of the separated * + H₂ system. Consequently, on metal surfaces, this step can be neglected in the overall thermochemistry of the reaction. The overpotential can therefore be described solely by ΔG_H, regardless of whether the reaction proceeds via the Volmer–Heyrovsky or Volmer–Tafel pathway. This principle also underpins the CHE model: on metal surfaces, only isolated H* intermediates are relevant, since once H₂ forms, it readily desorbs into the gas phase. Thus, ΔG_H is the key descriptor of the thermodynamic barrier of the process.

For SACs, however, the situation differs. SACs share important analogies with organometallic complexes. In their elegant study, Ye et al.²⁷ demonstrated that the ΔG_H values calculated for TM@4N-Gr pyridine and TM@4N-Gr pyrrole SACs show a linear correlation with those of the corresponding TM-phthalocyanine and TM-porphyrin complexes. This correlation arises from the similar nitrogen coordination environments present in both solid-state and molecular systems.

In coordination chemistry, the existence of stable dihydride, H*H*, and dihydrogen, (H₂)*, complexes, where two hydrogen atoms rather than just one are bound to the TM, is well established.^36,37 By tuning the ligands around the metal centre, the degree of H₂ activation can be modulated, yielding either a dihydride complex or the less conventional dihydrogen complex, in which the H–H bond remains only partially activated (typical H–H distance: 0.8–0.9 Å). In a recent study,³⁸ we demonstrated that analogous H*H* and (H₂)* intermediates can also form on SACs. Most of the SACs examined indeed stabilize dihydride or dihydrogen species, thereby introducing an additional reaction step in the HER that is not present in the classical CHE framework. The appearance of these new intermediates has direct relevant implications. Instead of a single step governed solely by ΔG_H, two distinct steps, the formation of H* and the formation of a hydrogen complex (either dihydrogen or dihydride), must now be considered before H₂ release. As a result, catalytic activity is more accurately represented by a three-dimensional volcano plot, as the number of intermediates increases from one to two.³⁸

Notably, the ability of SACs to bind more than one hydrogen atom has also been observed in other theoretical studies.³⁹ Moreover, as we discussed above, on certain supports the surface atoms can compete with the TM for binding hydrogen and other intermediates, giving rise to a complex interplay of adsorption and diffusion processes. This dynamic behaviour can critically influence the overall catalytic performance of SACs.

3.4 Is water just a solvent?

Electrochemical reactions can take place in water or in other solvents, including non-aqueous molecular solvents, ionic liquids, solid electrolytes, and molten salts. Among these, water is the most commonly used solvent, despite its relatively limited electrochemical potential window. In many computational studies of the HER, including high-throughput screening of potential new materials, water is often not explicitly included in the model. This simplification relies on the assumption that water primarily acts as a high-dielectric, polar medium that screens charges and stabilizes ions. To more accurately account for the role of water, several modeling approaches have been developed. Implicit solvent models treat water as a polarizable dielectric continuum characterized by a dielectric constant (ε ≈ 78),^40,41 but they do not capture its molecular structure or hydrogen-bonding network. Explicit solvent methods include a few or many water molecules near the electrode or adsorbates, sometimes combined with ab initio molecular dynamics (AIMD) to account for dynamic effects.⁴² Intermediate approaches, such as micro-solvation models,⁴³ have also been introduced to balance computational cost with a more realistic description of solvent–solute interactions. While screening a large number of novel catalysts without explicitly considering the role of water can help narrow down the list of potential candidates, it is clear that neglecting the solvent is a rather crude approximation. This simplification can have direct and significant effects on the calculated theoretical overpotentials. Here we discuss the role of water not merely as a solvent, but as a ligand capable of coordinating to TM centres and modifying their properties.

SACs embedded in graphene or N-doped graphene often show highly reversible redox switching under reaction conditions (e.g., cycling between TM^δ+ and TM^(δ+1) during the ORR, OER, CO₂ reduction, HER, etc.). Compared to nanoparticles, SACs exhibit more discrete and tunable redox states, largely due to strong interactions with the support. Accurately describing the relative stability of different oxidation states is therefore essential for reliable theoretical predictions. However, this task is often less straightforward than commonly assumed.

To illustrate this challenge, we consider Fe atoms stabilized in phthalocyanine (FePc).⁴⁴ FePc complexes are active in the reduction of CO₂ and NO₃⁻, as well as in the OER, which constitutes the anodic half-reaction of water splitting. In heterogeneous catalysis, Fe@N-Gr SACs are widely employed for the OER and are regarded as the solid-state counterparts of homogeneous FePc catalysts.⁴⁵ As noted above, the reactivity of Fe@N-Gr SACs correlates linearly with that of FePc molecular systems, which justifies focusing the analysis on FePc.²⁷ For this purpose, we performed calculations with the Gaussian molecular code,⁴⁶ using the PBE0 hybrid functional⁴⁷ to properly capture the magnetic character of Fe²⁺ and Fe³⁺ ions. Our aim is to determine the redox potential (V) of the Fe²⁺/Fe³⁺ couple.

In FePc, Fe donates two electrons to the phthalocyanine (Pc) ligand, resulting in Fe(II) (4s⁰3 d⁶). In the square-planar ligand field, the 3d orbitals split into the (3d_xy)²(3d_z²)²(3d_xz)¹(3d_yz)¹(3d_x²−y²)⁰ triplet configuration. Under oxidative conditions, such as at the OER equilibrium potential, the active species is Fe(III). This is consistent with experimental observations of the Fe²⁺/Fe³⁺ redox transition, which occurs around 0.6–0.8 V vs. SHE.^48–50 In the Fe(III) state, the electronic configuration becomes (3d_z²)²(3d_xy)¹(3d_xz)¹(3d_yz)¹(3d_x²−y²)⁰ (quartet spin state). The Fe(II) configuration can also be described as d⁶L, where L represents the ligand shell. Accordingly, Fe(III) in FePc is classified as d⁵L, since experimental data indicate that electron removal from the neutral system involves the Fe centre directly.

DFT modeling of this seemingly simple process reveals a more intricate picture.⁴⁴ Calculations predict that electron removal from FePc occurs primarily on the organic ligand, while the Fe centre remains in the +II oxidation state. This result contradicts experimental evidence and may lead to misleading interpretations of both the reaction mechanism and the associated energetics. The origin of this discrepancy lies in the electronic structure of FePc: the highest occupied molecular orbital (HOMO) is localized on the ligand (Fig. 4), a feature that persists regardless of the specific DFT functional applied. Consequently, removing one electron depletes the HOMO, producing a d⁶L⁻¹ configuration in which Fe remains in the +II oxidation state, at odds with experimental observations.

Fig. 4 Orbital energy diagram and [FePc] (a), [FePc(H₂O)₁] (b) and [FePc(H₂O)₂] (c) at the PBE0 level.

The root of the problem lies in the physical model itself, specifically, the coordination environment of the Fe atom in carbon-based SACs and Fe-based complexes. Water plays a crucial role in defining the redox properties of the Fe²⁺/Fe³⁺ couple, as it acts as a ligand that modifies the relative stability of Fe 3d orbitals with respect to ligand-based electronic levels. To illustrate this effect, we considered FePc complexes with explicit water coordination. In the ground state of FePc(H₂O)₁, with one H₂O molecule added, Fe is still in the +II oxidation state (triplet ground state). However, upon electron removal, the electron is now taken from the Fe centre, yielding Fe³⁺ in a d⁵L configuration. Thus, the coordination of even a single water molecule is sufficient to reconcile theory with the experimentally observed redox behavior. Water coordination destabilizes and raises the energy of the occupied frontier molecular orbital, which is dominated by Fe 3d contributions, making it nearly isoenergetic with the ligand-based orbital (Fig. 4). When two water molecules coordinate in a pseudo-octahedral geometry, [FePc(H₂O)₂], the same effect is observed, but even more pronounced.

Further evidence for the critical role of water as a ligand comes from the calculated oxidation potential:

V_{Fe²⁺/Fe³⁺} = ΔE_DFT − TΔS + ΔE_ZPE = (E_[FePc]⁺ − E_[FePc]) − T(S_[FePc]⁺ − S_[FePc]) + (ZPE_[FePc]⁺ − ZPE_[FePc]) (9)

where E_[FePc] and E_[FePc]⁺ are the DFT energies of the FePc systems, TS_[FePc] and TS_[FePc]⁺ are the entropic corrections and ZPE_[FePc] and ZPE_[FePc]⁺ are the corresponding zero-point energies of the complexes. Since the atomic coordinates of the complex and the water cluster have always been relaxed, the energetic terms related to the change in oxidation state of the metal atoms, and to some extent the reorganization energy of the water clusters, are included in the ΔE_DFT contribution. The DFT(PBE0)-calculated redox potential for FePc without water coordination is ∼1.2 V, significantly overestimating the experimental range of 0.6–0.8 V. Consequently, reaction profiles involving the Fe(II)/Fe(III) transition are strongly affected by this discrepancy. For FePc(H₂O)₁, the oxidation correctly involves the Fe atom and not the Pc ligand, but the calculated potential remains close to 1.2 V. However, upon coordination of a second water molecule, yielding FePc(H₂O)₂ with a coordinatively saturated Fe centre, the Fe²⁺/Fe³⁺ redox potential decreases to ∼0.7 V, in excellent agreement with experiment. The addition of further water molecules produces only minor fluctuations in this value (±0.2 V, Table 4).

Table 4 Calculated oxidation potential, V vs. RHE of the Fe²⁺/Fe³⁺ couple in FePc complexes (Gaussian code, PBE0+D3 functional).⁴⁴

System	V(Fe^2+/Fe^3+)/V
[FePc]/[FePc]^+	1.22
[FePc(H2O)]/[FePc(H2O)]^+	1.24
[FePc(H2O)2]/[FePc(H2O)2]^+	0.72
[FePc(H2O)4]/[FePc(H2O)4]^+	0.72
[FePc(H2O)6]/[FePc(H2O)6]^+	0.54
[FePc(H2O)8]/[FePc(H2O)8]^+	0.46
[FePc(H2O)10]/[FePc(H2O)10]^+	0.51
[FePc(H2O)12]/[FePc(H2O)12]^+	0.58
Exp	0.6–0.8 (ref. 48–50)

---

The key conclusion is that the catalytically active species is not the four-coordinate Fe centre in the bare FePc complex, but rather the five- or six-coordinate Fe present in FePc(H₂O)_1–2. Therefore, accurate quantum chemical modeling of FePc chemistry must explicitly include water coordination. This insight extends beyond FePc to Fe@N-Gr and, more broadly, to all TM atoms embedded in two-dimensional supports.

3.5 Evolution of the catalyst under operating conditions: (1) stability

“If you can identify it, it is probably not the catalyst”.⁵¹ With this remark, Jack Halpern, a chemist who devoted his career to unravelling the mechanisms of homogeneous catalysis, highlighted a central challenge: structural characterization efforts are often frustrated by the fact that catalysts may undergo dynamic changes during the reaction. Indeed, it is not uncommon for the active phase to emerge only under catalytic conditions. This realization has driven the development of operando spectroscopies, which provide crucial insights into the true nature of active sites and the species that form under working conditions.⁵² The dynamic evolution of catalytic sites is therefore of paramount importance in the modeling of catalytic reactions. In principle, sufficiently long ab initio molecular dynamics simulations could capture these transformations. However, the timescales involved remain prohibitive with current computational technologies. While significant advances are anticipated through the application of machine-learning-derived interatomic potentials, current approaches to catalyst evolution remain largely confined to the construction of specific, simplified models.

The discussion of catalyst changes under operating conditions can be divided into two key aspects. The first concerns the intrinsic stability of the catalyst itself. The second relates to the formation of species that bind to the active site, either poisoning or promoting it, but in all cases altering its activity. Both aspects are critical for establishing a meaningful comparison between theoretical predictions and experimental observations.

This paper focuses on SACs in the context of the HER. Accordingly, we restrict the discussion to catalyst stability under electrochemical conditions. In addition to atom diffusion and aggregation, two processes well known in thermal catalysis, electrocatalysis introduces further degradation pathways. The metal atom may undergo changes in oxidation state, interact with the solvent, or form stable complexes that ultimately lead to metal leaching and dispersion into the electrolyte (dissolution). Recently, we proposed a simple approach to estimate the relative stability of a SAC under specific electrochemical conditions of pH and applied potential.⁵³ The method, based on the Pourbaix diagram formalism, involves constructing a thermodynamic cycle from which the redox potential of a SAC, typically unknown, can be derived. This is achieved by combining tabulated experimental data with one theoretically accessible quantity: the binding energy of the TM atom to its support (Fig. 5).

Fig. 5 Thermodynamic cycle to determine the Gibbs free energy of the dissolution process (redox process) of a generic SAC, TM@S (TM = transition metal, S = support) in aqueous solution. For the sake of simplicity, the diagram shows only the formation of solvated ions and hydroxides, but in reality also complex oxoanions can form.

Depending on the electrochemical environment, the metal atom may dissolve by forming aquo-complexes, or, under certain pH and potential conditions, generate new oxide or hydroxide species that either precipitate or remain in solution (Fig. 5). Throughout these processes, the metal centre can undergo changes in oxidation state, yielding a wide variety of products. The associated free-energy changes, ΔG_1–4 (Fig. 5), represent the target quantities of interest and can be obtained through the thermodynamic cycle by decomposing the overall process into intermediate steps, each with Gibbs free energies either experimentally available or derivable from DFT.

The first step of the thermodynamic cycle is the detachment of the metal atom from the support, denoted as ΔG_1–2. This step represents the binding energy of a gas-phase TM atom to the support site S where it is stabilized. While this quantity is difficult to measure experimentally, it can be computed with high accuracy using quantum-chemical methods. Importantly, it is the only variable in the thermodynamic cycle that requires theoretical evaluation. The second step, ΔG_2–3, corresponds to the inverse of the bulk metal sublimation process and is equivalent to the cohesive energy of the metal. This quantity is well established and available from experimental data. Finally, the closing step of the cycle, ΔG_3–4, describes a classical redox process in which the bulk metal is oxidized to form aquo-complexes, oxides, hydroxides, or related species. The corresponding free-energy change is also known experimentally, as it coincides with the redox potential of the metal.

This information can be integrated with Pourbaix diagrams,⁵⁴ which indicate whether the active site corresponds to the as-prepared catalyst or is instead covered by species formed under reaction conditions. To do so, one must compute the free energy of the adsorbates as a function of pH and applied potential, using the CHE approach:

ΔG_x = ΔG⁰_x − nV_SHE − nk_bT ln(10)pH (10)

(n is the number of exchanged electrons). The stability diagrams are constructed by taking the most stable species for each value of pH and voltage, V. To illustrate this, we report the case of Ni atoms embedded in carbon nitride, C₃N₄, Ni@C₃N₄.

One of the most commonly reported models for TM atoms in C₃N₄ involves the coordination of the active metal species within heptazine cavities, bound to four non-equivalent nitrogen atoms (Fig. 6a).⁵⁵ This model has been extensively employed to simulate and predict the reactivity of TM@C₃N₄ systems.⁵⁶ However, stability analysis of Ni@C₃N₄ using the Pourbaix diagram approach described above indicates that the system is highly unstable (Fig. 6b). Under HER conditions (V ∼ 0 V at pH = 0), the catalyst in its clean state (no adsorbed species), indicated with *, is unstable and dissolves into Ni²⁺ ions in solution. Under a strongly reducing applied potential it becomes covered by H atoms, H*. Even at the equilibrium potential for the O₂/H₂O⁻ redox (V ∼ 1.23 V at pH = 0), Ni@C₃N₄ is unstable and dissolves or, at alkaline pH, it forms hydroxide species; see Fig. 6b. Therefore, Ni@C₃N₄ can exist only in a narrow window of low voltage and neutral pH, while it tends to form Ni²⁺ ions in solution or Ni(OH)₂ at moderate applied voltages. In short, Ni supported on the heptazine pore of a C₃N₄ nanosheet is expected to be very labile. This example highlights that performing calculations on systems that would not be stable under operando conditions offers limited value for understanding the true catalytic behavior.

Fig. 6 (a) Structure of a Ni atom in the heptazinic cavity of C₃N₄. (b) Pourbaix diagram of Ni@C₃N₄. The pH–voltage window where the catalyst is stable in the SAC state, either clean of covered by some intermediate, is reported with colors identifying H*, * (free catalyst), Ni²⁺ (aquo cation), and Ni(OH)₂ (hydroxide). Black dashed lines indicate the redox levels of H⁺/H₂ and O₂/H₂O couples.

In a recent study, we demonstrated that the actual structure of Ni@C₃N₄ is indeed different.⁵⁷ During the synthesis, Ni²⁺ ions become trapped by melem units, the precursor molecules for C₃N₄, during polymerization. As the melem units polymerize, they pair up to form filaments that incorporate Ni in a square-planar coordination environment. A combination of spectroscopic data and DFT calculations confirms the robustness of this structural assignment, showing that the as-prepared catalyst contains charged Ni²⁺ ions.⁵⁷

In general, the binding energy of an isolated metal atom to the support, E_b, serves as a simple yet reliable descriptor of SAC stability. Using the thermodynamic cycle illustrated in Fig. 5, one can calculate the minimum E_b required to ensure that a SAC remains stable under specific conditions.⁵³ More precisely, for given external conditions of pH and applied potential, it is possible to determine a threshold value of E_b that prevents dissolution of the metal centre. This approach provides a straightforward first estimate of SAC stability, which can be further refined by employing higher-accuracy DFT methods or more sophisticated models to account for the stability of reactive species. Importantly, we emphasize that evaluating the stability of a SAC is just as critical as analysing its reactivity.

3.6 Evolution of the catalyst under operating conditions: (2) changes of active site

Unlike the old “rigid site” view (where one imagined the catalyst as having a fixed set of sites), we now know that the active site can evolve during the reaction, adapting to the environment and reactants. For instance, reactants, intermediates, or even spectator molecules can adsorb onto the catalyst surface and modify the electronic or geometric structure of the active site. In electrocatalysis, adsorbed oxygen, hydrogen, or hydroxyl groups can “decorate” the active site, resulting either in poisoning or, more rarely, promoting the catalytic reaction.

Recently, we explored the formation of unconventional intermediates in the OER.^58,59 While on metal electrodes the reaction is typically assumed to proceed through OH*, O*, and OOH* intermediates, on SACs the chemistry is more complex, allowing for the formation of additional species (see Fig. 7).

Fig. 7 Possible OER intermediates on SACs embedded in N-doped graphene (4N-Gr). Reprinted from ref. 58, Copyright 2023, with permission from Elsevier.

Some of these adsorbed species form exceptionally stable bonds with the TM, resulting in deep minima on the potential energy surface. For example, on Mn@4N-Gr, the O* species binds with a ΔG of about −1.5 eV, corresponding to a very stable configuration. In this case, the adsorbed oxygen atom is so strongly bound that it suppresses further reactivity at the site.

Under oxidative potentials, water can split, with formation of an OH group adsorbed on the surface:

H₂O + * → H⁺ + OH* + e⁻ (11)

On Sc@4N-Gr, the OH group forms with a ΔG of −3 eV, giving rise to a species in which the Sc–OH bond is very difficult to break. A similar result was reported by Ye et al., who investigated the nature of the active site in SACs for the HER across the full pH range: they found that Sc@N-4Gr is always poisoned by OH groups regardless of pH.²⁷

However, reactive species formed under realistic conditions are not always detrimental to catalytic activity. We recently found that SACs based on TMs that bind reactants too strongly, thus forming overly stable intermediates that block the reaction, can regain activity by bonding to other molecular fragments.⁶⁰ In this case, the TM becomes coordinated with additional ligands, which decreases its binding strength toward reaction intermediates according to the bond-order conservation principle.⁶¹ Through this mechanism, a SAC can be shifted from the strong-binding left side to the apex of the volcano plot, where intermediates are bound neither too strongly nor too weakly. This results in labile complexes, a fundamental requirement for efficient catalysis.

This discussion highlights that predicting catalytic activity based solely on the hypothetical structure of the as-prepared catalyst can be highly misleading. For example, it has been suggested that the discrepancy between the high HER activity of V@4N-Gr predicted by DFT and the experimentally observed inert behaviour (Table 1) arises from the strong poisoning effect of OH* species.²⁷ Neglecting such dynamic effects can therefore lead to entirely erroneous conclusions.

3.7 Last but not least: role of the functional

At the conclusion of this analysis of factors affecting the predicted reactivity of SACs, it is important to address a fundamental limitation of DFT calculations: the dependence of chemical bond strengths on the choice of exchange–correlation (XC) functional. Among the most widely used XC functionals in the catalysis community are the local density approximation (LDA),⁶² representing the first rung of Perdew and Schmidt’s “Jacob’s ladder” of DFT,⁶³ and the generalized gradient approximation (GGA) in the Perdew–Burke–Ernzerhof (PBE) form,²⁰ which corresponds to the second rung. Hybrid functionals, such as PBE0,⁴⁷ are classified on the fourth rung.

An alternative strategy to mitigate the self-interaction errors inherent to LDA and GGA is the DFT + U approach, in which DFT is augmented with a Hubbard-like U parameter.⁶⁴ This correction, originally introduced in the Hubbard Hamiltonian to account for inter-electron repulsion, improves the description of localized electrons in systems where standard XC functionals are insufficient.

The accuracy of different computational approaches for studying TM chemistry, ranging from DFT to highly sophisticated wavefunction-based methods, such as coupled-cluster theory, has long been a central topic of discussion within the quantum chemistry community.^65,66 For example, Patel et al. investigated Cu-modified covalent triazine framework catalysts for the oxygen reduction reaction and benchmarked a variety of DFT methods against coupled-cluster calculations.⁶⁷ Their results showed unambiguously that hybrid functionals provide the most accurate description of ORR adsorbate adsorption energies on these complex systems.

Despite the rich literature on the subject, the impact of the choice of XC functional is often underestimated when predicting SAC reactivity in the HER using the ΔG_H = 0 criterion for optimal catalysts. We recently compared PBE, PBE+U, and PBE0 functionals for a series of TM@4N-Gr SACs.⁶⁸ For certain systems, the ΔG_H values obtained with these functionals differ substantially. For example, for Co@4N-Gr, ΔG_H = 0.13 eV at the PBE level compares with 0.62 eV (PBE+U) and 0.48 eV (PBE0) (Table 5). This is not merely a quantitative difference, it changes the qualitative assessment of the catalyst, from active (PBE) to inactive (PBE+U and PBE0). However, this behavior is system-dependent. For Ni@4N-Gr, all three functionals yield very similar ΔG_H values: 1.65 eV (PBE), 1.63 eV (PBE+U), and 1.65 eV (PBE0).⁶⁸ The difference arises from the electronic configuration: Co@4N-Gr has unpaired electrons in its ground state, whereas Ni@4N-Gr adopts a d⁸ closed-shell configuration. Self-interaction-corrected functionals, such as PBE+U or PBE0, better describe localized unpaired electrons, like those in the 3d orbitals of TM atoms, leading to significant variations in the calculated H adsorption energy.

Table 5 Gibbs free energy, ΔG_H, of adsorption of an H atom on TM@4N-Gr SACs [eV]

SAC	ΔGH (eV)
	r2SCAN	PBE^68	PBE+U^68	PBE0 ^68
Ti@4N-Gr	−0.63	−0.86	−0.58	−0.74
Cr@4N-Gr	0.54	0.31	0.65	0.67
Mn@4N-Gr	0.71	0.53	0.99	1.02
Fe@4N-Gr	0.67	0.37	1.04	1.12
Co@4N-Gr	0.45	0.13	0.62	0.48
Ni@4N-Gr	1.62	1.65	1.63	1.65
Cu@4N-Gr	1.64	1.72	1.98	2.14
Mo@4N-Gr	−0.30	−0.42	−0.14	−0.03
Ru@4N-Gr	−0.34	−0.46	−0.16	−0.26
Rh@4N-Gr	−0.30	−0.28	−0.29	−0.59
Pd@4N-Gr	2.02	1.97	2.00	2.17
Ag@4N-Gr	0.90	0.76	0.79	0.79
W@4N-Gr	−0.87	−1.05	−0.76	−0.84
Os@4N-Gr	−0.67	−0.70	−0.51	−0.71
Ir@4N-Gr	−0.38	−0.42	−0.40	−0.74
Pt@4N-Gr	1.76	1.63	1.58	1.98
Au@4N-Gr	2.19	2.10	2.19	2.31

---

In 2015, Sun, Ruzsinszky, and Perdew proposed the Strongly Constrained and Appropriately Normed (SCAN) meta-GGA,⁶⁹ on the third rung of Jacob’s ladder. To enhance its numerical efficiency, Furness et al. proposed the restored-regularized SCAN (r2SCAN).⁷⁰ Since their introduction, SCAN and r2SCAN have attracted considerable attention due to their improved performance and small additional computational cost compared to PBE.

We have recomputed a set of 17 SACs containing TM atoms from the first, second, and third rows of the periodic table using the r2SCAN functional and compared the new results with those obtained previously (Table 5). For the critical case of Co@4N-Gr discussed earlier, r2SCAN predicts ΔG_H = 0.45 eV, very close to the PBE0 value of 0.48 eV. For Ni@4N-Gr, ΔG_H is 1.62 eV at the r2SCAN level, in excellent agreement with the other functionals. However, this agreement is not universal. For Fe@4N-Gr, r2SCAN yields ΔG_H = 0.67 eV, which lies between the lower value from PBE (0.37 eV) and the higher value from PBE0 (1.12 eV). For second- and third-row TMs, r2SCAN generally provides results similar to those from other functionals, with a few exceptions, for example, Mo@4N-Gr and Ir@4N-Gr, where PBE0 predicts significantly different values. Overall, however, the uncertainties associated with the choice of XC functional are less pronounced for these heavier elements, making the functional choice less critical (Table 5).

Large variations are observed for first-row TMs, particularly Cr, Mn, Fe, and Co. As noted above, the presence of unpaired 3d electrons leads to functional-dependent differences in ΔG_H. In general, r2SCAN provides values intermediate between those obtained with PBE and PBE0. The only reliable way to determine which functional is most accurate is to perform high-level post-Hartree–Fock calculations, such as the coupled-cluster singles and doubles (CCSD) method, on molecular analogs of the SACs. Patel et al.⁶⁷ investigated the reactivity of a Cu complex in the ORR at the CCSD(T) level of theory. The Cu complex features a bonding environment similar to that of SACs on carbon supports. They found that binding energies of various intermediates predicted by the PBE0 hybrid functional were in close agreement with the CCSD(T) results, with a mean absolute error of 0.1 eV. In contrast, GGA functionals exhibited significantly larger errors (0.6–0.7 eV), as hybrid functionals provide a more accurate description of the electronic structure and bonding in these systems. Patel et al. concluded that “the blind use of GGA functionals to describe single-atom catalysts may produce inaccurate results,” a caution often overlooked by theoretical chemists in this field.⁶⁷ In the absence of high-level benchmarks, it is advisable to test the stability of results by comparing multiple functionals, particularly for first-row TM SACs.

4 Conclusions

The rational design of new catalysts from first principles is often cited as one of the ultimate goals of modern quantum chemistry. The enormous increase in computing power and the development of highly sophisticated electronic-structure codes have made this objective more attainable than ever before. Nevertheless, despite the undeniable progress achieved, the question of to what extent it is truly possible to predict the activity of a given catalyst from electronic-structure calculations remains open.

The field of single-atom catalysis (SAC), where the presence of an isolated metal atom reduces complexity compared with more traditional heterogeneous catalysts based on supported metal particles, offers new opportunities in this direction. The catalytic properties of nanoparticles depend strongly on factors such as particle size, shape, adsorption site, and the presence of defects, not to mention the dynamic effects that make metal clusters fluxional at finite temperatures, all of which are difficult to model accurately. In principle, SACs, comprising isolated metal atoms embedded in a solid matrix, offer more favourable conditions for predictive modelling. This is one of the reasons for the large number of studies devoted to predicting the properties of SACs supported on carbon-based materials (graphene, carbon nitride, carbon nanotubes, etc.). Today, with hundreds of such predictive studies reported in the literature and a smaller but significant number of systems tested experimentally, it is possible to attempt a critical assessment. However, as discussed in the Introduction, the results so far have been far from encouraging (Table 1).

In this work, we analysed the origin of this partial failure. The hydrogen evolution reaction (HER) in electrocatalytic processes serves as a prototypical case. The prediction of new catalysts is typically based on the computational hydrogen electrode (CHE), a powerful yet simplified approach to this problem. We discussed several reasons why comparisons between predicted and measured catalytic activities based on the CHE often lead to disappointing results.

First, even minor changes in the local environment of the active site can cause substantial differences in the thermodynamics of the process. This is a critical issue because, while the structure of a SAC is perfectly defined in theoretical models, its actual configuration is often unknown experimentally, or inferred only approximately, without the same level of structural detail. A second, though somewhat less critical, aspect is that theoretical models tend to focus exclusively on the metal site, even though atoms in the surrounding matrix may also participate in the reaction. In other words, all potentially reactive sites should be considered in the analysis.

The same reasoning applies to reaction intermediates. It is commonly assumed that the same intermediates formed on extended metal electrodes (e.g., H* in the HER; OH*, O*, and OOH* in the OER) also occur on SACs. However, SACs can bind a broader variety of chemical species and form a wider range of intermediates than metal surfaces. Neglecting these species can therefore limit the comparability between theory and experiment.

Another crucial aspect of electrochemical reactions is the role of the solvent. In aqueous environments, water can interact with the active site like any other ligand and even react with it. This exemplifies the dynamic evolution of a catalyst under operating conditions. For instance, under specific potentials, water can form OH groups that bind strongly to TM atoms, thereby blocking or altering the chemical properties of the active site. Predicting the activity of a SAC at a given state does not guarantee that it will remain unaltered throughout the reaction.

In many cases, the discrepancy between predicted and experimental activity arises simply because the newly designed catalysts are unstable under realistic pH and potential conditions, leading to dissolution, aggregation, precipitation, or other degradation phenomena that destroy the catalytic material.

Finally, one must bear in mind that DFT provides only an approximate solution to the Schrödinger equation. Different exchange–correlation functionals can yield noticeably different reaction free energies. This intrinsic uncertainty sometimes prevents reliable quantitative predictions of catalytic activity.

Considering all these factors (and the list is by no means exhaustive), it is not surprising that elegant but simplified approaches that omit key aspects of catalyst chemistry often fail to deliver reliable predictions for previously untested systems. At best, they can suggest general trends or narrow down the range of candidate materials, though even this requires caution. As we enter the promising new era of artificial intelligence applied to the design of materials with unprecedented properties, it becomes essential to identify clearly all the factors contributing to the overall chemical activity of SACs and their relative importance. Doing so will facilitate the identification of robust descriptors and simple correlations, helping to advance the field toward the long-sought goal of the rational design of new catalysts from first principles.

Author contributions

HB, AB and MS, contributed to this study with methodology, investigation, formal analysis of the computational results. GDL and GP contributed to this study with the conceptualisation, funding acquisition, project administration, supervision and writing (review and editing).

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting this article have been included as part of the supplementary information (SI). Supplementary information: adsorption properties and structures of H adsorption on 4N-Gr, Ti@4N-Gr, Rh@4N-Gr and Ce@4N-Gr. Input coordinates of the structures considered in the paper are available from the authors on request. See DOI: https://doi.org/10.1039/d5fd00112a.

Acknowledgements

The authors thanks Prof. Livia Giordano and Gianvito Vilè for useful discussions. GDL has received funding through the PRIN project “SACtoH2” (project code P2022AZETB) from the Italian Ministry for Universities and Research (MUR), in the context of the National Recovery and Resilience Plan and co-financed by the European Commission – Next Generation EU (Mission 4, Component 1). Access to the CINECA supercomputing resources was granted via ISCRAB and EuroHPC regular access projects.

References

1. A. Wang, J. Li and T. Zhang, Nat. Rev. Chem., 2024, 2, 65–81.
2. S. K. Kaiser, Z. Chen, D. Faust Akl, S. Mitchell and J. Pérez-Ramírez, Chem. Rev., 2020, 120, 11703–11809.
3. L. Gloag, S. V. Somerville, J. J. Gooding and R. D. Tilley, Nat. Rev. Mater., 2024, 9, 173–189.
4. C. Coperet, A. Comas-Vives, M. P. Conley, D. P. Estes, A. Fedorov and V. Mougel, et al., Chem. Rev., 2016, 116, 323–421.
5. E. Skúlason, V. Tripkovic, M. E. Björketun, S. Gudmundsdóttir, G. Karlberg and J. Rossmeisl, et al., J. Phys. Chem. C, 2010, 114, 18182–18197.
6. Web of Science (search: “DFT AND HER OR hydrogen evolution reaction”; “DFT AND HER OR hydrogen evolution reaction AND SAC”).
7. V. Fung, G. Hu, Z. Wu and D. E. Jiang, J. Phys. Chem. C, 2020, 124, 19571–19578.
8. H. Jin, M. Ha, M. G. Kim, J. H. Lee and K. S. Kim, Adv. Energy Mater., 2023, 13, 2204213.
9. M. D. Hossain, Z. Liu, M. Zhuang, X. Yan, G. L. Xu, C. A. Gadre, A. Tyagi, I. H. Abidi, C. J. Sun and H. Wong, et al., Adv. Energy Mater., 2019, 9, 1803689.
10. H. Fei, J. Dong and M. Arellano-Jiménez, et al., Nat. Commun., 2015, 6, 8668.
11. J. Guan, X. Wen, Q. Zhang and Z. Duan, Carbon, 2020, 164, 121–128.
12. Y. Liu, J. Wu, Y. Zhang, X. Jin, J. Li and X. Xi, et al., ACS Appl. Mater. Interfaces, 2023, 15, 14240–14249.
13. Z. Zhang, J. Cai, H. Zhu, Z. Zhuang, F. Xu and J. Hao, et al., Chem.–Eng. J., 2020, 392, 123655.
14. Q. Wu, H. Li, Y. Zhou, S. Lv, T. Chen and S. Liu, et al., Inorg. Chem., 2022, 61, 11011–11021.
15. N. Cheng, S. Stambula and D. Wang, Nat. Commun., 2016, 7, 13638.
16. H. Wu, Q. Zhang, S. Chu, H. Du, Y. Wang and P. Liu, Materials, 2024, 17, 5082.
17. S. Yadav, V. Dao, W. Wang, K. Chen, C. Kim, G. C. Kim and I. H. Lee, Mater. Adv., 2023, 4, 6498–6506.
18. G. Kresse and J. Hafner, Phys. Rev. B Condens. Matter, 1993, 47, 558–561.
19. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
20. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
21. P. E. Blöchl, Phys. Rev. B Condens. Matter, 1994, 50, 17953–17979.
22. H. J. Monkhorst and J. D. Pack, Phys. Rev. B Condens. Matter, 1976, 13, 5188.
23. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
24. M. Kottwitz, Y. Li, H. Wang, A. I. Frenkel and R. G. Nuzzo, Chem.:Methods, 2021, 1, 278–294.
25. H. V. Thang, G. Pacchioni, L. De Rita and P. Christopher, J. Catal., 2018, 367, 104–114.
26. T. Asset and P. Atanassov, Joule, 2020, 4, 33–44.
27. S. Ye, F. Liu, F. She, J. Chen, D. Zhang and A. Kumatani, et al., Angew. Chem., 2025, 137, e202425402.
28. G. Di Liberto, L. A. Cipriano and G. Pacchioni, ChemCatChem, 2022, e202200611.
29. L. Wang, H. Wang and J. Lu, Chem Catal., 2023, 3, 100492.
30. M. Spotti, G. Di Liberto, G. Pacchioni and G. Topics, Catal, 2025, 68, 1837–1847.
31. J. Yang, W. Li, D. Wang and Y. Li, Adv. Mater., 2020, 32, 2003300.
32. R. T. Hannagan, G. Giannakakis, M. Flytzani-Stephanopoulos and E. C. H. Sykes, Chem. Rev., 2020, 120, 12044–12088.
33. E. C. Manzi, M. Stamatakis, G. Di Liberto, G. Pacchioni and G. Surf, Science, 2025, 754, 122688.
34. C. Saetta, I. Barlocco, G. Di Liberto and G. Pacchioni, Small, 2024, 20, 2401058.
35. R. Prins, Chem. Rev., 2012, 112, 2714–2738.
36. G. J. Kubas, Acc. Chem. Res., 1988, 21, 120–128.
37. R. H. Crabtree, Chem. Rev., 2016, 116, 8750–8769.
38. G. Di Liberto, L. A. Cipriano and G. Pacchioni, J. Am. Chem. Soc., 2021, 143, 20431–20441.
39. H. Li, Z. Zhang and Z. Liu, Catalysts, 2019, 9, 84.
40. J. Tomasi and M. Persico, Chem. Rev., 1994, 94, 2027–2094.
41. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
42. A. Groß, F. Gossenberger, X. Lin, M. Naderian, S. Sakong and T. Roman, J. Electrochem. Soc., 2014, 161, E3015.
43. F. F. Calle-Vallejo, R. Morais, F. Illas, D. Loffreda and P. Sautet, J. Phys. Chem. C, 2014, 123, 5578–5582.
44. A. Bonardi, S. Xu, G. Di Liberto and G. Pacchioni, J. Phys. Chem. Lett., 2025, 16, 10049–10057.
45. Y. Zhou, G. Gao, Y. Li, W. Chu and I. W. Wang, PhysChemChemPhys, 2019, 21, 3024–3032.
46. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji et al., Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford CT, 2016.
47. C. Adamo and V. Barone, J. Chem.Phys., 1999, 110, 6158–6170.
48. C. S. Martin, P. Alessio, F. N. Crespilho, C. M. A. Brett and C. J. L. Constantino, J. Electroanal. Chem., 2019, 836, 7–15.
49. F. Van den Brink, W. Visscher and E. Barendrecht, J. Electroanal. Chem. Interfacial Electrochem., 1984, 175, 279–289.
50. C. Zúñiga Loyola, N. Troncoso, A. Gatica Caro and F. Tasca, ChemElectroChem, 2024, 11, e202400186.
51. J. E. Bercaw, Proc. Natl. Acad. Sci. U.S.A., 2018, 115, 5049–5050.
52. B. Roldán Cuenya and M. A. Bañares, Chem. Rev., 2024, 124, 8011–8013.
53. G. Di Liberto, L. Giordano and G. Pacchioni, ACS Catal., 2023, 14, 45–55.
54. E. McCafferty, Introduction to Corrosion Science, Springer, New York, 2009, pp. 95–117.
55. P. Suja, J. John, T. P. D. Rajan, G. M. Anilkumar, T. Yamaguchi and S. C. Pillai, et al., J. Mater. Chem. A, 2023, 11, 8599–8646.
56. G. Vilé, G. Di Liberto, S. Tosoni, A. Sivo, V. Ruta and M. Nachtegaal, et al., ACS Catal., 2022, 12, 2947–2958.
57. N. Allasia, S. Xu, S. F. Jafri, E. Borfecchia, L. A. Cipriano and G. Terraneo, et al., Small, 2025, 2408286.
58. I. Barlocco, L. A. Cipriano, G. Di Liberto and G. Pacchioni, J. Catal., 2023, 417, 351–359.
59. L. A. Cipriano, G. Di Liberto and G. Pacchioni, ACS Catal., 2022, 12, 11682–11691.
60. E. Inico, M. Spotti, G. Pacchioni and G. Di. Liberto, ACS Catal. , in press.
61. E. Shustorovich, in Advances in Catalysis, Academic Press, 1990, vol. 37, pp. 101–163.
62. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
63. J. P. Perdew and K. Schmidt, AIP Conf. Proc., 2001, 577, 1–20.
64. V. I. Anisimov, F. Aryasetiawan and A. I. Lichtenstein, J. Phys.: Condens. Matter, 1997, 9, 767.
65. X. Xu, W. Zhang, M. Tang and D. G. Truhlar, J. Chem. Theory Comput., 2015, 11, 2036–2052.
66. A. J. Cohen, P. Mori-Sánchez and W. Yang, Chem. Rev., 2012, 112, 289–320.
67. A. M. Patel, S. Ringe, S. Siahrostami, M. Bajdich, J. K. Nørskov, R. A and A. R. Kulkarni, J. Phys. Chem. C, 2018, 122, 29307–29318.
68. I. Barlocco, L. A. Cipriano, G. Di Liberto and G. Pacchioni, Adv. Theory Simul., 2023, 6, 2200513.
69. J. Sun, A. Ruzsinszky and J. P. Perdew, Phys. Rev. Lett., 2015, 115, 036402.
70. J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun, J. Phys. Chem. Lett., 2020, 11, 8208.

---