High-throughput screening and an interpretable machine learning model of single-atom hydrogen evolution catalysts with an asymmetric coordination environment constructed from heteroatom-doped graphdiyne

Abstract

Exploring high-activity and low-cost electrocatalysts for the hydrogen evolution reaction is the key to developing new energy sources, but it faces major challenges. Herein, a series of single atom catalysts with an asymmetric coordination environment constructed from heteroatom-doped graphdiyne are constructed to study HER activity. Based on DFT calculations, 120 configurations are screened and it is discovered that the ΔG_*H of Cr@N₂B1-GDY, Mn@N1B2-GDY, Fe@N1S2-GDY, Co@N1B4-GDY and Ni@N2-GDY was close to zero. Among them, Cr@N2B1-GDY has the lowest Gibbs free energy change (−0.0046 eV), exhibiting superior HER performance to known SACs of carbon-based supports. The constructed asymmetric coordination environment can reduce adsorption and enhance HER activity by adjusting the electron spin polarization of the central metal and the electronic structure of the D-band so that the e_g orbital is conducive to filling the antibonding orbital. The interpretable machine learning model's feature importance analysis further proves that the spin magnetic moment and the D-band center are highly correlated with the prediction accuracy of ΔG_*H. The hydrogen affinity of metal atoms also highly influences HER activity. RF is considered to be the best predictive machine learning model with R² of 96%. In this study, interpretive machine learning and DFT reveal that asymmetric doping of heteroatoms affects the D-band structure and spin states of the TM to regulate HER activity. It provides a novel way to establish the theory of HER SACs.

---

1 Introduction

Hydrogen is considered as a clean renewable energy source that can replace fossil fuels to address the increasingly serious problem of environmental pollution.^1–5 Currently, the hydrogen evolution reaction (HER), namely, electrocatalytic splitting of water, stands out among various hydrogen production technologies as a promising method for hydrogen production, because it generates high-purity hydrogen without any by-products during the reaction process.^6–8 Platinum (Pt) based catalysts have shown excellent performance for the HER;^9–11 however, the industrial implementation of these electrocatalysts has been thwarted by high cost and scarce resources.¹² Therefore, it is urgent that develop lower cost and more abundant electrocatalysts to enhance HER efficiency.^13–15

Since the successful stripping of graphene, two-dimensional (2D) materials have attracted much attention due to their unique layered structure and electronic properties,¹⁶ such as MXenes,^17,18 graphitic carbon nitride (g-C₃N₄),^19,20 transition metal dichalcogen compounds^21,22 and so on.^23–25 Recently, graphdiyne (GDY) has emerged as a novel electrocatalyst, possessing a distinctive structure of benzene rings with sp²-hybridized carbon atoms and diacetylenic linkages with sp-hybridized carbon atoms.²⁶ In experimental research, Jia et al. developed an in situ homocoupling reaction of hexaethynylbenzene on Cu foil for the preparation of large-area GDY films as early as 2010.²⁷ Besides, GDY also delivers abundant reactive sites for adsorption and catalysis, which shows high electronic conductivity and charge carrier mobility, rendering superior catalytic performance.^28–30 For instance, Xue et al. discovered that specific chemical interactions and electronic coupling between anchored single-atom Ni/Fe and GDY led to outstanding HER performance.³¹ In a recent study, Yu et al. synthesized a novel zero-valent Pd single-atom based on graphdiyne, and verified that the Pd⁰/GDY structure displayed remarkable HER activity with a smaller overpotential at 10 mA cm⁻².³² GDY offers significant advantages not only in HER catalytic performance, but also in various chemical reactions, such as the carbon dioxide reduction reaction (CO₂RR),²⁸ oxygen reduction reaction (ORR)³³ and nitrogen reduction reaction (NRR).³⁴

Moreover, adjusting the electronic properties of two-dimensional materials by doping atoms can effectively regulate the adsorption capacity between adsorbent atoms and the substrate. Single-atom catalysts (SACs) can improve the catalytic performance owing to their precise and well-defined active sites and extraordinary properties.^35–37 For example, Li et al. systematically calculated and investigated the catalytic properties of a series of single transition-metal atoms doped on a VSe₂ monolayer, and found that Pt@VSe₂ may be recognized as a triple-functional SAC (ΔG_*H = 0.11 eV, η^OER = 0.72 eV, η^ORR = 0.85 eV).³⁸ The coordination environments between the supported metal atoms and the substrate materials also affect the catalytic performance.³⁹ Researchers have extensively explored different non-metals (B, N, P and S) doped into various TM-SACs. Specifically, the catalytic activity of carbon-based iron–nitrogen systems (Fe@N_x) is affected by the coordination environment of single Fe atoms. Studied results identified that the Fe@N₂-opp configuration with two N atoms at the opposite sides of the Fe atom has high NRR catalytic performance, and the limiting potential is only −0.63 V.⁴⁰ Also, Umer et al. demonstrated that mono- or dual-type non-metals including B, N and P doped in diverse substrates influence the coordination environment and charge transfer behavior of the HER by means of the density functional theory (DFT) combined machine learning (ML) framework. Various types of descriptors of 364 candidate structures were selected to create the ML model to quickly predict catalyst stability and HER activity.⁴¹

Recently, ML methods have received prominent attention to provide an explanation for catalytic performance, and have significantly reduced the computational burden. All sorts of geometry and electronic descriptors have been proposed to interpret the design and screening of catalyst properties. For instance, Li et al. combined graph neural networks and crystal graph machine learning characterization with active learning to study Pt-based HER catalysts, without the help of expensive DFT geometry optimizations, and identified that Cu₃Pt and FeCuPt₂ show nearly identical catalytic properties to those of Pt (111).⁴² Moreover, relevant research adopted the active-learning method to predict the adsorption site of *CO, a pivotal CO₂ reduction reaction intermediate, and verified 6 promising catalysts with low limiting potential among 257 candidates.⁴³

In this study, we systematically studied the HER catalytic performance of GDY-based SACs with different central transition metal and surrounding atoms by DFT and ML. The C atoms in the GDY substrate were replaced by 1, 2 and 3 non-metallic atoms (B, N and S) respectively, and a transition metal atom (Cr, Mn, Fe, Co and Ni) embedded into the cavity, with a total of 120 potential components. This work employed high-throughput DFT calculations, through which the 5 best SACs were screened out, which exhibited good structural stability and suitable catalytic activity. Interpretable machine learning is used to construct predictions of HER activity and to explore structure–activity relationships. SHAP analysis is used to test the physicochemical principles of feature selection in ML models.

2 Computational methods

This subject used the Vienna Ab initio Simulation Package (VASP) to perform the first-principles calculations.⁴⁴ The generalized gradient approximation (GGA) was utilized to formulate the exchange-related potentials with Perdew–Burke–Ernzerhof (PBE).⁴⁵ The electron–ion interactions were described by the projector augmented wave (PAW) method.⁴⁶ The vacuum region was set at 20 Å to keep periodic images from interacting. The cutoff energy for the plane waves was set to 500 eV. The Brillouin zones were sampled using the Monkhorst–Pack method, in which a 3 × 3 × 1 mesh was used for structural relaxation and electronic properties. During the geometry optimization, the convergence tolerances for the electronic energy and the residual force were 1 × 10⁻⁴ eV and 0.02 eV Å⁻¹, respectively. In all the DFT calculations performed, all the energies take into account the kinetic energy of the electrons, the attraction of the electrons to the nucleus and the coulomb interaction energy between the electrons (potential energy), and the exchange correlation energy. In addition, the Ab initio molecular dynamics (AIMD) method was used to demonstrate the stability of the structure with the simulation time of 2 ps and the temperature of 500 K.

The formation energy of the GDY support structure doped by a heteroatom was calculated by E_{f_support} = E_support − E_HM −E_{defect_GDY}, where E_support, E_HM, and E_{defect_GDY} represent the energies of the GDY support doped with a heteroatom, the GDY with vacancy defects, heteroatom in the corresponding crystals, respectively. The formation energy of doping the GDY surface with a transition metal atom (E_{f, SAC}) was calculated by E_{f_SAC} = E_SAC − E_{defect_GDY} − E_TM − E_HM, where E_SAC, E_{defect_GDY}, E_TM, and E_HM represent the energies of the single-atom catalyst, the substrate of GDY with defects, transition metal atom and heteroatom in the corresponding crystals, respectively. A negative E_f value indicates that the metal atom could be stably loaded on the graphdiyne substrate surface.

For the HER, the first step of the reaction is the adsorption of a radical H* on the material (* + H⁺ + e⁻ → H*).⁴⁷ The Gibbs free energy change (ΔG) of the reaction step is based on the computational hydrogen electrode model (CHE) proposed by Nørskov et al., where the free energy of the (H⁺ + e⁻) pair is replaced by half of the chemical potential of the hydrogen molecule. The corresponding ΔG is calculated by ΔG = ΔE + ΔZPE − TΔS, where ΔE is the difference in self-consistent energy after H* adsorption, T is temperature (298.15 K), and ΔZPE and ΔS represent the change in zero-point energy and entropy, respectively. According to Nørskov's assumption,⁴⁸ the exchange current density (i₀) at pH = 0 is calculated by using the following equation: , where k₀ is the reaction rate constant, which was set to 1 under zero overpotential,⁴⁹ and K_B represents the Boltzmann constant.

Furthermore, ML is a branch of artificial intelligence where the analysis models are set up by the training data automatically, and decisions are made based on other data. Herein, to reveal the relationship between catalyst structure and catalytic activity, a machine learning (ML) model is constructed to investigate the effect of structural features on the catalytic activity of TM@HM-GDY using scikit-learn.

3 Results and discussion

3.1 Structures and stabilities of TM@HM-GDY

We first optimized the initial GDY configuration, and then, considering the symmetry of the GDY structure, carbon atoms at both ends of acetylene were selected as target sites to replace 1, 2 and 3 heteroatoms respectively, corresponding to positions C1, C2, C3 and C4 in Fig. 1a. We evaluate the stability by calculating the formation energy of HM doping, presented in Table S1.† It is obviously observed that the formation energy of the graphitylene substrates doped with heteroatoms is negative, which means that these structures are thermodynamically stable. These structures can be used as carriers to support metal atoms. Subsequently, 120 different TM@HM-GDY (shortened to TM@HM) structures were created, in which the transition metal atoms of Cr, Mn, Fe, Co and Ni were anchored on the GDY surface cavity. The formation energy (shown in Table S2†) is introduced to check the stability of these candidate catalysts. E_f for configurations such as Cr@B1S4 was not counted because of structural collapse during optimization. E_f of the remaining configurations is mostly negative except Cr@B₂, suggesting that this TM-HM atom-couple can be stably loaded on the GDY substrate. Table S2† shows that the coordination number of the central active metal atom with the surrounding atoms is 4 ∼ 6, thereby firmly anchoring the transition metal atom to GDY. After the stability check, the remaining 102 configurations are further screened for activity.

Fig. 1 (a) Top view of the structures of TM@NM-GDY (TM = Cr, Mn, Fe, Co and Ni; HM = N, B and S); (b) the optimized structure of Cr@N2B1-GDY as an example.

3.2 HER catalytic activity and stability

The hydrogen adsorption process is the first step of the HER, and appropriate chemisorption of *H is energetically favorable. Therefore, the adsorption energy of H on the candidate catalyst is calculated and shown in Table S3.† After the optimization *H static adsorption calculation is completed, the configuration of the TM with a change in the coordination environment is excluded. According to Table S3,† screening the structures with E_ads < 0 reveals that the TM–H bond lengths range from 1.445 ∼ 1.752 Å. In general, a shorter TM–H bond length generally suggests a stronger TM–H interaction. However, for the Cr@B2S1-GDY, the Cr–H bond length is 1.752 Å, yet it has the most negative adsorption energy among the BS heteroatom-doped combinations. It can be observed that the adsorption sites for H may not be limited to just the Cr atom, and the bond length between B and H is only 1.31 Å, suggesting that the doped B heteroatom facilitates the easier adsorption of H on the catalyst surface. Currently, the Volmer–Heyrovsky and Volmer–Tafel reaction mechanisms are popularly accepted to study the HER process based on the first-principles DFT calculations.⁵⁰ The appropriate reaction intermediate produces maximum catalytic activity. Excessively positive ΔG_*H values indicate that hydrogen is difficult to adsorb on the catalyst surface, while excessively negative ΔG_*H values reveal that the desorption of hydrogen from the substrate is arduous, both of which restrain HER performance. Therefore, the catalyst with a ΔG_*H value close to zero has been suggested to have superior catalytic properties.^51,52 All the ΔG_*H results of HER processing on TM@HM-GDY catalysts are presented in Table S4.†

It can be seen that there are eight ideal SACs whose ΔG_*H values are close to 0 eV (|ΔG_*H| < 0.05 eV), namely Cr@B2S3-GDY (0.0097 eV), Cr@B2S4-GDY (0.0437 eV), Cr@N2B1-GDY (−0.0046 eV), Mn@N1B2-GDY (0.0467 eV), Fe@N1S2-GDY (−0.0356 eV), Fe@N1S3-GDY (−0.0189 eV), Co@N1B4-GDY (0.0353 eV) and Ni@N2-GDY (−0.0209 eV), respectively. These values are all much lower than that of the ideal Pt-SACs for the HER (−0.07 eV).⁵³ In addition, we also calculated the HER catalytic activity of the initial GDY (ΔG_*H = 0.68 eV), evidencing that synergistic interaction of the doped TM atom with ambient heteroatoms improves the performance of the HER. In addition, for TM-loaded systems (Table S4†), it is observed that the doping of different heteroatoms at different sites had significant effects on HER efficiency. Specifically, for Cr systems, the immobilization of the N atom at the C1 site (Fig. 1a) exhibits a detrimental effect, manifested by an increased energy barrier for hydrogen adsorption (ΔG_*H = 0.81 eV). In contrast, loading at the C2 site facilitates the catalytic reaction, resulting in a lower energy barrier (ΔG_*H = 0.44 eV). Furthermore, doping both B and S positively impacts hydrogen adsorption. The doping of S exhibits a pronounced effect on hydrogen adsorption. S doping at the C1 position results in a reduction of ΔG_*H to 0.1074 eV; however, when anchored at the C2 site, it exhibits an excessively favorable effect, ultimately hindering the subsequent desorption process. The synergistic interaction of dual atoms notably enhances the activity of the hydrogen evolution reaction. For instance, B doping exerts a beneficial effect on the HER at the C1 site (hereafter referred to as B1), while the introduction of an additional heteroatom (S or N) significantly reduces the free energy, with both values being lower than that of the initial GDY catalyst. These configurations are denoted as B1S2 (ΔG_*H = −0.14 eV), B1S3 (ΔG_*H = 0.25 eV), and N2B1 (ΔG_*H = −0.0046 eV), respectively. However, the interaction in the N1S3 combination is excessively strong, resulting in an overly robust hydrogen adsorption characterized by ΔG_*H = −1.58 eV, and the excessive adsorption indicates that the desorption of *H from the catalyst surface is extremely challenging. For Mn-loaded systems, doping of a single heteroatom exerts a favorable influence on HER performance, and ΔG_*H indicates that the impact of S doping is the most pronounced. The co-introduction of B at C1 sites and S at C2 and C3 sites exhibits a positive impact on the activity of the HER (ΔG_*H = 0.21, 0.11 eV). However, the activity of the B2SX (B2S3, S4) type catalyst is weakened and ΔG_*H increased to 0.59 eV and 0.45 eV. When N and B were co-doped, N1B2 and N1B3 showed a favorable trend in HER activity, while other doping modes showed adverse effects. There is no catalyst with excellent HER activity in the co-doping mode of N and S. For Fe-loaded systems, all ∣ΔG_*H∣ values decreased significantly, which indicated that atom doping was beneficial to regulate the HER activity of Fe-based structures. Among them, the ΔG_*H of N1S2, N1S3 and N1S4 (0.04, 0.02 and 0.08 eV) configurations showed the desired HER activity. There is a positive effect on the activity of the HER for Co-based system doped with a heteroatom. A large proportion of ∣ΔG_*H∣ showed a decreasing trend; in particular, the ∣ΔG_*H∣ values of B1, N1B3, N1B4 and N1S3 are closest to ZERO. For Ni-loaded systems, modification of the coordination environment by single heteroatom doping has more positive effects than that by double heteroatom doping, because the partial co-doping mode greatly increased the ∣ΔG_*H∣. ∣ΔG_*H∣ of the single atom doping mode is generally decreased, while the ∣ΔG_*H∣ values of N1, N2, B2 and S1 are 0.14, 0.021, 0.20 and 0.11 eV, respectively.

Moreover, AIMD simulations are subsequently performed to verify the thermal stability of the screened structures at 500 K for 6 ps with a time step of 2 fs. As shown in Fig. S1,† the energies and temperatures of Cr@N2B1-GDY, Mn@N1B-GDY, Fe@N1S2-GDY, Co@N1B4-GDY and Ni@N2-GDY oscillate within small ranges and the structures remain stable without significant deformation. Therefore, these configurations are considered to be outstanding catalysts with both HER activity and stability. To verify the accuracy of the screening results, we reconsidered multiple possible loading TM sites, such as the C atom on the benzene ring and the substituted C atoms at both ends of ethyne (Fig. S2†). It can be seen from the calculated adsorption energy results that the configurations are most stable when TM atoms are adsorbed onto the surface cavity (Table S5†). Moreover, the adsorption of hydrogen on the doped heteroatoms in the five candidates was also explored, presented in Table S6.† In the Cr, Mn, and Co systems, the TM sites are more prone to the hydrogen evolution reaction, while the N heteroatom site is more favorable for H adsorption for the Ni configuration, though the difference is minimal. Considering the actual situation, solvation effects were used for further screening and validation, as shown in Table S7.† The calculated results show that the other catalysts still maintain excellent HER activity except Mn@N1B2-GDY which is affected by relatively large solvation, which means that our calculation and screening method is reasonable.

To validate the superior performance of the catalysts, we plotted the volcano curve as a function of the exchange current density (i₀) and ΔG_*H (shown in Fig. 2(a–e)). In principle, the best HER catalyst is located at the peak of the volcanic map,⁵⁴ that is, the closer the HER catalyst is to the top the better its HER activity. It is clear that the five screened SACs are located near the volcano peak with a pentagonal star in the diagrams; specifically, Cr@N2B1-GDY and Mn@N1B2 are situated exactly at the top, further confirming that the co-doping of TM and HM atoms is beneficial for the HER process.

Fig. 2 (a)–(e) The HER volcano curve relationship between ΔG_*H and theoretical exchange current i₀ for Cr-based, Mn-based, Fe-based, Co-based and Ni-based SACs, respectively; (f) free energy diagram of five TM@HM-GDY catalysts obtained by screening ΔG_*H close to 0.

3.3 Electronic properties

As mentioned earlier, the adsorption of hydrogen on the GDY surface is thermodynamically unfavorable, owing to the abundant sp²- and sp-hybridized carbon atoms and stable double and triple bonds, which makes it difficult for H atoms to adsorb and destabilize the structure. Therefore, additional metal-heteroatom combinations were introduced to activate the GDY surface to facilitate the HER process. To further elucidate the TM and HM co-doping modulation of HER activity, we performed partial density of states (PDOS) calculations on the pristine GDY and the single metal loaded GDY (abbreviated as SAC-GDY), which are presented in Fig. S3.† Meanwhile, the PDOSs of the heteroatom-doped GDY and the 5 screened doped SAC structures were also calculated, as detailed in Fig. 3. The band gap of pristine GDY (0.5 eV)⁵⁵ decreases with the intervention of heteroatoms (Fig. 3a, d, j, g, m). In addition, we can clearly see the strong interactions between the p-orbitals of the atoms in the substrate and the d-orbitals of the TM atoms (Fig. 3b, e, h, k, n), and the partially filled impurity states produced near the Fermi energy level, which facilitate catalytic reactions. In addition, PDOS also shows that N/B co-doping caused electron spin polarization of Cr, Mn and Co, implying that doping of heteroatoms at asymmetric sites on graphitylene may affect HER activity by regulating the spin states of electrons. To verify this conjecture, the PDOSs of SACs composed of several metals loaded with undoped atoms in GDY (Cr-GDY, Mn-GDY, Fe-GDY, Co-GDY and Ni-GDY) are calculated (Fig. S3 and S4†). Several metal centers are both in the high spin state before and after heteroatom doping. By analyzing the PDOS of metal d orbitals before and after doping, it is found that the PDOS of the d orbital shifts to the right away from the Fermi level. The d orbital PDOS components in different directions are used for further analysis. For Cr, Mn, and Co, the doping of HM intensifies the spin splitting of D-orbital electrons. For example, the energy level difference between the upward and downward spin components of d_z² in Cr atoms becomes larger. Similarly, a similar phenomenon can be observed in the PDOS of Co (d_z²) and Mn (d_xy, d_z², d_x²_–y²). Moreover, the spin asymmetry of these orbital components is further aggravated after HM doping compared with that before doping. Doping does not cause spin splitting of Fe and Ni atoms. However, the D-orbital components of the Ni (d_x²–y², d_z²) and Fe (d_z², d_x²–y²) atoms are observed to shift significantly to the right and away from the Fermi level. The components of d orbitals in various directions of several catalyst metals show that the doping of HM causes the components of e_g (d_z² and d_x²–y²) with higher energy to shift significantly to the right, which means that the electron filling of antibonding orbitals will increase when adsorption occurs. Table S7† also shows that the adsorption of H by TM@HM-GDY is significantly weakened. These results indicate that the hybridization of the doped p–d orbitals of HM reduces the adsorption strength of H through the high filling state of the anti-bonding orbital during adsorption and activates the reactive activity of the HER in the metal center.

Fig. 3 The partial density of states (PDOS) of heteroatom doping GDY, (a) N2B1-GDY, (d) N1B2-GDY, (g) N1S2-GDY, (j) N1B4-GDY and (m) N2-GDY; the PDOS of transition metal doped SACs, (b) Cr@N2B1-GDY, (e) Mn@N1B2-GDY, (h) Fe@N1S2-GDY, (k) Co@N1B4-GDY and (n) Ni@N2-GDY; top views of the three-dimensional charge difference plots for (c) Cr@N2B1-GDY, (f) Mn@N1B2-GDY, (i) Fe@N1S2-GDY, (l) Co@N1B4-GDY and (o) Ni@N2-GDY with an isovalue of 0.008 e ų.

Next, we calculated the charge transfer on 5 TM@HM-GDY catalysts, which is closely related to the active center, and the results are plotted in Fig. S5.† The electron-consuming and electron-accumulating zones are expressed as cyan and yellow isosurfaces, respectively. It is clear from the figure that electron redistribution occurs between the loaded metal atom and the heteroatoms and C atoms in the bonded substrate. Also, combined with Bader charge (Table S9†), we found that 0.52 ∼ 1.06e of TM atoms are transferred to the GDY to consolidate the structure and contribute to the adsorption of H.

After H adsorption, charge accumulation and depletion were clearly observed (Fig. 3c, f, i, l, o) between the active sites of TM atoms and the adsorbed H atom. In order to quantitatively estimate the charge transfer between H adsorption and the catalyst surface, we performed Bader charge analysis shown in Table S9.† For the Mn@N1B2 system, the maximum charge transfer (1.31e) to hydrogen guarantees the TM–H bond stabilization and the stable adsorption of *H intermediates. On the other hand, Co@N1B4 combination transports the least charge (1.05e) ensuring a weaker binding energy. In contrast, Cr@N2B1-GDY exhibits a moderate amount of charge transfer (1.28e) becoming the most superior HER catalyst (ΔG_*H = −0.0046 eV).

3.4 Construction of interpretable machine learning models

46 initial indicators are taken as the input variables including the atomic structure features, that is, atom radius (R_TM, R_C/N/B/S), distance between TM and HM (C/N/B/S-TM), distance between H and TM (dTM-H), TM and HA electronegativity (χ_TM, χ_C/N/B/S), magnetic moments of TM (magnetic), enthalpy of hydrogenation of TM (DH_(TM)_yH_x), formation energy of catalyst (E_f), D-band center of TM of catalysts and P-band center of all coordination atoms (D_band and P_band), the first ionization energy of TM and coordination atom (IE_C/N/B/S and IE_HA), band gap of catalysts (band_gap), relative atomic mass of TM and HA (M_TM and M_C/N/B/S), electron affinity of TM and HA (Eea_TM and Eea_C/N/B/S), Bard charge carried by TM (Q_bader_TM), number and type of coordination atoms (n_C/N/S/B), and number of valence electrons of TM and HA (V_TM and V_C/N/B/S). In order to avoid the impact of the initial feature input 46-dimension disaster and missing values, feature engineering is used for dimensionality reduction and to combine these features. For example, the product of C/N/B/S-TM, n_C/N/S/B, and Eea_TM, χ_TM, V_TM is Eea_corr, χ_corr, V_corr, sum of IE_HA (IE_HA), average bond length of N/B/S-TM (Bond_ave), total M_C/N/B/S (M_HA), sum of V_C/N/B/S(V_HA), and sum of R_C/N/B/S(R_HA). The principal component analysis technique is used for feature selection and the final 27 features are selected as the final inputs. ΔG_*H is the target predictive value. The data set is divided into 75% for training and the remaining 25% for testing. Six supervised ML regressors, Support Vector Machine (SVR), Random Forest (RF), eXtreme Gradient Boosting (XGBoost), Light Gradient Boosting Machine (LightGBM), k-nearest neighbor algorithm (kNN) and ElasticNet machine learning are trained as predictive models. Bayesian optimization algorithms and 10-fold cross validation are used to adjust hyperparameters. After hyperparameter optimization, the model is evaluated using the reserved test set and various evaluation indexes are calculated. After initial regression training, the test R² values of kNN, ElasticNet, SVR, RF, XGBoost and LightGBM models are 0.75, 0.51, 0.82, 0.85, 0.92 and 0.95, respectively (Fig. S6–S11†). As illustrated in Fig. 4a, RF shows the best predictability; however, RF achieving precise training can be attributed to limited data and specific training strategies. The model may contain redundant features that affect the predictive performance and interpretability and require further feature screening. Fig. 4b shows the SHAP value ranking of feature importance, with different color bars representing different types of parameters. The negligible proportion factors include V_HA, θd_TM, χ_H_ave, RM_HA, M_TM, Eea_TM and V_TM. The remaining 20 features were subsequently utilized as input for the retraining of RF. The results demonstrating the enhanced prediction accuracy are presented in Fig. 4c (R² = 0.96). It is evident from the figure that all test values are closely clustered around the fitted line, indicating a high degree of agreement between the predicted and actual values. SHAP summary visualization analysis is used to assess the importance of features in the RF model. The distribution of SHAP values in Fig. 4d shows that dTM-H, magnetic, and D_band are the three most important factors affecting the prediction. Meanwhile, the SHAP force plot shows (Fig. 4e) that there is a strong positive correlation between dTM-H, magnetic and ΔG_*H. This explains that the hydrogen affinity of metal atoms and the electronic structure regulated by the coordination environment play a decisive role in the activity of the HER. dTM-H symbolizes the affinity to hydrogen of TM, and the distribution of d electrons of M is closely related to the adsorption of H. The spin polarization of electrons affects the activity of the HER. The interpretability of the feature selection of RF models is reflected in physicochemical rules. Further, the two-dimensional partial dependencies of TDM-H with magnetic (Fig. 4f) and D_band (Fig. 4g), respectively, are rendered to represent the impact of their interaction on the RF model. It is worth noting that the values at medium sizes cause the strongest interactions, for example, the standard scaled dTM-H (0.48) and magnetic (0.6) give a dependency of 0.68, and the standard scaled dTM-H (0.46) and D_band (0.49) give a dependency of 0.67. This may be related to the “Sabatier principle” in catalytic reactions, and the quantitative measure of three features may enable the construction of volcanic descriptors that describe the activity of this series of catalysts. Finally, the feature selection of the RF model is consistent with the physicochemical principles.

Fig. 4 (a) R² of four machine learning models obtained through initial training based on XGBoost, RF, LGB and SVR; (b) feature importance of an initially trained RF machine learning model is demonstrated; (c) distribution of true and predicted values of the RF machine learning model after feature selection; (d) feature importance arrangement of the RF model based on SHAP analysis; (e) feature force plot analysis of the RF model based on the SHAP value; two-dimensional part dependence between (f) dTM-H and magnetic; and (g) dTM-H and D_band (all values are scaled to 0–1).

4 Conclusion

In summary, we investigated the HER performance of TM@NM-GDY (TM = Cr, Mn, Fe, Si, Ni and HM = N, B, S) as catalysts based on first-principles calculations. Among the 120 configurations screened computationally, 5 systems (Cr@N2B1-GDY, Mn@N1B2-GDY, Fe@N1S2-GDY, Co@N1B4-GDY and Ni@N2-GDY) with near-zero Gibbs free energy of H adsorption were certified, specifically, Cr@N2B1-GDY performed optimally with ΔG_*H = −0.0046 eV. The electron structure analysis shows that the asymmetric heteroatom-doping of graphdiyne can affect the adsorption of H by regulating the electron spin state of the TM or the D-band center. The results of AIMD calculations demonstrate that the five screened prototypes have better thermal stability. Interpretable machine learning based on SHAP analysis is used to construct these catalyst activity prediction models. Through feature engineering and hyper-parametric optimization it was determined that RF exhibited the best prediction accuracy, with an R² of 96%. SHAP analysis shows that the metal's hydrophilicity, electron spin state and D-band electron structure are the decisive factors affecting the prediction of the ML model, which is consistent with the electron structure analysis and fully explains the feature selection of the RF model. In this study, interpretive machine learning and density functional theory are used to reveal that asymmetric doping of heteroatoms on graphitic acetylene can regulate the electronic structure and spin state of TM and affect the activity of the HER. This provides a new way to establish the system theory of a single-atom HER catalyst.

Data availability

The data supporting this article have been included as part of the ESI.†

Author contributions

Conceptualization, methodology, formal analysis and investigation, writing – original draft preparation: Ying Zhao and Shuai–Shuai Gao; writing – review and editing: Ying Zhao; funding acquisition: Li-Shuang Ma; resources: Peng-Hui Ren; supervision: Ying Zhao and Xue-Bo Chen.

Conflicts of interest

The authors declare that there is no conflict of interest.

Acknowledgements

This work is funded by the Shandong Provincial Key Research and Development Program (Competitive Innovation Platform) (No. 2024CXPT036) and Distinguished Expert Program of Taishan scholars in Shandong Province (No. tstp20240844).

References

1. D. Er, H. Ye, N. C. Frey, H. Kumar, J. Lou and V. B. Shenoy, Prediction of Enhanced Catalytic Activity for Hydrogen Evolution Reaction in Janus Transition Metal Dichalcogenides, Nano Lett., 2018, 18(6), 3943–3949.
2. J. N. Tiwari, N. K. Dang, S. Sultan, P. Thangavel, H. Y. Jeong and K. S. Kim, Multi-heteroatom-doped carbon from waste-yeast biomass for sustained water splitting, Nat. Sustain., 2020, 3(7), 556–563.
3. C. Xiao, R. Sa, Z. Cui, S. Gao, W. Du, X. Sun, X. Zhang, Q. Li and Z. Ma, Enhancing the hydrogen evolution reaction by non-precious transition metal (Non-metal) atom doping in defective MoSi₂N₄ monolayer, Appl. Surf. Sci., 2021, 563, 150388.
4. Y. Zang, Q. Wu, W. Du, Y. Dai, B. Huang and Y. Ma, Activating electrocatalytic hydrogen evolution performance of two-dimensional MSi₂N₄ (M = Mo, W): a theoretical prediction, Phys. Rev. Mater., 2021, 5(4), 045801.
5. V. A. Panchenko, Y. V. Daus, A. A. Kovalev, I. V. Yudaev and Y. V. Litti, Prospects for the production of green hydrogen: review of countries with high potential, Int. J. Hydrogen Energy, 2023, 48(12), 4551–4571.
6. X. Xu, Q. Zhou and D. Yu, The future of hydrogen energy: bio-hydrogen production technology, Int. J. Hydrogen Energy, 2022, 47(79), 33677–33698.
7. X. Xiao, L. Yang, W. Sun, Y. Chen, H. Yu, K. Li, B. Jia, L. Zhang and T. Ma, Electrocatalytic Water Splitting: From Harsh and Mild Conditions to Natural Seawater, Small, 2022, 18(11), 2105830.
8. L. Zhao, H. Yuan, D. Sun, J. Jia, J. Yu, X. Zhang, X. Liu, H. Liu and W. Zhou, Active facet regulation of highly aligned molybdenum carbide porous octahedrons via crystal engineering for hydrogen evolution reaction, Nano Energy, 2020, 77, 105056.
9. J. N. Tiwari, S. Sultan, C. W. Myung, T. Yoon, N. Li, M. Ha, A. M. Harzandi, H. J. Park, D. Y. Kim, S. S. Chandrasekaran, W. G. Lee, V. Vij, H. Kang, T. J. Shin, H. S. Shin, G. Lee, Z. Lee and K. S. Kim, Multicomponent electrocatalyst with ultralow Pt loading and high hydrogen evolution activity, Nat. Energy, 2018, 3(9), 773–782.
10. L. Bu, N. Zhang, S. Guo, X. Zhang, J. Li, J. Yao, T. Wu, G. Lu, J.-Y. Ma, D. Su and X. Huang, Biaxially strained PtPb/Pt core/shell nanoplate boosts oxygen reduction catalysis, Science, 2016, 354(6318), 1410–1414.
11. Y. Shiraishi, Y. Kofuji, S. Kanazawa, H. Sakamoto, S. Ichikawa, S. Tanaka and T. Hirai, Platinum nanoparticles strongly associated with graphitic carbon nitride as efficient co-catalysts for photocatalytic hydrogen evolution under visible light, Chem. Commun., 2014, 50(96), 15255–15258.
12. Y. Zhang, H. Lei, D. Duan, E. Villota, C. Liu and R. Ruan, New Insight into the Mechanism of the Hydrogen Evolution Reaction on MoP(001) from First Principles, ACS Appl. Mater. Interfaces, 2018, 10(24), 20429–20439.
13. N. K. Dang, M. Umer, P. Thangavel, S. Sultan, J. N. Tiwari, J. H. Lee, M. G. Kim and K. S. Kim, Surface enrichment of iridium on IrCo alloys for boosting hydrogen production, J. Mater. Chem. A, 2021, 9(31), 16898–16905.
14. T. Wang, D. Gao, W. Xiao, P. Xi, D. Xue and J. Wang, Transition-metal-doped NiSe₂ nanosheets towards efficient hydrogen evolution reactions, Nano Res., 2018, 11(11), 6051–6061.
15. J. N. Tiwari, A. M. Harzandi, M. Ha, S. Sultan, C. W. Myung, H. J. Park, D. Y. Kim, P. Thangavel, A. N. Singh, P. Sharma, S. S. Chandrasekaran, F. Salehnia, J.-W. Jang, H. S. Shin, Z. Lee and K. S. Kim, High-Performance Hydrogen Evolution by Ru Single Atoms and Nitrided-Ru Nanoparticles Implanted on N-Doped Graphitic Sheet, Adv. Energy Mater., 2019, 9(26), 1900931.
16. A. Gupta, T. Sakthivel and S. Seal, Recent development in 2D materials beyond graphene, Prog. Mater. Sci., 2015, 73, 44–126.
17. V. Thirumal, R. Yuvakkumar, P. S. Kumar, G. Ravi, A. Arun, R. K. Guduru and D. Velauthapillai, Heterostructured two dimensional materials of MXene and graphene by hydrothermal method for efficient hydrogen production and HER activities, Int. J. Hydrogen Energy, 2023, 48(17), 6478–6487.
18. Y. Cheng, L. Wang, Y. Song and Y. Zhang, Deep insights into the exfoliation properties of MAX to MXenes and the hydrogen evolution performances of 2D MXenes, J. Mater. Chem. A, 2019, 7(26), 15862–15870.
19. S. Chandrasekaran, R. Sukanya, E. Arumugam, S. M. Chen and S. Vignesh, Effective sonochemical synthesis of titanium nitride nanoflakes decorated graphitic carbon nitride as an efficient bifunctional electrocatalyst for HER and OER performance, Colloids Surf., A, 2023, 665, 131190.
20. M. Gao, D. Liu, H. Yang, H. Huang, Q. Luo, Y. Huang, X.-F. Yu and P. K. Chu, Modification of Layered Graphitic Carbon Nitride by Nitrogen Plasma for Improved Electrocatalytic Hydrogen Evolution, Nanomaterials, 2019, 568.
21. J. Lee, S. Kang, K. Yim, K. Y. Kim, H. W. Jang, Y. Kang and S. Han, Hydrogen Evolution Reaction at Anion Vacancy of Two-Dimensional Transition-Metal Dichalcogenides: Ab Initio Computational Screening, J. Phys. Chem. Lett., 2018, 9(8), 2049–2055.
22. S. M. Tan, C. K. Chua, D. Sedmidubský, Z. Sofer and M. Pumera, Electrochemistry of layered GaSe and GeS: applications to ORR, OER and HER, Phys. Chem. Chem. Phys., 2016, 18(3), 1699–1711.
23. J. Kim, S. S. Baik, S. H. Ryu, Y. Sohn, S. Park, B.-G. Park, J. Denlinger, Y. Yi, H. J. Choi and K. S. Kim, Observation of tunable band gap and anisotropic Dirac semimetal state in black phosphorus, Science, 2015, 349(6249), 723–726.
24. J. Guan, T. Pal, K. Kamiya, N. Fukui, H. Maeda, T. Sato, H. Suzuki, O. Tomita, H. Nishihara, R. Abe and R. Sakamoto, Two-Dimensional Metal–Organic Framework Acts as a Hydrogen Evolution Co catalyst for Overall Photocatalytic Water Splitting, ACS Catal., 2022, 12(7), 3881–3889.
25. Y. Zhou, L. Sheng, Q. Luo, W. Zhang and J. Yang, Improving the Activity of Electrocatalysts toward the Hydrogen Evolution Reaction, the Oxygen Evolution Reaction, and the Oxygen Reduction Reaction via Modification of Metal and Ligand of Conductive Two-Dimensional Metal–Organic Frameworks, J. Phys. Chem. Lett., 2021, 12(48), 11652–11658.
26. Y. Zhao, J. Wan, H. Yao, L. Zhang, K. Lin, L. Wang, N. Yang, D. Liu, L. Song, J. Zhu, L. Gu, L. Liu, H. Zhao, Y. Li and D. Wang, Few-layer graphdiyne doped with sp-hybridized nitrogen atoms at acetylenic sites for oxygen reduction electrocatalysis, Nat. Chem., 2018, 10(9), 924–931.
27. Z. Jia, Y. Li, Z. Zuo, H. Liu, C. Huang and Y. Li, Synthesis and Properties of 2D Carbon—Graphdiyne, Acc. Chem. Res., 2017, 50(10), 2470–2478.
28. T. Liu, Q. Wang, G. Wang and X. Bao, Electrochemical CO₂ reduction on graphdiyne: a DFT study, Green Chem., 2021, 23(3), 1212–1219.
29. X. Gao, H. Liu, D. Wang and J. Zhang, Graphdiyne: synthesis, properties, and applications, Chem. Soc. Rev., 2019, 48(3), 908–936.
30. Z. Zuo, D. Wang, J. Zhang, F. Lu and Y. Li, Synthesis and Applications of Graphdiyne-Based Metal-Free Catalysts, Adv. Mater., 2019, 31(13), 1803762.
31. Y. Xue, B. Huang, Y. Yi, Y. Guo, Z. Zuo, Y. Li, Z. Jia, H. Liu and Y. Li, Anchoring zero valence single atoms of nickel and iron on graphdiyne for hydrogen evolution, Nat. Commun., 2018, 9(1), 1460.
32. H. Yu, Y. Xue, B. Huang, L. Hui, C. Zhang, Y. Fang, Y. Liu, Y. Zhao, Y. Li, H. Liu and Y. Li, Ultrathin Nanosheet of Graphdiyne-Supported Palladium Atom Catalyst for Efficient Hydrogen Production, iScience, 2019, 11, 31–41.
33. T. Liu, X. Hao, J. Liu, P. Zhang, J. Chang, H. Shang and X. Liu, Graphdiyne and Nitrogen-Doped Graphdiyne Nanotubes as Highly Efficient Electrocatalysts for Oxygen Reduction Reaction, Int. J. Mol. Sci., 2023, 16813.
34. L. Jasin Arachchige, Y. Xu, Z. Dai, X. L. Zhang, F. Wang and C. Sun, Double transition metal atoms anchored on graphdiyne as promising catalyst for electrochemical nitrogen reduction reaction, J. Mater. Sci. Technol., 2021, 77, 244–251.
35. Y. Shang, X. Duan, S. Wang, Q. Yue, B. Gao and X. Xu, Carbon-based single atom catalyst: synthesis, characterization, DFT calculations, Chin. Chem. Lett., 2022, 33(2), 663–673.
36. M. Younas, M. Yar, H. AlMohamadi, T. Mahmood, K. Ayub, A. L. Khan, M. Yasin and M. A. Gilani, A rational design of covalent organic framework supported single atom catalysts for hydrogen evolution reaction: a DFT study, Int. J. Hydrogen Energy, 2024, 51, 758–773.
37. Y. Lei, Y. Wang, Y. Liu, C. Song, Q. Li, D. Wang and Y. Li, Designing Atomic Active Centers for Hydrogen Evolution Electrocatalysts, Angew. Chem., Int. Ed., 2020, 59(47), 20794–20812.
38. Z. Li, X. Li, Y. Gu, X. Hu, L. Wang and P. Li, Theoretical screening of transition metal atoms doped two-dimensional VSe₂ monolayers as single-atom catalysts for HER, OER and ORR, Appl. Surf. Sci., 2024, 646, 158959.
39. X.-P. Yin, H.-J. Wang, S.-F. Tang, X.-L. Lu, M. Shu, R. Si and T.-B. Lu, Engineering the Coordination Environment of Single-Atom Platinum Anchored on Graphdiyne for Optimizing Electrocatalytic Hydrogen Evolution, Angew. Chem., Int. Ed., 2018, 57(30), 9382–9386.
40. X. Guo, J. Gu, X. Hu, S. Zhang, Z. Chen and S. Huang, Coordination tailoring towards efficient single-atom catalysts for N₂ fixation: a case study of iron–nitrogen–carbon (Fe@N–C) systems, Catal. Today, 2020, 350, 91–99.
41. M. Umer, S. Umer, M. Zafari, M. Ha, R. Anand, A. Hajibabaei, A. Abbas, G. Lee and K. S. Kim, Machine learning assisted high-throughput screening of transition metal single atom based superb hydrogen evolution electrocatalysts, J. Mater. Chem. A, 2022, 10(12), 6679–6689.
42. X. Li, R. Chiong, Z. Hu and A. J. Page, Low-Cost Pt Alloys for Heterogeneous Catalysis Predicted by Density Functional Theory and Active Learning, J. Phys. Chem. Lett., 2021, 12(30), 7305–7311.
43. X. Li, H. Li, Z. Zhang, J. Q. Shi, Y. Jiao and S.-Z. Qiao, Active-learning accelerated computational screening of A₂B@NG catalysts for CO₂ electrochemical reduction, Nano Energy, 2023, 115, 108695.
44. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47(1), 558–561.
45. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77(18), 3865–3868.
46. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50(24), 17953–17979.
47. Á. Valdés, Z. W. Qu, G. J. Kroes, J. Rossmeisl and J. K. Nørskov, Oxidation and Photo-Oxidation of Water on TiO₂ Surface, J. Phys. Chem. C, 2008, 112(26), 9872–9879.
48. J. K. Nørskov, T. Bligaard, A. Logadottir, J. R. Kitchin, J. G. Chen, S. Pandelov and U. Stimming, Trends in the Exchange Current for Hydrogen Evolution, J. Electrochem. Soc., 2005, 152(3), J23.
49. G. Gao, A. P. O'Mullane and A. Du, 2D MXenes: A New Family of Promising Catalysts for the Hydrogen Evolution Reaction, ACS Catal., 2017, 7(1), 494–500.
50. M. R. Gennero de Chialvo and A. C. Chialvo, Kinetics of hydrogen evolution reaction with Frumkin adsorption: re-examination of the Volmer–Heyrovsky and Volmer–Tafel routes, Electrochim. Acta, 1998, 44(5), 841–851.
51. F. Abdelghafar, X. Xu, S. P. Jiang and Z. Shao, Designing single-atom catalysts toward improved alkaline hydrogen evolution reaction, Mater. Rep.: Energy, 2022, 2(3), 100144.
52. T. Liang, J. Xie, A. Wang, D. Ma, Z. Mao, J. Wang and H. Li, Prediction of HER electrocatalyst with enhanced performance based on atoms-doped black phosphorene: a first-principles study, Appl. Surf. Sci., 2022, 604, 154508.
53. H. Xu, D. Cheng, D. Cao and X. C. Zeng, Retracted article: a universal principle for a rational design of single-atom electrocatalysts, Nat. Catal., 2018, 1(5), 339–348.
54. T. He, G. Gao, L. Kou, G. Will and A. Du, Endohedral metallofullerenes (M@C60) as efficient catalysts for highly active hydrogen evolution reaction, J. Catal., 2017, 354, 231–235.
55. G. Luo, Q. Zheng, W.-N. Mei, J. Lu and S. Nagase, Structural, Electronic, and Optical Properties of Bulk Graphdiyne, J. Phys. Chem. C, 2013, 117(25), 13072–13079.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ta08095e

‡ These authors contributed equally to this work and should be considered co-first authors.

---