Multifunctional electrocatalysis on transition metal-doped biphenylene: a first-principles and machine learning study of single-atom catalysts for HER, OER, and ORR

Abstract

Single-atom catalysts (SACs) anchored on two-dimensional substrates have emerged as a frontier in electrocatalysis due to their atomic-level efficiency and tunable reactivity. Herein, we present a comprehensive theoretical investigation of 3d, 4d, and 5d transition metal (TM) atoms embedded in single-vacancy biphenylene (SV-BPN), a recently synthesized carbon allotrope with a unique nonbenzenoid topology. By integrating density functional theory (DFT) and machine learning (ML), we evaluated the structural stability, electronic characteristics, and catalytic activity of TM@SV-BPN systems toward the hydrogen evolution reaction (HER), oxygen evolution reaction (OER), and oxygen reduction reaction (ORR). Our screening reveals several high-performing SACs with low overpotentials, including Mo@SV-BPN (η^HER = 0.006 V), Pd@SV-BPN (η^OER = 0.43 V), and Ag@SV-BPN (η^ORR = 0.67 V). Notably, Au@SV-BPN exhibits trifunctional catalytic activity across all three reactions. Electronic descriptors such as the d-band center and integrated crystal orbital Hamilton population (ICOHP) are correlated with intermediate adsorption energetics. A gradient boosting regression model accurately predicts adsorption energies (R² = 0.98), highlighting charge transfer as the most influential feature. This work not only identifies a novel class of trifunctional SACs but also establishes a data-driven paradigm for rational catalyst design based on biphenylene supports.

---

1. Introduction

Addressing environmental pollution and energy sustainability has emerged as an urgent global research focus in response to the pressing needs of our rapidly developing society. Currently, nonrenewable fossil fuels remain the predominant energy sources in both developed and developing regions, a reality that continues to pose significant environmental threats. The combustion of fossil fuels releases substantial greenhouse gases, exacerbating pollution, environmental degradation, and global climate change. Consequently, pursuing sustainable alternative energy sources has become a primary imperative of the 21st century.¹ Scientific advancements, such as water electrolysis, fuel cells, unitized regenerative fuel cells (URFCs), and metal–air battery systems, critically rely on electrochemical processes,^2,3 particularly the oxygen evolution reaction (OER),^4–7 oxygen reduction reaction (ORR)⁸ and hydrogen evolution reaction (HER).^9–13 However, HER, OER, and ORR are hindered by challenging proton-coupled electron transfer steps, resulting in high energy barriers and sluggish kinetics, thereby limiting efficient energy conversion.¹⁴ Currently, Pt-based alloys represent leading catalysts for HER and ORR,¹⁵ while IrO₂ and RuO₂ are widely recognized as the most effective for OER.¹⁶ Despite their remarkable performance, widespread utilization of these catalysts is significantly constrained by their scarcity and high economic cost.¹⁷ Additionally, employing different catalysts for HER and OER in water-splitting devices or for OER and ORR in metal–air batteries complicates cell design, inevitably increasing equipment and processing costs.¹⁸ Particularly, URFC systems, integrating both fuel cell and electrolyzer functionalities, require multifunctional electrocatalysts.¹⁹ Therefore, developing cost-effective, highly active, and durable electrocatalysts capable of performing dual roles (HER/OER or OER/ORR) or even triple roles (HER/OER/ORR) is of immense significance.

Single-atom catalysts (SACs) represent cutting-edge catalysis technologies,²⁰ demonstrating exceptional catalytic performance and outstanding atomic utilization efficiency, highlighting their critical role in catalytic advancements.^21,22 Among various catalyst supports, two-dimensional (2D) materials, notable for their unique geometric and electronic properties, have emerged as ideal platforms for stabilizing single atoms, thereby offering numerous advantages.²³ Recently, 2D materials have been extensively employed as substrates for SACs,²⁴ effectively catalyzing reactions such as HER,²⁵ OER,²⁶ ORR,²⁷ carbon-dioxide reduction reaction (CO₂RR),²⁸ nitrogen reduction reaction (NRR),^29,30 nitrate reduction reaction (NO₃RR)³¹ and CO-oxidation.³² Recent studies have reported high loadings of transition metal (TM) atoms stabilized on carbon-based materials.³³ For instance, Xia et al. achieved high single-atom loadings (approximately 40 wt%) of Ir-, Pt-, and Ni-based SACs utilizing graphene quantum dot (GQD)-assisted methodologies.³⁴ Although numerous reports on multifunctional catalysts exist, the development of bifunctional SACs remains limited, and trifunctional SACs are exceedingly rare.^35,36 Notably, purely carbon-based SACs without heteroatom doping are scarcely reported. Recently, Fan et al. synthesized a novel planar sp²-hybridized carbon allotrope, the biphenylene (BPN) network, via a dehydro-fluorination reaction (HF-zipping) on Au(111) surfaces.³⁷ This discovery has rapidly attracted significant research attention due to its promising applications.^38–42 Despite preliminary explorations, comprehensive studies on BPN-based SACs, particularly their bifunctional or trifunctional catalytic capabilities (HER, OER, and ORR), remain largely unexplored.

In this work, we systematically investigated a series of TM-doped single-vacancy BPN catalysts (TM@SV-BPN), to elucidate their electronic properties, including binding energy (E_bind), formation energy (E_form), d-band center (ε_d), and Bader charge analysis (e⁻). The electrocatalytic activities for HER, OER, and ORR were evaluated through Gibbs free energies (ΔG_*H, ΔG_*OH, ΔG_*O, and ΔG_*OOH), potential determining step (PDS), limiting potential (U_L) and overpotential (η). Our results indicate that nine catalysts (V, Mn, Fe, Mo, Tc, Ru, Ir, Pt and Au@SV-BPN) exhibit excellent HER activity, with η^HER less than 0.2 V; notably, Mo@ and Tc@SV-BPN demonstrate exceptionally low overpotentials (approximately 0.01 V). For OER, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ag, Ir, Pt, and Au@SV-BPN catalysts exhibit η^OER lower than 1.23 V with Pd@SV-BPN (0.43) and Ni@SV-BPN (0.55) displaying performances comparable to benchmark RuO₂ (0.42 V) and IrO₂ (0.56 V). Catalysts Fe@, Ru@, Ir@, Pt@, and Au@SV-BPN show promising bifunctional capabilities for HER/OER, while Cu@, Rh@, Pd@, Ag@, and Au@SV-BPN exhibit bifunctional activity for OER/ORR. Importantly, Au@SV-BPN is identified as a superior trifunctional catalyst, offering low overpotentials suitable for both water splitting (HER/OER) and metal–air battery applications (OER/ORR). Machine learning analysis further identifies charge transfer (Q_e) as the most critical descriptor for predicting the adsorption energies of OER and ORR intermediates on TM@SV-BPN catalysts. This work provides significant insights into designing cost-effective and high-performance bifunctional and trifunctional electrocatalysts.

2. Computational details

2.1 Methods

Density functional theory (DFT) calculations were performed using the Vienna Ab Initio Simulation Package (VASP).^43,44 The electron–ion interactions were described by projector-augmented wave (PAW) pseudopotentials⁴⁵ and the electron–electron exchange–correlation interactions were treated using the Perdew–Burke–Ernzerhof (PBE) functionals.⁴⁶ An energy cutoff of 500 eV was applied for the plane-wave expansion of wavefunctions. The Brillouin zone integration employed a 5 × 5 × 1 Monkhorst–Pack k-point grid for the structural optimizations, while a denser 10 × 10 × 1 k-point grid was utilized for electronic structure analyses.⁴⁷ All atomic positions were fully optimized until the residual forces on each atom were below 0.02 eV Å⁻¹. Spin-polarized calculations were systematically conducted for all investigated systems. To minimize artificial periodic interactions, a vacuum spacing of 20 Å was introduced perpendicular to the TM–BPN sheets. Dispersion corrections were accounted for using the Grimme DFT-D3 scheme to include long-range van der Waals interactions.⁴⁸ However, comparison of results obtained with and without DFT-D3 corrections revealed negligible differences, typically within ∼0.2 eV (see Table S1). The Gibbs free energy changes (ΔG) under acidic conditions (pH = 0) were calculated based on the computational hydrogen electrode (CHE) model (Note S1). Additional calculation details (Notes S2–S8) are provided in the SI. The Bader charge analysis was employed to determine charge transfer between TM atoms and the substrate.⁴⁹ The interaction strength between catalytically active centers and adsorbed intermediates was characterized via the crystal orbital Hamilton population (COHP) method.⁵⁰ To evaluate the kinetic stability, the climbing image nudged elastic band (CI-NEB) method was used to calculate the diffusion energy barrier of the TM atoms on the substrate.⁵¹ Furthermore, ab initio molecular dynamics (AIMD) simulations were performed to evaluate the thermodynamic stability of the catalysts.⁵² Finally, the post-processing analysis of thermodynamic quantities associated with catalytic reactions was carried out using the VASPKIT software package.⁵³

2.2 Model

Based on the HF-zipping reaction, a two-dimensional (2D) biphenylene network (BPN) was successfully synthesized by Fan et al.³⁷ The BPN structure comprises two distinct carbon atom types, labeled C₁ and C₂ (Fig. S1a). Unlike graphene, which consists exclusively of hexagonal rings, the BPN lattice incorporates tetragonal, hexagonal, and octagonal rings (Fig. S1a). Two possible single vacancies (SVs), V_C1 and V_C2, can form on pristine BPN, as illustrated by their initial and optimized structures in Fig. S1b and c, respectively. To evaluate the likelihood of forming these vacancies, their formation energies (E_vf) were calculated according to eqn (1):

ΔE_vf = E_v + µ_C − E_BPN (1)

Here, E_v represents the total energy of the defective BPN (either C₁ or C₂ vacancy). E_BPN is the total energy of BPN and µ_C denotes the chemical potential of the carbon atom, defined as the average total energy per carbon atom in graphene. Calculated formation energies for V_C1 and V_C2 in 2 × 2 and 3 × 3 supercells are provided in Table S2. The computed E_vf values for V_C1 and V_C2 are 5.029 and 6.543 eV, respectively. For comparison, the formation energy of a single vacancy in graphene was calculated to be 7.69 eV, significantly higher than those of SV-BPN substrates.⁵⁴ Consequently, the results indicate that the C₁ vacancy (V_C1) is energetically more favorable and stable than the C₂ vacancy (V_C2). Generally, lower vacancy formation energies imply higher feasibility for vacancy formation. Previous studies also suggest a preference for anchoring heteroatoms such as B/N atoms at the C₁ vacancy.⁵⁵ Therefore, the V_C1 substrate was selected for anchoring the transition metal (TM) atoms in this study.

The computational models employed supercells comprising 2 × 2 × 1 unit cells of the biphenylene network (BPN), as illustrated in Fig. 1a. Various transition metals (TMs), as listed in Fig. 1b, were introduced into the C1 vacancy site (V_C1) to form different TM@SV-BPN configurations. The formation energy (E_form) of each TM@SV-BPN was calculated using eqn (2).

ΔE_form = E_bind − E_coh (2)

ΔE_bind = E_TM@SV-BPN − E_SV-BPN − E_TM(g) (3)

ΔE_coh = E_TM(bulk)/n − E_TM(g) (4)

Here E_TM@SV-BPN and E_SV-BPN represent the total energy of the TM doped SV-BPN and the pristine defective BPN (C₁ vacancy), respectively. E_TM(bulk) denotes the energy of the TM in its bulk state, and E_TM(g) is the energy of an isolated TM atom in the gas phase. The parameter n refers to the number of metal atoms in the bulk structure.

Fig. 1 (a) Optimized structure of TM@SV-BPN and (b) transition metals considered in this work. (c) Binding energy (E_bind) and cohesive energy (E_coh) of the TM@SV-BPN. (d) Calculated formation energy (E_form) and (e) charge transfer for TM@SV-BPN (negative sign for charge refers to the electron-loss of TM atoms).

3. Results and discussion

3.1 Structural stability of TM@SV-BPN catalysts

The structural stability of TM@SV-BPN was investigated, and representative atomic structures are shown in Fig. 1a. The vacancy formation energy of the C₁ site in BPN was calculated to be 5.029 eV. Compared to the monovacancy formation energy in graphene (7.69 eV), the lower energy in BPN suggests that C₁ vacancies are more likely to form, highlighting the feasibility of vacancy engineering in BPN. To evaluate the potential aggregation tendency of TM atoms on the SV-BPN substrate, both binding energy (E_bind) and cohesive energy (E_coh) were computed. The bond lengths between TM atoms and adjacent carbon atoms in TM@SV-BPN are summarized in Table S3 (Fig. S2). These bond lengths are consistent with values reported for TM-doped SV-graphene, confirming the structural reliability of our models. E_bind and E_coh were calculated using eqn (3) and (4), and the results are presented in Fig. 1c and Table S4. TMs with a d¹⁰ configuration (Zn, Cd, and Hg) exhibit weak binding energies (less than −1 eV), indicating poor anchoring ability; these were excluded from further analysis. Cohesive energy reflects the intrinsic aggregation tendency of the metal atoms, while formation energy (E_form = E_bind − E_coh) serves as an indicator of thermodynamic stability. A negative E_form suggests a preference for isolated TM atoms over clustering. However, a slightly positive E_form does not necessarily preclude stability, as kinetic barriers may still prevent aggregation. Our calculated E_form values span from −1.14 to 2.20 eV. After excluding the d¹⁰ metals, 26 SAC candidates remain. Among them, 13 TM@SV-BPN systems (Sc, V, Cr, Mn, Fe, Co, Ni, Y, Ru, Rh, Pd, Ir, and Pt) exhibit negative E_form values (−0.14 to −1.14 eV), indicating strong anchoring. The other 13 systems (Ti, Cu, Zr, Nb, Mo, Tc, Ag, Hf, Ta, W, Re, Os, and Au) show moderately positive E_form values (0.06 to 2.20 eV). Notably, the highest E_form (2.20 eV for Nb@SV-BPN) is still below the formation energy of experimentally realized Pd/g-C₃N₄ (2.4 eV), suggesting good prospects for experimental synthesis. In order to evaluate the metal migration on TM@SV-BPN, we performed nudged elastic band calculations along the metal atom diffusion pathway to the nearest hexagonal hole on the catalysts. We chose Mn@SV-BPN (E_form = 0.78 eV), Nb@SV-BPN (E_form = 2.20 eV) and Mo@SV-BPN (E_form = 0.41) catalysts for further study. The energy barriers are as high as 1.00, 0.62 and 2.19 eV for Mn, Nb and Mo (Fig. S3). The high energy barriers suggest that these single TM atoms can remain in the vacancy site of biphenylene (SV-BPN). These results support the feasibility of TM@SV-BPN structures and justify further investigation of their electrocatalytic properties.

To investigate the interaction between transition metal (TM) atoms and the SV-BPN substrate, Bader charge analysis was performed. As illustrated in Fig. 1e and summarized in Table S4, charge accumulation at the TM–substrate interface confirms the formation of TM–C bonds in TM@SV-BPN, consistent with the previously discussed negative binding energies. The amount of charge transfer decreases across each period from left to right, correlating well with the increasing electronegativity of the TM atoms. Additionally, the net positive charge on the TM atoms is advantageous for the adsorption of electrochemical intermediates, as it facilitates electronic interactions with negatively charged species. To further evaluate the electronic characteristics, partial density of states (PDOS) calculations were carried out to determine the d-band center (ε_d) of the TM atoms. As shown in Fig. S4 and detailed in Table S4, the ε_d shifts toward lower energies relative to the Fermi level with increasing d-orbital occupation across the TM series. In general, a more negative ε_d indicates weaker binding between the TM and the adsorbates, while a less negative or positive ε_d corresponds to stronger adsorption strength. This relationship provides valuable insight into the catalytic activity trends of TM@SV-BPN systems.

3.2 HER performance on TM@SV-BPN

The HER and the OER are the two half-cell processes of electrochemical water splitting, both of which require efficient catalysts to reduce the associated electrochemical overpotentials. Among them, the HER serving as the cathodic half-reaction has been extensively investigated in recent years. Here, we focus on evaluating the HER electrocatalytic performance of the TM@SV-BPN catalysts under acidic conditions. In an acidic electrolyte, protons (H⁺) are adsorbed onto the catalyst surface, forming *H, where * represents the active site. This initial adsorption step is referred to as the Volmer step. Subsequently, the adsorbed hydrogen (*H) combines with another proton and an electron to form a hydrogen molecule (H₂), which is known as the Heyrovsky step (the Tafel pathway is not considered in this study). The overall HER mechanism, based on the Volmer–Heyrovsky pathway, is summarized by eqn (S4) and (S5) in the SI.

The computational evaluation of HER activity typically relies on a three-state model consisting of the initial state (H⁺ + e⁻), the adsorbed state (intermediate *H), and the final state (1/2H₂). At 298 K, the energy of the H⁺ + e⁻ pair can be approximated by half the energy of a hydrogen molecule (½H₂).⁵⁶ The Gibbs free energy of hydrogen adsorption (ΔG_*H) serves as a key descriptor for HER performance. According to the Sabatier principle, an optimal HER electrocatalyst should exhibit thermoneutral adsorption, i.e., ΔG_*H ≈ 0,⁵⁷ as excessively strong or weak hydrogen binding would hinder either adsorption or desorption, thus reducing the overall catalytic efficiency. In this study, we computed the ΔG_*H values (as defined in eqn (S7)) for hydrogen atoms adsorbed on TM sites of the SV-BPN monolayers to assess their HER activity. The resulting ΔG_*H values span a range of −0.507 to 0.693 eV (see Table S5), as illustrated in Fig. 2a and S5. Notably, Mo@SV-BPN and Tc@SV-BPN exhibit |ΔG_*H| values lower than that of commercial Pt (0.09 eV),⁵⁸ suggesting their potential as highly efficient, non-precious HER catalysts. Several TM@SV-BPN systems demonstrated |ΔG_*H| values close to the optimal 0.2 eV threshold, indicating favorable HER kinetics. These include Mo (0.006 eV), Tc (−0.013 eV), Mn (−0.11 eV), V (0.14 eV), Au (0.15 eV), Pt (−0.16 eV), Fe (0.16 eV), Ir (−0.17 eV), and Ru (0.17 eV). These results highlight their promise as stable and effective HER electrocatalysts.

Fig. 2 (a) The Gibbs free energy of HER and (b) exchange current density (i₀) of TM@SV-BPN.

To further identify high-performance HER catalysts, a volcano plot was constructed by plotting the exchange current density (i₀) as a function of the hydrogen adsorption free energy (ΔG_*H), as illustrated in Fig. 2b. In this relationship, catalysts situated near the apex of the volcano correspond to ΔG_*H values approaching zero, indicating an optimal balance between hydrogen adsorption and desorption. In contrast, catalysts on the left or right slopes of the volcano exhibit overly strong or weak hydrogen binding, respectively, leading to lower catalytic activity. Among the evaluated systems, Mo@SV-BPN and Tc@SV-BPN are positioned close to the summit of the volcano curve, consistent with their near-zero ΔG_*H values. This placement signifies their ability to achieve high exchange current densities and supports their classification as highly efficient HER electrocatalysts. These findings further validate the use of ΔG_*H as a reliable descriptor for HER activity and highlight the potential of TM@SV-BPN materials as viable alternatives to platinum-based catalysts.

3.3 OER performance on TM@SV-BPN catalysts

To systematically evaluate the OER performance of TM@SV-BPN single-atom catalysts, we conducted a comprehensive computational screening to identify the most promising OER-active systems. The OER mechanism involves four sequential proton-coupled electron transfer (PCET) steps, as detailed in Note S4. Initially, a water molecule adsorbs on the catalyst surface, releasing a proton–electron pair to form *OH (eqn (S11)). This is followed by the formation of *O through deprotonation of *OH (eqn (S12)). A subsequent water molecule reacts with *O to generate *OOH (eqn (S13)), which ultimately decomposes to release O₂ and regenerate the active site (eqn (S14)). Structural optimizations reveal that all oxygen-containing intermediates preferentially bind atop the TM atom, forming TM–O bonds analogous to those observed in metallic⁵⁹ and metal oxide systems.⁶⁰ The Gibbs free energies of the adsorbed intermediates (ΔG_*OH, ΔG_*O, and ΔG_*OOH) were calculated as described in Note S6 and summarized in Table S6. The rate-determining step (RDS) was identified by comparing the Gibbs free energy change associated with each step. Ideally, the energy difference between any two intermediates would be 1.23 eV, corresponding to the theoretical thermodynamic potential for OER under standard conditions.⁵⁶ Free energy diagrams at U = 0 and U = 1.23 V were constructed for all TM@SV-BPN systems (Fig. 3 and S6), and the potential-determining step (PDS) was identified for each catalyst. The associated free energy changes (ΔG₁ to ΔG₄) are listed in Table S7, with overpotentials (η^OER) determined according to Note S7. The η^OER values are plotted in Fig. 3b and tabulated in Table S8. Analysis of the reaction energetics reveals distinct trends in intermediate adsorption and PDS assignment. In the first step (H₂O → OH + H⁺), 15 catalysts namely V, Cr, Mn, Fe, Co, Ni, Nb, Mo, Tc, Ru, W, Re, Os, Ir, and Pt@SV-BPN, exhibit negative ΔG_*OH values (−0.04 to −1.93 eV), indicating strong interactions between the TM center and OH species. In contrast, Cu, Rh, Pd, Ag, and Au@SV-BPN show positive ΔG_*OH values, suggesting weaker *OH binding. In the second step (*OH → *O + H⁺), positive ΔG_*O values were observed for Cr, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ag, Os, Ir, Pt, and Au@SV-BPN (0.07 to 2.34 eV), whereas V, Mn, Nb, Mo, Tc, W, and Re@SV-BPN exhibited negative values. Notably, this step was identified as the PDS for Cu, Pd, Ag, and Au@SV-BPN. The third step, the formation of *OOH from *O, presents the highest energy barrier in several systems including Cr, Co, Mo, Tc, Ru, Rh, Re, Ir, and Pt@SV-BPN, highlighting the difficulty of O–O bond formation in these cases. For catalysts such as V, Mn, Fe, Ni, W, and Os@SV-BPN, the final step (*OOH → O₂ + H⁺) was found to be the PDS, primarily due to the challenge in *OOH desorption. The calculated η^OER values range from 0.43 to 1.87 V across the series: V (1.65 V), Cr (1.23 V), Mn (1.49 V), Fe (0.91 V), Co (0.70 V), Ni (0.56 V), Cu (0.91 V), Mo (1.54 V), Tc (1.74 V), Ru (0.82 V), Rh (0.65 V), Pd (0.43 V), Ag (1.11 V), W (1.87 V), Re (1.83 V), Os (1.33 V), Ir (1.04 V), Pt (0.89 V), and Au (0.81 V). Remarkably, Pd@SV-BPN and Ni@SV-BPN display overpotentials of 0.43 V and 0.56 V, respectively—comparable to or even lower than that of IrO₂(110) (0.56 V), as reported in previous theoretical studies,⁶⁰ indicating their exceptional OER catalytic performance.

Fig. 3 (a) The elementary reaction steps of the OER/ORR pathway on TM@SV-BPN. (b) OER and (d) ORR overpotential of TM@SV-BPN. The OER/ORR Gibbs free energy diagrams for (c) Pd@SV-BPN and (e) Ag@SV-BPN catalysts. The OER is the counter-clock wise, and the reverse direction refers to the ORR.

3.4 ORR performance on TM@SV-BPN catalysts

To further assess the electrocatalytic properties of TM@SV-BPN systems, we investigated their ORR activity and reaction pathways, as detailed in Note S5. The ORR typically proceeds via two mechanistic pathways: (1) a four-electron (4e⁻) pathway in which O₂ is reduced to H₂O through either an associative or dissociative mechanism, and (2) a two-electron (2e⁻) pathway leading to the formation of H₂O₂. According to previous studies, a ΔG_*O value less than 3.52 eV is indicative of thermodynamically favorable O₂ reduction via the 4e⁻ pathway to produce H₂O.⁶¹ As summarized in Table S6, all studied TM@SV-BPN catalysts exhibit ΔG_*O values below this threshold, suggesting a strong thermodynamic preference for the 4e⁻ pathway. These studied catalysts are more inclined to undergo the 4e⁻ ORR pathway to generate H₂O, this result can also be verified experimentally by in situ techniques.⁶² Given the absence of adjacent active sites on single-atom catalysts, the dissociative mechanism is unlikely to occur,⁶¹ making the associative 4e⁻ pathway the most plausible under acidic conditions (pH = 0). This reaction can be considered the reverse of the OER via the associative mechanism. The corresponding overpotential for ORR (η^ORR) was calculated according to eqn S39 (Note S7), which considers four key adsorption and desorption steps. Among the investigated systems, TM@SV-BPN catalysts such as V, Cr, Mn, Fe, Co, Ni, Mo, Tc, Ru, W, Re, Os, Ir, and Pt exhibit negative ΔG_*OH values and relatively large η^ORR values (Fig. 3d). These high overpotentials reflect strong binding of *OH intermediates, which in turn hinders their subsequent desorption, thus lowering the overall catalytic efficiency. The specific overpotential values are as follows: V (2.28 V), Cr (1.90 V), Mn (2.13 V), Fe (1.61 V), Co (1.40 V), Ni (1.28 V), Mo (2.16 V), Tc (1.93 V), Ru (1.45 V), W (2.55 V), Re (2.29 V), Os (1.97 V), Ir (1.64 V), and Pt (1.33 V). These results indicate that these catalysts require significant driving energy and are therefore less favorable for ORR applications. In contrast, TM@SV-BPN systems with positive ΔG_*OH values, namely Cu, Rh, Pd, and Au, exhibit lower overpotentials, indicating weaker *OH adsorption and thus enhanced ORR activity. The corresponding η^ORR values are 1.13 V (Cu), 1.12 V (Rh), 1.08 V (Pd), and 0.92 V (Au), as shown in Fig. 3e. Among all candidates, Ag@SV-BPN exhibits the lowest ORR overpotential of 0.68 V, suggesting it to be the most promising ORR catalyst in this study. For Au@SV-BPN, the potential-determining step was identified as the third step (*OOH → *O), which imposes a moderate energy barrier on the overall reaction. Furthermore, we examined the influence of solvation effects and DFT + U corrections on the computed OER and ORR energetics (Tables S9–S11). The resulting changes in adsorption energies and overpotentials were minimal, demonstrating that these corrections do not substantially alter the reaction energetics or the qualitative catalytic trends.

3.5 Selectivity of bifunctional/multifunctional catalysts

We further evaluated the bifunctional and multifunctional catalytic performance of TM@SV-BPN systems, with overpotentials for HER, OER, and ORR summarized in Table S12. In the context of water splitting, catalysts capable of efficiently catalyzing both HER and OER are highly desirable.⁶³ Similarly, for electrochemical energy conversion systems such as fuel cells and metal–air batteries, effective bifunctional activity toward OER and ORR is essential, as both reactions govern the cathodic processes.⁶⁴ The development of catalysts with dual or multiple functionalities offers a promising strategy to simplify device architecture and reduce overall system cost by eliminating the need for separate catalysts for each half-reaction.⁶⁵ Thus, multifunctional catalysts represent a more efficient, cost-effective, and sustainable approach for next-generation clean energy technologies. Fig. S7a and b illustrate the bifunctional activities of TM@SV-BPN catalysts, with dotted reference lines denoting the performance benchmarks: 0.2 V for HER and 1.23 V for OER/ORR.^66,67 Our results show that Fe, Ru, Ir, Pt, and Au@SV-BPN exhibit bifunctional capability for water splitting (HER/OER), while Cu, Rh, Pd, Ag, Pt, and Au@SV-BPN demonstrate bifunctional behavior suitable for metal–air batteries (OER/ORR). To quantify the overall performance, we define the combined overpotential for water splitting as W–η^Bi = |η^HER| + |η^OER|, W–η^Bi ≤ 1.43 V, with a target threshold of ≤1.43 V, and for metal–air batteries as M–η^Bi = |η^OER| + |η^ORR|, M–η^Bi ≤ 2.46 V, with a target of ≤2.46 V. Based on these criteria, ten TM@SV-BPN catalysts, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ir, Pt, and Au, meet the W–η^Bi requirement for water splitting, while Co, Ni, Cu, Ru, Rh, Pd, Ag, Pt and Au@SV-BPN catalysts fulfill the M–η^Bi criterion for metal–air battery applications (Table S12). Among the candidates, Pd@SV-BPN, Ni@SV-BPN, Au@SV-BPN, and Ru@SV-BPN stand out as the most promising multifunctional catalysts, exhibiting low overpotentials across all three key reactions (HER, OER, and ORR), thus offering strong potential for integrated energy conversion and storage applications.

3.6 AIMD calculations

Structural stability is the basic prerequisite of electrocatalytic materials. Especially, the active center TM atom can be stably embedded in the substrate and resist aggregation, which is closely related to the experimental feasibility and the long-term catalytic activity of the SAC. For structural stability and practical application of the TM@SV-BPN catalysts, particularly those exhibiting promising HER, OER, and ORR activity (TM = V, Cr, Mn, Fe, Co, Ni Cu, Mo, Tc, Ru, Rh, Pd, Ag, W, Re, Os, Ir, Pt, and Au) under actual electrocatalytic conditions, the thermal stability at room temperature or higher is very important. We performed ab initio molecular dynamics (AIMD) simulations at 500 K to evaluate the thermal stability of TM@SV-BPN catalysts (TM = V, Cr, Mn, Fe, Co, Ni Cu, Mo, Tc, Ru, Rh, Pd, Ag, W, Re, Os, Ir, Pt, and Au). The results of AIMD simulations over a 4 ps time period at 500 K are given in Fig. S8. It is evident that the skeleton and structural integrity of the systems remains intact throughout the 4 ps simulation, indicating effective bonding between TM and biphenylene (@SVBPN). The high thermal stability findings confirm that TM@SV-BPN (V, Cr, Mn, Fe, Co, Ni Cu, Mo, Tc, Ru, Rh, Pd, Ag, W, Re, Os, Ir, Pt, and Au) catalysts maintain their structural integrity at elevated temperatures, thereby reinforcing their potential for practical application in electrochemical energy conversion and storage devices.

3.7 Origin of TM@SV-BPN activity

To elucidate the fundamental trends governing oxygen electrocatalysis, we systematically analyzed the scaling relationships among the adsorption free energies of key intermediates (OH, O, and OOH) involved in the OER and ORR processes. Scaling behavior is a well-established phenomenon in oxygen-based catalysis, where the adsorption energies of intermediates tend to vary synchronously across different catalyst surfaces. Accordingly, we examined the correlations between ΔG_*OH, ΔG_*O and ΔG_*OOH to gain insight into the intrinsic limitations of the TM@SV-BPN catalysts. As shown in Fig. 4a, strong linear correlations were observed. In particular, ΔG_*OOH exhibited an excellent linear relationship with ΔG_*OH, with a coefficient of determination (R²) of 0.98. This relationship can be expressed as: ΔG_*OOH = 0.93ΔG_*OH + 3.15. This high correlation suggests that ΔG_*OH can serve as a reliable predictor for ΔG_*OOH in this catalyst class. A secondary but still significant linear relationship was also observed between ΔG_*OH and ΔG_*O (R² = 0.89), indicating coupled behavior in oxygen intermediate binding. To further assess catalytic performance for OER, we explored the thermodynamic descriptor ΔG_*O–ΔG_*OH, which reflects the relative binding strength of oxygen intermediates and has been identified as a critical factor in determining overpotential. When plotting ΔG_*OH–ΔG_*O against the negative of the potential-determining step (PDS) for each catalyst, a volcano-type relationship was observed (Fig. 4b). Catalysts positioned near the volcano apex exhibit a balanced interaction with oxygen intermediates, neither too strong nor too weak, resulting in optimal OER activity. Among all TM@SV-BPN systems studied, Pd@SV-BPN stands out, exhibiting the lowest η^OER value (0.42 V), placing it closest to the volcano peak and indicating its superior catalytic performance. This volcano behavior is in line with the Sabatier principle, reinforcing the importance of fine-tuning intermediate adsorption strengths to minimize overpotentials. Thus, Pd@SV-BPN emerges as a highly promising electrocatalyst for the OER. In the case of the ORR, we evaluated the correlation between ΔG_*OH and η^ORR (Fig. 4c). An exceptionally strong linear relationship was observed, with a coefficient of determination R² = 1.00. This result confirms that ΔG_*OH is an excellent descriptor for predicting ORR activity on TM@SV-BPN catalysts, further supporting its use in rational catalyst design.

Fig. 4 (a) The scaling relations between the adsorption energy of intermediates (*OH, *O and *OOH) on TM@SV-BPN. (b) Volcanic plot of OER on TM@SV-BPN with η^OER vs. ΔG_*O–G_*OH. (c) ΔG_*OH vs. η^ORR. (d) Relationship between the d-band center and adsorption energy of intermediates (ΔG_*OH, ΔG_*O, and ΔG_*OOH).

To further rationalize the observed trends in intermediate adsorption, we examined the relationship between the electronic structure of the catalysts and the adsorption free energies of oxygen intermediates. Specifically, the d-band center (ε_d) of the active TM site was correlated with the adsorption free energies of *OH, *O, and *OOH, as shown in Fig. 4d. Moderate linear relationships were observed for all three intermediates, with coefficients of determination of R² = 0.76 for ΔG_*OH, R² = 0.60 for ΔG_*O, and R² = 0.76 for ΔG_*OOH. These findings are consistent with the d-band theory, which posits that a higher d-band center facilitates stronger adsorption due to increased overlap between the metal d orbitals and the adsorbate molecular orbitals. This analysis suggests that the position of the d-band center serves as an effective electronic descriptor for interpreting and predicting the catalytic behavior of TM@SV-BPN systems. To gain further insights into the bonding characteristics between the metal center and adsorbed intermediates, we also performed crystal orbital Hamilton population (COHP) analysis. This approach provides a quantitative evaluation of the bonding interactions, offering a complementary perspective to the electronic descriptors derived from density of states.^68,69 To further elucidate the bonding interactions governing catalytic performance, we conducted a crystal orbital Hamilton population (COHP) analysis focused on the *OH intermediate adsorption. Three representative TM@SV-BPN catalysts were selected: Pd@SV-BPN, which exhibited the best OER performance (η^OER = 0.43 V; η^ORR = 1.07 V); Ag@SV-BPN, the most efficient ORR catalyst η^OER = 1.11 V; η^ORR = 0.67 V); and Au@SV-BPN, which demonstrated trifunctional activity (η^HER = 0.15 V; η^OER = 0.81 V, η^ORR = 0.91 V). The projected COHP profiles describe energy-resolved bonding (negative values) and antibonding (positive values) interactions between the TM site and the OH adsorbate. As shown in Fig. S9a, the Pd–OH system displays a broad and deep bonding region below the Fermi level, indicative of strong metal–adsorbate interaction. This is quantitatively reflected by the most negative integrated COHP (ICOHP) value of −2.93 eV and a correspondingly low ΔG_*OH value of 0.15 eV, suggesting strong and favorable adsorption. For Au@SV-BPN (Fig. S9c), the Au–OH interaction also shows substantial bonding below the Fermi level, albeit slightly weaker than that of Pd. The ICOHP value of −2.84 eV and the moderate ΔG_*OH value of 0.31 eV indicate a balanced interaction conducive to multifunctional catalysis. In contrast, Ag@SV-BPN (Fig. S9b) exhibits a noticeably shallower and narrower bonding region. Its ICOHP value is significantly less negative at −2.15 eV, which correlates with the weakest OH adsorption energy among the three (ΔG_*OH = 0.76 eV). We further analyzed the quantitative relationship between ΔG_*OH and the corresponding ICOHP values across the catalyst series. As shown in Fig. S9d, a strong linear correlation was observed, with a coefficient of determination R² = 0.91. This result confirms that ICOHP serves as a reliable electronic descriptor for OH adsorption strength. Since ICOHP quantifies the integrated bonding interaction between the metal center and the adsorbate, more negative values signify stronger binding, which directly corresponds to more negative ΔG_*OH values.⁷⁰ These electronic structure insights provide a mechanistic understanding of the observed catalytic activity and offer a valuable basis for the rational design of high-performance single-atom catalysts through fine-tuning of metal–adsorbate bonding characteristics.

3.8 Machine learning

In recent years, machine learning (ML) approaches have been increasingly employed to accelerate the design of electrocatalysts and to provide mechanistic insights into electrocatalytic processes. In this study, ML was utilized to establish the relationship between the intrinsic properties of TM@SV-BPN catalysts and their catalytic performance. A set of descriptors encompassing both electronic and structural features was selected, as summarized in Table S13. These include atomic number (Z), atomic mass (M), electron affinity (EA), first ionization energy (E_i), atomic electronegativity (X), d-orbital electron count (N_e), atomic radius (r_d), charge transfer (Q_e), d-band center (ε_d), binding energy (E_bind), and formation energy (E_form). Among them, Z, M, EA, E_i, X, N_e, and r_d represent intrinsic atomic properties, while Q_e, ε_d, E_bind, and E_form were derived from DFT calculations. To evaluate the intercorrelation among these features, Pearson correlation matrices were constructed (Fig. S10). Subsequently, a gradient boosting regression (GBR) model was developed to predict the adsorption free energies of key oxygen intermediates, ΔG_*OH, ΔG_*O and ΔG_*OOH, traditionally obtained via computationally demanding DFT simulations. To evaluate numerical stability and minimize stochastic deviations inherent to ensemble-based learning algorithms, the GBR model was trained 1000 times under identical conditions. Each repetition used the same training and test sets but differed only in internal stochastic elements of the boosting framework. For each run, the coefficient of determination (R²), mean absolute error (MAE), and root mean square error (RMSE) were recorded to quantify predictive accuracy. The GBR model exhibited consistently strong performance, with the initial optimized model achieving R² = 0.986, MAE = 0.097, and RMSE = 0.115 for ΔG_*O. Across the 1000 independent repetitions, the average R² was 0.983 with a standard deviation of only 0.003, demonstrating excellent numerical stability. The mean MAE and RMSE values across repetitions were 0.1070 and 0.1169, respectively, both exhibiting narrow variances. These results confirm that model predictions are highly reproducible and that stochastic fluctuations have a negligible effect on predictive trends.

Feature selection was refined using the recursive feature elimination (RFE) algorithm, which iteratively removes the least informative descriptors based on model performance. The feature ranking process is illustrated in Fig. 5a for ΔG_*O, Fig. S11a for ΔG_*OH, and Fig. S12a for ΔG_*OOH, with corresponding Pearson correlations shown in Fig. 5b, S11b and S12b, respectively. The dataset was randomly split into a training set (90%) and a test set (10%). The GBR algorithm was chosen due to its capacity to capture complex, nonlinear relationships between input features and target values. A well-performing model is characterized by a coefficient of determination (R²) close to 1, and low values of root mean square error (RMSE) and mean absolute error (MAE). The ML-predicted adsorption energies exhibit excellent agreement with DFT-calculated values, as shown in Fig. 5c (ΔG_*O), Fig. S11c (ΔG_*OH), and Fig. S12c (ΔG_*OOH). Performance metrics, including R², RMSE, and MAE, are summarized in Table S14. Notably, the model predicting ΔG_*O achieved outstanding accuracy, with R² values of 1.00 and 0.99 for the test and training sets, respectively. To further interpret feature importance, the contribution of each descriptor to model performance was evaluated. Among all descriptors, charge transfer (Q_e) emerged as the most influential feature, contributing 81% for ΔG_*O, 80% for ΔG_*OH, and 55% for ΔG_*OOH predictions. These findings highlight the significance of charge redistribution at the active site in governing oxygen intermediate adsorption and catalytic performance. In summary, the integration of ML techniques provides a powerful framework for elucidating structure–activity relationships and enables accurate prediction of catalytic behavior, thereby facilitating the rational design of advanced SACs such as TM@SV-BPN.

Fig. 5 (a) The number of selected features vs. features ranking. (b) Pearson correlation between selected features. (c) Calculated ΔG_*O vs. ML predicted ΔG_*O by the GBR algorithm. (d) Feature importance of the GBR model for ΔG_*O.

4. Conclusions

In summary, we have demonstrated that biphenylene-based single-atom catalysts (TM@SV-BPN) provide a highly versatile and thermodynamically stable platform for efficient electrocatalysis across HER, OER, and ORR. The combination of the unique carbon framework geometry, strong metal anchoring ability, and tunable electronic structure enables a diverse range of catalytic behaviors. Among the 39 evaluated TM@SV-BPN systems, several exhibit excellent single- and bifunctional performance, while Au@SV-BPN emerges as a rare trifunctional catalyst. Our in-depth analysis using d-band theory, COHP, and machine learning models reveals robust correlations between structural/electronic descriptors and catalytic activity, enabling predictive screening beyond traditional trial-and-error approaches. The identification of charge transfer as a dominant descriptor in ML models further provides a mechanistic foundation for future high-throughput catalyst discovery. This work not only expands the fundamental understanding of SACs on nonbenzenoid carbon supports but also introduces a predictive, data-guided strategy for developing cost-effective multifunctional electrocatalysts. These findings hold broad implications for the rational design of next-generation energy conversion systems, including electrolyzers, fuel cells, and metal–air batteries.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are available within the article and its supplementary information (SI). Supplementary information: additional computational details; computational hydrogen electrode (CHE) model (Note S1), HER (Note S2), theoretical exchange current i₀ (Note S3), OER (Note S4), ORR (Note S5), adsorption energy calculations for OER/ORR (Note S6), theoretical overpotential calculations for OER/ORR (Note S7), and d-band center (Note S8). Optimized biphenylene with and without single vacancy (Fig. S1 and S2). The diffusion barriers on @SV-BPN (Fig. S3). The ε_d of 3d TM, 4d TM and 5d TM (Fig. S4). The free energy diagram of HER (Fig. S5). The Gibbs free energy diagram of the OER/ORR on TM@SVBPN catalysts (Fig. S6). Overpotential for water-splitting (η^HER/η^OER) and for metal–air batteries (η^OER/η^ORR) (Fig. S7). AIMD simulations at temperature of 500 K for TM@SVBPN catalysts (Fig. S8). pCOHP of Pd@, Ag@, and Au@SV-BPN with reaction intermediate *OH and the relationship between ΔG_*OH and ICOHP (Fig. S9). Pearson correlation of 11 feature variables (Fig. S10). ML predicted results for ΔG_*OH and ΔG_*OOH (Fig. S11 and S12). Adsorption energy of reaction intermediates with and without dispersion correction on the Pd@SV-BPN catalyst. (Table S1); defective C₁ and/or C₂ biphenylene vacancy formation energies (Table S2), bond length between the TM atom and SV-BPN substrate (Table S3); E_form and E_bind of TM atoms doped on SV@BPN, E_coh of TM atoms in bulk, charge transfer (e⁻) between TM atoms and SV@BPN and ε_d (Table S4); adsorption energy of ΔG_*H on TM@SV-BPN catalysts (Table S5), adsorption energies of *OH, *O and *OOH (ΔG_*OH, ΔG_*O and ΔG_*OOH) on TM@SV-BPN catalysts (Table S6), free energy change of each OER elementary step on TM@SV-BPN (Table S7); theoretical over-potentials for OER (η^OER) and ORR (η^ORR) on TM@SV-BPN catalysts (Table S8); solvent effect (Table S9), magnetic moment DFT vs. DFT + U (Table S10), DFT + U results (Table S11); bifunctional over-potential for water-splitting (W–η^Bi) and metal–air batteries (M–η^Bi) (Table S12); features details for ML approaches (Table S13); evaluated R², MAE, and RMSE metrics to predict ΔG_*OH, ΔG_*O and ΔG_*OOH for test and train data (Table S14). See DOI: https://doi.org/10.1039/d5ta10121b.

Additional raw data are available from the corresponding author upon reasonable request.

Acknowledgements

H.-T. Chen gratefully thanks Chung Yuan Christian University (CYCU), the National Center for Theoretical Sciences (NCTS), and the National Science and Technology Council (NSTC), Taiwan, under Grant Numbers, NSTC 114-2113-M-033-002, 114-2113-M-007-026, 113-2113-M-033-007, and 113-2113-M-007-027 for the financial support and the National Center for High-performance Computing (NCHC), Taiwan, for providing computational and storage resources.

References

1. S. Lin, H. Xu, Y. Wang, X. C. Zeng and Z. Chen, J. Mater. Chem. A, 2020, 8, 5663–5670.
2. Z. P. Cano, D. Banham, S. Ye, A. Hintennach, J. Lu, M. Fowler and Z. Chen, Nat. Energy, 2018, 3, 279–289.
3. D. Larcher and J.-M. Tarascon, Nat. Chem., 2015, 7, 19–29.
4. X. Lv, Z. Xiao, H. Wang, X. Wang, L. Shan, F. Wang, C. Wei, X. Tang and Y. Chen, J. Energy Chem., 2021, 54, 626–638.
5. F. Dawood, M. Anda and G. Shafiullah, Int. J. Hydrogen Energy, 2020, 45, 3847–3869.
6. E. Vijayakumar, S. Ramakrishnan, C. Sathiskumar, D. J. Yoo, J. Balamurugan, H. S. Noh, D. Kwon, Y. H. Kim and H. Lee, Chem. Eng. J., 2022, 428, 131115.
7. N. Ma, C. Gong, H. Xie, C. Shi, J. Sha, C. He, F. He, N. Zhao and E. Liu, Int. J. Hydrogen Energy, 2022, 47, 29762–29770.
8. J. Quílez-Bermejo, S. García-Dalí, A. Daouli, A. Zitolo, R. L. Canevesi, M. Emo, M. T. Izquierdo, M. Badawi, A. Celzard and V. Fierro, Adv. Funct. Mater., 2023, 33, 2300405.
9. J. Wang, Y. Fan, S. Qi, W. Li and M. Zhao, J. Phys. Chem. C, 2020, 124, 9350–9359.
10. Q. Wu, Y. Ma, H. Wang, S. Zhang, B. Huang and Y. Dai, ACS Appl. Mater. Interfaces, 2020, 12, 24066–24073.
11. J. Yang, L. Chang, H. Guo, J. Sun, J. Xu, F. Xiang, Y. Zhang, Z. Wang, L. Wang and F. Hao, J. Mater. Chem. A, 2020, 8, 1229–1237.
12. Y. Zhou, G. Gao, J. Kang, W. Chu and L.-W. Wang, Nanoscale, 2019, 11, 18169–18175.
13. G. Wang, C. Chen, B. S. Teketel, B. Xu and B. Lin, Int. J. Hydrogen Energy, 2021, 46, 39183–39194.
14. Z. W. Seh, J. Kibsgaard, C. F. Dickens, I. Chorkendorff, J. K. Nørskov and T. F. Jaramillo, Science, 2017, 355, eaad4998.
15. M. D. Bhatt and J. Y. Lee, Energy Fuels, 2020, 34, 6634–6695.
16. R. Qin, G. Chen, X. Feng, J. Weng and Y. Han, Adv. Sci., 2024, 11, 2309364.
17. Y.-J. Wang, N. Zhao, B. Fang, H. Li, X. T. Bi and H. Wang, Chem. Rev., 2015, 115, 3433–3467.
18. Y. Yang, J. Zhou, Z. Zhao, G. Sun, S. Moniri, C. Ophus, Y. Yang, Z. Wei, Y. Yuan and C. Zhu, Nat. Catal., 2024, 7, 796–806.
19. T. Sadhasivam, G. Palanisamy, S.-H. Roh, M. D. Kurkuri, S. C. Kim and H.-Y. Jung, Int. J. Hydrogen Energy, 2018, 43, 18169–18184.
20. B. Qiao, A. Wang, X. Yang, L. F. Allard, Z. Jiang, Y. Cui, J. Liu, J. Li and T. Zhang, Nat. Chem., 2011, 3, 634–641.
21. C. Jia, Q. Sun, R. Liu, G. Mao, T. Maschmeyer, J. J. Gooding, T. Zhang, L. Dai and C. Zhao, Adv. Mater., 2024, 36, 2404659.
22. R. Lang, X. Du, Y. Huang, X. Jiang, Q. Zhang, Y. Guo, K. Liu, B. Qiao, A. Wang and T. Zhang, Chem. Rev., 2020, 120, 11986–12043.
23. B. Hinnemann, P. G. Moses, J. Bonde, K. P. Jørgensen, J. H. Nielsen, S. Horch, I. Chorkendorff and J. K. Nørskov, J. Am. Chem. Soc., 2005, 127, 5308–5309.
24. H.-Y. Zhuo, X. Zhang, J.-X. Liang, Q. Yu, H. Xiao and J. Li, Chem. Rev., 2020, 120, 12315–12341.
25. H.-C. Huang, Y. Zhao, J. Wang, J. Li, J. Chen, Q. Fu, Y.-X. Bu and S.-B. Cheng, J. Mater. Chem. A, 2020, 8, 9202–9208.
26. X. Ma, J. Wen, S. Zhang, H. Yuan, K. Li, F. Yan, X. Zhang and Y. Chen, ACS Sustain. Chem. Eng., 2017, 5, 10266–10274.
27. H.-T. Chen, Y.-T. Chiou, T.-H. Chen and H.-L. Chen, Chem. Phys., 2025, 591, 112540.
28. K.-Y. Chiu, C.-W. Chan and H.-T. Chen, Electrochim. Acta, 2025, 511, 145389.
29. N. Sathishkumar and H.-T. Chen, ACS Appl. Mater. Interfaces, 2023, 15, 15545–15560.
30. D. K. D. Chittibabu, N. Sathishkumar, S.-Y. Wu and H.-T. Chen, ACS Appl. Energy Mater., 2023, 6, 6636–6645.
31. N. Sathishkumar, S.-Y. Wu and H.-T. Chen, Appl. Surf. Sci., 2022, 598, 153829.
32. R. Krishnan, S.-Y. Wu and H.-T. Chen, Phys. Chem. Chem. Phys., 2019, 21, 12201–12208.
33. Y. Wang, R. Hu, Y. Li, F. Wang, J. Shang and J. Shui, Nano Res., 2022, 15, 1054–1060.
34. C. Xia, Y. Qiu, Y. Xia, P. Zhu, G. King, X. Zhang, Z. Wu, J. Y. Kim, D. A. Cullen and D. Zheng, Nat. Chem., 2021, 13, 887–894.
35. N. Ullah, R. Ullah, S. Khan and Y. Xu, Front. Mater. Sci., 2021, 1–10.
36. A. Raza, J. Z. Hassan, U. Qumar, A. Zaheer, Z. U. D. Babar, V. Iannotti and A. Cassinese, Mater. Today Adv., 2024, 22, 100488.
37. Q. Fan, L. Yan, M. W. Tripp, O. Krejčí, S. Dimosthenous, S. R. Kachel, M. Chen, A. S. Foster, U. Koert and P. Liljeroth, Science, 2021, 372, 852–856.
38. K. Pandiyan, D. K. D. Chittibabu and H.-T. Chen, ACS Appl. Energy Mater., 2024, 7, 10758–10769.
39. D. K. D. Chittibabu and H.-T. Chen, Electrochim. Acta, 2024, 497, 144578.
40. J. Hao, Z. Zhao, C. Chen, C. Zhang, L. Li, S. Gao, B. Jia and P. Lu, Sustain. Energy Fuels, 2022, 6, 3446–3452.
41. T. Liu, Y. Jing and Y. Li, J. Phys. Chem. Lett., 2021, 12, 12230–12234.
42. Z. Wang, Z. Yin, Y. Gao, H. Wang, J. Gao and J. Zhao, J. Mater. Chem. A, 2024, 12, 2748–2759.
43. J. Hafner, J. Comput. Chem., 2008, 29, 2044–2078.
44. G. Sun, J. Kürti, P. Rajczy, M. Kertesz, J. Hafner and G. Kresse, J. Mol. Struct.:THEOCHEM, 2003, 624, 37–45.
45. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953.
46. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
47. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188.
48. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
49. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
50. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, LOBSTER: a tool to extract chemical bonding from plane‐wave based DFT, Wiley Online Library, 2016.
51. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
52. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558.
53. V. Wang, N. Xu, J.-C. Liu, G. Tang and W.-T. Geng, Comput. Phys. Commun., 2021, 267, 108033.
54. Y. Ma, P. Lehtinen, A. S. Foster and R. M. Nieminen, New J. Phys., 2004, 6, 68.
55. Z. Feng, B. Zhang, R. Li, F. Li, Z. Guo, S. Zheng, G. Su, Y. Ma, Y. Tang and X. Dai, J. Power Sources, 2023, 558, 232613.
56. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jonsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
57. P. Sabatier, Ber. Dtsch. Chem. Ges., 1911, 44, 1984–2001.
58. H. Xu, D. Cheng, D. Cao and X. C. Zeng, Nat. Catal., 2024, 7, 207–218.
59. J. Rossmeisl, A. Logadottir and J. K. Nørskov, Chem. Phys., 2005, 319, 178–184.
60. J. Rossmeisl, Z.-W. Qu, H. Zhu, G.-J. Kroes and J. K. Nørskov, J. Electroanal. Chem., 2007, 607, 83–89.
61. X. Guo, S. Lin, J. Gu, S. Zhang, Z. Chen and S. Huang, ACS Catal., 2019, 9, 11042–11054.
62. K. Dong, J. Liang, Y. Wang, L. Zhang, Z. Xu, S. Sun, Y. Luo, T. Li, Q. Liu and N. Li, ACS Catal., 2022, 12, 6092–6099.
63. M. Chatenet, B. G. Pollet, D. R. Dekel, F. Dionigi, J. Deseure, P. Millet, R. D. Braatz, M. Z. Bazant, M. Eikerling and I. Staffell, Chem. Soc. Rev., 2022, 51, 4583–4762.
64. L. Du, L. Xing, G. Zhang, M. Dubois and S. Sun, Small Methods, 2020, 4, 2000016.
65. Y. Jiao, Y. Zheng, M. Jaroniec and S. Z. Qiao, Chem. Soc. Rev., 2015, 44, 2060–2086.
66. S. Ram, A. S. Lee, S.-C. Lee and S. Bhattacharjee, Chem. Mater., 2025, 37, 3608–3621.
67. B. Wei, Z. Fu, D. Legut, T. C. Germann, S. Du, H. Zhang, J. S. Francisco and R. Zhang, Adv. Mater., 2021, 33, 2102595.
68. R. Dronskowski and P. E. Bloechl, J. Phys. Chem., 1993, 97, 8617–8624.
69. F. Li, H. Ai, D. Liu, K. H. Lo and H. Pan, J. Mater. Chem. A, 2021, 9, 17749–17759.
70. R. Hu, Y. Li, Q. Zeng and J. Shang, Appl. Surf. Sci., 2020, 525, 146588.

---