MBenes-supported single-atom catalysts for oxygen reduction and oxygen evolution reactions by first-principles study and machine learning

: Oxygen reduction reaction (ORR) and oxygen evolution reaction (OER) are key catalytic processes in various renewable energy conversion and energy storage technologies. Herein, we systematically investigated the ORR and OER catalytic activity of the single-atom catalysts (SACs) composed of 4d/5d period transition metal (TM) atoms embedded on MBene substrates (TM-M 2 B 2 O 2 , M = Ti, Mo, and W). We found that TM dominates the catalytic activity compared to the MBene substrates. The SACs embedded with Rh, Pd, Au, and Ir exhibit excellent ORR or OER catalytic activity. Specifically, Rh-Mo 2 B 2 O 2 and Rh-W 2 B 2 O 2 are promising bifunctional catalysts with ultra-low ORR/OER overpotentials of 0.39/0.21 V and 0.19/0.32 V, respectively, lower than that of Pt/RuO 2 (0.45/0.42 V). Importantly, through machine learning, the models containing 10 element features of SACs were developed to quickly and accurately identify the superior ORR and OER electro-catalysts. Our findings provide several promising SACs for ORR and OER, and offer effective models for catalyst

Finally, combining data from our previous work [1,19], an ML dataset was constructed, which contained 117 SAC systems and 24 material elemental features that can be quickly obtained without complex computations.Subsequently, the feature engineering and decision tree method were used to filter out the top 10 features with lower correlation coefficients.Based on these features, two models were established to rapidly and accurately predict ORR and OER catalytic activity, respectively.Furthermore, based on symbolic regression and genetic algorithms, two symbolic expressions were developed to describe the ORR and OER catalytic activity of SACs, respectively.Our work provides several ORR and OER electrocatalysts with superior catalytic activity and provides simple and effective methods for screening and designing new catalysts.

COMPUTATIONAL METHODS
The spin-polarized DFT calculations were employed in Vienna ab initio simulation package (VASP) code with the projector augmented wave (PAW) method [35].The exchange-correlation interaction was described by the Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA) [36,37].The kinetic energy cutoff was set to 500 eV.A 3 × 3 × 1 Monkhorst-Pack k-point grid and a denser 7 × 7 × 1 Monkhorst-Pack k-point grid was chosen for structural optimization and electronic structure calculations, respectively.The interlayer interactions were eliminated efficiently by setting a vacuum of 25 Å along the z-axis.The convergence thresholds for energy and force were set to 10 −5 eV/atom and 0.01 eV/Å, respectively.Given that the studied systems involve surface adsorption and desorption processes, the impact of the van der Waals interactions was treated by using the empirical correction DFT-D3 approach [38].Microkinetic simulations were carried out using the Catalysis Microkinetic Analysis Package (CATMAP) software package [39].The thermal stability of SAC structures under aqueous solution was evaluated by performing ab initio molecular dynamics (AIMD) simulations at 300 K for 10 ps.It is worth noting that the selection of the above parameters is roughly consistent with our previous work to facilitate subsequent ML [1,19].
The ORR process of MBenes-supported SACs is typically performed by the widely reported four-electron dissociation mechanism, which includes four elementary steps and can be described by the following equations [40]: where *, H + , and e − represent the active sites of SAC surface, proton, and transferred electron, respectively.g and l denote the gas phase and liquid phase, respectively.OER can be regarded as the reverse reaction process of ORR, which is composed of the following four corresponding elementary steps [41]: In order to evaluate the ORR catalytic activity, it is necessary to calculate the Gibbs free energy change (∆G) of each elementary step.According to the computational hydrogen electrode (CHE) model proposed by Nørskov et al. [9], the general formula of ∆G can be described as the following equation for each elementary step.
where ∆E is the total energy change before and after the elementary reaction, which can be directly obtained from DFT calculations.∆ZPE and T∆S represent the changes in zero-point energy and entropic energy, represents the Gibbs free energy correction at the specified pH, where k B is the Boltzmann constant.The last term is the correction of electrical potential to Gibbs free energy, denoted as ∆G U = −eU, where U is the applied electrode potential with respect to the standard hydrogen electrode and was set to 0 V, and e is the elementary charge.Herein, all reactions were assumed to occur under standard conditions (T = 298.15K, P = 1 bar and pH = 0).Correspondingly, the free energy changes of the four elementary steps of ORR can be specifically expressed as Equations ( 10)- (13).
For ORR, the elementary step with the maximum energy change would be determined as the potential determining step.To evaluate the catalytic activity of ORR, the well-accepted theoretical overpotential η ORR was defined as the equation below: According to the definition of overpotential, for both ORR and OER, the lower the overpotential, the higher the catalytic activity.

Structural stability
The long-lasting stability of the catalysts' structure is indispensable for practical applications.In order to access the feasibility of SAC formation, we first investigated all possible embedding sites of TM in 54 SACs composed of 18 kinds of 4d/5d period TM atoms anchored on three kinds of MBene substrates (donated as TM-M 2 B 2 O 2 , TM = 4d/5d period atoms, M = Ti, Mo and W), and selected the structure with optimal embedding sites for subsequent research, as illustrated in Figure 1.Then, the binding strength of single atoms anchored on different substrates was calculated, which can be quantified in the form of binding energy, defined as Natl Sci Open, 2024, Vol.3, 20230043 where E TM-MBene , E MBene , and E TM are the total energies of the TM atom anchored SAC systems, the free MBene substrates, and isolated TM atoms, respectively.The more negative binding energy E b , the more favorable the load of TM atoms.The binding energies of all SAC systems are negative, indicating that all TM atoms can be successfully anchored (Table S1).
As is well known, maintaining uniform and isolated distribution of single atom on substrates is a prerequisite for the advantages of SACs.Thus, to evaluate the aggregation possibility of TM atoms, the metal cluster energy (E cluster = E b − E coh ) was calculated, where E coh is the metal cohesive energy.Herein, E coh can be calculated by the energy difference between a metal atom in the most stable bulk phase and in the vacuum, expressed as E coh = E TM-bulk /n − E TM-single , where n is the number of atoms in the most stable bulk phase.A value of E cluster less than 0 eV is usually considered a powerful criterion for determining the thermodynamic stable of SAC systems.In addition, the dissolution potentials (U diss = U diss-bulk − E cluster / (eNe)) of TM atoms were further calculated to explore the stability of the studied SAC systems in electrochemical environments.U diss- bulk represents the standard dissolution potential of TM bulk that can be directly obtained from the handbook of chemistry and physics, and Ne is the number of transferred electrons [33,42].Theoretically, a positive value of U diss indicates that the SACs are electrochemically stable.The E cluster and U diss of all SACs are depicted in Figure 2, where the SACs located in the gray background area are thermodynamic and electrochemical stable (U diss > 0 V and E cluster < 0 eV).In addition, we found that the substrate has a significant impact on the thermodynamic stability and electrochemical stability of SACs, with Mo-based being the optimal substrate materials, which may benefit from the more charge transfer between the Mo-based substrate and metal atoms (Figure S1).Finally, in order to verify the stability of SACs under aqueous solution, the stability of the Rh-Mo 2 B 2 O 2 under aqueous solution composed of 32 water molecules was further performed through AIMD simulations at 300 K.As shown in Figure S2, the AIMD simulations confirm that Rh-Mo 2 B 2 O 2 is thermodynamically stable under aqueous solution.

ORR/OER catalytic activity
We systematically evaluated the ORR and OER catalytic activities of all the 54 SAC systems in this section.Generally, ORR occurring on the surface of MBene-supported SACs tends to undergo a four-electron association mechanism rather than the dissociation mechanism.This is because that the whole ORR process only involves a single active site, and the fracture of O-O bond needs to overcome a large energy barrier [1,19,43].Therefore, we only discussed the association mechanism in this work, which is consistent with our previous work [19].Firstly, the free energies of all adsorbed oxygen intermediates (O*, OH*, and OOH*) in the ORR process were calculated according to Equation ( 9) and summarized in Tables S2-S4.Furthermore, based on Equations ( 14) and ( 15), the overpotentials of all studied SAC systems for ORR and OER were obtained and summarized in Figure 3C and Tables S2 and Pt-W 2 B 2 O 2 exhibit superior ORR catalytic activity with overpotentials of 0.41, 0.39, 0.42, 0.19, and 0.30 V, respectively, comparable to or even lower than that of Pt (0.45 V) [9].Moreover, Rh-

ORR/OER activity trend
The adsorption strength between TM in SACs and the three intermediates (O*, OH*, and OOH*) generally presents a scaling relationship, which has been reported in many related studies [33,44,45].Then, we organized the adsorption free energies of three intermediates (ΔG O* , ΔG OH* , and ΔG OOH* ) in all SACs and explored their relationships.As shown in Figure 4 Based on the scaling relationship, the volcanic curves of ΔG OH* versus −η ORR , as well as ΔG O* − ΔG OH* versus −η OER were plotted in Figure 5A and B, respectively.As shown in Figure 5A, the PDS of the ORR process in all studied SACs is mostly concentrated in the elementary step of OH* removal to form H 2 O (OH* + H + + e − → * + H 2 O), which is reflected in the highly linear fitting part of the left branch in the volcanic curve, corresponding to the strong adsorption of OH* intermediate.On the contrary, the right branch part exhibits weak OH* binding strength, corresponding to the formation of OOH* as the PDS in these SACs (* + O 2 + H + + e − → OOH*).According to Sabatier principle [14], both too strong and too weak adsorption strength between the catalysts and the adsorbent are detrimental to catalytic performance.These catalysts located at the top of the volcano plot demonstrate moderate adsorption strength with the adsorbent and  10)-( 13), the four free-energy changes can be further expressed as the following Equations ( 17)-( 20 Subsequently, the contour map of ORR catalytic activity regarding ΔG OH* and ΔG O* − ΔG OH* is plotted in Figure 5C.Obviously, this contour map is divided into four regions by dashed lines, representing the four elementary reactions of the ORR process.If a catalyst falls in a certain region, it means that the corresponding elementary reaction in that region is the PDS of this catalyst.For early TMs of 4d and 5d periods (such as Y, Zr, Nb, Mo, Tc, Hf, Ta, W, and Re), MBenes-supported SACs exhibit strong OH* adsorption (ΔG OH* < −0.5 eV), making the PDS solely determined by ΔG OH* , that is, the desorption process of water in the fourth elementary step.For the later TM (such as Ru, Rh, Pd, Ag, Cd, Os, Ir, Pt, and Au), the further Finally, the microkinetic simulations were performed to evaluate the respective turnover frequency (TOF, molecules/sites) of H 2 O generation for the selected catalysts at 0.9 V based on possible ORR pathways (the specific simulation details are provided in the Supplementary information), as shown in Figure 6A.Due to the strong scaling relationships between ΔG OH* and ΔG OOH* , the SACs are dispersed at both sides of the fitting

OER/ORR activity origin
To reveal the essential mechanism of high catalytic activity of the studied systems, the electronic structures of all SACs have been calculated.Since the single transition atoms acting as active sites in all SACs are directly involved in OER and ORR processes, the electronic structure characteristics of the single atoms were given special attention.Previous studies have proved that the d-band center (ε d ) is a reliable descriptor to characterize the catalytic activity of catalysts [5,[47][48][49][50], typically defined as  In this work, the derivative of the Fermi-Dirac function, f E E ( ) T F , was chosen as the weight function, which has been validated in many previous studies (Figure 7B) [51,52].Then, the improved d-band center (ε id ) is defined as where E F is the energy of the Fermi level and f E E ( ) T F is the Fermi-Dirac function, expressed as where k is the Boltzmann's constant, T represents the temperature, and kT represents the nominal electron temperature.Obviously, the spreading of the weight function can be adjusted through the kT parameter.Through systematic testing of kT values ranging from 0.1 to 1 eV, it was confirmed that the best fitting effect was achieved when the kT was set to 0.5 eV.As shown in Figure 7C, the relationships between the adsorption free energies of intermediate adsorbents (O*, OH*, and OOH*) and ε id were established.Compared to the adsorption free energies versus ε d (Figure 7A), the adsorption free energies versus ε id exhibit a higher linear correlation.According to the d-band center theory [53], the interaction between the sp orbitals of adsorbed molecules and the d orbitals of metals would cause hybridization of orbitals and form binding and antibonding states.With the decrease of the d-band center, the anti-bonding energy band generated by coupling Natl Sci Open, 2024, Vol.3, 20230043 decreases accordingly and is then filled with electrons, reducing the bond strength and further leading to the decrease of the adsorption free energy.This also indicates that the results of this work follow the d-band center criterion, that is, the d-band center is negatively correlated with the adsorption free energy.Nevertheless, the fact is that the improved d-band center still cannot effectively describe the catalytic activity of the studied SACs.
To provide a clearer explanation of the adsorption behavior mechanism of adsorbents, the bonding strength between OH* intermediate and metal atoms for all SACs was first investigated by the bond length and charge transfer, as shown in Figure S4.It can be seen that the bond length and charge transfer from the metal atoms to the substrates cannot accurately describe the ΔG OH* due to the different radii and valence electron structures for different metal atoms.Then, the crystal orbital Hamiltonian population (COHP) of OH* intermediate for all SACs was further calculated.Furthermore, by integrating COHP, the parameter ICOHP can be obtained to quantify the bond strength between OH* and metal atoms.A more intuitive relationship between ICOHP and ΔG OH* is depicted in Figure 8A.It can be observed that there is a linear correlation between ICOHP and ΔG OH* .The linear relationship between ICOHP and adsorption free energy has also been reflected in our previous work (Figure S5) [19].Besides, ICOHP presents a similar pattern as the increase of the atomic number of TM atoms for a given period, that is, it decreases first and then increases, except for Pt-M 2 B 2 O 2 (M = Ti, Mo, and W), as shown in Figure 8B.The more negative ICOHP, the smaller ΔG OH* , and the stronger the bond strength between OH* and metal atoms.The SACs embedded with Rh, Ir, Pd, and Au atoms exhibit suitable ICOHP values, thus displaying superior ORR or OER catalytic activity, which also explains the reason why single metal atoms are the main factor affecting catalytic activity.

Machine learning methods
Based on the above findings, it is known that the catalytic activity of SACs is closely related to their adsorption surface structures and adsorbed TM atoms.In order to obtain a more explicit and intuitive relationship expression, the models containing 10 intrinsic elements features of SACs were established, and the quantitative expressions of ORR and OER catalytic activity were further obtained through symbolic regression.In this work, the construction of the models and the symbolic regression expressions for the OER and ORR performance mainly encompasses the following three core aspects:

Dataset for machine learning
The dataset used for ML consists of a total of 117 initial data, sourced from 54 SAC systems in this work, 27 SAC systems of MBenes-supported TM atoms (TM-M 2 B 2 O 2 , TM = Ti, V, Cr, Mn, Fe, Co, Ni, Cu, and Zn, M = Ti, Mo, and W) [19], 24 SAC systems of Ti 2 CO 2 -supported TM atoms (TM = Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zr, Nb, Mo, Tc, Ru, Rh, Pd, Ag, Hf, Ta, W, Re, Os, Ir, Pt, and Au) and 12 SAC systems of Ti 3 C 2 O 2 -supported TM atoms (TM = Mn, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ag, Ir, Pt, and Au) [1].Furthermore, based on the crossvalidation method, the aforementioned dataset was divided into a training set containing 75% data and a validation set containing 25% data.Within the training set, 70% data were allocated for training, and the remaining 30% was used for the testing dataset.

Descriptors
Material descriptors, often termed as "features" in the ML community, are typically used to transform material composition or structural information into ML input datasets.In this study, 37 features (descriptors) were assigned for each given data-point, which can be directly extracted from element characteristics.These features include both numerical and categorical variables (the configuration of extra-nuclear electrons).Through feature engineering, 13 features with similar maximum and minimum values are cleared, and the remaining 24 features are listed in Table S5.Their correlation coefficients are showcased in Figure 9A.A notable number of these features demonstrate a high correlation coefficient, implying that these features are not independent variables (Figure S6).Therefore, the features with an absolute value of correlation coefficients exceeding 0.85 have been removed, and the subsequent 10 independent features are represented in   S6.

Machine learning model and symbolic regression
With the advancements in computational capability and GPU acceleration, we have separately established two ML models for ORR and OER performances using the Random Forest algorithm.These models consist of 100 decision tree sub-models, and use the mean-square error (MSE) as the accuracy measurement standard of the model.The minimum number of samples required for node splitting in the decision tree is 2, and the minimum number of samples for leaf nodes is 1.A maximum of one feature is considered during each decision tree branching.Overfitting is prevented by setting the maximum search depth of the decision tree to 10 and incorporating cross-validation.As a result, we have constructed two models for predicting ORR and OER performances, with accuracies of 92.17% and 91.72%, respectively, as depicted in Figure 10A and B. These models do not rely on any computational parameters and can rapidly and accurately predict the ORR and OER performances of materials based on ten element-related features.
Furthermore, to quantitatively describe the mathematical relationship between material features and target performances, we have integrated the Random Forest method with genetic algorithms to construct a symbolic regression approach, as detailed in Github link (https://github.com/ALKEMIE-Lab/sraes).In symbolic regression, the cross-validation method is used to ensure that the data does not overfit.Initially, we obtained the importance of each feature through the decision tree method.The importance of each feature for ORR and OER is depicted in Figure S7.It can be found that N e in TM atoms is the dominant factor associated with ORR and OER performance, followed by the covalent bond of the TM element (rcov TM ), the atomic number of TM (Z TM ), and the number of valence electrons in A and B atoms (V A and V B ), etc.Furthermore, targeting the top 10 significant feature descriptors, we constructed quantitative expressions to describe the catalytic activity of ORR and OER, as illustrated in Equations ( 24) and (25).The mean square errors for the aforementioned equations were 0.095 and 0.094, respectively, corresponding to the prediction accuracy of 79.91% and 75.36% for ORR and OER, respectively, as depicted in Figure 10C and D. The accuracy of symbolic regression models is comparable to that of the relevant work [54,55].The aforementioned equations have proposed new ideas for exploring the physicochemical laws between material features by using compact mathematical expression.

Figure 1
Figure 1 The schematic diagram of SAC structures, the considered TM atoms as well as four possible embedding sites of single atoms on the O-terminated MBenes surface, i.e., Hollow, Bridge1, Bridge2, and Top.(A) Top view; (B) side view.
and Au-W 2 B 2 O 2 exhibit superior OER catalytic activity with overpotentials of 0.54, 0.51, 0.20, 0.34, 0.31, 0.45, and 0.55 V, comparable to or even lower than that of IrO 2 /RuO 2 (0.56/0.42 V)[10][11][12].After a comprehensive comparison, Rh-Mo 2 B 2 O 2 and Rh-W 2 B 2 O 2 exhibit the highest comprehensive performance due to the combination of ultra-low ORR/ OER overpotential.It is worth noting that these SAC catalysts with high ORR and/or OER catalytic activity are mainly concentrated in several specific single atoms such as Rh, Ir, Pd, and Au, indicating that the essential regulatory role of single atoms dominates the catalytic activity in SAC systems.In order to further understand the detailed processes of catalysis, the Gibbs free energy diagrams were drawn and shown in Figure3A and Band Figure S3, from which the energy change of each elementary step and the potential determining steps (PDS) can be clearly observed.For Rh-Mo 2 B 2 O 2 and Rh-W 2 B 2 O 2 , the PDS of ORR are both the process of O 2 protonation to OOH*, and the PDS of OER is both the process of OH* deprotonation to O*.It is worth noting that the absolute value of Gibbs free energy change of the four elementary steps of Rh-Mo 2 B 2 O 2 and Rh-W 2 B 2 O 2 are close to the equilibrium potential of 1.23 V for ORR and OER, and thus

Figure 2
Figure 2The cluster energies and the dissolution potentials of TM atoms.The SACs in gray background area are thermodynamically and electrochemically stable, that is, U diss > 0 V and E cluster < 0 eV.

Figure 3
Figure 3 The Gibbs free energy diagrams for ORR and OER of (A) Rh-Mo 2 B 2 O 2 , and (B) Rh-W 2 B 2 O 2 .(C) Histograms of ORR and OER overpotentials for all SAC systems.

Figure 4
Figure 4 Scaling relations between the adsorption free energy of intermediate O*, OH*, and OOH* on all TM-M 2 B 2 O 2 systems (TM = 4d/5d period atoms; M = Ti, Mo, and W).
E) is the density of states (DOS) of the d orbitals of the TM atoms.Subsequently, the relationship between the adsorption free energy of the intermediates (O*, OH*, and OOH*) and the d-band center of the SACs are depicted in Figure7A.It can be observed that the ΔG O* , ΔG OH* , and ΔG OOH* exhibit a similar trend as the value of the d-band center increases.However, the correlation of the ΔG O* , ΔG OH* , and ΔG OOH* versus ε d is relatively poor, making ε d unsuitable to describe the OER/ORR catalytic activity for SACs.In fact, the closer the transition metal atom's electronic states to the Fermi level, the greater its contribution to the interaction with the adsorbents.Therefore, it is necessary to introduce a weight function to quantitatively consider the unequal weight contribution of electronic states at different distances from Fermi energy level.

Figure 6 (
Figure 6 (A) 2D maps of H 2 O production at 0.9 V using ΔG OH* and ΔG OOH* as variables.(B) The potential-dependent TOFs of H 2 O production on the pure Pt, Rh-Mo 2 B 2 O 2 SACs and Rh-W 2 B 2 O 2 SACs.

Figure 7 (
Figure 7 (A) The Gibbs free energy of intermediate adsorbents (O*, OH*, and OOH*) in TM-M 2 B 2 O 2 as a function of d-band center.(B) Schematic diagram of weight function and improved d-band center.(C) The Gibbs free energy of intermediate adsorbents (O*, OH*, and OOH*) in TM-M 2 B 2 O 2 as a function of the improved d-band center using kT = 0.5 eV (TM = 4d/5d period atoms; M = Ti, Mo, and W).

Figure 8 (
Figure 8 (A) The linear relationship between ICOHP and the Gibbs free energy of OH* intermediate in all TM-M 2 B 2 O 2 .(B) Histograms of ICOHP for all TM-M 2 B 2 O 2 systems (TM = 4d/5d period atoms; M = Ti, Mo, and W).

Figure 9
Figure 9The correlation coefficients of (A) 24 descriptors and (B) final 10 features.

Figure 9B and Table
Figure 9B and TableS6.

CONCLUSION
In summary, a series of SACs composed of 4d and 5d period transition atoms anchored onto three MBenes (Ti 2 B 2 O 2 , Mo 2 B 2 O 2 , and W 2 B 2 O 2 ) were designed to comprehensively study their ORR and OER catalytic activity based on first principles calculations and ML.
The results show that Rh-Mo 2 B 2 O 2 , Ir-Mo 2 B 2 O 2 , Rh-W 2 B 2 O 2 , and Pt-W 2 B 2 O 2 have superior ORR catalytic activity, with overpotentials of 0.41, 0.39, 0.42, 0.19, and 0.30 V, respectively, comparable to or even lower than that of Pt (0.45 V).In addition, Rh-Ti 2 B 2 O 2 , Ir-Ti 2 B 2 O 2 , Rh-Mo 2 B 2 O 2 , Au-Mo 2 B 2 O 2 , Rh-W 2 B 2 O 2 , Ir-W 2 B 2 O 2 , and Au-W 2 B 2 O 2 have superior OER catalytic activity, with overpotentials of 0.54, 0.51, 0.20, 0.34, 0.31, 0.45, and 0.55 V, respectively, comparable to or even lower than that of IrO 2 /RuO 2 (0.56/0.42 V).Specifically, Rh-Mo 2 B 2 O 2 and Rh-W 2 B 2 O 2 exhibit both ultra-low overpotentials of ORR and OER, making them promising bifunctional catalysts.Moreover, we found that the influence of transition metal atoms on the catalytic activity of SACs is dominant compared to the influence of the MBene substrates.Combining the scaling relationships of intermediate species, volcano maps and contour maps were established to explore the adsorption behavior of the SACs.Furthermore, the analysis of d-band center and ICOHP explained that the moderate interaction strength between the SACs and intermediate species is a source of excellent catalytic activity.Finally, through ML, the models containing 10 intrinsic element features of SACs were established, and the symbolic expressions describing the catalytic activity were provided to quickly and accurately identify catalysts with excellent ORR and OER catalytic activity.Our work provides several novel and promising ORR and OER electrocatalysts, and develops efficient models to reasonably guide and design new catalysts.