Beyond scaling relations for the description of catalytic materials

TL;DR

This paper introduces a novel compressed sensing approach that predicts adsorption energies on complex catalyst surfaces using simple surface properties, outperforming traditional scaling relations.

Contribution

It develops a general, accurate method for predicting adsorption energies across various catalysts using non-linear descriptors derived from DFT calculations.

Findings

Outperforms previous scaling relations in accuracy.

Works well across pure and alloy metal surfaces.

Predicts stable surface geometries from a single DFT calculation.

Abstract

Computational screening for new and improved catalyst materials relies on accurate and low-cost predictions of key parameters such as adsorption energies. Here, we use recently developed compressed sensing methods to identify descriptors whose predictive power extends over a wide range of adsorbates, multi-metallic transition metal surfaces and facets. The descriptors are expressed as non-linear functions of intrinsic properties of the clean catalyst surface, e.g. coordination numbers, d-band moments and density of states at the Fermi level. From a single density-functional theory calculation of these properties, we predict adsorption energies at all potential surface sites, and thereby also the most stable geometry. Compared to previous approaches such as scaling relations, we find our approach to be both more general and more accurate for the prediction of adsorption energies on alloys with mixed-metal surfaces, already when based on training data including only pure metals. This accuracy can be systematically improved by adding also alloy adsorption energies to the training data.

Figures

Tables

Table 1. Table 1: Primary features used for the feature construction.

Class	Name	Abbreviation
Atomic	Pauling electronegativity	PE
	Ionization potential	IP
	Electron affinity	EA
Bulk	fcc nearest neighbor distance	bulknnd
	Radius of $d$-orbitals	$r_d$
	Coupling matrix element squared	$V_{\mathrm{ad}}^2$
Surface	Work function	$W$
Site	Number of atoms in ensemble	siteno
	Coordination number	CN
	Nearest neighbor distance	sitennd
	$d$-band center	${\epsilon }_d$
	$d$-band width	$W_d$
	$d$-band skewness	$S_d$
	$d$-band kurtosis	$K_d$
	$d$-band filling	$f_d$
	$sp$-band filling	$f_{sp}$
	Density of $d$-states at Fermi level	DOSd
	Density of $sp$-states at Fermi level	DOSsp

Equations

$$\begin{array}{rl}{P}^j=& c_1^j\cdot \left(\frac{V_{\mathrm{ad}}^2}{\mathrm{PE}}-\sin ({\mathrm{site}}_{\mathrm{no}})+\frac{\mathrm{IP}\cdot \sin ({\epsilon }_d)}{{\epsilon }_d}\right)\\ +& c_2^j\cdot \left(\mathrm{PE}\cdot f_d-\frac{f_{sp}\cdot \log ({\mathrm{DOS}}_{sp})\cdot ({\epsilon }_d-W)}{V_{\mathrm{ad}}^2}\right)\\ +& c_0^j,\end{array}$$

Full text

{tocentry}

1 Introduction

Surface catalysis has an enormous impact on the environment and our society’s health and prosperity 1. However, the reliable description of catalytic properties and the prediction of what materials may be even better catalysts than what we know today is still weak. This is due to the highly nonlinear and intricate relationship between the ”catalyst material” (the static material that is introduced before the catalytic process is running) and the strongly kinetically controlled surface reactions at realistic conditions 2, 3. Simply speaking, the basic understanding of heterogeneous catalysis is given by the Sabatier principle and the Brønsted-Evans-Polanyi (BEP) relation 4. The first states that there is an optimum adsorption strength for which the reactants bind strong enough to allow for adsorption and dissociation into reaction intermediates, but weak enough to allow for consecutive desorption of products. In turn, the BEP concept tells that the energy barriers of the chemical reactions scale approximately linearly with the adsorption energies of the molecules.

In consequence, the reliable prediction of adsorption energies is a key element of any theoretical description and search for new catalyst materials. For this, we here present a data-driven approach that does not start from a specific physical model, e.g. the tight-binding description of chemical bonding, but accepts that the intricacy of processes that cooperate or compete in materials properties may not necessarily be describable by a closed physical equation. This has been described as the fourth paradigm of materials science 5. Previous data-driven approaches to the prediction of adsorption energies 6, 7, 8, 9, 10 have exploited an approximately linear correlation between the adsorption energies of certain adsorbates on pure transition metal (TM) surfaces (scaling relations 11) to extend the data-driven predictions of one or two species to other adsorbates involved in the reaction. We instead directly learn the adsorption energies of a whole range of atoms and molecules at all potential adsorption sites (thereby also the most stable site) only from properties of the clean surface. This allows us to go beyond scaling relations and the often unfulfilled assumptions tied to this particular physical model. It furthermore opens the perspective of directly searching for outliers to scaling relations, which can be highly interesting catalyst materials missed by the standard approach 12, 13, 14.

The method for identifying the key descriptive parameters is the recently developed compressed sensing method SISSO 15 (sure independence screening and sparsifying operator), which enables us to identify the best multidimensional descriptor out of an immensity of candidates (billions). Our descriptors are more general and less costly to use than previous approaches and allow for making predictions for a huge number of surfaces including both multi-metallics and various facets. Through BEP relations our approach can also describe chemical reactions, as well as the diffusion of atoms and molecules at the surface visiting metastable adsorption sites.

2 Computational details

2.1 Density Functional Theory

The data sets employed in the present work were obtained from plane-wave density-functional theory (DFT) calculations (Quantum ESPRESSO code 16) using the van der Waals-corrected BEEF-vdW exchange-correlation functional 17. The larger data set consists of adsorption energies of atomic and molecular adsorbates (C, CH, CO, H, O, OH) on the stepped fcc(211) facets of nine TMs (Ni, Cu, Ru, Rh, Pd, Ag, Ir, Pt and Au) and selected single-atom (SA) and AB bimetallic alloys. The metals were modeled in fcc stacking using a $(1\times 3)$ ($(1\times 2)$ for AB alloys) supercell with fifteen metal layers. The considered adsorption sites are illustrated in Fig. 1(a) and cover both high-symmetry terrace and step sites. For each adsorbate, all adsorption sites that correspond to local minima on the potential energy surface were included. The SA alloys 18, 13, 14 were constructed by replacing one metal atom at the step with a different metal (see Fig. 1(b)). Specifically, we considered Ag@Cu (Ag atom in Cu surface), Pt@Rh, Pd@Ir and Au@Ni. For the AB L10 alloys 19 (AgPd, IrRu, PtRh and AgAu) the considered surface termination is depicted in Fig. 1(c). The total numbers of adsorption energies in the data set are 344 (metals), 281 (SA alloys), and 259 (AB alloys). The smaller facets data set consists of one adsorption site for each adsorbate on the fcc(111), (110) and (100) facets of the nine TMs, leading to a total of 54 adsorption energies on each facet. All data sets are compiled in Supplementary Sec. S1 together with further computational details.

2.2 Compressed sensing

The SISSO method 15 employed for descriptor identification makes the ansatz that the properties of interest $P^{j}_{1},...,P^{j}_{N}\in\mathbb{R}$, (in this case a vector of $N$ adsorption energies of adsorbate $j$) can be expressed as linear functions of candidate features ${\mathbf{d}}_{1},...,{\mathbf{d}}_{M}\in\mathbb{R}^{N}$, where the features are constructed as non-linear functions of user-defined primary features (see below). SISSO identifies the few best features (the number of which corresponds to the dimensionality of the descriptor) out of immense feature spaces by use of the sparsifying $\ell_{0}$ constraint. This is carried out in a smaller feature subspace selected by a screening procedure (sure independence screening (SIS)). The size of the subspace is equal to a user-defined SIS value times the dimension of the descriptor.

In this work we make use of multi-task learning 20 to identify common descriptors for the adsorption energies of several different adsorbates simultaneously. That is, the identified features are constrained to be identical for every adsorbate, while the fitting coefficients are allowed to vary between the adsorbates. We find that multi-task learning gives a better predictive performance compared to the identification of separate descriptors for each adsorbate (see Supplementary Fig. S2). We consider two hyperparameters in the SISSO method: the dimension of the descriptor and the feature space rung (see below) as well as the SIS value (see Supplementary Sec. S2), that we fix for the current application through a validation data set (see below).

2.3 Primary features

The decision which primary features to use as input for the feature construction is crucial for the predictive performance of the resulting descriptors. Inspired by previous studies 21, 22, 23, 24, 25, 26, 27, 6, 7 we consider four classes of primary features (see Table 1) related to the metal atom, the metal bulk, the metal surface and the metal adsorption site. For pure metals the primary features of the site class were calculated as averages over the metal atoms making up the site ensemble, while for alloys this was the case for the primary features of all classes. We note that the consideration of fixed adsorption sites as well as the averaging over the site ensemble is an approximation. It may break down in case of surface reconstruction or any other appearance of new adsorption motifs that were not accounted for in the calculation of the primary features. Further details regarding the primary features and all data are given in Supplementary Sec. S1.

2.4 Feature construction

As discussed above, candidate features are constructed as non-linear functions of the primary features. In the SISSO method this is achieved in practice by applying algebraic/functional operators such as addition, multiplication, exponentials, powers, roots etc. to the features 15. A full list of the used operators can be found in Supplementary Sec. S2. Arbitrarily large feature spaces can be constructed by iteratively applying these operators to the already generated features. The starting point $\mathbf{\Phi}_{0}$ corresponds to the 18 primary features listed in Table 1. We consider up to three iterations, generating thereby the feature spaces $\mathbf{\Phi}_{1}$, $\mathbf{\Phi}_{2}$ and $\mathbf{\Phi}_{3}$. Note that a given feature space $\mathbf{\Phi}_{n}$ contains also all of the lower rung feature spaces. The $\mathbf{\Phi}_{1}$ and $\mathbf{\Phi}_{2}$ feature spaces are still comparatively small. They consist of 783 and about 106 features, respectively. For the third iteration generating $\mathbf{\Phi}_{3}$ the approach we chose consisted in carrying out two rounds of feature construction and descriptor identification, each for a subset of only 16 out of the 18 primary features, in order to limit the $\mathbf{\Phi}_{3}$ feature space of each round to a tractable value of about $10^{11}$. In the first round the skewness and kurtosis of the $d$-band were excluded, since the higher order $d$-band moments are expected to be less important than the lower order moments. Among the identified best descriptors (i.e. with lowest validation errors, see below) of the first round, two primary features of the site class never appeared, namely the nearest neighbor distance and the density of $d$-states at the Fermi level. In the second round these two primary features were then excluded, while the skewness and kurtosis of the $d$-band were re-included. At every dimension the best performing $\mathbf{\Phi}_{3}$ descriptor originated from the first round and therefore only the results of the first round are presented below.

3 Results and discussion

3.1 Scaling relations

We begin by evaluating the performance of prevalent scaling relations for predicting adsorption energies on SA and AB alloys. In Fig. 2(a) and (b) we show two examples of scaling relations constructed by linear fits to the DFT-calculated adsorption energies on the pure TMs (black stars). Corresponding explicitly calculated adsorption energies on SA and AB alloys are also indicated by colored stars. While many bimetallics are well described by the linear scaling relations, there are also a number of serious outliers. Some systems with particularly large prediction errors of the order of 1 eV are highlighted. They typically contain mixed-metal sites made up of metals with very different reactivity towards O (e.g. Cu and Ag) or C (e.g. Ag and Pd). This poor performance of scaling relations derives from their calculation of the descriptors at one specific site, which fails to account for the variation in metal composition of the many other sites on the alloy surface. This issue likely occurs most severely for the considered thermochemical scaling relations. BEP relations for activation energies, in contrast, are more local in the sense that often both the transition state and the initial and final reaction intermediates coordinate (or can be chosen to coordinate) to the same metal atoms at the considered site ensemble. Correspondingly, BEP relations are typically found to exhibit significantly lower errors than thermochemical scaling relations even for the pure metals 28.

We note that an alternative scaling-relation-based approach is to calculate all potential adsorption sites for the descriptors on a mixed-metal surface and then to consider only the most stable adsorption sites 29. However, at concomitantly increased screening costs this still does not alleviate the problem since different adsorbates (e.g. O and OH) generally adsorb to different site types (e.g. O typically prefers higher coordinated sites than OH). As a consequence, the metal composition of the preferred sites could be different. We will come back to this point in the discussion of the compressed sensing results. In addition, not only the most stable sites, but also metastable sites missed by this approach, can get populated at higher coverages and then play an important role in the catalytic pathway 28.

3.2 Descriptor identification

The demonstrated failure of scaling relations to predict accurate adsorption energies on alloys with mixed-metal surfaces, as well as the high cost associated with the calculation of two or more adsorption energies on each alloy to be screened, emphasizes the need for new, accurate and low-cost descriptors for computational screening. In Fig. 3(a) we compare the performance of scaling relations to new descriptors identified by SISSO in terms of the root-mean-square error (RMSE) on training and validation data sets. We define the best descriptor as the descriptor that achieves the lowest RMSE on the validation data set. In the calculation of the RMSE the same weight is given to every adsorbate considered in the multi-task learning irrespective of how many data points exist for the adsorbate (see Supplementary Table S1). The training data consists exclusively of the adsorption energies on the pure metals and the validation data consists of 50% of each of the SA and AB alloys data. SISSO data are shown for 1D-8D descriptors identified from each considered feature space $\mathbf{\Phi}_{n}$. For each case a number of SIS values have also been tested (see Supplementary Fig. S4) and the best descriptor (i.e. the descriptor with the lowest RMSE on the validation data set) is shown. As expected, the SISSO training errors systematically decrease when increasing either the complexity and size of the feature space (larger $n$) or the dimensionality of the descriptor. The validation errors show the same trend, but the errors level out around the 5D to 8D descriptor depending on the rung of the used feature space. We would expect the validation errors to increase again at even higher dimensions. However, such higher dimensions are outside of the scope of the present study since the leveling out of the validation errors suggests that going beyond 8D is unlikely to result in descriptors with lower validation errors.

The 5D to 8D descriptors of $\mathbf{\Phi}_{3}$ all have very low validation errors, differing from each other by only about ten meV. Likely, there is no statistical significant difference in their performance. A detailed statistical analysis to derive error bars is outside of the scope of this work, which aims at a simple comparison to scaling relations. The latter are usually derived based on only one fixed training data set considering only the pure metals and we therefore follow the same approach here. In the absence of error bars we choose the descriptor with the lowest observed validation RMSE (of 0.15 eV) as our best descriptor, i.e. the 8D, $\mathbf{\Phi}_{3}$ descriptor. This descriptor is significantly better than scaling relations, for which the validation RMSE is 0.28 eV (horizontal dashed line). In fact, already the 2D descriptor of $\mathbf{\Phi}_{3}$ (with a validation RMSE of 0.22 eV) performs better than scaling relations. The best descriptor among the primary features (the SISSO 1D, $\mathbf{\Phi}_{0}$ descriptor) is found to be the $d$-band center, i.e. SISSO identifies the physics that has already been discovered in form of the $d$-band model more than twenty years ago 23.

A comparison of the performance of the $d$-band center, scaling relations and the best (8D, $\mathbf{\Phi}_{3}$) SISSO descriptor on the test data set (the remaining 50% of the SA and AB alloys data) is shown in Fig. 3(b). The RMSEs are: $d$-band center: 0.37 eV, scaling relations: 0.28 eV, SISSO: 0.15 eV, and the maximum absolute errors (maxAEs) are: $d$-band center: 1.31 eV, scaling relations: 1.43 eV, SISSO: 0.61 eV. Here, the maxAE is the maximal error observed for any of the adsorbates considered in the multi-task learning. As already observed for the validation data, the test results thus demonstrate the great improvement of the new SISSO descriptor compared to previous approaches. In Supplementary Table S7 we provide additional information about the largest deviations between calculated and SISSO-predicted adsorption energies for the alloys. In general, the C adsorption energies are the most difficult to predict since they vary much more over the TM series than e.g. the H adsorption energies. In addition, the most difficult alloys to predict adsorption energies for are those that combine a more noble and a more reactive metal such as Au@Ni or Ag@Cu. The maximum absolute error of 0.61 eV is found for the OH adsorption energy on top of the Au atom in the Au@Ni alloy. The embedding of the larger and more noble Au atom in the smaller lattice constant and more reactive Ni surface (see Fig. 1(b)) probably provides an adsorption site that is both geometrically and electronically very different from everything else in the data set and thus harder to predict.

3.3 Transferability of descriptors

To further test the predictive performance of the best (8D, $\mathbf{\Phi}_{3}$) SISSO descriptor we show in Fig. 3(c) the error distribution for the prediction of adsorption energies on three new facets, the (111), (110), and (100) facets. Some sites, in particular those found on the (111) facet, are very similar to sites found on the (211) facet in the training data. However, the (110) and (100) facets contain sites that are very different, albeit of similar coordination numbers (7 – 9.5), compared to the (211) facet sites. It is seen that the descriptor with the best performance for (211) facets performs very well for the (111) facet (RMSE of 0.17 eV), but significantly worse for the (110) and (100) facets (RMSEs of 0.29 eV and 0.30 eV, respectively). This shows a well-known limitation of compressed sensing (and machine learning) methodologies, namely that a good transferability cannot be expected a priori for cases not previously encountered in the training or validation data. However, to put this in perspective, we note that this ”poor” RMSE for the other facets is essentially of the same level as the RMSEs obtained for the widely used scaling relations in the first place.

In order to identify a descriptor that has a good compromise between accuracy for alloys and facets, we include 50% of the facets data in the validation data (see Supplementary Fig. S5). The new best descriptor is found for the hyper parameters 8D, $\mathbf{\Phi}_{1}$. It is interesting to note that since this is a $\mathbf{\Phi}_{1}$ descriptor, the functional form of the features is much less complex than for the $\mathbf{\Phi}_{3}$ descriptor optimized for the alloys alone. This suggests that a less complex mathematical form is required for a descriptor that is transferable across both alloys and active site motifs.

The RMSEs (maxAEs) of the (8D, $\mathbf{\Phi}_{1}$) descriptor for the remaining 50% of the facets data are found to be: (111): 0.17 eV (0.44 eV), (110): 0.21 eV (0.59 eV), (100): 0.24 eV (0.48 eV). These very moderate errors show that it is possible to identify a descriptor with a good predictive performance for a wide range of structural motifs as exemplified by the low-index fcc facets. The improved performance on the facets data sets comes at the very moderate expense of increasing the RMSEs for the alloys to 0.18 eV compared to 0.15 eV before. We therefore suggest the (8D, $\mathbf{\Phi}_{1}$) SISSO descriptor for cases where simultaneous screening of alloys and a wide range of active site motifs is desired. It should be emphasized though that in general a good performance can only be expected for active sites types that resemble to some extent those types that the descriptor was optimized for. Likely, more varied training and validation data will be required to identify a descriptor that would work for very different active site motifs such as kinks, vacancies and adatoms.

3.4 Composition of descriptors

Having confirmed the predictive performance of the identified descriptors, we now move on to discuss their composition. As an example of a descriptor identified by SISSO, we give in eq. 1 the best 2D descriptor of $\mathbf{\Phi}_{3}$ from Fig. 3(a) (alloy validation data set):

[TABLE]

where $P^{j}$ is an adsorption energy of adsorbate $j$ and the primary features entering the descriptor are evaluated for the material/site combination relevant for $P^{j}$. The $d$-band center and $sp$-band properties such as the Pauling electronegativity are identified as highly important primary features. This was also found in previous studies employing artificial neural networks 6, 7, but is here expressed in an explicit non-linear functional form owing to the compressed sensing methodology. Among the remaining primary features entering the descriptor we especially highlight the DOS at the Fermi level (here of the $sp$-band), which is a feature not considered in these previous studies, even though its importance for the reactivity of TM surfaces was discussed already more than 30 years ago by Yang and Parr 22.

An overview of the identified descriptors of each dimension and rung for both the alloy validation data set and the combined alloy and facets validation data set together with the fitting coefficients $c_{i}^{j}$ is given in Supplementary Table S8 and S9.

3.5 Enlarging the training data set

The predictive performance of the identified descriptors for alloy screening is already impressive, given that no explicit information on alloys was given in the training data. However, a further advantage of data-driven approaches is that the learning can be systematically improved by enlarging the training data set. In contrast, the rigid format of linear scaling relations does not allow for significant improvements, even if fitting also to the alloys data, as evidenced by the scattering of the alloy data points around the fitted line in Fig. 2(a) and (b).

To provide a simple estimate of the learning improvement possible when including also alloys in the training data, we identify a new SISSO descriptor (see Supplementary Table S10) based on the pooled metal and alloys data sets, but excluding the 23 DFT-calculated alloy data points (colored stars) shown in Fig. 2(a) and (b). For this we use the hyperparameters that were found to be best for alloys (8D, $\mathbf{\Phi}_{3}$). The colored dots in Fig. 2(a) and (b) show the SISSO predictions for the data points left out in the training. Already a visual inspection reveals that the agreement is very good. The maxAE (of 0.52 eV) is found for the OH adsorption energy of the dark green point in Fig. 2(a), which corresponds to OH adsorption on top of the Au atom in the Au@Ni alloy. Note that exactly this data point is also the maxAE (with the slightly larger value of 0.61 eV) for the descriptor trained only on the pure metals. The RMSE over the 23 predicted alloy data points in Fig. 2(a) and (b) decreases from 0.23 eV (training on pure metals only) to 0.18 eV (training on pooled data set).

Overall, the good agreement between the DFT-calculated values and the SISSO predictions shows that our models have the required accuracy to systematically search for outliers to scaling relations. In addition, our approach is computationally cheap enough to allow for the screening of immense alloy spaces. For this, we will next present a simple example.

3.6 A first screening example

In the following we make use of a SISSO descriptor (see Supplementary Table S11) identified using the hyperparameters (8D, $\mathbf{\Phi}_{3}$) and the entire pooled metals and alloys data sets for training. We predict the adsorption energies for the adsorbates and sites considered in Fig. 2(a) and (b) on the additional 32 possible AB alloys (those that were not explicitly calculated by DFT) as shown with the gray dots. Similar to the explicitly DFT-calculated alloy data points, there is a considerable scatter around the scaling relation lines. This shows that there indeed exist many materials with potentially interesting catalytic properties, which would be missed by a scaling-relation-based screening approach.

A particularly interesting perspective for catalyst screening is to be able to search directly for candidate materials that break scaling relations. For many reactions, the incentive would be to break scaling relations in a desired way, since it has been suggested that scaling relations impose an upper limit to the possible catalyst activity. For example, it is known that for the oxygen evolution reaction it would be desirable to find a material where O is destabilized relative to OOH 32, and for electrochemical CO2 reduction it is desirable to destabilize CO relative to CHO 33. For other reactions, where optimum catalytic activity has hitherto been exclusively formulated in terms of singular descriptors, the interest would be to evaluate the effect of scatter in the binding of other important intermediates that hitherto has been assumed as fixed through scaling relations. For the oxygen reduction reaction (ORR), optimum catalytic activity has for instance been associated with an optimum oxygen adsorption energy 34, while for selective ethanol synthesis both the carbon and the oxygen adsorption energy must be simultaneously optimized 31. The activity at the top of this theoretical volcano curve is then independent of e.g. the OH (ORR) or CH (ethanol) adsorption energy, as the latter are connected to the optimum O or C adsorption energy through a scaling relation, respectively.

In Fig. 2(c) and (d) we specifically check on the scatter by showing histograms of the predicted adsorption energies for (c) O and OH and (d) C and CH at all potential adsorption sites of all 36 AB alloys. There are more predicted points for OH (around 1000) than for the other adsorbates (around 400), since OH adsorption can take place at more site types, i.e. also at top and bridge sites. The black shaded areas highlight in Fig. 2(c) the ORR ”volcano optimal” O adsorption energy on (111) facets 30 and in Fig. 2(d) the optimal C adsorption energy on (211) facets for selective ethanol synthesis 31. For this simple screening example, we will assume that only the most stable adsorption site of a given adsorbate plays a catalytic role, keeping in mind that in reality also less stable (meta-stable) sites could get populated at higher coverages. The SISSO approach directly gives us the energetics for the most stable and all meta-stable sites, so in general we are not limited to considering only most stable sites. Since the ORR volcano was developed for (111) facets, we search for materials for which the most stable O adsorption energy among the (111)-like (terrace) sites of the (211) facet falls within the desired range. This results in three candidate materials: PdPt, AgPd, and AgPt. The latter material is highlighted in Fig. 2(c), since it has an OH adsorption energy on its most stable bridge1-$t$ site that is 0.23 eV lower than the value that would be predicted for this site from scaling relations (indicated by the black arrow) based on the O adsorption energy on its most stable hcp-$t$ site. The opposite behavior (a lower O adsorption energy relative to OH) is seen for the material shown with the green dot (PdPt). The cause of this breaking of scaling relations is thereby the slightly different oxygenate adsorption energy for Pt and Ag, and the fact that the preferred adsorption sites for OH and O have a different composition of Pt and Ag. A similar breaking of scaling relations is observed for the ethanol synthesis example. Here, SISSO recovers the scatter in the most stable CH adsorption energies for materials (RuAg, RhAu, RuAu, IrAu, and CuIr) that all have about the same most stable ”optimum” C adsorption energy. This scatter shows the extent to which it is possible to tune the CH adsorption energy independently of the C adsorption energy and thereby further tailor the catalytic activity.

3.7 A high-throughput screening perspective

It should be noted that the moderate breaking of scaling relations observed in the simple examples from the previous section is related to the consideration of only a handful of ”near-optimal” materials (out of a total of only 36 considered materials) and in particular the consideration of only most stable adsorption sites. A full assessment of the extent to which scaling relations can be broken on alloys with mixed-metal surfaces would only be revealed by a full high-throughput screening of hundreds of thousands of materials that is beyond the scope of the present study. Such a high-throughput screening could also involve the evaluation of a microkinetic model for each catalyst material that takes into account all possible adsorption sites of every adsorbate, as well as kinetic barriers for reaction and diffusion steps through BEP relations. If such a microkinetic model was initially carried out within the simplifying mean-field approximation 35, the cost of its evaluation would still only be a negligible fraction of the (already small) cost of carrying out a DFT calculation of the primary features of the clean catalyst surface for the descriptor evaluation. Once a selection of promising catalyst materials had been identified, a next step could then be the evaluation of a more thorough microkinetic model from e.g. a kinetic Monte Carlo simulation 28, possibly taking into account also lateral interactions between the adsorbates through cluster expansion methods 36. A full assessment of identified promising catalyst materials would ultimately also need to take into account other aspects such as bulk and surface segregation stability under realistic surface coverages for the chemical reaction and reaction conditions of interest, as well as stability against metal stripping in electrocatalysis applications.

4 Conclusions

In summary, we have used compressed sensing to identify new and better descriptors that allow to predict adsorption energies for a whole range of atoms and molecules at all potential surface sites of TMs and bimetallics formed of TMs. The descriptors can be obtained from a single DFT calculation of the clean surface and their predictive power extends over both multi-metallics and various surface facets. Importantly, this enables low-cost catalyst screening not only in materials, but also in active site space 37 with unprecedented accuracy. With respect to materials, the thereby enabled systematic identification and analysis of outliers to traditional scaling relation energetics seems particularly promising. With respect to active sites, the availability of energetic data for a wide range of site types paves the way to actively embrace the uncertainty in surface structure and composition of working catalysts.

{suppinfo}

Supporting Information. Additional DFT and SISSO computational details, multi-task learning versus learning of separate descriptors, overview of largest prediction errors, all scaling relation plots, tested SIS values, hyperparameter testing with facets data in validation data set, and all identified descriptors. This material is available free of charge via the Internet at http://pubs.acs.org.

{acknowledgement}

This project has received funding from the European Unions Horizon 2020 research and innovation program under grant agreement No. 676580, The NOMAD Laboratory, a European Center of Excellence. The authors also gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). We also thank Luca Ghiringhelli for helpful discussions and a careful proofreading of the manuscript. Runhai Ouyang has written the multi-task SISSO code 20 used in the present work and provided technical assistance.