Predicting Segregation Energy in Single Atom Alloys Using Physics and Machine Learning.

Single atom alloys (SAAs) show great promise as catalysts for a wide variety of reactions due to their tunable properties, which can enhance the catalytic activity and selectivity. To design SAAs, it is imperative for the heterometal dopant to be stable on the surface as an active catalytic site. One main approach to probe SAA stability is to calculate surface segregation energy. Density functional theory (DFT) can be applied to investigate the surface segregation energy in SAAs. However, DFT is computationally expensive and time-consuming; hence, there is a need for accelerated frameworks to screen metal segregation for new SAA catalysts across combinations of metal hosts and dopants. To this end, we developed a model that predicts surface segregation energy using machine learning for a series of SAA periodic slabs. The model leverages elemental descriptors and features inspired by the previously developed bond-centric model. The initial model accurately captures surface segregation energy across a diverse series of FCC-based SAAs with various surface facets and metal-host pairs. Following our machine learning methodology, we expanded our analysis to develop a new model for SAAs formed from FCC hosts with FCC, BCC, and HCP dopants. Our final, five-feature model utilizes second-order polynomial kernel ridge regression. The model is able to predict segregation energies with a high degree of accuracy, which is due to its physically motivated features. We then expanded our data set to test the accuracy of the five features used. We find that the retrained model can accurately capture E seg trends across different metal hosts and facets, confirming the significance of the features used in our final model. Finally, we apply our pretrained model to a series of Ir- and Pd-based SAA cuboctahedron nanoparticles (NPs), ranging in size and FCC dopants. Remarkably, our model (trained on periodic slabs) accurately predicts the DFT segregation energies of the SAA NPs. The results provide further evidence supporting the use of our model as a general tool for the rapid prediction of SAA segregation energies. By creating a framework to predict the metal segregation from bulk surfaces to NPs, we can accelerate the SAA catalyst design while simultaneously unraveling key physicochemical properties driving thermodynamic stabilization of SAAs.

Introduction

Over the years, metal alloys have received great attention due to their broad range of applications; these include among others, drug delivery,[1] catalysis,[2,3] and optics.[4] These systems exhibit electronic and chemical properties that differ from their monometallic counterparts, which can yield enhanced catalytic activity, selectivity, and overall stability.[2] Moreover, metal alloys have tunable properties (by controlling size, shape, and metal composition) that can be tailored for a wide range of catalytic reactions, including oxygen reduction (ORR),[5] CO oxidation,[6] and alkane dehydrogenation reactions.[7]

Single atom alloys (SAAs) makeup a class of bimetallic materials that are experiencing rapid growth in interest, especially for catalysis due to their efficiency[8] and high selectivity.[9] Many SAA catalysts are synthesized by doping a single highly reactive precious metal into a less reactive (and typically cheaper) host metal, resulting in structurally uniform surface active sites.[10] These active sites provide a unique and well-defined electronic structure that can be tuned precisely to achieve high catalytic activity[8,11] and selectivity[12] for a range of reactions. For instance, Liu et al. demonstrated that doping a single Pt on a Cu(111) surface resulted in reduced CO binding strength (compared to pure Pt systems), which can prevent CO poisoning during catalysis.[13] In other cases, the host metal can provide a second type of active site that contributes to the enhanced performance. Kyriakou et al. doped Pd into the less active Cu(111) to examine the catalytic activity for hydrogenation reactions. It was found that the addition of Pd led to a bifunctional surface that could selectively hydrogenate styrene and acetylene with high activity.[9] Under such catalytic environments, surface reconstruction and metal diffusion may occur depending on the mixing behavior of the two metals, facet, and the morphology (i.e., size and shape).[14,15] Moreover, sintering can lead to decreased catalytic activity.[16] Previously, Tan et al. developed a statistical model to predict the adsorption energy of different SAAs supported on metal oxides to test for sintering.[17] Therefore, understanding the thermal stability of SAAs, which is dictated by the interactions between the host and dopant metals,[10] is critical to the design of stable catalysts that maintain active surface sites.

A critical measurement of SAA stability is surface segregation energy (Eseg),[10] which can be defined as how thermodynamically favorable it is for a heterometal dopant to occupy a surface site compared to a bulk site in a metal host.[18] Since the surface environment is crucial to processes like adsorption[19] —and thus catalysis—Eseg plays a key role in the SAA catalyst design as it defines whether a dopant will preferentially reside on the surface of the catalyst. Understanding dopant segregation tendencies enables the proper selection of host–dopant pairs for targeted catalytic reactions. For example, it was previously found that Pd atoms segregate to the surface of inert coinage metals (i.e., Au and Ag), preventing CO poisoning and enhancing the CO oxidation catalytic activity.[20] Conversely, in ORR, the segregation of X (X = Cu, Ni, and Co) to the Pt surface may destabilize the alloy in an acidic setting.[21] Therefore, understanding metal surface segregation tendencies of the catalyst without any adsorbates, is the first, but crucial step for effective catalyst design.

Various methods have been implemented such as the tight-binding theory[22,23] and density functional theory (DFT)[24,25] to gain insights into the stability of the SAA. Back in 1999, Ruban et al. applied Friedel’s rectangular state density model[26,27] to predict Eseg in FCC(111), BCC(110), and HCP(0001) surfaces.[28] Although the model predicts Eseg trends well, it is unable to correctly predict all SAA cases due to its lack of ability in capturing structural effects. Despite being effective, utilizing ab initio methods like DFT can be computationally expensive and time consuming. Thus, new approaches are needed to accurately predict the metal segregation and reveal the physicochemical reasons driving segregation across different metals. Previous work leveraged machine learning techniques with DFT calculations to develop models that capture Eseg for the constrained subsets of SAAs. In 2019, Farsi and Deskins developed a statistical model, based on atomic properties and surface energies, that predicts DFT Eseg on FCC (100), (110), (111), and (210) surfaces in Pt, Ir, Pd, and Rh systems.[18] Ologunagba and Kattel developed a non-linear model using gradient boosting regression that predicts DFT Eseg across different FCC, BCC, and HCP surfaces.[29] Rao et al. developed a model (based on bond-counting theory) to predict formation energies on FCC (111), BCC (110), and HCP (0001) of different SAAs with surface, subsurface, dimer, and adatom dopants.[30] Although these models can capture Eseg trends in certain systems, there is still a need for a general framework capable of predicting Eseg across the complete metal material space, spanning from bulk SAA systems to nanoparticles (NPs). Such a model requires descriptors that capture the underlying physics dictating segregation in SAAs and consider different coordination environments on the catalyst surface. In 2018, Yan et al. developed the bond-centric model (BCM), a coordination number (CN)-based model that captures the stability of bimetallic NPs of any chemical ordering, morphology, and metal composition.[31] Since the BCM captures explicit bond energies, it can be used to probe Eseg in NPs. In its original development, the authors employed the BCM to investigate the metal segregation in Cu54Au NP and found qualitative agreement with DFT, suggesting its potential to capture the physics of metal segregation.

We herein develop a framework to predict DFT Eseg on a range of FCC surfaces in Pt-, Ir-, Pd-, and Rh-based SAAs (tabulated from the extensive work of Farsi and Deskins[18]) by leveraging the BCM,[31] elemental properties, and machine learning techniques. After exploring the limitations of predicting Eseg on bulk surfaces with the BCM (which was originally developed on metal NPs[31]), we introduce additional elemental properties and machine learning to develop a model with improved accuracy. By dissecting physical relationships found in the BCM, we enrich our feature set and train a model that accurately predicts Eseg,DFT for SAAs with FCC dopants. We then extend our approach to build a single model that captures Eseg of SAAs made from FCC, BCC, and HCP dopants with little loss in performance. Finally, we employ our model (trained on SAA bulk surfaces) on a sample of SAA NPs, ranging in size and dopant coordination to explore the model’s extrapolation capabilities beyond bulk systems. Our model achieves similar accuracy on the SAA NPs, revealing its ability to capture size effects on Eseg.

Results and Discussion

The surface segregation energy on FCC (100), (110), (111), and (210) surfaces in Ir-, Pt-, Rh-, and Pd-based SAAs doped with FCC metals was calculated using the BCM and eq (see the methodology section for additional details on these calculations). To test the accuracy of the Eseg,BCM, it is compared against the Farsi and Deskins data set,[18] as shown in Figure . It is important to note again that the FCC metal combinations reported in the data set were not all used in the comparison due to BCM’s metal pair limitations.[31] The results reveal that the BCM cannot predict Eseg,DFT across the data set, resulting in a mean absolute error (MAE) of 0.513 eV. Although the BCM captures Eseg,DFT trends accurately for Pd metal host and somewhat for Pt and Ir, it fails completely for Rh-based systems. Moreover, there is a clear systematic error arising in Pt-based systems (purple points in Figure ). This could potentially be due to the fact that the BCM does not account for any strain effects between the dopant and host metals. In an effort to reduce errors from our Eseg,BCM predictions, we can implement machine learning techniques to train a model that captures the missing physics.

Figure 1

Parity plot of Eseg,BCM vs Eseg,DFT of SAAs (FCC host metals and FCC dopants). Different host metals are indicated in different colors.

Machine learning provides a powerful means to unravel complex patterns within a high-dimensional feature space. To leverage machine learning with our problem, we first generated a data set of physically relevant features. Our initial list of features includes elemental properties of hosts and dopant metals. Importantly, these properties do not require expensive first principles calculations to determine. They cover physical (e.g., atomic radii), electronic (e.g., atomic electronegativity), and bulk [e.g., bulk cohesive energy (CE)] properties. Since Eseg is dependent on host–dopant interactions,[10] we also incorporate features based on the difference in properties (e.g., Δradius = radiushost – radiusdopant). We take into consideration the different contributions of each metal to ensure a systematic and unbiased approach in the feature selection process. Additionally, expanding each term allows us to understand which representation of these descriptors has a more dominant effect on the Eseg. Combining these descriptors with Eseg,BCM, we achieve an initial set of 29 tabulated features (see Table S1) for each host–dopant pair.

After tabulating our initial set of elemental features, we conducted feature selection through recursive feature elimination with cross-validation (RFECV) analysis. The RFECV results, as shown in Figure a, reveal that five features provide a reasonable tradeoff between low MAE and model complexity (i.e., degrees of freedom). For example, shifting down to four features leads to a significant decrease in accuracy, while moving up to six features yields minimal gains. Therefore, the top five features were selected for further analysis, which are (in order of importance), the difference in the bulk CE of the host and dopant (ΔCEbulk), Gordy electronegativity of the host (χgordy,host), the bulk CE of the host (CEbulk,host), the difference in the lattice constant (Δa), and Eseg,BCM. Next, we constructed Pearson’s correlation table to determine whether selected features strongly correlate and thus convey the same information. The results in Figure b show that there is a strong correlation between ΔCEbulk, CEbulk,host, and Eseg,BCM, which can be expected since the BCM is a function of host and dopant CEbulk values.

Figure 2

Applying (a) LASSO-based RFECV to the initial feature set (Table S1). Dark blue points indicate LOOCV-MAEs of the best-performing model at the given number of features. Shaded blue represents the standard deviation of MAEs found during LOOCV of each top-performing model. (b) Pearson’s correlation table of the five most significant features (based on RFECV).

Our overall findings suggest that the initial feature set is insufficient for predicting Eseg,DFT of FCC metal pairs. Even with five features, we still only achieve a leave one out cross-validation (LOOCV)-MAE of ∼0.23 eV (Figure a). Furthermore, based on RFECV rankings, there are four features found to be more important than Eseg,BCM. Not only does ΔCEbulk contribute the most but also CEbulk appears in multiple top features, indicating its significance in predicting Eseg,DFT. Therefore, before attempting to train additional machine learning models, we sought to expand our potential feature set. Rather than constrain our solution to improving the predictions from Eseg,BCM, we instead chose to remove Eseg,BCM and replace it with three additional descriptors inspired directly from the physics that BCM is based upon. The switch results in a new extended feature set that contains nine new features (as highlighted in Table S1). By removing Eseg,BCM, we completely avoid the need for any DFT calculations to generate the features. However, we note that the DFT calculations for parameterizing the BCM are minimal (i.e., require a single dimer calculation between each metal pair). Nevertheless, a Pearson correlation table, as shown in Figure S2, reveals strong correlations between Eseg,BCM and some of the new features, indicating that we still capture most of the information within the BCM-calculated value. Thus, by dissecting the BCM, we achieve new, physically motivated features to improve our Eseg,DFT predictions.

Leveraging our extended feature set, we again employed RFECV to explore feature importance (Figure S3a). By incorporating these new features into our data set, we can now achieve comparable MAEs to the previously studied five-feature model with only three features (∼0.23 eV). Once again, we find that five features give the best tradeoff between accuracy (low MAE) and complexity (fewer terms). These include, in order of importance, ΔCEbulk/CNdopant, χgordy,host, atomic radius of the dopant (rdopantatomic), difference in the electron affinity (ΔEA), and CEbulk,host/√CNdopant. Importantly, the results verify that Eseg,BCM is not required to predict Eseg,DFT, and it is instead better to include terms inspired by the BCM. Further analysis of the five features using Pearson’s correlation table (Figure S3b) reveals that the first and last terms are highly correlated (Pearson’s correlation = 0.75); hence, CEbulk,host/√CNdopant is dropped from our selected features. Previously, it was found that the CE, atomic radius, and electronegativity influence the segregation energy of metal alloys.[2,23] Of note, these are the same terms—with the addition of the difference in the electron affinity—selected by the RFECV analysis. Using the selected features, different regression models were trained and compared (see the methodology section for additional training details). We note that model training was carried out with three and four features to compare the performance of less complex models and investigate any potential overfitting. The results for three-feature and four-feature models are shown as parity plots in Figures S4 and S7, respectively. Kernel ridge regression (KRR) was found to produce the lowest MAE for both three- and four-feature models (Figures S4 and S7). The results show that there is no significant difference between the MAEs when the models are simplified to three features (∼0.05 eV for a second-order polynomial KRR and ∼0.02 eV for a third-order polynomial KRR). Additionally, we checked for overfitting in the different models using bootstrapping analysis[32] (Figures S5 and S8). We note that the larger the difference between the bootstrapped training and test errors, the more likely the model is overfitting. Although both KRR models showed similar differences in the MAE between the test and training set, a better MAE distribution is observed using three features (Figures S6 and S9). In other words, there is a greater overlap between training and testing error distributions for the three-feature KRR model compared to the four-feature case. Based on our overall analysis, we selected the second-order polynomial KRR (α = 0.0015 and γ = 0.03), which uses three features (ΔCEbulk/CNdopant, χgordy,host, and rdopantatomic) to predict Eseg,DFT of SAAs formed from FCC metal pairs.

The trained KRR model exhibits a LOOCV-MAE of 0.18 eV, a root-mean-squared error (RMSE) of 0.24 eV, and an R2 of 0.94. Compared to the initial analysis performed (Figure ), the model uses fewer features without having to comprise accuracy, which reiterates the importance of extending our feature set. As shown in Figure a, the KRR model captures the Eseg,DFT trends across all four FCC hosts. However, the model systematically underpredicts Eseg of Pt-based SAAs. We found a similar result with the Eseg,BCM data plotted in Figure , though the new model has much smaller deviation. Further separating the results reveals a relatively small range of accuracies (0.077–0.23 eV) between host metals, as shown in Figure b. The order of host-based MAE follows: Pd < Ir ≈ Rh < Pt. In addition, Figure b colors the data by dopant type, revealing no systematic errors arising from dopants. This is good indication that the selected features accurately describe the interactions between the host and dopant metals. Overall, employing our extended feature set and machine learning techniques enabled the development of a robust model to predict Eseg,DFT of FCC-based SAAs.

Figure 3

Parity plots of Eseg,model (second-order polynomial KRR using three features) vs Eseg,DFT of FCC metal hosts doped with FCC metal dopants. Data are presented as (a) a single plot of all the results (where color represents host metal) and (b) multiple subplots (where color indicates the dopant) separated by the individual host metal: (i) Ir, (ii) Pt, (iii) Rh, and (iv) Pd. MAE for the model (a) is calculated from LOOCV, whereas the subplot MAEs (b) are determined post-training from the segmented data.

Based on our model’s predictions, we find that when Pd is doped with FCC metal, segregation of the dopant rarely occurs. Similarly, when Pt is doped with FCC metal, there is either no segregation or weak segregation of the dopant, as opposed to Ir and Rh metal hosts. Additionally, Au dopants almost always segregate, regardless of the facet and metal host. On the other hand, Ir dopants are less likely to segregate, as shown in Figure . The Ir-host has a wider Eseg distribution (i.e., range of values) followed by Rh and Pt and Pd. If we further compare the host metals to each other, Ir host has a wider distribution compared to Rh, and Pt host has a wider distribution than Pd. Based on this observation, we can conclude that the Eseg distribution increases as we go from period 5 to 6. These distributions are crucial in determining the most stable SAAs.

In previous analyses, our focus was on predicting Eseg of FCC SAAs. The question now becomes: can we extend our methodology to capture SAAs formed from FCC host metals with FCC, BCC, and HCP dopants? To answer this question, we used the extended feature set with the complete Farsi and Deskins SAA data set[18] and employed RFECV analysis. As shown in Figure S10a, the lowest MAE is observed with just five features. Based on the correlation table in Figure S10b, there are no repetitive features used; thus, we select five features for subsequent model training. The five selected features for the final model are, in order of importance, ΔCEbulk/CN, χgordy,host, ΔEA, rdopantatomic, and first ionization potential of the dopant (IPdopant). We note that the top four features obtained from RFECV analysis on predicting FCC SAAs (Figure S3a) are the same top four features found in this generalized analysis. With our top features selected, we trained a series of different regression models. Compared to the other models investigated, KRR again produced the lowest MAE (Figure S11). Although the MAE resulting from the third-order polynomial KRR is lower compared to the second-order polynomial KRR, it is observed that the third-order polynomial KRR model is overfitting, which is confirmed through bootstrapping analysis. As shown in Figure S13, the difference between the bootstrapped training and test errors for third-order and second-order polynomial KRR is ∼0.1 and ∼0.03 eV, respectively. Hence, we selected the second-order polynomial KRR (α = 0.002 and γ = 0.02) as our final model to predict Eseg,DFT. Importantly, we have provided all data and Python code to generate the model on our GitHub. The DFT data set along with our descriptors and model can be found here https://github.com/mpourmpakis/BCSeg.

Model performance on the Farsi and Deskins data set[18] is illustrated in Figure . The model captures the Eseg trends of the bulk systems with a high degree of accuracy, producing a LOOCV-MAE of 0.22 eV, a RMSE of 0.28 eV, and an R2 of 0.91. The Farsi and Deskins model (which also uses five features) exhibits a RMSE of 0.43 eV and an R2 of 0.77. Furthermore, although the Farsi and Deskins model also captures FCC, BCC, and HCP Eseg,DFT trends, our model eliminates the need to use DFT.[18]

Figure 4

Parity plots of Eseg,model (second-order polynomial KRR using five features) vs Eseg,DFT of FCC metal hosts doped with FCC, BCC, and HCP metal dopants. Data are presented as (a) single plot of all the results (where color represents host metal) and (b) multiple subplots (where color indicates the surface facet) separated by (i) FCC, (ii) BCC, and (iii) HCP dopants. The MAE for the model (a) is calculated from LOOCV, while the subplot MAEs (b) are determined post-training from the segmented data.

We note that our model predicts a reversed segregation behavior for a subset of SAAs that exhibit weak segregation preference (i.e., all fall within the range of |Eseg,DFT| ≤ 0.38 eV). Importantly, the opposite segregation trend accounts for only 3.75% of the entire data set present in Figure . Moreover, despite the opposite predictions, the model predicts the Eseg of these cases to be within 0.5 eV from thermoneutral, thus capturing their weak segregation behavior relative to the wide Eseg range of the data set. Thus, despite the few incorrect predictions of segregation preference, our final model still correctly predicts the weak segregation tendencies of these systems.

To gain more insights on the accuracy of our model, we compared it to another model reported in the literature trained on the same data set. Farsi and Deskins refitted their data set to Yu et al.’s model[25] and found that it produced an R2 of 0.61 and a RMSE of 0.60 eV.[18] Based on these metrics, we can conclude that our model shows enhanced performance. Additionally, Yu et al.’s model does not differentiate between the different facets due to its use of metal-dependent surface energy, which is blind to facet type. Our model, on the other hand, can capture Eseg trends across (111), (100), (110), and (210) surfaces because it importantly incorporates CN of the dopant. Indeed, extending the features set to incorporate BCM-inspired descriptors enabled us to generate a powerful model for rapid Eseg,DFT prediction. To break down the accuracy of the model further, we divided the data into three plots based on the dopant metal type (Figure b). The model captures the trends across each dopant metal type with about the same accuracy (∼0.22 eV). The subplots also demonstrate that there are no systematic errors arising from different surface facets. For FCC and HCP dopants, there is a wider Eseg distribution compared to BCC dopants, as illustrated in Figures b and S12. Additionally, BCC dopants have fewer cases where the dopant is promoted to the surface and opposed to FCC and HCP dopants, regardless of the different facets. This could be due to the difference in the bulk CN of BCC dopants (CB = 8), as compared to FCC and HCP dopants (CB = 12). From the BCM, we can verify that the CB is an important descriptor in determining the stability of the NPs. Moreover, compared to the other dopants, the segregation of BCC dopants in FCC metal hosts is thermodynamically unfavorable; hence, the segregation of BCC dopants is less likely to be observed in SAA synthesis. Along with Figure S12, this indicates the great influences on the Eseg behavior. These trends can aid in determining which of the SAA combinations are experimentally feasible. From our previous observations, it was found that the metal dopant type significantly affects not only the Eseg distribution but also the promotion of the dopant itself. We also show that the model captures Eseg trends across the surface facets with similar accuracy (i.e., the variation of MAEs across different facets is within 0.05 eV). These trends are particularly critical as the position of the dopant plays a key role in determining the stability of the SAA alloys, which in turn will affect its application.

The results highlight the ability of our model to capture the underlying physics that govern Eseg,DFT across a wide range of SAAs. This stems from the features present in the final model, which reflect the physics behind segregation energy. For instance, prior work has revealed that CEbulk of metals contributes to alloy mixing behavior (e.g., core–shell tendencies).[23] Furthermore, it was previously proven that CEbulk and CN are key features in capturing the stability of a local site on a metal surface because they directly contribute to the bond energies of the metals.[33] Since the stability of the metal alloy depends on the segregation tendencies of the metal core, the first (and most significant) feature in the model is essential in describing the stability of the SAA. Of note, it is also the only feature that differentiates between surface facets. In the work by Ruban et al., charge-transfer effects were taken into consideration to compute surface segregation energy.[28] Importantly, our developed model captures charge-transfer effects between the host and the dopant by incorporating electron affinity, first ionization potential, and the electronegativity of the metals. Finally, the atomic radius of the dopant accounts for strain effects present in the system, which is an important structural property that can influence metal segregation.[34]

To further test the accuracy of KRR and the descriptors used in capturing different facets and metal hosts, we expanded our Eseg search to incorporate more DFT data from the literature. Our expanded data set includes Farsi and Deskins’ data set plus the Eseg of Ni(111)[25] and d9 host metals doped with platinum group metals on (211),[20] (100),[20] and (111)[35] surfaces. We next retrained our KRR model using the expanded data set. The results are shown in Figure S15 of the Supporting Information. Comparing Figures a to S15, the difference in the overall accuracy is ∼0.04 eV, revealing the ability of our model to capture a broader range of SAA systems. To further understand the performance of our model, we compared MAEs based on the different host metals. We observe from Figure S16a that the model captures the trends of the different metal hosts with similar accuracy (∼0.23 eV) with an exception of Ni-host cases, which exhibits significantly worse accuracy (∼0.45 eV). The poor accuracy can be qualitatively observed in Figure S15, though we note that the overall Eseg trend between Ni-host samples is largely captured. We next analyzed the results of our outlier cases (i.e., Ni host) to understand where the discrepancy between DFT and the model is arising. By splitting the Ni-host cases based on dopant packing type (FCC, BCC, HCP; Figure S16b), we found that the model captures the Eseg trends of most Ni-host cases but fails to predict the Eseg when Ni is doped with HCP metals. This trend is not observed across the dopant packing in other host metals, making it difficult to explain why the model fails to predict Ni-host HCP-dopant SAAs. Overall, the results suggest that our model does not reliably predict Eseg for Ni-host HCP-dopant SAAs, which is an important limitation to specify. Moreover, future work should focus on overcoming this limitation toward a complete, generalized model.

To analyze the model’s predictions further of the entire data set, we investigated the model’s accuracy in capturing the different dopant types. Similar to our previous analysis, we find that there are no systematic errors arising from the dopant type as we observe a similar model performance on FCC, BCC, and HCP dopants (Figure S15b). It is not a coincidence that the retrained model performs well in predicting Eseg trends across different metal hosts and facets as the descriptors utilized in our final model capture the driving forces that govern surface segregation. Additionally, this finding is crucial as we confirm that the ability of our model in capturing the Eseg trends is not just due to the complexity of the model but importantly the descriptors used in our final model.

Up until this point, the focus has been to predict Eseg,DFT in bulk systems (periodic surfaces). In catalysis, NPs are widely investigated and used due to their high surface-to-volume ratio and tunable catalytic behavior.[3] NP properties can be tuned by adjusting morphology (size and shape), which can in-turn be controlled through different synthesis methods.[36,37] Similar to bulk systems, Eseg is a critical property to synthesize SAA NPs.[10] Previous work has revealed size-dependent behavior of Eseg in SAA NPs.[34] Thus, by developing a model that can capture Eseg in NPs, we can provide crucial insight into the SAA NP catalyst design. Toward this goal, we performed DFT calculations to determine Eseg,DFT for a series of Pd-, Ir-based, and CuIr SAA cuboctahedron NPs. NP sizes of 147, 309, and 561 are investigated due to their unique coordination sites. We focus on cuboctahedra since it is a highly stable shape within this size regime.[38] Instead of reapplying our machine learning methodology, we directly use our model trained on the Farsi and Deskins data to test its efficacy in capturing Eseg for SAA NPs. Since the model was previously trained on periodic surfaces, we were able to use tabulated CEbulk values for our first feature. To extend it to SAA NPs and avoid DFT, the BCM is used to calculate CE of the systems. The ΔCEbulk term is replaced by ΔCENP since this term can be considered CE of the monometallic reference systems (CEbulk is representative of periodic systems, while CENP represents NPs). We note that this does not change the form of our model since the BCM scales to CEbulk for infinitely large systems (i.e., periodic slabs). The results are presented in Figure as a parity plot. Remarkably, the pretrained model gave an MAE of 0.24 eV, which indicates that it is capturing the underlying physics behind Eseg at the NP size regime. Furthermore, the model captures the Eseg,DFT of Pd-based NPs more accurately as compared to Ir-based NPs, which is the same behavior observed within the periodic slab predictions (Figure a). Importantly, we note that the CuIr case utilizes a different host metal (Cu) that was not included in the original slab training data. It should also be noted that this metal host (Cu) is at a different period and group from the ones studied in this work. This could explain why we find the other metal combinations, which were used in model training, closer to the parity line, as compared to CuIr. Nevertheless, the model is still able to capture the size effects across all the systems, revealing its ability to generalize to new sizes (NPs) and metal hosts. From Figure , we also observe that the model appears to perform better for 147-NP size, as compared to 309 and 561. It is crucial to understand because the ΔCEbulk term is replaced by ΔCENP, we cannot conclude that the model will perform better for larger NP sizes. Still, the results are extremely accurate, considering that the model has never seen SAA NP data during its training. In this way, we can conclude that the second-order polynomial KRR model exhibits promise as a general predictive tool for capturing Eseg,DFT across all size regimes—from NP to bulk SAAs.

Figure 5

(a) Parity plot between Eseg,model and Eseg,DFT for a set of Pd- and Ir-based SAA NPs. Marker type indicates NP size (147-, 309-, or 561-atom NPs), and color indicates host–dopant pair (e.g., yellow represents Ir-based NPs with the Ni dopant). (b) Example NPs of different sizes and CNs (i) 147-atom NP doped on CN = 7 (blue), (ii) 309-atom NP doped on CN = 8 (purple), and (iii) 561-atom NP doped on CN = 9 (green).

We note that our final results reinforce the importance of using the BCM in our analysis. Although using the original BCM expression failed to accurately capture segregation behavior, by dissecting the BCM, we formulated new terms to extend our feature set, which caused significant improvement in the performance of our final model. Additionally, as we shifted from bulk systems to NPs (Figure ), we leveraged the BCM again to calculate the CENP in our final model. Doing so enabled us to capture Eseg in NPs while still avoiding DFT calculations—proving the significance of the BCM in developing the final model.

Despite the model’s ability in predicting Eseg in periodic slabs and NPs, we do acknowledge that the model could have limitations. Since machine learning models are known to struggle with extrapolation,[39] we must be critical when predicting how our model will perform on new systems. For example, we would expect our final model to struggle to predict metal combinations and facets not included in its original training set. Although we note the remarkable agreement with SAA NP results, the model could still have difficulties predicting Eseg of CN sites that it has not previously seen (e.g., CN = 5). For future work, we plan to improve upon the final model by investigating Eseg across a broader range of metal types and surfaces, including more complex surfaces and NP shapes. Nevertheless, our model is an important step toward a universal predictor that captures segregation behavior across the spectrum of the SAA configuration space.

Conclusions

In this work, we utilized the fundamental principles behind the BCM, elemental properties, and machine learning to capture DFT surface segregation energy in FCC metal SAAs, ranging in host metal type (Ir, Pd, Pt, and Rh), surface facets [(100), (110), (111), and (210)], and metal dopants (FCC, BCC, and HCP). Initially focusing on FCC-based SAAs (FCC hosts and dopants), we first employed the BCM, which was previously shown to capture bimetallic NP stability, to investigate the Eseg in bulk metal systems. The results illustrated that Eseg,BCM failed to accurately capture Eseg,DFT of FCC-based SAAs. After exploring machine learning approaches to improve the Eseg,BCM predictions, we found that replacing Eseg,BCM with new descriptors inspired by the BCM allowed us to generate a rich set of features that completely avoid the need for DFT parameterization. With our extended feature set, we redeployed our machine learning methodology to develop a model capable of predicting Eseg,DFT of FCC-based SAAs with a high degree of accuracy (LOOCV-MAE = 0.18 eV). The model, a second-order polynomial KRR that uses three features showed a significant improvement in the performance relative to our analyses with Eseg,BCM. Next, we repeated our machine learning approach to develop a model that captures Eseg,DFT of SAAs formed from FCC, BCC, and HCP metal dopants. Again, our results yielded a second-order polynomial KRR but this time with five features. Importantly, the same three features used in the first model (FCC dopants) were found in the final model (FCC + BCC + HCP dopants). The overall model exhibited robust performance across all host, dopant, and surface types, suggesting that the selected features capture the underlying physics that govern SAA Eseg. Indeed, the five features describe critical host–dopant interactions, including thermodynamic stability (metal CEbulk), structural effects (atomic radius and CN), and electronic effects (responsible for charge transfer) such as the electron affinity, first ionization potential, and electronegativity. We then incorporated more data from the literature into our data set to further test the accuracy of these descriptors. We found that the retrained model accurately predicted the Eseg across different metal hosts and facets, confirming the importance of the features utilized in our final model. To increase the applicability and impact of our model, we conducted an extrapolation test and applied it to a series of SAA NPs. The periodic slab-trained model was able to capture Eseg,DFT of the NPs with remarkable accuracy (MAE = 0.24 eV). The results highlight the physical nature of the features that parameterize our model, suggesting its potential use as a generalized tool to predict Eseg,DFT. Overall, our accurate model enables the rapid screening of the broad materials space toward the improved SAA catalyst design.

Methodology

Density Functional Theory

DFT Eseg of bulk FCC surfaces was obtained from a previously published work, which we refer to as the Farsi and Deskins data set.[18] The data were used to probe the performance of BCM for predicting Eseg as well as train and validate all machine learning models reported herein.

The BCM differentiates bond types (homo- vs heteroatomic) through a γ term, which captures relative contributions of different metal types to a bond.[31] These γ values present in the BCM are computed from DFT dimer bond dissociation energies (BDEs). Therefore, we performed DFT calculations to compute BDEs (and thus γ) for a sample of FCC metal pairs in the Farsi and Deskins data set.[18] We note that since the BCM is limited to specific metal pairs (∼85% of potential transition metal alloys),[31] some metal combinations were excluded from the γ calculations. The calculations were performed using PBE exchange–correlation functional[40] with def2-SV(P)[41] basis set, accelerated with the resolution of identities approximation,[42] as implemented in the Turbomole package.[43] PBE was used to maintain consistency with the functional used on the Farsi and Deskins data set. All geometries were relaxed utilizing a self-consistent field (SCF) energy optimization of each step (with an SCF convergence criterion of 10–6 Ha) until the interatomic forces were no greater than 10–3 Ha/Bohr.

Eseg,DFT of 147-, 309-, and 561-atom SAA cuboctahedron NPs was calculated in the CP2K package[44] using the PBE exchange–correlation functional.[40] A computational box size of 34 × 34 × 34 Å was used for all calculations. Double-zeta valence-polarized basis set with a cut-off of 500 Ry in combination with Goedecker–Teter–Hutter pseudopotentials was used.[45] All geometries were relaxed until the interatomic forces converged to 0.025 eV/Å with an SCF convergence criterion of 10–7 Ha for each step. All NP structures were constructed from the Atomic Simulation Environment Python package.[46] These DFT calculations were performed to address the applicability of our developed models on NPs (i.e., beyond periodic surfaces).

We note that both our dimer and NP calculations used the same functional and similar basis sets to ensure that our results are at the same level of theory. Our package choice varied due to the sizes of our systems. Turbomole is more commonly used to calculate the energetics of small, molecular systems. On the other hand, CP2K is often used to compute the energetics of larger systems, such as periodic surfaces and NPs, due to its scaling capabilities. Hence, we used Turbomole for our metal dimer calculations (i.e., computing BDE to parameterize the BCM) and CP2K to calculate Eseg in SAA NPs.

Machine Learning

We employed a supervised machine learning approach to develop a generalized segregation energy model (i.e., predict Eseg,DFT) for SAAs. Along with BCM-calculated segregation energy (Eseg,BCM), we first tabulated several elemental properties of the host and the dopant. These include the lattice constant, electronegativity, electron affinity, and first ionization potential, which were obtained from the Mendeleev Python package.[47] Other elemental features include atomic radius[48] and bulk CE.[49] Before feature selection and model training, all features were standardized by transforming the data such that each feature distribution has a mean of 0 and a standard deviation of 1. This allows for equal contribution of the different features used in analyses. A full list of the elemental properties can be found in Table S1 in the Supporting Information.

RFECV[50] was employed with the Scikit-Learn Python package[51] as a means of feature selection, or to determine which of the tabulated features contribute most to predicting Eseg,DFT. The method begins by training a LASSO model using all given features. Next, it eliminates the least important feature (based on model coefficients) and retrains a new LASSO model. The process continues by recursively eliminating the least important feature until a one-feature model is generated, containing the most important feature. Model accuracy is measured using MAE from LOOCV where the MAEs are calculated usingwhere y is the actual output value, ŷ is the predicted output value, and x is the total number of data points. In order to ensure an optimal model performance during the process, the regularization parameter (hyperparameter, α) in LASSO varied from 0.001 to 1 (by factors of 10) to minimize MAE. Once completed, RFECV provided a range of top-performing models for all possible number of features. By comparing the model performance versus number of features, we can rank feature importance (e.g., most important feature is selected for the one-feature model) while simultaneously determine an optimal number of features to select (based on error vs complexity) for predicting Eseg,DFT.

After feature selection based on RFECV analysis, we trained a series of regression models to compare their performance for predicting Eseg,DFT. These include LASSO,[52] linear regression (OLS),[53] support vector regression (SVR),[54] and KRR.[55] Hyperparameters found in KRR, LASSO, and SVR were tuned using GridSearchCV[56] by minimizing MAE calculated using the LOOCV method. All optimized hyperparameters are reported in Tables S2–S4 in the Supporting Information. In addition to the MAE, the RMSE is also computed (eq ), using LOOCV, to directly compare the model performance to other studies.

All model implementations, including training and validation, were performed using the Scikit-Learn Python package.[51] To aid in the final model selection and further ensure that there is no overfitting, we carried out bootstrap analysis[32] (also implemented in Scikit-Learn). Bootstrapping is a sampling technique to estimate population errors from a sample of data. Since our Eseg,DFT data set only covers a sample of the potential SAA material space, bootstrapping provides a more robust measure of the model performance and degree of overfitting. Upon further analysis of the model performance, we complete our process by selecting our final, top-performing model to predict Eseg,DFT.

Calculating Eseg with the BCM

The use of the BCM has been shown to be effective for predicting the stability of metal NP sites. For instance, Dean et al. utilized the BCM to derive a descriptor (CElocal) that captures the local stability of a metal site. The descriptor was included with other descriptors and machine learning to accurately predict adsorption energy in periodic systems and NPs.[33] Following a similar idea of extending to periodic systems, we leverage the BCM to calculate the Eseg of SAAs in metal slabs. Eseg of SAAs can be calculated by taking the difference in energy of the system when the dopant is on the surface and when the dopant is in the bulk. However, by reframing Eseg as a function of CE, we can employ the BCM to calculate CEs and thus determine Eseg,BCM, as shown in eq ,where n is the total number of atoms in a system and CE is the CE of the SAA when the dopant is on the surface (CN = 9, 8, 7, and 6) or bulk (CN = 12). Under this formulation, a positive Eseg value would indicate that the dopant would prefer to stay in the bulk, while a negative Eseg value indicates that the dopant would segregate to the host metal surface. Additional details on the derivation of eq and the DFT-based CE equations are provided in Section S1 of the Supporting Information. We note that computing CE terms within eq is carried out using the BCM, which computes CE of a metal NP (equation provided in prior work[31]). It is important to note that the BCM was developed to calculate CE of monometallic and alloy NPs and not periodic surfaces.[31] However, the BCM is a coordination-based model that computes explicit bond energies based on local coordination environments (i.e., metal types and CNs). For large NP SAAs, the local coordination of a dopant on the surface (i.e., CN of first and second neighbors) becomes constant. Since Eseg is calculated from a difference in CEs, the consistent coordination environment causes Eseg,BCM to converge at large NP sizes, as shown in Figure S1. Thus, by leveraging large enough SAA NP models with the correct surface facets, we can use the BCM to determine Eseg,BCM for bulk systems (periodic surfaces). Details regarding the morphology effects on Eseg,BCM can be found in Section S2 of the Supporting Information.