Deciphering the structure–activity–selectivity relationship of high-entropy alloys for CO₂ reduction via interpretable machine learning

Abstract

High-entropy alloys (HEAs) have emerged as a promising class of multisite catalysts that exhibit high levels of performance due to a diversity of active sites; however, establishing their structure–performance relationships remains a grand challenge. Herein, we systematically explored the structure–activity–selectivity relationship of HEAs for the CO₂ reduction reaction (CO₂RR) with the assistance of a machine learning framework and density functional theory computations. Statistical analysis of hundreds of thousands of binding energies of *CO, *CHO, and *H on (FeCoNiCuMo)₅₅ clusters revealed that HEAs can break the well-established scaling relationship of pure metal catalysts, but they also face an activity–selectivity tradeoff. This originates from the positive role of the unpaired d electron number in enhancing the binding strength of *CHO and *H and limits the overall performance. Moreover, an activity and a selectivity descriptor were constructed, giving accurate predictions for the performance variations of the reported experiments. On this basis, rapid screening among 26 334 types of HEAs was performed, and 10 promising candidates that balanced activity and selectivity were selected. Our workflow not only provides quantitative criteria to accelerate the rational design of HEA catalysts for the CO₂RR, but it also offers a systematic approach to unraveling the intricate structure–performance relationship in complex systems.

---

Introduction

The interaction between reaction intermediates and catalysts usually plays a dominant role in catalysis.^1–3 Due to the existence of an intrinsic correlation between different binding strengths of different intermediates, pure-metal or single-site catalysts usually exhibit functioning that is far below that of an ideal performance.^4–6 This situation is more pronounced in multistep reactions that require balancing the interaction of a series of intermediates. It is expected that multisite catalysts will overcome this limitation due to the presence of diverse active sites that will play different roles during the entire process.^7–9 Therefore, as typical multisite catalysts, high-entropy alloys (HEAs) have attracted great attention and shown superior catalytic performance in widespread reactions, such as the hydrogen evolution reaction (HER),^10–13 oxygen evolution reaction (OER),^14–16 and oxygen reduction reaction (ORR).^16,17 Despite these achievements, successful applications of HEAs in a typical multistep reaction, i.e., the CO₂ reduction reaction (CO₂RR), are quite rare,^18–20 and how to design and optimize their structure to realize superior performance remains a grand challenge.

The structure–performance relationship unveils how material structure influences its properties and ultimately performance.^21–24 and has been widely adopted to guide catalyst design and optimization. However, the diverse active sites of HEAs also give rise to an extremely complicated structure–performance relationship,^25,26 in addition to the great promise for realizing unattainable performance. Because of the large number of possible active sites in HEAs, it is computationally costly to comprehensively assess the binding strength of reaction intermediates.^27,28 However, diversity in structure undoubtedly leads to the wide difference in binding strength across different active sites,^29,30 and therefore, it is also very challenging to determine how to computationally evaluate the overall performance. As a result, the structure–performance of HEAs remains underexplored, although there have been great efforts in the study of HEA catalysts.

Aiming at this challenge, a machine learning (ML) workflow, coupled with density functional theory (DFT) calculations and statistical approaches, was adopted to explore the structure–activity–selectivity relationship of HEAs for the CO₂RR. Using (FeCoNiCuMo)₅₅ HEA clusters as prototypes, the binding energies of three key intermediates, *CO, *CHO, and *H (ΔE_*CO, ΔE_*CHO and ΔE_*H), on an appropriate portion of possible sites were calculated using DFT computations. On this basis, three ML regression models were trained, realizing the rapid and accurate prediction of ΔE_*CO, ΔE_*CHO, and ΔE_*H on all the possible sites of (FeCoNiCuMo)₅₅ HEA clusters (81 900, 488 250, and 488 250 for *CO, *CHO, and *H, respectively).

Statistical analysis of the obtained hundreds of thousands of binding energies revealed an activity–selectivity tradeoff that leads to a low overall performance of HEAs for the CO₂RR. This rationalizes the aforementioned phenomenon that experimental reports on HEAs for the CO₂RR are very rare, despite the abundance of research in electrocatalysis. Additionally, classification models for activity and selectivity evaluations were established to determine the key factors regarding catalytic performance, and the structure–activity–selectivity relations were explored based on the SHapley Additive exPlanations (SHAP) analyses of these features. Furthermore, we constructed activity and selectivity descriptors using SHAP values, and related experimental reports were referenced to validate the reliability. Finally, these descriptors were employed for the rapid screening of 26 334 types of unknown HEAs, where 10 promising candidates with significantly enhanced performance as compared to reported HEAs were determined.

Computational methods

DFT calculation

The first-principles calculations were implemented using the projector augmented wave method as employed in the Vienna Ab initio Simulation Package (VASP).^31,32 The Perdew–Burke–Ernzerhof (PBE) functional was used to treat the exchange–correlation interactions.^33,34 The van der Waals interactions were incorporated using Grimme's D3-type of semi-empirical method.³⁵ All calculations were performed in a 20 × 20 × 20 ų periodic cubic simulation box, sufficiently large to exclude interactions between adjacent nanoparticles, and the Brillouin zone was sampled by the gamma point. All the geometries were fully optimized, where the convergence criteria for total energy and forces were set at 10⁻⁵ eV and 0.02 eV Å⁻¹, respectively. Spin polarization was considered, and the energy cutoff for the plane-wave basis set was 450 eV. LOBSTER software was used for the crystal orbital Hamilton population (COHP) analysis.³⁶

In the calculation of binding energies, to compare the scaling relations of HEA clusters with those of transition metals, we used the same calculation equations reported by Nørskov et al.^37,38 The binding energies of *CO, *CHO, and *H are calculated through:

ΔE_*CO = E_*CO − E_* − E_CO

where ΔE_*CO, ΔE_*CHO, and ΔE_*H are the binding energies of *CO, *CHO, and *H, respectively. E_*CO, E_*CHO, E_*H, and E_* represent the total energies of *CO, *CHO, *H adsorbed on HEA clusters, and isolated HEA clusters, respectively. E_CO and E_H2 represent the energy of gas-phase species for CO and H₂, respectively.

ML methods

The workflow employing ML methods, including ML regression, ML classification, and SHAP analyses, to investigate the structure–performance relationship in catalysis is presented in Fig. S1. For ML regression, the XGBRegressor algorithm was utilized to build prediction models for the binding energies of *CO/*CHO/*H, and the mean square error (MSE) was adopted to evaluate the model performance.^39,40 For ML classification, the XGBClassifier algorithm was utilized to build models to distinguish the activity/selectivity of HEAs, and the area under the receiver operating characteristic (ROC) curve was utilized to evaluate model performance.^39,41 During the training process, 5-fold cross-validation was employed to mitigate the bias resulting from data splitting and obtain more reliable estimates of model performance, and the optimal hyper-parameters for each ML model were selected to obtain the optimal model (Table S1).⁴⁰ To break the ‘black box’ of ML models for interpretability, SHAP values were computed with the SHAP packages to quantify the marginal contribution of features.⁴²

Results and discussion

To explore the structure–performance for CO₂RR, a Cu-based HEA cluster, namely, (FeCoNiCuMo)₅₅, consisting of 55 metal atoms (M₅₅) with random occupation of five commonly used metal elements, including Fe, Co, Ni, Cu, and Mo, was used as the prototype.^43,44 Due to the random distribution of multiple atoms, the total number of the possible structures for this (FeCoNiCuMo)₅₅ HEA is huge (approximately 16 562 000). Considering the fact that the coordination environment of the active site usually dominates its catalytic performance,^27,29,30 further evaluations were carried out based on the consideration of the local microstructures (MS-HEAs) of the active site, where two distinct local microstructures (penta-coordinate sites (MS_penta-HEAs) and hexa-coordinate sites (MS_hex-HEAs), as shown in Fig. 1a) were involved.

Fig. 1 (a) Structure of FeCoNiCuMo HEA clusters, with top and side views of microstructures (MS_penta-HEAs and MS_hex-HEAs). The different regions of the microstructures (R₁, R₂, and R₃) have been outlined. (b) Coordination features extracted from top (*CO), bridge (*CHO), and hollow (*H) sites. (c) Element vectors extracted from 6 regions (R₁, R₂, R₃, B₁, B₂, and H₁) on MS_penta-HEAs and MS_hex-HEAs. (d) Composition features converted from element vectors by group counting of different elements in different regions. (e) Atomic property features generated from the composition features by grouping the average for atomic properties in different regions. Plots of DFT-calculated vs. ML-predicted values of (f) ΔE_*CO, (g) ΔE_*CHO, and (h) ΔE_*H, and the mean square errors (MSEs) of training and testing sets for model evaluation.

For each structure, the active center (region 1, R₁), along with its nearest neighboring coordination atoms on the surface (region 2, R₂) and the subsurface (region 3, R₃), were considered. As a result, over 80 000 (3150 for MS_penta-HEAs and 78 750 for MS_hex-HEAs, Fig. S2) distinct active centers in M₅₅ HEA clusters were identified, which significantly decreased as compared to the total number of the clusters, but is a sufficiently large number for examination of the activity–performance relationship for the CO₂RR. Therefore, the diverse local structure of active sites in M₅₅ clusters enables comprehensive exploration of the performance trend, although it is relatively small.

The binding energies of three key intermediates, i.e., *CO, *CHO, and *H (ΔE_*CO, ΔE_*CHO, and ΔE_*H), were adopted for the assessment of electrocatalytic performance in the CO₂RR, according to numerous experimental and computational studies.^30,45–47 These intermediates (*CO, *CHO, and *H) prefer to adsorb on the top, bridge, and hollow site of the clusters (Fig. S3), respectively, leading to 81 900 (3150 for MS_penta-HEAs and 78 750 for MS_hex-HEAs), 488 250 (3150 × 5 + 78 750 × 6), and 488 250 (3150 × 5 + 78 750 × 6) distinct adsorption sites (Fig. S4) for *CO, *CHO, and *H, respectively. The difference in the adsorption configuration was characterized using the coordination features X_1, X_2, and X_3, which denote the surface coordination number of the atoms at the adsorption sites bonded to the intermediates (Fig. 1b). Furthermore, apart from the R₁, R₂, and R₃ regions, the *O-connected atoms of bridge sites (bridge region 1, B₁), the nearest neighboring coordination atoms of bridge sites on the surface (bridge region 2, B₂), and the *H-connected atoms of hollow sites (hollow region 1, H₁) were extracted as active regions (Fig. 1c).

On this basis, composition features and atomic property features were introduced to describe the variation of the active sites. Specifically, element vectors corresponding to element names in different regions (R₁, R₂, R₃, B₁, B₂, and H₁) of MS_penta-HEAs and MS_hex-HEAs were extracted and converted into composition features through counting the number of distinct elements within different regions (Fig. 1d). In this way, 15, 25, and 20 features for top, bridge, and hollow sites, respectively, are generated. In addition, the number of unpaired d electrons (N_d-up), the first ionization energy (I₁), and electronegativity (χ) were used to describe the physical and chemical properties of different metal atoms (Fig. 1e). By incorporating the composition feature and atomic property feature into the following function based on the grouping average, 9, 15, and 12 features for top, bridge, and hollow sites, respectively, were generated:

where F_j^k represents the number of j-th elements in the k-th region, F_i^j represents the value of i-th elemental properties of the j-th element, and represents the mean feature value of i-th elemental properties in the k-th region. After feature engineering, the total number of features utilized for predicting ΔE_*CO, ΔE_*CHO, and ΔE_*H is 25, 42, and 35, respectively.

Based on the designed features and the moderate-scale DFT calculations (Fig. S5–S7), three datasets for ΔE_*CO, ΔE_*CHO, and ΔE_*H were compiled, and three adsorption prediction models designated as XGBR_*CO, XGBR_*CHO, and XGBR_*H were trained by employing the Extreme Gradient Boost Regression (XGBR) algorithm.^39,40 The training and testing datasets were randomly selected and divided into an 80% and 20% ratio. To evaluate the predictive accuracy of these models, the mean square errors (MSE) between DFT-calculated and ML-predicted ΔE_*CO, ΔE_*CHO, ΔE_*H were calculated for training and testing datasets.

Our results showed that the XGBR_*CO, XGBR_*CHO, and XGBR_*H models exhibited excellent prediction accuracy when the data scale reached 322, 368, and 831, respectively, where the corresponding MSE values are 0.009/0.011/0.005 eV and 0.032/0.057/0.028 eV for the training and testing sets (Fig. 1f–h). Using the trained models, the ΔE_*CO, ΔE_*CHO, and ΔE_*H on all the possible sites of the (FeCoNiCuMo)₅₅ HEA cluster (81 900, 488 250, and 488 250 for *CO, *CHO, and *H, respectively) were predicted (Fig. S8–S10). Some general tendencies can be obtained from predicted energies: (i) the binding strengths of *CO on the active sites with Cu atoms in the R₁ region are relatively weak (Fig. S8); (ii) the binding strengths of *CHO on the active sites with Mo atoms in the B₁ region are usually strong (Fig. S9); (iii) additional Cu atoms in the H₁ region usually lead to weaker adsorption of *H (Fig. S10). The adsorption energies on most sites deviate from the well-established scaling relations between ΔE_*CO and ΔE_*CHO on the (211) and (111) surfaces of mono-metals (Fig. 2a).^37,38,48 This deviation leads to a low reaction free energy for *CO reduction into *CHO, which is usually the potential-determining step for the CO₂RR into deep-reduced products. Consequently, HEAs demonstrate promising catalytic activity for the CO₂RR by facilitating this crucial reaction pathway.

Fig. 2 (a) ML-predicted ΔE_*CO and ΔE_*CHO values for all MS-HEAs. The solid lines denote the scaling relations on metal (211) substrate (black) and metal (111) substrate (red). (b) ROC curve and normalized confusion matrix of XGBC_activity. (c) SHAP summary plots of every feature in XGBC_activity. The bidirectional arrows at the top of the figure show the relationship between the SHAP values and activity. (d) The relationship between the feature value of N_d-up-B₁ and SHAP values of N_d-up-B₁ in XGBR_*CHO and XGBC_activity. (e) ICOHP and (f) COHP between the O atoms of *CHO and metal atoms in the B₁ region (with different values of N_d-up-B₁). (g) The relationship between the feature value of N_d-up-R₂ and SHAP values of N_d-up-R₂ in XGBR_*CO and XGBC_activity. The (h) mean ICOHP and (i) COHP between the Cu atoms in the R₁ region and metal atoms in the R₂ region (with different values of N_d-up-R₂).

To uncover the driving force behind this phenomenon, an activity classification model (XGBC_activity) based on the XGBoost algorithm was developed, where a dataset combining ΔE_*CO and ΔE_*CHO over 488 250 sites with activity labels (where the sites breaking the scaling relationship were recognized to be highly active in the CO₂RR and were labeled as 1; otherwise were labeled as 0) and 7 key structural features were used. Herein, the strategy for selecting 7 key structural features considered the correlation with intermediate adsorption and the transferability of the model. Specifically, the top 10 features from XGBR_*CO and XGBR_*CHO were selected by taking their union to ensure a high correlation with *CO and *CHO adsorption (Fig. S11). Then, features with atomic properties (such as N_d-up-B₁, N_d-up-R₂, and I₁-R₁) were retained, while features with specific elemental composition information (such as Cu-R₁, Co-R₁, and Mo-B₁) were removed.

The accuracy of the trained classification model was very high, where the area under the curve (AUC) value of the operating characteristic (ROC) curve was 0.987, and the main diagonal values of the normalized confusion matrix reached 0.95 and 0.98 (Fig. 2b). On this basis, SHAP analysis was further carried out,^41,42 and it showed that the highest importance to the activity was attributed to the number of unpaired d electrons at the B₁ region (N_d-up-B₁) (Fig. 2c). This feature also dominated the binding strength of *CHO (Fig. S12). Therefore, the relationship between the feature values of N_d-up-B₁, SHAP values of N_d-up-B₁ in XGBR_*CHO, and SHAP values of N_d-up-B₁ in XGBC_activity were plotted to understand the structure–property–activity relationship from N_d-up-B₁.

As shown in Fig. 2d, highly active sites (the points in the blue dashed box with positive SHAP values of XGBC_activity) generally possess positive SHAP values in XGBR_*CHO (corresponding to the enhancement of *CHO binding strength) and relatively large values (3–5) of N_d-up-B₁. This can be ascribed to the adsorption configuration of *CHO. Specifically, a higher N_d-up of B₁ atoms usually corresponds to a higher oxygen affinity of the metal atom, leading to enhancement of the binding strength of the O–M bond (Fig. 2e, f and S13). Interestingly, a large N_d-up of the coordination atoms also resulted in the relatively weak binding strength of *CO (negative SHAP values of N_d-up-R₂ in XGBR_*CO, Fig. 2g), due to the enhanced bonding strength between the active center and coordination atoms (Fig. 2h, i and S15). Specifically, the mean ICOHP values between the active center and coordination atoms significantly decrease when N_d-up-R₂ values are 3–5, while they remain almost unchanged when N_d-up-R₂ is in the range of 0–3. Therefore, the large number of unpaired d electrons of coordination atoms stabilizes the *CHO and simultaneously weakens the binding strength of *CO, resulting in the breaking of the scaling relationship between ΔE_*CO and ΔE_*CHO on pure metal surfaces and potentially high activity of MS-HEAs for the CO₂RR.

Nevertheless, these active sites exhibit low selectivity towards the CO₂RR, due to the preference of *H adsorption. Our results show that most of the studied active sites possess a more negative *H adsorption energy than that on a pure Cu₅₅ cluster (Fig. 3a). The stronger *H binding strength indicates that these sites will be more easily covered by *H intermediates as compared to Cu.^27,30,49 This facilitates the competing HER or blocks the related active sites, both of which lowers the reaction rate of the CO₂RR. Using a dataset of over 488 250 sites with selectivity labels (where the sites with *H binding strength weaker than Cu were labeled as 1, and otherwise were labeled as 0) and 5 key structural features, the trained classification model exhibited high accuracy for AUC values reaching 0.988 (Fig. 3b). Similar to the binding strength of *CHO, the mean N_d-up of the hollow site (N_d-up-H₁) is also the determining factor for the ΔE_*H (Fig. S16), as well as the selectivity (Fig. 3c), where a larger value of N_d-up-H₁ usually leads to stronger *H-binding strength and lower CO₂RR selectivity (Fig. 3d). Therefore, a high number of unpaired d electrons facilitates the activity but lowers the selectivity of the CO₂RR, suggesting an activity–selectivity tradeoff in the HEAs that results in the low overall performance for the CO₂RR.

Fig. 3 (a) Distributions of the difference between the ML-predicted ΔE_*H values of all MS-HEAs and ΔE_*H values on a pure Cu₅₅ cluster. (b) Receiver operating characteristic (ROC) curve and normalized confusion matrix of XGBC_selectivity for model evaluation. (c) SHAP summary plots of every feature in XGBC_selectivity. The color of each dot represents the feature value. The bidirectional arrows at the top of the figure show the relationship between the SHAP values and selectivity. (d) The relationship between the feature value of N_d-up-H₁ and SHAP values of N_d-up-H₁ in XGBR_*H and XGBC_selectivity. The ratios of active sites (P_activity, P_selectivity, and P_performance) in FeCoNiCuMo MS-HEAs with various key feature values of (e) N_d-up-B₁ or (f) N_d-up-H₁.

To provide a more intuitive picture of the above phenomenon, we calculated the ratios of active sites with high-activity (P_activity), high-selectivity (P_selectivity), and high-performance (with high-activity and high-selectivity, P_performance) with the variation of the key features (N_d-up-B₁ and N_d-up-H₁ for the activity and selectivity, respectively). The ratio is defined as where N_high denotes the number of active sites with high activity, selectivity, or performance and N_total denotes the total number of active sites under specific situations. A general tendency can be observed from Fig. 3e and f, where P_activity increases, while P_selectivity decreases with the increase in N_d-up-B₁ and N_d-up-H₁. As a result, P_performance constantly remains at a low level. These findings highlight the activity–selectivity tradeoff in HEAs, which likely explains the current lack of experimental reports on HEAs for the CO₂RR, despite extensive research on their potential for electrocatalysis.

The structure–activity–selectivity relationship was further explored, to determine which sites of (FeCoNiCuMo)₅₅ HEAs can achieve high activity and selectivity for the CO₂RR. To this end, the sum of the SHAP values of the key features was adopted as the descriptor. Accordingly, two descriptors, labeled as SHAP_activity and SHAP_selectivity, were constructed for the activity and selectivity evaluation, respectively, where the positive value corresponds to high activity or selectivity. Our results show that these two descriptors can provide accurate classification of the activity and selectivity of the active sites, where the values on the main diagonal of the confusion matrix are as high as 0.93/0.97 and 0.99/0.94 for the activity and selectivity, respectively (Fig. 4a and b). By combining SHAP_selectivity and SHAP_activity, the probability density distribution of different types of sites is obtained, where the sites with high-activity and low-selectivity were found to be dominant (Fig. 4c).

Fig. 4 Density distribution of (a) SHAP_activity and (b) SHAP_selectivity for high-activity (high-selectivity)/low-activity (low-selectivity) MS-HEAs with a true label. The red dashed line denotes the boundary between high-activity (high-selectivity) and low-activity (low-selectivity) from the descriptors. (c) The density distribution of SHAP_activity and SHAP_selectivity for all FeCoNiCuMo MS-HEAs. (d) The distribution of P_activity, P_selectivity, and P_performance for different local environments ((M₁M₂M₃)_penta/hex). (e) Violin diagrams of SHAP_activity and SHAP_selectivity for the top 8 local environments ((M₁M₂M₃)_penta/hex).

On the contrary, active sites with high activity and selectivity are very rare, supporting the activity–selectivity tradeoff and poor performance for the CO₂RR of (FeCoNiCuMo)₅₅ HEAs once again. Furthermore, the local environments with excellent performance were extracted, and were defined as (M₁M₂M₃)_penta or (M₁M₂M₃)_hex. Specifically, M₁, M₂, and M₃ are the three metal atoms of a hollow site in MS-HEA, while other atoms are random (Fig. 4d). The 8 local environments with the highest P_performance were selected due to their ability to maintain high-performance while the surrounding atoms change. Our results show that (CuMoCu)_penta, (MoMoNi)_penta, (CuMoCu)_hex, (MoMoNi)_hex, (NiMoMo)_hex, (CuFeCu)_hex, (CuFeCu)_penta, and (MoMoCu)_penta are promising local structures for the CO₂RR (Fig. 4e). In particular, 75% of the active sites that contain a (CuMoCu)_penta center are expected to exhibit excellent performance. This can be well understood from the aforementioned feature analysis. Specifically, the large N_d-up of Mo enhances the *CHO binding strength, and the small N_d-up of Cu atoms leads to the low mean N_d-up of the active center to weaken the *H binding strength, leading to the excellent performance.

The feasibility of the constructed structure–activity–selectivity relationship in other HEAs was further evaluated. Specifically, the performance of two representative HEAs, i.e., AgCuAuPdPt¹⁸ and PtMoPdRhNi¹⁰ that have been experimentally proven to be promising catalysts for the CO₂RR and HER, respectively, were systematically explored. Using (AgCuAuPdPt)₅₅ and (PtMoPdRhNi)₅₅ clusters as models, and SHAP_activity and SHAP_selectivity as the descriptors, the probability density distribution of activity and selectivity was predicted (Fig. 5a and b). The aforementioned activity–selectivity tradeoff still exists in these two systems, where a large portion of active sites exhibits high activity but low selectivity for the CO₂RR. Nevertheless, (AgCuAuPdPt)₅₅ possesses a considerable number of sites with high activity and selectivity (Fig. 5a), in agreement with the experimentally reported satisfactory performance in the CO₂RR.¹⁸ On the contrary, almost no active site lies in the region that meets the requirement of high catalytic performance for (PtMoPdRhNi)₅₅ (Fig. 5b), in accordance with its high HER activity.¹⁰

Fig. 5 Density distribution of SHAP_activity and SHAP_selectivity for all MS-HEAs in (a) AgCuAuPdPt and (b) PtMoPdRhNi HEAs. (c) The ratios of active sites with high-activity (P_activity), high-selectivity (P_selectivity), and high-performance (P_performance) in reported HEAs. (d) All possible five-element combinations of 22 metal atoms. (e) The distribution of P_activity, P_selectivity, and P_performance for all possible five-element combinations. (f) The ratios of active sites for the top 10 high-performance HEAs.

The structure–activity–selectivity relationship was further adopted to other reported HEAs, including FeCoNiCuAl,¹⁹ CoMnNiCuZn,²⁰ CoCuNiZnSn,²⁷ CoCuGaNiZn,³⁰ FeCoNiRuMn,¹² RuRhPdPtIr,⁵⁰ NiFeCoMoW,⁵¹ PtCoNiMoRh,⁵² PtRuFeCoNi,⁵³ PtPdRuMoNi,⁵⁴ IrRuRhMoW,¹⁴ PtCoNiRuIr,⁵⁵ NiFeCrVTi,⁵⁶ and NiCoFePtRh.¹¹ Surprisingly, our model provides a satisfactory description of the performance variation of these systems, supporting the reliability of the proposed models (Fig. S17 and S18). Specifically, for the reported HEAs with a high CO₂RR performance, including FeCoNiCuAl,¹⁹ CoMnNiCuZn,²⁰ CoCuNiZnSn,²⁷ and CoCuGaNiZn,³⁰ there is a relatively high ratio of active sites with high-performance (P_performance = 0.12–0.18) as compared to others with a high HER activity (P_performance = 0–0.05) (Fig. 5c). Therefore, the constructed structure–activity–selectivity relationship deciphers the universal rule of the CO₂RR over HEAs.

With the assistance of the constructed structure–activity–selectivity relation, a high-throughput screening of promising HEAs for the CO₂RR was performed. The candidates were 26 334 HEAs via the possible combination of 22 elements (Fig. 5d). Although the combinations are complex, the descriptors are transferable because the features are based on fundamental atom properties rather than specific elemental composition information. With the assistance of the SHAP_activity and SHAP_selectivity, all the high-performance sites of the 26 334 HEAs were efficiently evaluated, and the top 10 HEAs with the highest P_performance were selected as excellent catalysts for the CO₂RR (Fig. 5e), including AgCuZnCdTa, AgCuZnMnCd, AgCuZnCrCd, AgCuZnVCd, AgCuZnNbCd, AgCuZnCdW, AgCuNiZnCd, CuPdZnCdTa, AgPdZnCdTa, and AgCuPdZnCd (Fig. S19). Clearly, the high-activity ratios (P_activity = 0.72–0.96) and the high-selectivity ratios (P_selectivity = 0.52–0.77) of these promising HEAs reached a balance, leading to a significantly higher P_performance (0.51–0.72) than the reported HEAs (0.12–0.18) (Fig. 5f). Therefore, these HEAs are expected to exhibit excellent overall performance for the CO₂RR.

Conclusions

We have deciphered the structure–activity–selectivity relationship of attractive but complicated HEA catalysts for the CO₂RR, with the assistance of a ML workflow and DFT computations. A systematic and deep analysis of the adsorption properties of *CO, *CHO, and *H across 488 250 different sites of (FeCoNiCuMo)₅₅ HEA clusters uncovered the activity–selectivity tradeoff of HEAs for the CO₂RR. This originates from the positive effect of unpaired d electrons of coordination atoms on enhancing the binding strength of *CHO and *H intermediates. Therefore, most sites are expected to possess high activity by breaking the scaling relation, but low selectivity exists due to the preference for *H adsorption. This degrades the overall performance of HEAs for the CO₂RR and rationalizes the current situation in this field, i.e., there are rare experimental reports on HEAs for the CO₂RR despite the abundance of research into electrocatalysis.

Using the sum of the SHAP values of key features as the descriptors for activity and selectivity, the activity–selectivity tradeoff in (FeCoNiCuMo)₅₅ HEAs is intuitively presented. Moreover, the diversity in the structure and components of the active sites enables the proposed descriptors to be generally applicable to other HEAs, where the reliability is supported by the accurate description of the performance variation of reported HEAs. On this basis, rapid screening among 26 334 types of HEAs was performed, where AgCuZnCdTa, AgCuZnMnCd, AgCuZnCrCd, AgCuZnVCd, AgCuZnNbCd, AgCuZnCdW, AgCuNiZnCd, CuPdZnCdTa, AgPdZnCdTa, and AgCuPdZnCd were selected as the top 10 promising catalysts to achieve a balance between activity and selectivity for the CO₂RR. These predictions are expected to guide experimental research to realize superior performance of HEAs in the CO₂RR. Overall, our research, based on quantitative descriptors, deciphers the structure–activity–selectivity of HEAs for the CO₂RR, provides insights into related phenomenon and guidelines for the design of promising HEAs, and also provides a useful method for the exploration of the structure–performance relationship of complicated systems.

Author contributions

The manuscript was jointly written with contributions from all authors, and all have approved the final version. X. Z., J. W. and C. L. conceived the idea, supervised the project. J. S. performed the machine learning training and DFT calculations, and drafted the manuscript. X. X., Y. M., A. C. and S. W. contributed to selected DFT calculations and the associated data analysis. X. Z., J. W., C. L. and L. S. critically revised the manuscript.

Conflicts of interest

The authors declare no competing financial interests.

Data availability

The datasets and source code utilized for training the models are accessible from GitHub: https://github.com/sjx0327/HEA_CO2RR.

Supplementary information: details of the ML workflow to investigate the structure–performance relation; the number of MS-HEAs and all adsorption sites; ML-predicted ΔE_*CO/ΔE_*CHO/ΔE_*H on all adsorption sites; COHPs and ICOHPs of O–M/M₁–M_2−n bonds; and the density distribution of SHAP_activity and SHAP_selectivity for reported/predicted HEAs. See DOI: https://doi.org/10.1039/d5sc05762k.

Acknowledgements

This work is supported by the National Key Research and Development Program of China (2022YFA1503103), the National Natural Science Foundation of China (22422303, 22033002, 92261112, 22303011, 22203046), the Six Talent Peaks Project in Jiangsu Province (XCL-104), the Basic Research Program of Jiangsu Province (BK20220800), the Open Research Fund of the Key Laboratory of Quantum Materials and Devices (Southeast University), the Natural Science Research Start-up Foundation of Recruiting Talents of Nanjing University of Posts and Telecommunications (NY221128), and the Project of the State Key Laboratory of Organic Electronics and Information Displays, Nanjing University of Posts and Telecommunications (GZR2023010003). The authors are thankful for computational resources from the BigData Computing Center of Southeast University and the Center for Fundamental and Interdisciplinary Sciences of Southeast University.

References

1. X. Tian, X. Zhao, Y.-Q. Su, L. Wang, H. Wang, D. Dang, B. Chi, H. Liu, E. J. M. Hensen and X. W. Lou, et al., Engineering bunched Pt-Ni alloy nanocages for efficient oxygen reduction in practical fuel cells, Science, 2019, 366, 850–856.
2. S. Popović, M. Smiljanić, P. Jovanovič, J. Vavra, R. Buonsanti and N. Hodnik, Stability and Degradation Mechanisms of Copper-Based Catalysts for Electrochemical CO2 Reduction, Angew. Chem., Int. Ed., 2020, 59, 14736–14746.
3. C. Ren, S. Lu, Y. Wu, Y. Ouyang, Y. Zhang, Q. Li, C. Ling and J. Wang, A Universal Descriptor for Complicated Interfacial Effects on Electrochemical Reduction Reactions, J. Am. Chem. Soc., 2022, 144, 12874–12883.
4. M. B. Gawande, A. Goswami, F.-X. Felpin, T. Asefa, X. Huang, R. Silva, X. Zou, R. Zboril and R. S. Varma, Cu and Cu-Based Nanoparticles: Synthesis and Applications in Catalysis, Chem. Rev., 2016, 116, 3722–3811.
5. A. Wang, J. Li and T. Zhang, Heterogeneous single-atom catalysis, Nat. Rev. Chem., 2018, 2, 65–81.
6. Y. Shang, X. Xu, B. Gao, S. Wang and X. Duan, Single-atom catalysis in advanced oxidation processes for environmental remediation, Chem. Soc. Rev., 2021, 50, 5281–5322.
7. Z. Cao, Q. Chen, J. Zhang, H. Li, Y. Jiang, S. Shen, G. Fu, B.-a. Lu, Z. Xie and L. Zheng, Platinum-nickel alloy excavated nano-multipods with hexagonal close-packed structure and superior activity towards hydrogen evolution reaction, Nat. Commun., 2017, 8, 15131.
8. B. Wang, K. Zhao, Z. Yu, C. Sun, Z. Wang, N. Feng, L. Mai, Y. Wang and Y. Xia, In situ structural evolution of the multi-site alloy electrocatalyst to manipulate the intermediate for enhanced water oxidation reaction, Energy Environ. Sci., 2020, 13, 2200–2208.
9. J. Li, Y. Chen, B. Yao, W. Yang, X. Cui, H. Liu, S. Dai, S. Xi, Z. Sun and W. Chen, et al., Cascade Dual Sites Modulate Local CO Coverage and Hydrogen-Binding Strength to Boost CO2 Electroreduction to Ethylene, J. Am. Chem. Soc., 2024, 146, 5693–5701.
10. M. Wei, Y. Sun, J. Zhang, F. Ai, S. Xi and J. Wang, High-entropy alloy nanocrystal assembled by nanosheets with d–d electron interaction for hydrogen evolution reaction, Energy Environ. Sci., 2023, 16, 4009–4019.
11. G. Feng, F. Ning, J. Song, H. Shang, K. Zhang, Z. Ding, P. Gao, W. Chu and D. Xia, Sub-2 nm Ultrasmall High-Entropy Alloy Nanoparticles for Extremely Superior Electrocatalytic Hydrogen Evolution, J. Am. Chem. Soc., 2021, 143, 17117–17127.
12. J. Hao, Z. Zhuang, K. Cao, G. Gao, C. Wang, F. Lai, S. Lu, P. Ma, W. Dong and T. Liu, et al., Unraveling the electronegativity-dominated intermediate adsorption on high-entropy alloy electrocatalysts, Nat. Commun., 2022, 13, 2662.
13. H. Minamihara, K. Kusada, D. Wu, T. Yamamoto, T. Toriyama, S. Matsumura, L. S. R. Kumara, K. Ohara, O. Sakata and S. Kawaguchi, et al., Continuous-Flow Reactor Synthesis for Homogeneous 1 nm-Sized Extremely Small High-Entropy Alloy Nanoparticles, J. Am. Chem. Soc., 2022, 144, 11525–11529.
14. H. Luo, L. Li, F. Lin, Q. Zhang, K. Wang, D. Wang, L. Gu, M. Luo, F. Lv and S. Guo, Sub-2 nm Microstrained High-Entropy-Alloy Nanoparticles Boost Hydrogen Electrocatalysis, Adv. Mater., 2024, 36, 2403674.
15. Y. Mei, Y. Feng, C. Zhang, Y. Zhang, Q. Qi and J. Hu, High-Entropy Alloy with Mo-Coordination as Efficient Electrocatalyst for Oxygen Evolution Reaction, ACS Catal., 2022, 12, 10808–10817.
16. R. He, L. Yang, Y. Zhang, D. Jiang, S. Lee, S. Horta, Z. Liang, X. Lu, A. Ostovari Moghaddam and J. Li, et al., A 3d-4d-5d High Entropy Alloy as a Bifunctional Oxygen Catalyst for Robust Aqueous Zinc–Air Batteries, Adv. Mater., 2023, 35, 2303719.
17. X. Zhao, H. Cheng, X. Chen, Q. Zhang, C. Li, J. Xie, N. Marinkovic, L. Ma, J.-C. Zheng and K. Sasaki, Multiple Metal–Nitrogen Bonds Synergistically Boosting the Activity and Durability of High-Entropy Alloy Electrocatalysts, J. Am. Chem. Soc., 2024, 146, 3010–3022.
18. S. Nellaiappan, N. K. Katiyar, R. Kumar, A. Parui, K. D. Malviya, K. G. Pradeep, A. K. Singh, S. Sharma, C. S. Tiwary and K. Biswas, High-Entropy Alloys as Catalysts for the CO2 and CO Reduction Reactions: Experimental Realization, ACS Catal., 2020, 10, 3658–3663.
19. W. L. Ng, A. C. M. Loy, C. Kundu, S. Ebrahimian, Y. Hora, T. Williams and S. Bhattacharya, Customizable Quinary FeCoNiCuAl-Based High-Entropy Alloy Nanoparticles Supported on 3D-Printed Monolith for CO2 Hydrogenation, ACS Appl. Nano Mater., 2024, 7, 25314–25321.
20. J. Zhao, J. Bao, S. Yang, Q. Niu, R. Xie, Q. Zhang, M. Chen, P. Zhang and S. Dai, Exsolution–Dissolution of Supported Metals on High-Entropy Co3MnNiCuZnOx: Toward Sintering-Resistant Catalysis, ACS Catal., 2021, 11, 12247–12257.
21. C. Ren, Q. Li, C. Ling and J. Wang, Mechanism-Guided Design of Photocatalysts for CO2 Reduction toward Multicarbon Products, J. Am. Chem. Soc., 2023, 145, 28276–28283.
22. Y. Yuan, K. He and J. Lu, Structure–Property Interplay Within Microporous Manganese Dioxide Tunnels For Sustainable Energy Storage, Angew. Chem., Int. Ed., 2024, 63, e202316055.
23. C. Wu, N. Corrigan, C.-H. Lim, W. Liu, G. Miyake and C. Boyer, Rational Design of Photocatalysts for Controlled Polymerization: Effect of Structures on Photocatalytic Activities, Chem. Rev., 2022, 122, 5476–5518.
24. X. Lin, S. Zhen, X. Wang, L. V. Moskaleva, P. Zhang, Z.-J. Zhao and J. Gong, Data-Driven Design of Single-Atom Electrocatalysts with Intrinsic Descriptors for Carbon Dioxide Reduction Reaction, Trans. Tianjin Univ., 2024, 30, 459–469.
25. H. Li, J. Lai, Z. Li and L. Wang, Multi-Sites Electrocatalysis in High-Entropy Alloys, Adv. Funct. Mater., 2021, 31, 2106715.
26. J.-T. Ren, L. Chen, H.-Y. Wang and Z.-Y. Yuan, High-entropy alloys in electrocatalysis: from fundamentals to applications, Chem. Soc. Rev., 2023, 52, 8319–8373.
27. D. Roy, S. C. Mandal and B. Pathak, Machine Learning Assisted Exploration of High Entropy Alloy-Based Catalysts for Selective CO2 Reduction to Methanol, J. Phys. Chem. Lett., 2022, 13, 5991–6002.
28. T. A. A. Batchelor, J. K. Pedersen, S. H. Winther, I. E. Castelli, K. W. Jacobsen and J. Rossmeisl, High-Entropy Alloys as a Discovery Platform for Electrocatalysis, Joule, 2019, 3, 834–845.
29. X. Wan, Z. Zhang, W. Yu, H. Niu, X. Wang and Y. Guo, Machine-learning-assisted discovery of highly efficient high-entropy alloy catalysts for the oxygen reduction reaction, Patterns, 2022, 3, 100553.
30. J. K. Pedersen, T. A. A. Batchelor, A. Bagger and J. Rossmeisl, High-Entropy Alloys as Catalysts for the CO2 and CO Reduction Reactions, ACS Catal., 2020, 10, 2169–2176.
31. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558–561.
32. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
33. E. J. Baerends, Perspective on “Self-consistent equations including exchange and correlation effects”, Theor. Chem. Acc., 2000, 103, 265–269.
34. G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15–50.
35. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
36. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, Analytic projection from plane-wave and PAW wavefunctions and application to chemical-bonding analysis in solids, J. Comput. Chem., 2013, 34, 2557–2567.
37. A. A. Peterson and J. K. Nørskov, Activity Descriptors for CO2 Electroreduction to Methane on Transition-Metal Catalysts, J. Phys. Chem. Lett., 2012, 3, 251–258.
38. Y. Ouyang, L. Shi, X. Bai, Q. Li and J. Wang, Breaking scaling relations for efficient CO2 electrochemical reduction through dual-atom catalysts, Chem. Sci., 2020, 11, 1807–1813.
39. J. Gong, S. Chu, R. K. Mehta and A. J. H. McGaughey, XGBoost model for electrocaloric temperature change prediction in ceramics, npj Comput. Mater., 2022, 8, 140.
40. S. Lu, Q. Zhou, Y. Guo and J. Wang, On-the-fly interpretable machine learning for rapid discovery of two-dimensional ferromagnets with high Curie temperature, Chem, 2022, 8, 769–783.
41. J. Sun, A. Chen, J. Guan, Y. Han, Y. Liu, X. Niu, M. He, L. Shi, J. Wang and X. Zhang, Interpretable Machine Learning-Assisted High-Throughput Screening for Understanding NRR Electrocatalyst Performance Modulation between Active Center and C-N Coordination, Energy & Environ. Mater., 2024, 7, e12693.
42. S. M. Lundberg, G. Erion, H. Chen, A. DeGrave, J. M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal and S.-I. Lee, From local explanations to global understanding with explainable AI for trees, Nat. Mach. Intell., 2020, 2, 56–67.
43. W. A. Saidi, W. Shadid and G. Veser, Optimization of High-Entropy Alloy Catalyst for Ammonia Decomposition and Ammonia Synthesis, J. Phys. Chem. Lett., 2021, 12, 5185–5192.
44. P. Xie, Y. Yao, Z. Huang, Z. Liu, J. Zhang, T. Li, G. Wang, R. Shahbazian-Yassar, L. Hu and C. Wang, Highly efficient decomposition of ammonia using high-entropy alloy catalysts, Nat. Commun., 2019, 10, 4011.
45. Y. Feng, W. An, Z. Wang, Y. Wang, Y. Men and Y. Du, Electrochemical CO2 Reduction Reaction on M@Cu(211) Bimetallic Single-Atom Surface Alloys: Mechanism, Kinetics, and Catalyst Screening, ACS Sustainable Chem. Eng., 2020, 8, 210–222.
46. W. Luo, X. Nie, M. J. Janik and A. Asthagiri, Facet Dependence of CO2 Reduction Paths on Cu Electrodes, ACS Catal., 2016, 6, 219–229.
47. L. Ou, W. Long, J. Huang, Y. Chen and J. Jin, Theoretical insight into effect of doping of transition metal M (M = Ni, Pd and Pt) on CO2 reduction pathways on Cu(111) and understanding of origin of electrocatalytic activity, RSC Adv., 2017, 7, 11938–11950.
48. S. Zhen, X. Lin, G. Zhang, D. Cheng, C. Jiang, X. Shi, S. Wu, Z.-J. Zhao and J. Gong, Accessing the Nature of Active Sites and Particle Size Effect for Reduction of Carbon Dioxide over Copper-Based Catalysts, J. Phys. Chem. C, 2023, 127, 4975–4983.
49. D. Roy, S. C. Mandal and B. Pathak, Machine Learning-Driven High-Throughput Screening of Alloy-Based Catalysts for Selective CO2 Hydrogenation to Methanol, ACS Appl. Mater. Interfaces, 2021, 13, 56151–56163.
50. T.-H. Hu, C.-Y. Wu, Z. Y. He, Y. Chen, L.-C. Hsu, C.-W. Pao, J.-T. Lin, C.-W. Chang, S.-C. Lin and R. Osmundsen, et al., Unconventional Hexagonal Close-Packed High-Entropy Alloy Surfaces Synergistically Accelerate Alkaline Hydrogen Evolution, Adv. Sci., 2025, 12, 2409023.
51. Z. L. Carroll, M. J. R. Haché, B. Wang, L. Chen, S. Wu, U. Erb, S. Thorpe and Y. Zou, Electrodeposited NiFeCoMoW High-Entropy Alloys with Nanoscale Amorphous Structure as Effective Hydrogen Evolution Electrocatalysts, ACS Appl. Energy Mater., 2024, 7, 8412–8422.
52. M. Wei, Y. Sun, J. Zhang, J. Zeng and J. Wang, Coupling High-Entropy Core with Rh Shell for Efficient pH-Universal Hydrogen Evolution, Small, 2024, 20, 2403353.
53. C.-Y. Wu, Y.-C. Hsiao, Y. Chen, K.-H. Lin, T.-J. Lee, C.-C. Chi, J.-T. Lin, L.-C. Hsu, H.-J. Tsai and J.-Q. Gao, et al., A catalyst family of high-entropy alloy atomic layers with square atomic arrangements comprising iron- and platinum-group metals, Sci. Adv., 2024, 10, eadl3693.
54. R. Nandan, H. Nara, H. N. Nam, Q. M. Phung, Q. P. Ngo, J. Na, J. Henzie and Y. Yamauchi, Tailored Design of Mesoporous Nanospheres with High Entropic Alloy Sites for Efficient Redox Electrocatalysis, Adv. Sci., 2024, 11, 2402518.
55. X. Cui, Y. Liu, X. Wang, X. Tian, Y. Wang, G. Zhang, T. Liu, J. Ding, W. Hu and Y. Chen, Rapid High-Temperature Liquid Shock Synthesis of High-Entropy Alloys for Hydrogen Evolution Reaction, ACS Nano, 2024, 18, 2948–2957.
56. Z. Liu, H. Li, C. Yang, M. Jiang, J. Zhang and C. Fu, High-Throughput Designed and Laser-Etched NiFeCrVTi High-Entropy Alloys with High Catalytic Activities and Corrosion Resistance for Hydrogen Evolution in Seawater, Small, 2024, 20, 2309078.

---