
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38478684
202320232
10.1073/pnas.2320232121
research-articleResearch ArticlechemChemistry410
Physical Sciences
Chemistry
Interpreting chemisorption strength with AutoML-based feature deletion experiments
Li Zhuo a b https://orcid.org/0000-0002-8830-9891

Zhao Changquan c 1
Wang Haikun d 1
Ding Yanqing e https://orcid.org/0009-0008-4503-1851

Chen Yechao d
Schwaller Philippe f g https://orcid.org/0000-0003-3046-6576

Yang Ke h
Hua Cheng cheng.hua@sjtu.edu.cn
d 2 https://orcid.org/0000-0002-1662-2424

He Yulian yulian.he@sjtu.edu.cn
a b 2
aUniversity of Michigan-Shanghai Jiao Tong University Joint Institute, Shanghai Jiao Tong University, Shanghai 200240, China
bSchool of Chemistry and Chemical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
cSchool of Mathematical Science, Shanghai Jiao Tong University, Shanghai 200240, China
dAntai College of Economics and Management, Shanghai Jiao Tong University, Shanghai 200240, China
eFu Foundation School of Engineering and Applied Science, Columbia University, New York, NY 10027
fLaboratory of Artificial Chemical Intelligence, Institut des Sciences et Ingénierie Chimiques, École Polytechnique Fédérale de Lausanne, Lausanne 1015, Switzerland
gNational Centre of Competence in Research Catalysis, École Polytechnique Fédérale de Lausanne, Lausanne 1015, Switzerland
hKey Laboratory of Advanced Energy Materials Chemistry, Nankai University, Tianjin 300071, China
2To whom correspondence may be addressed. Email: cheng.hua@sjtu.edu.cn or yulian.he@sjtu.edu.cn.
Edited by Catherine Murphy, University of Illinois at Urbana-Champaign, Urbana, IL; received November 17, 2023; accepted February 10, 2024

1C.Z. and H.W. contributed equally to this work.

13 3 2024
19 3 2024
13 9 2024
121 12 e232023212117 11 2023
10 2 2024
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

Explainable AI (XAI) algorithms are considered potentially paradigm-changing for the research of physical sciences, yet they often meet problems of inconsistency and stochasticity in practice, hence are limited in complicated problems. This work provides a different approach for knowledge extraction from a statistical perspective showcasing with a crucial question in catalytic science: determining which factor dominates the strength of chemisorption. With AutoML-based feature deletion experiments, the local structure of the adsorption site was revealed as the dominating factor of chemisorption strength on bimetallic alloy surface, showing AutoML as a reliable tool for extracting theoretically meaningful knowledge from large datasets.

The chemisorption energy of reactants on a catalyst surface, Eads, is among the most informative characteristics of understanding and pinpointing the optimal catalyst. The intrinsic complexity of catalyst surfaces and chemisorption reactions presents significant difficulties in identifying the pivotal physical quantities determining Eads. In response to this, the study proposes a methodology, the feature deletion experiment, based on Automatic Machine Learning (AutoML) for knowledge extraction from a high-throughput density functional theory (DFT) database. The study reveals that, for binary alloy surfaces, the local adsorption site geometric information is the primary physical quantity determining Eads, compared to the electronic and physiochemical properties of the catalyst alloys. By integrating the feature deletion experiment with instance-wise variable selection (INVASE), a neural network-based explainable AI (XAI) tool, we established the best-performing feature set containing 21 intrinsic, non-DFT computed properties, achieving an MAE of 0.23 eV across a periodic table-wide chemical space involving more than 1,600 types of alloys surfaces and 8,400 chemisorption reactions. This study demonstrates the stability, consistency, and potential of AutoML-based feature deletion experiment in developing concise, predictive, and theoretically meaningful models for complex chemical problems with minimal human intervention.

automatic machine learning
chemisorption
surface science
explainable AI
MOST | National Natural Science Foundation of China (NSFC) 501100001809 22202131 Cheng HuaYulian He Shanghai Science and Technology Development Foundation (Shanghai Science and Technology Foundation) 100012543 22YF1419400 Cheng HuaYulian He MOST | National Natural Science Foundation of China (NSFC) 501100001809 72301172 Cheng HuaYulian He Shanghai Science and Technology Development Foundation (Shanghai Science and Technology Foundation) 100012543 21YF1420200 Cheng HuaYulian He Shanghai Education Development Foundation (SEDF) 501100003024 22CGA12 Cheng Hua
==== Body
pmcThe chemisorption energy of an adsorbate onto a catalytic surface, Eads, is a universally important quantity that is extensively used to pinpoint the optimal catalyst recipe for a specific catalytic reaction, following the empirical Sabatier–Balandin volcano plot analysis established in the mid-twentieth century (1, 2) However, the ongoing challenge remains in identifying the pivotal properties of catalytic systems that dictate chemisorption strength (3). This issue arises from constraints in approximating realistic adsorption energies and the intricate nature of the chemisorption process itself. Over recent decades, quantum mechanical (QM) calculations, particularly density functional theory (DFT) methods, have showcased their potential in expediting catalyst design and mechanistic explorations through high-throughput screening (4–6). In contrast to extracting chemisorption behaviors solely from experimental catalytic processes, these QM calculations construct chemisorption reactions from the ground up. This approach directly ensures the uniformity and clarity of the Eads data, thereby serving as a potential playground for theoretical investigation into fundamental questions and offering valuable insights for catalyst design (7, 8). However, QM methods are limited by computational costs and often rely on chemists’ hard-coded instructions. As a result, high-throughput QM methods prioritize providing information over interpretation.

Due to the aforementioned challenges, researchers are increasingly turning to machine learning (ML) methods, particularly explainable artificial intelligence (XAI), to uncover new insights into catalysis. XAI models aim at replacing “black-box” models with “glass-box” models (Fig. 1A), from which knowledge is derived from understanding the interactions between features within the mathematical framework (3, 9). Yet, we want to point out several drawbacks of the “glass-box” approach. First, the “transparency” in XAI is typically achieved by restricting the algorithm to a relatively simple mathematical structure, such as a generalized additive model (GAM, see SI Appendix, Note 1) (10). This limitation on structural complexity may reduce the predictive capabilities of ML models for challenging problems. Notably, the predictive ability of ML methods to handle complex and scientifically unclear systems is often the primary reason for their adoption. Second, XAI explanations are model-based, which means that both the model and its corresponding explanations are subject to the stochastic nature of the ML algorithm. This includes factors such as the selection of training/testing. datasets and the stopping point of the training process (9). While stochasticity may be less of a concern for purely predictive tasks, it can lead to fundamentally different explanations from similarly predictive models, causing confusion when choosing between ambiguous results. Finally, XAI explanations are also feature-specific; they are inherently tied to the manner in which the descriptor is constructed. However, every chemical problem involves at least one real-world system, and there are infinite ways to create physically equivalent descriptors. Thus, in the rapidly evolving field of AI for chemistry, the model- and feature-specific explanations provided by XAI may struggle to meet the level of clarity and certainty that chemists require.

Fig. 1. Overview of this research. (A) The typical XAI approach is considered a “glass box.” Descriptors of the investigated phenomenon are input into the model, with knowledge being extracted from the interior of the “box.” The system’s behavior is predicted through an individual ML model. (B) This work aggregates numerous ML models within an AutoML framework and extracts knowledge from the statistical results of these AutoML outputs. Based on experimental design, sets of physical quantities (feature sets, F^) are grouped and processed by the AutoML frame in parallel. The system’s behavior is predicted via an optimal feature set Fopt, which is independent of the model used to establish it. An analogy for the difference between this approach and typical XAI is that of thermodynamic/kinetics vs. quantum mechanical tools such as DFT. The latter focuses on explaining and establishing physical insight into each microscopic system (molecule or mathematical structure), while the former focuses on measuring and explaining the overall behavior of macroscopic systems (thermodynamic functions or statistics of ML models).

Thus, we propose an alternative, Automatic Machine Learning (AutoML)-guided approach for knowledge extraction (Fig. 1B). Instead of delving into the inner workings of the ML algorithm, we bundle numerous comparable ML models together for collective analysis. Specifically, we establish physical insights based on a plain and fundamental principle not exclusive to ML methods. We posit that “critical” physical quantities should significantly influence the predictability of a physical model; removing these quantities would thus degrade model effectiveness, and vice versa. An initial baseline feature set (Ftotal, F for feature set) is constructed and verified to ensure its descriptiveness, i.e., models using this feature set should exhibit acceptable predictive performance. Subsequently, internally related features are removed from Ftotal to examine any changes in model predictability. Three benefits can be found in this approach. First, physical insights are gleaned by comparing the performance of different feature sets, thereby explicitly incorporating physical considerations. With a carefully designed experimental setup, changes in predictability can be linked to physical assumptions. Second, model stochasticity is mitigated by analyzing the statistics of comparable models. Last, this approach obviates the need to comprehend the detailed mathematical structure of the ML algorithm during the knowledge extraction process, thereby avoiding the trade-off between model complexity and explainability. This strategy shifts the source of scientific insight from elucidating model behavior to assessing feature set performance, aligning more closely with the traditional methodologies used to establish physical insights in experimental science.

Such an approach also raises the issue of training numerous comparable ML models on different datasets, a task that is impractical to achieve through manual parameterization of each model. Fortunately, recent advancements in AutoML have made the implementation of this method cost-effective in terms of computational resources. In this work, AutoML tools such as Autogluon have enabled the training of over 50,000 linked ML models without human intervention (11). Autogluon automatically parameterizes each model with Hyperband algorithms until optimal model performance is reached. In this work, Autogluon conducts exhaustive testing across a diverse range of models and configurations, encompassing neural networks, LightGBM and CatBoost boosted trees, Random Forests, and Extremely Randomized Trees, leading to the creation of thousands of distinct models with varying structures and hyperparameter settings. Furthermore, it leverages the parallel processing capabilities of multi-core CPUs and GPUs. This efficient use of hardware accelerates the evaluation and training of a wide array of models, significantly speeding up the machine-learning workflow. Each ML model is logically connected to a series of models by a deletion sequence, which implies a physical assumption of how different features affect the target property (detailed in Section 1.3 of Results and Discussion). In the following sections, we will demonstrate the construction of a descriptive baseline feature set (Ftotal) for the dissociative chemisorption process of diatomic molecules. Additionally, we will show how AutoML-based feature experiments reveal geometric features, specifically the local structure of the adsorption site, as crucial physical information governing the strength of the diatomic chemisorption process. Finally, by integrating AutoML-based feature experiments with Instance-wise Variable Selection (INVASE), a neural network-based XAI algorithm, we demonstrated the feasibility of simplifying a complex and realistic problem like the diatomic chemisorption process into a predictive model reliant on merely 11 non-DFT features. This showcases the potential of AutoML as a stable, consistent, and theoretically insightful tool for knowledge extraction.

1. Results and Discussion

1.1. Dataset Construction.

A high-throughput DFT-calculated dissociative chemisorption energy dataset was selected as the ground truth for this study. The data quality was verified by carefully reproducing the adsorption energy using the same DFT protocol suggested by Mamun et al. (12) The database contains DFT-calculated Eads values of various adsorbates on binary alloy surfaces, formed by 37 different metal elements from the d, ds, and p regions. The 37 selected metals form alloys with the L12 and L10 Strukturbericht designations, corresponding to ordered, homogeneous binary alloys with A3B and AB stoichiometries, respectively. Along with 37 pure metal surfaces in the A1 face-centered cubic structure, the stoichiometric B fraction of 0%, 25%, and 50% were sampled. The dataset excluded reactions that caused complete reconstructions of alloy surfaces, such as horizontal sliding or top layer dissociation, to ensure that the energy change reflects the dissociative adsorption of diatomic molecules on bimetallic surfaces. We subsequently refined the dataset from 88,587 entries, involving tens of different adsorbates in chemisorption reactions, to include only five diatomic molecular adsorbates (H2, O2, N2, CO, and NO, Table 1), resulting in 8,418 entries. The main reason of limiting the adsorbate to diatomic molecules was to minimize the complexity arising from adsorbate structure and to unify the adsorbate descriptors, allowing the ML models to focus on the surface behavior of the involved alloys (i.e., the catalysts).

Table 1. Diatomic chemisorption reactions investigated in this work

#	Reactions	
1	CO(g) + 2* →C* + O*	
2	NO(g) + 2* →N* + O*	
3	H2(g) + 2* → 2H*	
4	N2(g) + 2* → 2N*	
5	O2(g) + 2* → 2O*	

To identify which set of physical quantities is crucial for chemisorption strength, we first established a baseline feature set (Ftotal) using a “feature flooding” strategy, i.e., listing as many related features as possible. The overall workflow of this study is illustrated schematically in Fig. 2. Instead of extracting features from the DFT process, we chose to use the intrinsic properties of the involved elements and adsorbates, i.e., diatomic molecules and alloy surfaces (comprising a major alloy component A and a minor alloy component or ligand B). Intrinsic properties such as the electronic configuration of alloy elements and polarizability of diatomic adsorbates can be easily sourced from publicly available databases, thereby greatly reducing data collection costs. Using intrinsic features also significantly enhances the practicality of the optimized chemisorption model, as DFT-based features may not always be readily available for new systems. In catalytic science, three categories of physical properties are generally considered relevant to chemisorption strength, including physiochemical (describing the ensemble behavior of the involved species), electronic (electronic configuration of the involved element and d-band information), and geometric (detailing the adsorption site) information (13, 14). Fig. 2 also shows the structure of Ftotal. The construction of Ftotal can be visualized using the analogy of weaving a fabric. The three involved species are considered the “warps”(vertical threads of a fabric), and the three categories of physical properties are the “wefts”(horizontal threads of a fabric). Each feature is situated at the intersection of a “warp” and a “weft.” The following process of using feature deletion experiments to determine the relative importance of these physical properties can be envisioned as trimming the “threads” to assess whether the “fabric” remains intact.

Fig. 2. Workflow of this research. Ftotal, a predictive feature set for the diatomic chemisorption process, serves as the testing ground for identifying the most prominent physical quantities. Using feature deletion experiments, geometric information emerges as the critical features for predicting the chemisorption process. By combining greedy deletion experiments with INVASE, a neural network-based XAI method, a leaner and physically meaningful feature set, F21 was selected from Ftotal, which outperforms Ftotal across various ML algorithms. Examination of F21 through feature deletion experiments also highlights geometric information, specifically regarding the local adsorption site, as a prominent descriptor of the chemisorption process.

A total of 139 features were woven into the “fabric,” including 80 electronic features, 32 physiochemical features, and 27 geometric features. A detailed feature list is included in SI Appendix, Table S1. For electronic features, two main types of information are included: d-band information and electronic configuration. Please note that the d-band information of the pure metal was used instead of the d-band of the binary alloy’s surface, which would require costly DFT calculations, and as such, it is considered an intrinsic property of the element. d-band information, such as d-band center, d-band skewness, and d-band filling, are known to be crucial descriptors for metal adsorption behavior and, in some simpler cases, can exclusively determine adsorption strength (15–18). The spatial behaviors of d orbitals are also taken into consideration, characterized by the spatial extent of the d orbital (rd) and the coupling Hamiltonian matrix element Vad2. Vad2 is a function of rd and d-band information, and is an intrinsic property of a given element based on Muffin-Tin-Orbital theory (19). The electronic configurations of the alloy elements, including information on both inner and outer electron shells, are detailed and expressed as dummy variables (binary values indicating the absence or presence of categorical variables, see SI Appendix, Note 4) to fulfill the prerequisites of our INVASE algorithm (Section 4).

The site specificity of adsorption on heterogeneous surface is captured through geometric features. These features are formulated by counting the numbers of neighboring atoms within certain cutoff radii (20, 21). Importantly, all geometric features are extracted from the initial crystal structures by positioning the adsorbate’s central atom 1 Å away from the surface to simulate the adsorption process (SI Appendix, Fig. S1A), eliminating the need for costly DFT-relaxation calculations. As shown in SI Appendix, Fig. S1B, various cutoff radii were considered to account for both short- and long-range interactions at the adsorption site. These include 2 Å, 3 Å, 4 Å, 5 Å, and 6 Å, as well as a natural cutoff based on the adsorbate-metal bond length (the sum of the atomic radii of the adsorbate’s central atom and the adsorbent metal atom, typically between 2∼3 Å) (22). Generalized local coordination number and local electron affinity are also included in the geometric features, as they provide information about the bonding behavior around the adsorption site (Method for detail) (18, 23–25) The bulk nearest-neighbor distance of the element is also included in the geometric features, as it has been proposed to relate to the adsorption energy on bimetallic alloys (26).

Finally, intrinsic features regarding the physiochemical information of the system at both the atomic and molecular levels were also collected. These features include atomic number, atomic radius, Van Der Waals radius, electron affinity, Pauling electron negativity, dipole polarizability, gas phase heat of formation, first ionization energy, work function, thermal conductivity, and periodic information, such as atomic number and element class. These contribute additional information to the model. All data sources used can be found in the Materials and Methods.

1.2. Primary Model Selection.

The reliability of Ftotal was first evaluated using the AutoML method. To be considered valid, Ftotal must achieve an acceptable mean absolute error (MAE) to confirm that these features provide sufficient information for predicting the energy changes in the diatomic chemisorption process. In the meantime, an appropriate ML algorithm balancing model predictability and computational cost must be identified for subsequent data experiments. SI Appendix, Fig. S2 shows the results of model comparisons using Autogluon (version 0.7.0) on CPU (11). Among the models offered by Autogluon, three tree-based models, XGBoost, LightGBM, and CatBoost excelled, with average testing set MAEs of 0.29, 0.28, and 0.23 eV across five different training-testing sets, respectively. The rationale for averaging the MAE across various training-testing sets is to minimize uncertainty caused by different training/testing combinations and to better assess the predictability of the investigated feature set. Given that the uncertainty associated with the DFT method typically falls within a range of 0.1∼0.2 eV, it is reasonable to conclude that Ftotal contains sufficient information for constructing a predictive model for chemisorption energy (27). Among these tree-based models, CatBoost (MAE = 0.23 eV) and LightGBM (MAE = 0.28eV) outperform XGBoost (MAE = 0.29 eV). However, their training times are significantly longer (20 times and 2.5 times that of XGBoost, respectively). We observed that LightGBM models usually train faster than XGBoost models (especially on GPU). The time difference noted in this case may be related to the automatic hyperparameter tuning process of Autogluon. Given that the subsequent feature deletion experiments will involve training thousands of comparable models, XGBoost was chosen as the balanced option between model predictability and computational cost.

1.3. Examining Feature Importance with Feature Deletion Experiments.

With the as-constructed Ftotal and chosen ML algorithm, we proceed to examine the importance of different physical quantities. The examination is conducted using a custom-designed method called feature deletion experiment. Based on the fundamental assumption that removing crucial information will significantly impact model predictability, an idea inspired by Occam’s razor, we gradually remove selected feature sets from Ftotal. Each new feature set is evaluated by automatically running the XGBoost model 5 times, each time with a different training and testing set (SI Appendix, Fig. S3A). If a significant change is observed in the average model MAE on test set (i.e., impacting model predictability), the deleted features should be considered crucial for predicting chemisorption strength.

We began by validating the methodology of the proposed feature deletion experiment. It is evident that information from all three components in the chemisorption reaction, i.e., the adsorbate and alloy elements A and B, is crucial for predicting chemisorption strength. Removing this information will inevitably impact model predictability. SI Appendix, Fig. S4 shows the result of gradually removing the three “reactants” from Ftotal. For alloy components A and B (SI Appendix, Fig. S4 A and B), removing some of the features has a minor impact on predictability, likely due to the “flooding” strategy employed during the construction of Ftotal. Eliminating more alloy component descriptors from Ftotal results in a loss of model predictability and eventually leads to an unacceptable MAE (>0.4 eV) when one of the alloy components is entirely removed. For the adsorbate (SI Appendix, Fig. S4C), randomly removing 9 out of 10 adsorbate features has a minor impact on model predictability, while removing the last one causes the model to collapse entirely. The reason of only one adsorbate feature is needed is very likely due to the structure of Ftotal. As shown in Table 1, the reaction matrix of Ftotal consists of 5 adsorbate × 1,684 binary alloy, meaning each alloy of elements A and B corresponds to five different reactions. The huge difference between the dimension of adsorbates and alloys causes the ML algorithm to mainly recognize patterns from the alloy dimension instead of adsorbate dimension. In other words, adsorbate features served only as labels of different reactions on the same alloy surface, thus only one adsorbate feature is required in an optimized feature set, irrelevant of which adsorbate feature is passed. These results are in alignment with the fundamental assumption that gradually removing crucial information from feature sets will cause a collapse in model predictability, thereby validating the effectiveness of the feature deletion experiment for chemisorption problems.

Next, the feature deletion experiment is used to evaluate the three categories of physical quantities: namely, electronic, physiochemical, and geometric features. The results are shown in Fig. 3. Initially, it is observed that the removal of electronic features from Ftotal leads to only a slight decrease in predictability (ΔMAE ≈ 0.005 eV, Fig. 3A). Removing physiochemical features has a stronger, yet still minor, impact on model predictability (ΔMAE ≈ 0.01 eV, Fig. 3B). Interestingly, removing all geometric features from Ftotal results in a dramatic impact, causing the model to completely lose its predictability (ΔMAE ≈ 0.4 eV, Fig. 3C). To rule out the possibility that the models might confuse alloys with the same elements but different ratios after most geometric features are deleted, the A-to-B atomic ratio information was added explicitly to the baseline set. The result is shown in SI Appendix, Fig. S5 (SI Appendix, Note 2). Despite the ratio information mitigating the collapse of model predictability, it was found that ΔMAE still reached 0.22 eV, suggesting that geometric information holds special importance over the other two sets of physical quantities. It was also found that in both experiments (SI Appendix, Fig. S4 and Fig. 3), the trend of MAE change is independent from the deletion sequence, confirming the high stability of feature deletion experiment.

Fig. 3. Results of feature deletion and feature addition experiments. In the feature deletion experiments (A–C), the algorithm first generates a random sequence of the investigated features. It then sequentially removes each feature from Ftotal and assesses the predictability of the resulting feature set by automatically segregating training/testing sets and the training model five times. To mitigate the influence of any specific feature deletion sequence, the deletion experiment was repeated five times with different sequences to exclude the impact of deletion sequence. The lightly colored region represents the SD of the feature set MAE, arising from the random selection of training/testing sets. The result of the feature deletion experiment is summarized in (D) by averaging the results in Ftotal, showing that geometric features are crucial for accurately predicting chemisorption strength.

Integrating these results, we find that only geometric information is vital for accurately predicting chemisorption strength, as the omission of electronic or physiochemical features has a minor impact on model predictability. However, the model also fails to accurately predict chemisorption strength when only geometric information is provided, yielding an MAE ≈ 0.6 eV (SI Appendix, Fig. S6). Taken together, these findings suggest that an effective prediction of chemisorption strength requires geometric information in conjunction with at least one set of either electronic or physiochemical information. Identifying the governing properties of a catalyst surface related to chemisorption strength has remained one of the most important fundamental questions to be solved for decades. As summarized in Fig. 3D, the results strongly indicate that geometric information, specifically regarding the adsorption site, is overwhelmingly more important for the prediction of chemisorption strength than the other two categories. Electronic and physiochemical properties offer similar insights but from different perspectives, so only one set is needed to complement geometric information for a sufficient model. However, the value of physiochemical and electronic features shows no apparent correlation, as the Pearson correlation coefficient shows little interdependence between the two sets (SI Appendix, Fig. S7). Both electronic and physiochemical information describe the electronic energy level of the involved metal species, thus are implicitly related and somewhat interchangeable. This is consistent with previous statistical analyses using Kullback–Leibler divergence and DFT-derived features, which have also revealed the importance of site geometric information for predicting chemisorption strength in bimetallic alloy systems (20).

Since the chemisorption dataset is based solely on high-throughput calculation, it is important to confirm that the knowledge generation algorithm is robust while polluted data are involved. Given that data pollution can take many forms, we propose that it can be described by a 2 × 2 matrix (Table 2). The (instance-systematic) dimension indicates that the pollution could occur in individual instances and/or across the entire dataset, and the (feature-label) dimension means that the pollution could affect both the feature and the label. Following this scheme, we deliberately corrupted the dataset by garbling the features of some instances (Fig. 4A), randomizing the Eads of some instances (Fig. 4B), adding random error to each continuous features of the whole dataset (Fig. 4C) and adding random error to Eads of the whole dataset (Fig. 4D), then applied the feature deletion experiment to examine whether the knowledge remains stable. It is observed that the predictability of Ftotal declines as more pollution is added to the dataset, yet the importance rank remains Geometric > Physiochemical ≈ Electronic in all four cases, showing that feature deletion experiment is stable against significant perturbations in the dataset. This stability makes it suitable for future applications in experiment-based datasets, in which both occasional mislabeling and experimental systematic error must be considered.

Table 2. Forms of data pollution

	Feature (x)	Label (y)	
Instance	Feature error	Label error	
Label	Feature distribution	Label distribution	

Fig. 4. Robustness of the feature deletion experiment for knowledge generation. Average MAE of 5 deletion sequence is presented in (A) 1%/3% instances have wrong features, such as mistakenly marking NO → N* + O* on Ti3Fe alloy as CO → N* + O* on Ti3Fe, or as NO → N* + O* on Ti3Ni; (B) 1%/3% instances have randomized Eads label, where the incorrect Eads is kept within the range of original Eads to avoid significant influence on the average Eads of the dataset; (C) 0.5%/1.5% random error, relative to the range of each feature, added to every continuous feature in the dataset; (D) 0.5%/1.5% random error, relative to the range of Eads, added to every Eads in the dataset.

As geometric features (FGeo) appear to be the most crucial category, an in-depth investigation into these features is naturally warranted. The investigation started by separating the short-range features (Fs, features with ≤ 3 Å cutoff) from the long-range features (Fl, features with > 3 Å cutoff). The relative importance ranking of features within Fi (the investigated feature set) can be determined using a recursive greedy algorithm (SI Appendix, Fig. S3B). This algorithm tentatively removes each feature from Fi within Ftotal and identifies the feature whose removal results in either the most improvement or the least detriment. The identified feature is then removed from Ftotal, and the algorithm recurs with Ftotal−1 until only one feature from Fi remains. This reverse deletion sequence can be interpreted as the relative importance ranking of features within Fi.

The same greedy deletion algorithm was applied to both Fs and Fl, and the results are shown in Table 3 and SI Appendix, Fig. S8. For Fs, the relative importance ranking strongly favors natural cutoff properties based on covalent radii. We posit that such a natural cutoff is more effective for describing the first-nearest-neighbor interactions across vastly different alloy structures compared to an arbitrary fixed distance (such as 2Å and 3 Å). In addition, the information related to host element A seems to weigh more than that of ligand element B for Fs, while the opposite is true for Fl. A thorough inspection of FGeo revealed that short-range cutoff radii cannot effectively capture the coordination information about neighbor B, which has relatively large atomic radii (i.e., most neighbor_B_x with x < 4 are equal to 0). Consequently, neighbor B information on a longer-range scale (> 3 Å) is required to provide valid information on local B coordination. This analysis further suggested 4 Å as an optimal range for including effective B coordination information among the long-range interactions considered.

Table 3. Relative importance ranking of features in Fs and Fl

Importance ranking	Fs (short range)	Fl (long range)	
1	neighbor_A_n	neighbor_B_4	
2	neighbor_B_n	neighbor_A_6	
3	local_electron_affinity_n	coordination_number_4	
4	bulknn	local_electron_affinity_6	
5	neighbor_B_3	neighbor_B_5	
6	local_electron_affinity_3	neighbor_B_6	
7	neighbor_A_2	neighbor_A_5	
8	neighbor_A_3	coordination_number_5	
9	local_electron_affinity_2	coordination_number_6	
10	coordination_number_2	neighbor_A_4	
11	coordination_number_n	local_electron_affinity_4	
12	neighbor_B_2	local_electron_affinity_5	
13	coordination_number_3	

Overall, the results from the local geometric feature analysis on short- and long-range scales suggest that the coordination states of elements A and B are crucial for accurately predicting chemisorption strength. Specifically, the first-nearest neighbor interactions between the adsorbate central atom and the host element A, as captured by neighbor_A_n, together with the neighbor information of ligand B on a longer-range scale (neighbor_B_4), are identified as the most effective descriptors compared to others defined by various cutoff radii.

1.4. Establishing Chemically Meaningful Feature Set.

Feature deletion experiment not only revealed the importance of geometric information over the other two sets, but also suggested that not all Ftotal features are needed for a predictive chemisorption model. It is intuitive to try to establish a lean, predictive, and chemically meaningful model from Ftotal using the deletion strategy. However, to down select such a chemically meaningful model, one must first examine the uniformity of the problem, i.e., if all included instances can be described by the same set of physical quantities. Thus, an XAI model with the ability of identify feature importance on an instance-by-instance basis is needed. A neural-network-based XAI technique, Instance-wise Variable Selection (INVASE), was first introduced into catalytic research in this work. INVASE is designed to discern important features from irrelevant ones on an instance-wise basis (28). A schematic illustration of the INVASE algorithm used in this work is shown in Fig. 5A. Compared with vanilla INVASE algorithm, we modified the loss function to better fit the requirements of chemistry problems, and the detailed modifications are described in SI Appendix, Note 3.

Fig. 5. INVASE results. (A) Schematic illustration of the modified INVASE algorithm used in this work. xi are the components of the feature vector, pi are the components of the importance vector, Si are the components of the selection vector generated by the random sampler, y~P is the prediction output from the predictor network, and y~BL is the prediction output from the baseline network; (B) Importance matrix of Ftotal, generated by the modified INVASE algorithm. All reaction instances can be predicted by the same set of features, suggesting high uniformity of the problem. This uniformity paves the way for subsequent deletion experiments to develop a more streamlined physical description.

Using the modified INVASE algorithm, irrelevant features in Ftotal were also identified. The baseline network was first established with MAE = 0.34 eV. After 20,000 training iterations, both the selector and predictor networks individually achieved MAE = 0.31 eV, confirming that removing some of the features can indeed enhance predictability. The resulting importance matrix is presented in Fig. 5B. The first characteristic of the matrix is its instance-wise uniformity across reaction instances; a feature that is useful for predicting any single instance is also useful for the entire dataset, and vice versa. Such uniformity is rarely observed in ML datasets. For comparison, SI Appendix, Fig. S10 shows the importance matrix for three medical and economic datasets using the modified INVASE algorithm; none display such uniformity. Nonuniform importance matrixes are also typically found in previous works using vanilla INVASE (29). The instance-wise uniformity confirms that all reactions included in this study can be predicted with the same feature set. The importance of using high-throughput DFT data for establishing physically meaningful explanations for the chemisorption process is also highlighted by INVASE, as mentioned in the introduction. Fig. 5B shows that at least 28 features can be excluded from Ftotal without adversely affecting the predictability of the neural network, with 111 features remaining (F111, subscript stands for number of features in the set). The full list of selected/excluded features is given in SI Appendix, Table S2. Using the percentage of features retained in F111 as a qualitative measurement of the relative importance of the three categories of quantities, geometric features (88.9%) again prove most important, followed by electronic features (81.2%) and physiochemical features (68.8%). Note that all non-ordered features were converted into dummy variables due to the feature formal rationality in INVASE (SI Appendix, Note 4).

The above results reaffirm that a streamlined feature set, with non-essential features removed, can indeed improve the model predictability due to the variance-bias tradeoff (30). This tradeoff implies that each feature added to the set contributes both new information and high-dimensional noise. By eliminating non-critical information, the model’s generalizability may improve through noise reduction. Moreover, the removal of less informative features allows for clearer scientific insights.

Identifying the globally optimal feature set among all possible combinations of the 139 features (Σi=1139C139i=2139−1≈7×1041) is practically impossible. However, the same recursive greedy deletion algorithm, used to establish the relative importance ranking for geometric features, can prune less important features from the INVASE-selected feature set F111. The results are shown in Fig. 6. Removing the first 70 features only slightly improves the model MAE from Ftotal. However, starting from F35, the MAE significantly decreases as more features are removed. At F21, a global minimum is reached, with MAEXGB = 0.259 eV (SI Appendix, Fig. S11A), which is a considerable improvement from Ftotal of MAE = 0.29 eV. To validate the sufficiency of information in F21, this feature set was tested with CatBoost and LightGBM models. The results revealed that F21 has comparable or even better performance than Ftotal (MAECAT = 0.23 eV for both sets in SI Appendix, Fig. S11B, MAELGBM = 0.26 eV for F21 and 0.28 eV for Ftotal). This confirms that F21 retains all essential information from Ftotal while substantially reducing high dimensional noise. The Bayesian information criterion (BIC, Fig. 6B) for XGBoost models applied to sets F111 to F21 also suggests that F21 is among the most balanced feature sets between maximum likelihood and overfitting. This shows that the feature deletion experiment is able to generate a more efficient and optimized feature set.

Fig. 6. Greedy scanning of chemically meaningful feature set. (A) Average and SD of the performance of feature sets. (B) Bayesian information criterion (BIC) of the models on the corresponding feature set. BIC awards the model for accuracy and punishes it for increased complexity, and is defined as BIC=nln(SSE)+ln(n)p, where n is the number of data points in the training set, SSE is sum of squared error, and p is the number of features used in the model. (C) Use XGBoost to measure the predictability of 5,000 randomly selected feature sets, each containing 20 alloy features + 1 adsorbate feature.

To validate the efficacy of feature deletion experiment in identifying the optimal feature set, and to demonstrate that F21 is fundamentally different from a random selection, 5,000 feature sets, each containing 20 alloy features plus one adsorbate feature, were randomly generated, and their predictability was evaluated using XGBoost. The results are shown in Fig. 6C. Using the performance of F21 (0.259 eV) and Ftotal (0.29 eV) as benchmarks, none out of the 5,000 sets outperform F21, and only 3.1% of the randomly selected feature sets outperform Ftotal. These results prove that the predictability improvement of F21 is not only due to the reduction of high dimensional noise but also the right combination of physical quantities. The following analysis will show that F21 provides clear chemical indications.

Table 4 shows the details of F21, containing 1 adsorbate feature, 3 geometric features, seven physiochemical features, and 10 electronic features. Only three geometric features that describe the number of A atoms in the natural cutoff range (neighbor_A_n), the number of B atoms in 4 Å cutoff range (neighbor_B_4), and the number of A atoms in 5 Å cutoff range (neighbor_A_5) remain in F21. Two of the selected features in F21 (neighbor_A_n and neighbor_B_4) are also ranked first in their own subsets Fs and Fl, respectively, as discussed above, suggesting that the feature deletion method has high internal consistency.

Table 4. Details of the best performing feature set from greedy scanning

Physiochemical	Geometric	Electronic	
Polarizability	neighbor_A_n	rd_A	
atomic_number_A	neighbor_B_4	db_filling_B	
electron_affinity_A	neighbor_A_5	inner_4d_A_10	
work_function_A		inner_4f_A_0	
273K_thermal_conductivity_A		inner_4f_A_0	
conductivity_A		outer_d_A_10	
atomic_number_B		outer_d_A_8	
atomic_radius_B		outer_d_A_5	
class_transition_B		inner_5s_B_0	
		inner_5s_B_2	
		inner_6s_B_0	

Using the same feature deletion method applied to Ftotal, the relative importance of geometric, physiochemical, and electronic features was determined for F21. The results are shown in Fig. 7. As shown in Fig. 7A, removing electronic features from F21 led to ΔMAE ≈ 0.04 eV, resulting in MAE = 0.30 eV, which is comparable to that of Ftotal. Given that the range of chemisorption energy involved in this study is (−12.9, 11.2) eV, this suggests that an effective minimum feature set could consist of as few as 11 intrinsic features, including one adsorbate feature, three adsorption site geometric features, and seven physiochemical features, without any electronic quantum mechanical information.

Fig. 7. Detailed analysis of F21. (A)–(C) Feature deletion experiments of electronic, geometric, and alloy physiochemical features from F21. (D) comparison between the three categories of physical quantities.

Similar to Ftotal, geometric information also plays the most critical role in F21, despite only 3 being selected, as Fig. 7B shows a ΔMAE of ≈0.4 eV. Fig. 7C indicates that removing alloy physiochemical information from F21 has a stronger impact (ΔMAE ≈ 0.15 eV) than that of electronic features. A specific feature of alloy component B, atomic_radius_B, was found to be especially important. A ΔMAE ≈ 0.1 eV was observed when atomic_radius_B was removed, irrespective of the deletion sequence. The importance of atomic_radius_B is likely associated with the well-known “ligand” and/or “strain” effects in bimetallic nanocrystals. Introducing a second metal B into the host metal matrix A could induce notable variations in electronic states (i.e.,d-band center, d-band width) and/or lattice strain (compressive or tensile), thereby affecting chemisorption strength (31). The size of the ligand atom will not only determine the bond distance between the ligand and host atom, but also induce lattice strain, resulting in different degrees of d-orbital overlapping and thus different d-band width and d-band center. Similar findings were reported in a recent XAI work by Linic et al. using generalized additive models (iGAM), where the size of ligand metal was also identified as one of the most crucial features for describing the chemisorption strength (10). As summarized in Fig. 7D, the relative importance ranking was found to be geometric > physiochemical > electronic on F21, which is consistent with the findings from Ftotal.

The eight physiochemical features in F21 highlight the importance of the ensemble behavior of electrons (electron_affinity_A, work_function_A, 273K_thermal_conductivity_A), which correlates well with the electronic state distributions, as well as the periodic identity of the alloy elements (atomic_number_A, atomic_number_B) that provide effective and unique labels for each element. Interestingly, the physiochemical properties of host metal A are found to be more important than those of ligand metal B. This is reasonable, as only L12 and L10 bimetallic crystal structures (corresponding A3B and AB stoichiometries, respectively) were included in the investigated dataset. The atomic concentration of host metal A is always greater than, or equal to, that of the ligand. Only one adsorbate information, polarizability, remains in F21, ensuring the model distinguishes between different reactions, as discussed previously in SI Appendix, Fig. S4C.

A total of 10 electronic features were selected in F21. As previously mentioned, these electronic features only had a minor impact on the overall prediction of chemisorption strength. In this sense, we speculate that the primary purpose of involving these electronic features is to effectively categorize elements into several subgroups, thereby providing additional labels to the catalytic systems, as discussed further below.

Two d-band behaviors, rd_A and db_filling_B are retained in F21, and they are believed to coordinate with other quantities from the physiochemical category to set the energy level of the alloy surface. The rest of the electronic features chosen by F21 are ground-state electronic configuration of A and B, and was found quite interesting. F21 only contains information about d, f, and s electrons, which are also the most influential electronic shells in catalytic chemistry. Furthermore, F21 specifically asks for information about fully filled orbitals (inner_4d_A_10, outer_d_A_10, inner_5s_B_2) and completely empty orbitals (inner_4f_A_0, inner_5s_B_0, inner_6s_B_0). Since the dataset used in this work contains no f1∼13 elements, these features can be expressed as lines partitioning the periodic table (Fig. 8). A particularly intriguing selection of electronic features can be found in Table 4 (outer_d_A_5, outer_d_A_8, outer_d_A_10). These features mark three special d electron configurations of A, in specific, d5 (half-filled), d8, and d10 (fully filled), respectively. Compared with d5 and d10, which are known to have special importance in the periodic table, the selection of outer_d_A_8 seems to be awkward. We suggest that the selection of outer_d_A_8 feature is strongly related to Pt family, as d8 configuration is very special within the periodic table and is a special “character” of the Pt family elements. The “electron filling order,” i.e., electron energy diagram for ground-state atoms (SI Appendix, Fig. S12) is a fundamental knowledge of chemistry. The d-related exceptions only happen at d5, d8, and d10, which were unexpectedly chosen by F21. The “relativistic effect,” i.e., abnormality of electron configuration among heavy atoms, is clearly centered around d8 configuration, or Pt family (Rh, Pd, and Pt). Thus, the inclusion of outer_d_A_8 may be interpreted as the algorithm recognizing Pt family elements, of which their special catalytic behavior has been intensively investigated for centuries.

Fig. 8. Chemical implication of F21. Each electronic dummy variable (SI Appendix, Note 3) included in F21 can be used to partition involved section of the periodic table. inner_5s_B_0 and inner_5s_B_2 separate the 4th and 5th/6th period (only 5th/6th period elements have 5s electrons); inner_4d_A_10 separate the 6th d-block elements and ds-block/p-block (only 6th elements have 4d electrons as inner d electrons); inner_4f_A_0 and inner_6s_B_2 separate the 5th and 6th period (since La is the only La family element included in this work, f electron is either 0 or 14 in Ftotal, thus f electron information is related to 6th period elements). The outer_d_A_5,8,10 features specifically ask for elements with d5, d8, and d10 configurations, which are related to VIB elements, IB∼VA elements, and Pt family elements. These d electron numbers are directly linked to the exception of the electron filling order (SI Appendix, Fig. S12 and inlet).

The results shown above again explicitly highlight the prime importance of adsorption site geometric information over both electronic and physiochemical properties. This also explains why tree-based models such as XGBoost, CatBoost, LightGBM, etc. outperform more complex Neural Network models in this context, as the best feature set F21 implies a decision tree to group elements based on electronic configuration information to refine the energy level calculated from geometric and physiochemical features. Collectively, these findings strongly suggest that the feature deletion experiment is not only a method of identifying the optimal feature set for predicting physical quantities but also a powerful tool to identify the most chemically relevant features to construct simplified models with theoretical significance.

2. Conclusion

In conclusion, through customized AutoML-based feature deletion experiments, this study identified local geometric information of the adsorption site as the primary physical quantity for determining chemisorption strength, in comparison with electronic and physiochemical properties of the alloy surface. Specifically, factors like the number of major alloy component atoms within the natural cutoff radius, the number of minor alloy component atoms within the 4 Å range, and the size of the minor alloy component atom (ligand metal) emerged as the most potent predictors of chemisorption strength. Although less pronounced, both electronic and physiochemical properties still contribute to fine-tuning the overall prediction, likely by partitioning the elements into finer subgroups. These results again pointed out the importance of critical descriptors in the prediction of catalyst surface behaviors, as previous works have demonstrated that quantitative structure-adsorption property relationship between adsorption energy and generalized coordination number can be used to directly guide the design of Pt nanoparticle structures in electrocatalysis (32, 33).

More importantly, we introduced a methodology for machine learning-assisted knowledge extraction. As shown in Fig. 1, the explanations or knowledge generate based on the statistical performance of automatically created, comparable ML models. This ensures that the resulting insights are independent of the specific mathematical structures of the employed ML algorithms. No human intervention, such as hyperparameter tuning, is required during this knowledge extraction process. The effectiveness of this methodology is further evidenced by the fact that F21 outperforms Ftotal not just in XGBoost but also in CatBoost and LightGBM algorithms, even though these latter two were not involved in the selection of F21.

Compared to traditional explainable AI approaches, which seek to replace the “Black Box” with a “Glass Box” and elucidate what happens within the Box, as seen in Fig. 1A, the AutoML-based feature deletion experiment avoids delving into the internal workings of the ML algorithm. Instead, it extracts physically meaningful knowledge from the statistical outputs of the Box. This approach parallels real-world chemical experiments, depicted in Fig. 1B, where the actual microscopic processes within a reactor remain elusive, yet the statistical outcomes, i.e., kinetic and thermodynamic phenomena, can be easily observed and understood.

Furthermore, the consistency in the results obtained from investigating Ftotal, FGeo, and F21 through the feature deletion experiment underscores the robustness of the proposed method, irrespective of the specific ML algorithm within the AutoML framework. We acknowledge that not every feature selection can be perfectly rationalized, and it may never be, given the natural limitations of the dataset and the intrinsic complexity of the problem (SI Appendix, Note 5). However, the proposed AutoML-based feature analysis methodology is a robust and adaptable tool for revealing the statistical feature importance in complex physical sciences, even beyond catalysis.

3. Materials and Methods

3.1. Dataset Construction.

The adsorption energy dataset used in this work was obtained from the open-source dataset MamunHighT2019 at Catalysis-Hub.org via Graphql API (https://www.catalysis-hub.org/) (12). The dataset contains 88,587 optimized adsorption entries calculated by DFT for small molecule adsorption on 1,998 bimetallic alloys with L12 and L10 as Strukturbericht designation, as well as 37 pure metal surfaces in the A1 face-centered cubic structure, with sampled stoichiometric A:B ratios of 0%, 25%, 50%, 75%, and 100%. The dataset also considered all unique adsorption sites, including top, bridge, and hollows. The data quality was carefully verified before the experimentation, following the same DFT calculation protocols as suggested by Mamun et al., and the absolute errors between our evaluations and the dataset were found on the order of 10−3 eV.

The full list of Ftotal is provided in SI Appendix. Electronic, geometric, and physiochemical features of alloy surfaces are collected from open-source databases and references (15–26, 34). The adsorbate features were collected from Computational Chemistry Comparison and Benchmark Database (CCCBDB) (34). Local coordination numbers (local CN) and local electron affinities (local EA) of the adsorption sites are calculated using local connectivity matrix from the ASE library:(24)[1] local  CN=ΣiΣjIi,jΣi,

where i is the total number of atoms within the cutoff radius and j is the number of bonded atoms, Ii,j indicates whether a bond is formed between atoms i and j (1 for bonding and 0 for nonbonding).[2] local  EA=∏jEea,j1/N,

where j is the index of atoms within the cutoff radius, N is the total number of atoms within the cutoff radius and Eea,j is the electron affinity of atom j.

Note that there exist other definitions of local electron affinity. In this work, we used the geometric average of the atoms involved at the adsorption site for its simplicity and clear physical implications (35).

3.2. AutoML Feature Experiments.

AutoML feature experiments were carried out using Autogluon-Tabular (version 0.7.0, distributed Mar. 2023) (11). Autogluon-Tabular is an open-source AutoML framework for tabular datasets, which automates the process of applying ML models to real-world tasks with minimal need for fine-tuning. ML models are automatically parameterized using Hyperband algorithms until an internal criterion is met so that no human intervention is needed while training a large number of models. A 7:3 training/testing split was used for all AutoML training in this work. Validation set (1/10 of the training set) was automatically reserved from training set by Autogluon during each training process. Model MAE on the test set was used to evaluate the performance of the feature set. For neural network algorithms involved in Autogluon experiments, three hidden layers with 200 neurons were used.

Two data experiment algorithms were used in this work, named feature deletion experiment and greedy deletion experiment, respectively. The process of the feature deletion experiment is described in SI Appendix, Fig. S3A. The algorithm takes in two feature sets, Fbase (baseline feature set for the following analysis) and Fi (investigated feature set), where Fi is a subset of Fbase. The algorithm continues until Fi is empty. If not, a randomly selected feature f∈Fi is removed from Fbase and Fi to form Fbase−1 and Fi−1. The Autogluon module of the algorithm then measures and records the predictability of Fbase−1 by averaging the MAEs of Fbase−1 with 5 different training/testing (7:3) sets. Then Fbase−1 and Fi−1 are set as the new Fbase and Fi, and the algorithm recurs with the updated Fbase and Fi.

The process of the greedy deletion experiment is described in SI Appendix, Fig. S3B. The algorithm takes in two feature sets, Fbase (baseline feature set) and Fi (investigated feature set), where Fi is a subset of Fbase. The algorithm breaks if Fi contains only one element. If not, the algorithm tentatively removes each fi∈Fi and records the predictability of (Fbase - fi) for each i. Then the algorithm finds and records the least damaging/most improved fmax using the same Autogluon module with feature deletion experiment and removes it from Fbase and Fi, and set (Fbase - fmax) as new Fbase and (Fi - fmax) as new Fi. The algorithm then recurs with the updated Fbase and Fi.

3.3. INVASE Experiments.

INVASE experiments were carried out using a modified INVASE algorithm (11). The selector, predictor, and baseline network each has three hidden layers, with 200 neurons in each layer. The expression of the loss function of the selector network, Lselector, was modified for chemistry-related problems as follows:[3] Lselector=R∑iSipi+λ1∑ipi+λ2∑ipi(1−pi),

where R is the difference in loss between the predictor and baseline networks (Lpredictor−Lbaseline), Si is the i-th component of the selection vector, pi is the ith component of the importance vector, and λ1 and λ2 are learning rates. The selector loss function penalizes the selection network for 1) causing the predictor network to lose predictability, 2) selecting too many features, and 3) having an importance component pi close to 0.5. Thus, the importance of each feature will converge to either 0 or 1 during the training process if possible.

Supplementary Material

Appendix 01 (PDF)

Y.H. would like to thank the National Natural Science Foundation of China (22202131) and the Shanghai Science and Technology Development Funds of the “Rising Star” Sailing Program (22YF1419400) for the financial support. C.H. would like to thank the National Natural Science Foundation of China (72301172 and 72394370:72394375), the Shanghai Science and Technology Development Funds of the Sailing Program (21YF1420200), and the Shanghai Education Commission and Shanghai Education Development Foundation Chenguang Program (22CGA12) for the financial support.

Author contributions

Z.L. and Y.H. designed research; Z.L., C.Z., H.W., Y.D., and C.H. performed research; Z.L., H.W., Y.C., P.S., K.Y., and C.H. contributed new reagents/analytic tools; Z.L. and C.Z. analyzed data; Y.H. supervision of the research; and Z.L. wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

All study data are included in the article and/or SI Appendix. The chemisorption dataset, algorithms and computer code used in this work are available for download from a public repository (36).

Supporting Information

This article is a PNAS Direct Submission.
==== Refs
1 J. K. Nørskov , Universality in heterogeneous catalysis. J. Catal. 209 , 275–278 (2002).
2 A. Vijh, Volcano relationships in catalytic reactions on oxides. J. Catal. 33 , 385–391 (1974).
3 J. A. Esterhuizen, B. R. Goldsmith, S. Linic, Interpretable machine learning for knowledge generation in heterogeneous catalysis. Nat. Catal. 5 , 175–184 (2022).
4 B. Hammer, J. K. Nørskov, Theoretical Surface Science and Catalysis-Calculations and Concepts (Academic Press, 2000), vol. 45, pp. 71–129.
5 A. H. Motagamwala, J. A. Dumesic, Microkinetic modeling: A tool for rational catalyst design. Chem. Rev. 121 , 1049–1076 (2021).33205961
6 B. W. J. Chen, L. Xu, M. Mavrikakis, Computational methods in heterogeneous catalysis. Chem. Rev. 121 , 1007–1048 (2021).33350813
7 A. A. Latimer, A. Kakekhani, A. R. Kulkarni, J. K. Nørskov, Direct methane to methanol: The selectivity-conversion limit and design strategies. ACS Catal. 8 , 6894–6907 (2018).
8 T. Z. Gani, H. J. Kulik, Understanding and breaking scaling relations in single-site catalysis: Methane to methanol conversion by FeIV=O. ACS Catal. 8 , 975–986 (2018).
9 H. Wang , Scientific discovery in the age of artificial intelligence. Nature 620 , 47–60 (2023).37532811
10 J. A. Esterhuizen, B. R. Goldsmith, S. Linic, Theory-guided machine learning finds geometric structure-property relationships for chemisorption on subsurface alloys. Chem 6 , 3100–3117 (2020).
11 N. Erickson et al., Autogluon-tabular: Robust and accurate automl for structured data. arXiv [Preprint] (2020). https://arxiv.org/abs/2003.06505 (Accessed 22 August 2023).
12 O. Mamun, K. T. Winther, J. R. Boes, T. Bligaard, High-throughput calculations of catalytic properties of bimetallic alloy surfaces. Sci. Data 6 , 76 (2019).31138814
13 B. R. Goldsmith, J. Esterhuizen, J. Liu, C. J. Bartel, C. Sutton, Machine learning for heterogeneous catalyst design and discovery. AIChE J. 64 , 2311–2323 (2018).
14 H. Xin, A. Holewinski, S. Linic, Predictive structure-reactivity models for rapid screening of pt-based multimetallic electrocatalysts for the oxygen reduction reaction. ACS Catal. 2 , 12–16 (2011).
15 B. Hammer, J. K. Norskov, Theoretical surface science and catalysis - calculations and concepts. Adv. Catal. 45 , 71–129 (2000).
16 I. Takigawa, K. Shimizu, K. Tsuda, S. Takakusagi, Machine-learning prediction of the d-band center for metals and bimetals. RSC Adv. 6, 52587–52595 (2016).
17 A. Vojvodic, J. K. Nørskov, F. Abild-Pedersen, Electronic structure effects in transition metal surface chemistry. Top. Catal. 57 , 25–32 (2013).
18 M. H. Hansen et al., An atomistic machine learning package for surface science and catalysis. arXiv [Preprint] (2019). https://arxiv.org/abs/1904.00904 (Accessed 22 August 2023).
19 W. Gao , Determining the adsorption energies of small molecules with the intrinsic properties of adsorbates and substrates. Nat. Commun. 11 , 1196 (2020).32139675
20 X. Liu, C. Cai, W. Zhao, H. J. Peng, T. Wang, Machine learning-assisted screening of stepped alloy surfaces for C1 catalysis. ACS Catal. 12 , 4252–4260 (2022).
21 Z. W. Ulissi , Machine-learning methods enable exhaustive searches for active bimetallic facets and reveal active site motifs for CO2 reduction. ACS Catal. 7 , 6600–6608 (2017).
22 A. Obeidat, A. Jaradat, B. Hamdan, H. Abu-Ghazleh, Effect of cutoff radius, long range interaction and temperature controller on thermodynamic properties of fluids: Methanol as an example. Phys. A: Stat. Mech. Appl. 496 , 243–248 (2018).
23 O. Mamun, K. T. Winther, J. R. Boes, T. Bligaard, A Bayesian framework for adsorption energy prediction on bimetallic alloy catalysts. npj Comput. Mater. 6 , 177 (2020).
24 A. Hjorth Larsen , The atomic simulation environment-a Python library for working with atoms. J. Phys.: Condens. Matter 29 , 273002 (2017).28323250
25 Z. Li, X. Ma, H. Xin, Feature engineering of machine-learning chemisorption models for catalyst design. Catal. Today 280 , 232–238 (2017).
26 N. Wu, Feature analysis: Revealing the importance of bulk nearest-neighbor distance feature for CO adsorption energy on bimetallic alloys. J. Phys. Chem. C 125 , 19268–19275 (2021).
27 J. A. Keith , Combining machine learning and computational chemistry for predictive insights into chemical systems. Chem. Rev. 121 , 9816–9872 (2021).34232033
28 J. Yoon, J. Jordon, M. van der Schaar, “INVASE: Instance-wise variable selection using neural networks” in 7th International Conference on Learning Representations, ICLR 2019 (Open Review, New Orleans, LA, USA, 2019).
29 S. Letzgus , Toward explainable artificial intelligence for regression models: A methodological perspective. IEEE Sig. Process. Magaz. 39 , 40–58 (2022).
30 G. James, D. Witten, T. Hastie, R. Tibshirani, An Introduction to Statistical Learning (Springer, New York, 2013).
31 C. Li, S. Yan, J. Fang, Construction of lattice strain in bimetallic nanostructures and its effectiveness in electrochemical applications. Small 17 , 2102244 (2021).
32 F. Calle’Vallejo, J. I. Martínez, J. M. García’Lastra, P. Sautet, D. Loffreda, Fast prediction of adsorption properties for platinum nanocatalysts with generalized coordination numbers. Angewandte Chemie Int. Edn. 53 , 8316–8319 (2014).
33 F. Calle-Vallejo , Finding optimal surface sites on heterogeneous catalysts by counting nearest neighbors. Science 350 , 185–189 (2015).26450207
34 NIST computational chemistry comparison and benchmark database, NIST standard reference database number 101, release 22, May 2022, editor: Russell D. Johnson III (2022). http://cccbdb.nist.gov/, 10.18434/t47c7z.
35 B. Ehresmann, B. Martin, A. H. C. Horn, T. Clark, Local molecular properties and their use in predicting reactivity. J. Mol. Mod. 9 , 342–347 (2003).
36 Z. Li , AI-Catalysis. GitHub. https://github.com/knoevenagel/AI-Catalysis. Deposited 30 September 2023.
