Issue
Natl Sci Open
Volume 3, Number 2, 2024
Special Topic: AI for Chemistry
Article Number 20230043
Number of page(s) 18
Section Chemistry
DOI https://doi.org/10.1360/nso/20230043
Published online 01 February 2024

© The Author(s) 2023. Published by Science Press and EDP Sciences.

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

INTRODUCTION

The oxygen evolution reaction (OER) and oxygen reduction reaction (ORR) are important half-reactions in metal-air batteries and proton-exchange membrane full cells (PEMFC), respectively, which play crucial roles in green energy storage and conversion [18]. However, their four-electron reaction process results in overall slow kinetics, requiring efficient catalysts to accelerate the reaction. Conventionally, noble metal Pt and noble metal oxides RuO2 or IrO2 are still the most commonly used catalysts to promote ORR or OER, with overpotentials of 0.45 V (Pt for ORR) [9] and 0.42/0.56 eV (RuO2/IrO2 for OER) [1012], but the high cost, scarcity, poor robustness, and low reaction selectivity severely restricted their large-scale practical applications [1315]. Therefore, it is necessary to design stable, efficient, and non-precious metal ORR and OER electrocatalysts.

Due to the large specific surface area, unusual electronic properties, and high stability, two-dimensional (2D) MBenes have received widespread attention since their proposal by Guo et al. [16] in 2017, especially in the field of electrocatalysts [1721]. Experimentally, MBenes can be obtained by selectively etching the A atomic layer from the MAB parent phase using different synthesis means (where M represents transition metal, A is the main group element of IIIA or IVA, and B is boron element) [2226]. The outermost exposed M metal layer of MBenes is easily functionalized by O, OH, and F terminal groups in solution, which greatly reduces the adsorptions of ORR and OER reactants on the surface [19,24,27]. Therefore, increasing the adsorption capacity of ORR and OER reactants by regulating the surface electronic structure of functionalized MBenes is the key to the practical applications of MBenes in electrocatalytic ORR and OER.

Introducing single metal atoms on the surface of 2D materials is one of the commonly used methods to adjust the surface electronic structure [2832]. In fact, the introduced single metal atoms can directly participate in ORR and OER as the catalytic activity site. On the one hand, due to the specific single active site, the formed single-atom catalysts (SACs) have high reaction selectivity, maximum atom utilization, and high intrinsic catalytic activity. On the other hand, 2D materials as the substrate can not only change the electron distribution of a single atom to activate suitable reactants but also firmly fix the single atom to prevent clusters caused by excessive surface energy of SACs, thus maintaining the overall high stability. However, finding efficient SACs for ORR and OER from numerous combinations of single-atom and 2D substrate materials is a huge challenge, and the mechanism of metal-support interactions also urgently needs to be further revealed.

Recently, driven by artificial intelligence technology, machine learning (ML) has become a powerful assistant for quickly screening superior ORR and OER catalysts. Liu et al. [33] proposed that the number of d electrons in transition metal (TM) atoms (Ne), the radius of TM atoms (rd), and the charge transfer of TM atoms (Qe) are the three most important features affecting the catalytic behavior by implementing the ML algorithm. Performing ML, Anand et al. [34] found that the embedding energy (Ee) of doped TM atoms, the distance between adsorbed OH and doped TM atoms (d-OH), and the distance between adsorbed O and doped TM atoms (d-O), the differences between Ee and cohesive energy (E-C) are main features determining the catalytic activity of ORR and OER. These extracted features are based on the density functional theory (DFT) and the architecture of systems, which are time-consuming and complex. Thus, proposing a more convenient descriptor based solely on element features is of great significance.

Herein, based on our previous work [19], the most stable O-terminal MBenes (Ti2B2O2, Mo2B2O2, and W2B2O2) are selected as substrates. Using the first principles calculations, a systematic study was conducted on the ORR and OER catalytic activity of 54 SACs composed of 18 kinds of 4d, 5d period TM atoms anchored onto the three substrates mentioned above. By calculating the Gibbs free energy and overpotential of ORR and OER, Pd-Ti2B2O2, Rh-Mo2B2O2, Ir-Mo2B2O2, Rh-W2B2O2, and Pt-W2B2O2 were selected as superior ORR catalysts, 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). Moreover, Rh-Ti2B2O2, Ir-Ti2B2O2, Rh-Mo2B2O2, Au-Mo2B2O2, Rh-W2B2O2, Ir-W2B2O2, and Au-W2B2O2 were selected as superior OER catalysts, 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 IrO2/RuO2(0.56/0.42 V). Specifically, Rh-Mo2B2O2 and Rh-W2B2O2 exhibit high ORR and OER catalytic activity at the same time, and are considered promising bifunctional catalysts. In addition, we found that the adjusting single atom has a more significant effect on catalytic activity compared to adjusting the substrates for SACs. 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]:

(1)

(2)

(3)

(4)

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]:

(5)

(6)

(7)

(8)

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.

(9)

where ∆E is the total energy change before and after the elementary reaction, which can be directly obtained from DFT calculations. ∆ZPE and TS represent the changes in zero-point energy and entropic energy, respectively. represents the Gibbs free energy correction at the specified pH, where kB is the Boltzmann constant. The last term is the correction of electrical potential to Gibbs free energy, denoted as ∆GU = −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.15 K,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).

(10)

(11)

(12)

(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:

(14)

Similarly, the theoretical overpotential ηOER was defined as

(15)

According to the definition of overpotential, for both ORR and OER, the lower the overpotential, the higher the catalytic activity.

RESULTS AND DISCUSSION

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-M2B2O2, 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

thumbnail 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.

(16)

where ETM-MBene, EMBene, and ETM 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 Eb, 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 (Ecluster = EbEcoh) was calculated, where Ecoh is the metal cohesive energy. Herein, Ecoh can be calculated by the energy difference between a metal atom in the most stable bulk phase and in the vacuum, expressed as Ecoh = ETM-bulk/nETM-single, where n is the number of atoms in the most stable bulk phase. A value of Ecluster less than 0 eV is usually considered a powerful criterion for determining the thermodynamic stable of SAC systems. In addition, the dissolution potentials (Udiss = Udiss-bulkEcluster / (eNe)) of TM atoms were further calculated to explore the stability of the studied SAC systems in electrochemical environments. Udiss-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 Udiss indicates that the SACs are electrochemically stable. The Ecluster and Udiss of all SACs are depicted in Figure 2, where the SACs located in the gray background area are thermodynamic and electrochemical stable (Udiss > 0 V and Ecluster < 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-Mo2B2O2 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-Mo2B2O2 is thermodynamically stable under aqueous solution.

thumbnail Figure 2

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

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–S4. Obviously, Pd-Ti2B2O2, Rh-Mo2B2O2, Ir-Mo2B2O2, Rh-W2B2O2, and Pt-W2B2O2 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-Ti2B2O2, Ir-Ti2B2O2, Rh-Mo2B2O2, Au-Mo2B2O2, Rh-W2B2O2, Ir-W2B2O2, and Au-W2B2O2 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 IrO2/RuO2(0.56/0.42 V) [1012]. After a comprehensive comparison, Rh-Mo2B2O2 and Rh-W2B2O2 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 Figure 3A and B and Figure S3, from which the energy change of each elementary step and the potential determining steps (PDS) can be clearly observed. For Rh-Mo2B2O2 and Rh-W2B2O2, the PDS of ORR are both the process of O2 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-Mo2B2O2 and Rh-W2B2O2 are close to the equilibrium potential of 1.23 V for ORR and OER, and thus they have ultra-low ηORR and ηOER, exhibiting excellent ORR and OER bifunctional catalytic activity.

thumbnail Figure 3

The Gibbs free energy diagrams for ORR and OER of (A) Rh-Mo2B2O2, and (B) Rh-W2B2O2. (C) Histograms of ORR and OER overpotentials for all SAC systems.

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 (ΔGO*, ΔGOH*, and ΔGOOH*) in all SACs and explored their relationships. As shown in Figure 4, both ΔGOOH* and ΔGO* show a positive correlation with the variation of ΔGOH*. Their correlation equations can be linearly fitted as ΔGOOH* = 0.79ΔGOH* + 3.19 (R2 = 0.97) and ΔGO* = 1.4ΔGOH* + 1.35 (R2 = 0.76), respectively.

thumbnail Figure 4

Scaling relations between the adsorption free energy of intermediate O*, OH*, and OOH* on all TM-M2B2O2 systems (TM = 4d/5d period atoms; M = Ti, Mo, and W).

Based on the scaling relationship, the volcanic curves of ΔGOH* versus −ηORR, as well as ΔGO* − ΔGOH* 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 H2O (OH* + H+ + e → * + H2O), 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 (* + O2 + 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 exhibit optimal catalytic activity, such as Pd-Ti2B2O2, Rh-Mo2B2O2, and Rh-W2B2O2. Besides, a small amount of catalyst deviates from the volcanic curves, such as Ag-Ti2B2O2, Cd-Mo2B2O2, and Cd-W2B2O2, mainly due to the formation of O* intermediate being their PDS (OOH* + H+ + e → O* + H2O). Similarly, for ΔGO* −ΔGOH* versus −ηOER volcanic curves, the PDS of the catalysts located on the left and right branches are the formation of OOH* intermediate (O* + H2O → OOH* + H+ + e) and O* intermediate (OH* → O* + H+ + e), respectively. The SACs with the best OER catalytic activity were located at the top of the volcanic plots, occupied by Rh-Ti2B2O2, Rh-Mo2B2O2, Au-Mo2B2O2, Rh-W2B2O2, and Ir-W2B2O2. Moreover, the PDS of these SACs located off the volcanic curves is the desorption of O2 (OOH* + H2O → O2 + H+ + e), such as Hf-Ti2B2O2, Zr-Ti2B2O2, Hf-Mo2B2O2, Zr-Mo2B2O2, Hf-W2B2O2, Zr-W2B2O2, and Ta-W2B2O2. Therefore, the combination of superior OER and ORR catalytic activity renders Rh-Mo2B2O2 and Rh-W2B2O2 SACs promising comprehensive bifunctional catalysts.

thumbnail Figure 5

(A) The volcano plot of ΔGOH* versus −ηORR as well as (B) ΔGO* − ΔGOH* versus −ηOER on TM-M2B2O2. The contour maps of (C) ORR and (D) OER catalytic activity on TM-M2B2O2 with ΔGO* − ΔGOH* and ΔGOH* as variables (TM = 4d/5d period atoms; M = Ti, Mo, and W). Here the potential-determining steps are labeled.

Furthermore, in order to more intuitively display the PDS and overpotential distribution regions of all SACs, variables ΔGOH* and ΔGO* −ΔGOH* were selected to depict the catalytic activity of ORR and OER. Specifically, by replacing ΔGOOH* with ΔGOOH* = 0.79ΔGOH* + 3.19 (R2 = 0.97) in Equations (10)–(13), the four free-energy changes can be further expressed as the following Equations (17)–(20), containing only ΔGOH* and ΔGO* − ΔGOH* variables.

(17)

(18)

(19)

(20)

Subsequently, the contour map of ORR catalytic activity regarding ΔGOH* and ΔGO* − ΔGOH* 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 (ΔGOH* < −0.5 eV), making the PDS solely determined by ΔGOH*, 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 weakening of the adsorption strength of MBenes-supported SACs for OH* intermediate results in the PDS jointly dominated by ΔGOH* and ΔGO* − ΔGOH*, and the region where the SACs located correspondingly shifts from Zone 4 (OH* → H2O) to Zone 2 (OOH* → O*) and Zone 1 (O2 → OOH*). Due to the limitations of the scaling relationship, the idea ORR SACs should be located at the boundaries of the four zones. Among all SACs, Pd-Ti2B2O2, Rh-Mo2B2O2, Ir-Mo2B2O2, Rh-W2B2O2, and Pt-W2B2O2 are closest to the boundaries of the four zones and exhibit excellent ORR catalytic activity with an ultra-low overpotential of 0.41, 0.39, 0.42, 0.19 and 0.30 V, respectively.

In the same way, the OER contour map diagram is presented in Figure 5D, which can also be divided into four regions by dashed lines. Accompanied by an increase in ΔGO* − ΔGOH*, the PDS of the SACs shifts from the third step (O* → OOH*) to the fourth step (OOH* → O2) and the second step (OH* → O*). Rh-Mo2B2O2, Au-Mo2B2O2, and Rh-W2B2O2 are closest to the boundary of the four zones and exhibit optimal OER catalytic activity with an ultra-low overpotential of 0.21, 0.34, and 0.32 V, respectively.

Finally, the microkinetic simulations were performed to evaluate the respective turnover frequency (TOF, molecules/sites) of H2O 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 ΔGOH* and ΔGOOH*, the SACs are dispersed at both sides of the fitting line, where Rh-W2B2O2 SACs are closest to the optimal region, exhibiting the highest TOF of H2O production. At the same time, it is not difficult to see that it is worth looking forward to breaking the original scaling relationship through further surface modification to improve the ORR catalytic activity of SACs, that is, the SACs fall into the optimal region. The potential-dependent TOFs of H2O production are plotted and depicted in Figure 6B. Compared with the simulation results of the half-wave potentials on pure Pt for ORR [46], the corresponding half-wave potentials of Rh-Mo2B2O2 SACs and Rh-W2B2O2 SACs exhibit positive shifts of 40 and 140 mV, respectively, further demonstrating the superior ORR catalytic activity of Rh-Mo2B2O2 SACs and Rh-W2B2O2 SACs.

thumbnail Figure 6

(A) 2D maps of H2O production at 0.9 V using ΔGOH* and ΔGOOH* as variables. (B) The potential-dependent TOFs of H2O production on the pure Pt, Rh-Mo2B2O2 SACs and Rh-W2B2O2 SACs.

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,4750], typically defined as

(21)

where D(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 Figure 7A. It can be observed that the ΔGO*, ΔGOH*, and ΔGOOH* exhibit a similar trend as the value of the d-band center increases. However, the correlation of the ΔGO*, ΔGOH*, and ΔGOOH* 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. In this work, the derivative of the Fermi-Dirac function, , 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

thumbnail Figure 7

(A) The Gibbs free energy of intermediate adsorbents (O*, OH*, and OOH*) in TM-M2B2O2 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-M2B2O2 as a function of the improved d-band center using kT = 0.5 eV (TM = 4d/5d period atoms; M = Ti, Mo, and W).

(22)

where EF is the energy of the Fermi level and is the Fermi-Dirac function, expressed as

(23)

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 anti-bonding states. With the decrease of the d-band center, the anti-bonding energy band generated by coupling 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 ΔGOH* 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 ΔGOH* is depicted in Figure 8A. It can be observed that there is a linear correlation between ICOHP and ΔGOH*. 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-M2B2O2 (M = Ti, Mo, and W), as shown in Figure 8B. The more negative ICOHP, the smaller ΔGOH*, 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.

thumbnail Figure 8

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

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-M2B2O2, TM = Ti, V, Cr, Mn, Fe, Co, Ni, Cu, and Zn, M = Ti, Mo, and W) [19], 24 SAC systems of Ti2CO2-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 Ti3C2O2-supported TM atoms (TM = Mn, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ag, Ir, Pt, and Au) [1]. Furthermore, based on the cross-validation 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 Figure 9B and Table S6.

thumbnail Figure 9

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

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.

thumbnail Figure 10

Comparison between ML-predicted and DFT-calculated (A) ηORR, (B) ηOER using models and (C) ηORR, (D) ηOER using symbolic regression models.

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 Ne in TM atoms is the dominant factor associated with ORR and OER performance, followed by the covalent bond of the TM element (rcovTM), the atomic number of TM (ZTM), and the number of valence electrons in A and B atoms (VA and VB), 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.

(24)

(25)

CONCLUSION

In summary, a series of SACs composed of 4d and 5d period transition atoms anchored onto three MBenes (Ti2B2O2, Mo2B2O2, and W2B2O2) were designed to comprehensively study their ORR and OER catalytic activity based on first principles calculations and ML. The results show that Rh-Mo2B2O2, Ir-Mo2B2O2, Rh-W2B2O2, and Pt-W2B2O2 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-Ti2B2O2, Ir-Ti2B2O2, Rh-Mo2B2O2, Au-Mo2B2O2, Rh-W2B2O2, Ir-W2B2O2, and Au-W2B2O2 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 IrO2/RuO2(0.56/0.42 V). Specifically, Rh-Mo2B2O2 and Rh-W2B2O2 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.

Acknowledgments

We are very grateful for the support of the high-performance computing (HPC) resource platform of Beihang University.

Funding

This work was supported by the National Key Research and Development Program of China (2022YFB3807200).

Author contributions

J.Z. and Z.S. proposed and supervised the project. E.W. and G.W. conceived and wrote the manuscript. All authors read and approved the submission of the manuscript.

Conflict of interest

The authors declare no conflict of interest.

Supplementary information

The supporting information is available online at https://doi.org/10.1360/nso/20230043. The supporting materials are published as submitted, without typesetting or editing. The responsibility for scientific accuracy and content remains entirely with the authors.

References

  • Peng Q, Zhou J, Chen J, et al. Cu single atoms on Ti2CO2 as a highly efficient oxygen reduction catalyst in a proton exchange membrane fuel cell. J Mater Chem A 2019; 7: 26062-26070. [Article] [Google Scholar]
  • Chen Y, Ji S, Wang Y, et al. Isolated single iron atoms anchored on N-doped porous carbon as an efficient electrocatalyst for the oxygen reduction reaction. Angew Chem Int Ed 2017; 56: 6937-6941. [Article] [Google Scholar]
  • Fu Z, Ling C, Wang J. A Ti3C2O2 supported single atom, trifunctional catalyst for electrochemical reactions. J Mater Chem A 2020; 8: 7801-7807. [Article] [Google Scholar]
  • Gao G, Bottle S, Du A. Understanding the activity and selectivity of single atom catalysts for hydrogen and oxygen evolution via ab initial study. Catal Sci Technol 2018; 8: 996-1001. [Article] [Google Scholar]
  • Gao G, Waclawik ER, Du A. Computational screening of two-dimensional coordination polymers as efficient catalysts for oxygen evolution and reduction reaction. J Catal 2017; 352: 579-585. [Article] [Google Scholar]
  • Ahsan MA, He T, Eid K, et al. Tuning the intermolecular electron transfer of low-dimensional and metal-free BCN/C60 electrocatalysts via interfacial defects for efficient hydrogen and oxygen electrochemistry. J Am Chem Soc 2021; 143: 1203-1215. [Article] [Google Scholar]
  • Ahsan MA, He T, Eid K, et al. Controlling the interfacial charge polarization of MOF-derived 0D-2D vdW architectures as a unique strategy for bifunctional oxygen electrocatalysis. ACS Appl Mater Interfaces 2022; 14: 3919-3929. [Article] [Google Scholar]
  • Bai L, Hsu CS, Alexander DTL, et al. Double-atom catalysts as a molecular platform for heterogeneous oxygen evolution electrocatalysis. Nat Energy 2021; 6: 1054-1066. [Article] [Google Scholar]
  • Nørskov JK, Rossmeisl J, Logadottir A, et al. Origin of the overpotential for oxygen reduction at a fuel-cell cathode. J Phys Chem B 2004; 108: 17886-17892. [Article] [Google Scholar]
  • Man IC, Su H, Calle-Vallejo F, et al. Universality in oxygen evolution electrocatalysis on oxide surfaces. ChemCatChem 2011; 3: 1159-1165. [Article] [Google Scholar]
  • Cao L, Luo Q, Chen J, et al. Dynamic oxygen adsorption on single-atomic ruthenium catalyst with high performance for acidic oxygen evolution reaction. Nat Commun 2019; 10: 4849. [Article] [Google Scholar]
  • Rossmeisl J, Qu ZW, Zhu H, et al. Electrolysis of water on oxide surfaces. J Electroanal Chem 2007; 607: 83-89. [Article] [Google Scholar]
  • Yu Y, Zhou J, Sun Z. Novel 2D transition-metal carbides: Ultrahigh performance electrocatalysts for overall water splitting and oxygen reduction. Adv Funct Mater 2020; 30: 2000570. [Article] [Google Scholar]
  • Yang X, Zhang Y, Fu Z, et al. Tailoring the electronic structure of transition metals by the V2C MXene support: Excellent oxygen reduction performance triggered by metal-support interactions. ACS Appl Mater Interfaces 2020; 12: 28206-28216. [Article] [Google Scholar]
  • Xu H, Wang D, Yang P, et al. Atomically dispersed M–N–C catalysts for the oxygen reduction reaction. J Mater Chem A 2020; 8: 23187-23201. [Article] [Google Scholar]
  • Guo Z, Zhou J, Sun Z. New two-dimensional transition metal borides for Li ion batteries and electrocatalysis. J Mater Chem A 2017; 5: 23530-23535. [Article] [Google Scholar]
  • Li F, Tang Q. First-principles calculations of TiB mbene monolayers for hydrogen evolution. ACS Appl Nano Mater 2019; 2: 7220-7229. [Article] [Google Scholar]
  • Zhang B, Zhou J, Guo Z, et al. Two-dimensional chromium boride MBenes with high HER catalytic activity. Appl Surf Sci 2020; 500: 144248. [Article] [Google Scholar]
  • Wang E, Zhang B, Zhou J, et al. High catalytic activity of MBenes-supported single atom catalysts for oxygen reduction and oxygen evolution reaction. Appl Surf Sci 2022; 604: 154522. [Article] [Google Scholar]
  • Zhu X, Zhou X, Jing Y, et al. Electrochemical synthesis of urea on MBenes. Nat Commun 2021; 12: 4080. [Article] [Google Scholar]
  • Wu H, Gao Z, Ma F, et al. Surface functionalization of two-dimensional boridene family: Enhanced stability, tunable electronic property, and high catalytic activity. Appl Surf Sci 2022; 602: 154374. [Article] [Google Scholar]
  • Zhang B, Zhou J, Sun Z. MBenes: Progress, challenges and future. J Mater Chem A 2022; 10: 15865-15880. [Article] [Google Scholar]
  • Alameda LT, Lord RW, Barr JA, et al. Multi-step topochemical pathway to metastable Mo2AlB2 and related two-dimensional nanosheet heterostructures. J Am Chem Soc 2019; 141: 10852-10861. [Article] [Google Scholar]
  • Zhou J, Palisaitis J, Halim J, et al. Boridene: Two-dimensional Mo4/3B2−x with ordered metal vacancies obtained by chemical exfoliation. Science 2021; 373: 801-805. [Article] [Google Scholar]
  • Wang J, Ye TN, Gong Y, et al. Discovery of hexagonal ternary phase Ti2InB2 and its evolution to layered boride TiB. Nat Commun 2019; 10: 2284. [Article] [Google Scholar]
  • Alameda LT, Moradifar P, Metzger ZP, et al. Topochemical deintercalation of Al from MoAlB: Stepwise etching pathway, layered intergrowth structures, and two-dimensional MBene. J Am Chem Soc 2018; 140: 8833-8840. [Article] [Google Scholar]
  • Li B, Wu Y, Li N, et al. Single-metal atoms supported on MBenes for robust electrochemical hydrogen evolution. ACS Appl Mater Interfaces 2020; 12: 9261-9267. [Article] [Google Scholar]
  • Liu JC, Ma XL, Li Y, et al. Heterogeneous Fe3 single-cluster catalyst for ammonia synthesis via an associative mechanism. Nat Commun 2018; 9: 1610. [Article] [Google Scholar]
  • Liu C, Li Q, Zhang J, et al. Conversion of dinitrogen to ammonia on Ru atoms supported on boron sheets: A DFT study. J Mater Chem A 2019; 7: 4771-4776. [Article] [Google Scholar]
  • Wang S, Li L, San Hui K, et al. Computational screening of single atoms anchored on defective Mo2CO2 MXene nanosheet as efficient electrocatalysts for the synthesis of ammonia. Adv Eng Mater 2021; 23: 2100405. [Article] [Google Scholar]
  • Zhang M, Du J, Chen Y. Single Cu atom supported on modified h-BN monolayer as n-p codoped catalyst for CO oxidation: A computational study. Catal Today 2021; 368: 148-160. [Article] [Google Scholar]
  • Qi R, Zhu B, Han Z, et al. High-throughput screening of stable single-atom catalysts in CO2 reduction reactions. ACS Catal 2022; 12: 8269-8278. [Article] [Google Scholar]
  • Liu X, Zhang Y, Wang W, et al. Transition metal and N doping on AlP monolayers for bifunctional oxygen electrocatalysts: Density functional theory study assisted by machine learning description. ACS Appl Mater Interfaces 2022; 14: 1249-1259. [Article] [Google Scholar]
  • Anand R, Ram B, Umer M, et al. Doped MXene combinations as highly efficient bifunctional and multifunctional catalysts for water splitting and metal-air batteries. J Mater Chem A 2022; 10: 22500-22511. [Article] [Google Scholar]
  • Kresse G, Furthmüller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B 1996; 54: 11169-11186. [Article] [Google Scholar]
  • Perdew JP, Wang Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys Rev B 1992; 45: 13244-13249. [Article] [Google Scholar]
  • Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Phys Rev Lett 1996; 77: 3865-3868. [Article] [Google Scholar]
  • Moellmann J, Grimme S. DFT-D3 study of some molecular crystals. J Phys Chem C 2014; 118: 7615-7621. [Article] [Google Scholar]
  • Medford AJ, Shi C, Hoffmann MJ, et al. CatMAP: A software package for descriptor-based microkinetic mapping of catalytic trends. Catal Lett 2015; 145: 794-807. [Article] [Google Scholar]
  • Ying Y, Fan K, Luo X, et al. Unravelling the origin of bifunctional OER/ORR activity for single-atom catalysts supported on C2N by DFT and machine learning. J Mater Chem A 2021; 9: 16860-16867. [Article] [Google Scholar]
  • Wang E, Guo M, Zhou J, et al. Reasonable design of MXene-supported dual-atom catalysts with high catalytic activity for hydrogen evolution and oxygen evolution reaction: A first-principles investigation. Materials 2023; 16: 1457. [Article] [Google Scholar]
  • Guo X, Gu J, Lin S, et al. Tackling the activity and selectivity challenges of electrocatalysts toward the nitrogen reduction reaction via atomically dispersed biatom catalysts. J Am Chem Soc 2020; 142: 5709-5721. [Article] [Google Scholar]
  • Kan D, Lian R, Wang D, et al. Screening effective single-atom ORR and OER electrocatalysts from Pt decorated MXenes by first-principles calculations. J Mater Chem A 2020; 8: 17065-17077. [Article] [Google Scholar]
  • Hwang J, Noh SH, Han B. Design of active bifunctional electrocatalysts using single atom doped transition metal dichalcogenides. Appl Surf Sci 2019; 471: 545-552. [Article] [Google Scholar]
  • Zhong L, Li S. Unconventional oxygen reduction reaction mechanism and scaling relation on single-atom catalysts. ACS Catal 2020; 10: 4313-4318. [Article] [Google Scholar]
  • Wang YY, Chen DJ, Allison TC, et al. Effect of surface-bound sulfide on oxygen reduction reaction on Pt: Breaking the scaling relationship and mechanistic insights. J Chem Phys 2019; 150: 041728. [Article] [Google Scholar]
  • Zhang Y, Xu J, Ding Y, et al. Tuning the d-band center enables nickel-iron phosphide nanoprisms as efficient electrocatalyst towards oxygen evolution. Int J Hydrogen Energy 2020; 45: 17388-17397. [Article] [Google Scholar]
  • Lu C, Lee IC, Masel RI, et al. Correlations between the heat of adsorption and the position of the center of the D-band: Differences between computation and experiment. J Phys Chem A 2002; 106: 3084-3091. [Article] [Google Scholar]
  • Nørskov JK, Abild-Pedersen F, Studt F, et al. Density functional theory in surface chemistry and catalysis. Proc Natl Acad Sci USA 2011; 108: 937-943. [Article] [Google Scholar]
  • Toyoda E, Jinnouchi R, Hatanaka T, et al. The d-band structure of pt nanoclusters correlated with the catalytic activity for an oxygen reduction reaction. J Phys Chem C 2011; 115: 21236-21240. [Article] [Google Scholar]
  • Huang B, Xiao L, Lu J, et al. Spatially resolved quantification of the surface reactivity of solid catalysts. Angew Chem Int Ed 2016; 55: 6239-6243. [Article] [Google Scholar]
  • Zhang Y, Chen X, Huang Y, et al. The role of intrinsic defects in electrocatalytic activity of monolayer VS2 basal planes for the hydrogen evolution reaction. J Phys Chem C 2017; 121: 1530-1536. [Article] [Google Scholar]
  • Hammer B, Norskov JK. Why gold is the noblest of all the metals. Nature 1995; 376: 238-240. [Article] [Google Scholar]
  • O’Connor NJ, Jonayat ASM, Janik MJ, et al. Interaction trends between single metal atoms and oxide supports identified with density functional theory and statistical learning. Nat Catal 2018; 1: 531-539. [Article] [Google Scholar]
  • Sun X, Zheng J, Gao Y, et al. Machine-learning-accelerated screening of hydrogen evolution catalysts in MBenes materials. Appl Surf Sci 2020; 526: 146522. [Article] [Google Scholar]

All Figures

thumbnail 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.

In the text
thumbnail Figure 2

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

In the text
thumbnail Figure 3

The Gibbs free energy diagrams for ORR and OER of (A) Rh-Mo2B2O2, and (B) Rh-W2B2O2. (C) Histograms of ORR and OER overpotentials for all SAC systems.

In the text
thumbnail Figure 4

Scaling relations between the adsorption free energy of intermediate O*, OH*, and OOH* on all TM-M2B2O2 systems (TM = 4d/5d period atoms; M = Ti, Mo, and W).

In the text
thumbnail Figure 5

(A) The volcano plot of ΔGOH* versus −ηORR as well as (B) ΔGO* − ΔGOH* versus −ηOER on TM-M2B2O2. The contour maps of (C) ORR and (D) OER catalytic activity on TM-M2B2O2 with ΔGO* − ΔGOH* and ΔGOH* as variables (TM = 4d/5d period atoms; M = Ti, Mo, and W). Here the potential-determining steps are labeled.

In the text
thumbnail Figure 6

(A) 2D maps of H2O production at 0.9 V using ΔGOH* and ΔGOOH* as variables. (B) The potential-dependent TOFs of H2O production on the pure Pt, Rh-Mo2B2O2 SACs and Rh-W2B2O2 SACs.

In the text
thumbnail Figure 7

(A) The Gibbs free energy of intermediate adsorbents (O*, OH*, and OOH*) in TM-M2B2O2 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-M2B2O2 as a function of the improved d-band center using kT = 0.5 eV (TM = 4d/5d period atoms; M = Ti, Mo, and W).

In the text
thumbnail Figure 8

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

In the text
thumbnail Figure 9

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

In the text
thumbnail Figure 10

Comparison between ML-predicted and DFT-calculated (A) ηORR, (B) ηOER using models and (C) ηORR, (D) ηOER using symbolic regression models.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.