Introduction

The chemisorption of molecular and atomic species on solid-state material surfaces is a central concept in chemistry, physics, and material science. The ability to identify the key surface and adsorbate properties that govern the chemisorption strength is crucial in understanding chemical processes in surface science1,2,3,4,5,6,7,8,9. In heterogeneous catalysis, it is recognized that the bond strength between intermediates in a reaction and the surface provide decisive information about the catalytic activity and selectivity10,11. Efficient methods that can accurately predict adsorption energies on complex catalyst formulations based on limited input data are essential to push the field forward. Despite decades of focused research, interpretable modeling methods capable of predicting adsorption energies with the desired accuracy on complex catalyst structures and compositions with active site resolution are still lacking.

Elemental transition metals, transition metal alloys, and intermetallics, constitute an important class of materials in catalysis. Significant research has been performed on these systems and several simplifying models have been proposed to describe interactions on these surfaces. The d-band model by Hammer & Nørskov1 is a particularly successful model, because it effectively correlates electronic-structure features of the system obtained before interaction with the chemisorption strength. The model is based on observations from tight-binding models like the Newns-Anderson model12,13, where certain systematic perturbations in the adsorbate electronic structure, and hence the bond strength, are identified as important for describing chemisorption trends and then linked to relevant electronic features in the surface. Transition metals are defined by their electronic similarities, thus allowing for a specific distinction between interactions with metal sp-electrons and metal d-electrons. Due to the delocalized nature of metallic sp-states, it is reasonable to assume that the interaction between a given adsorbate and these states in the metal is essentially constant, and therefore any observed variation in the bond strength is assigned entirely to changes in the metal d-electronic states.

Most commonly, the position of the d-band center with respect to the Fermi-level has been applied to explain adsorption trends on pure transition metals and some alloys10,12,14,15,16,17,18,19. However, especially for alloys and intermetallics, the d-band model is known to have shortcomings20,21. The shortcomings mainly arise because the center of the d-band carries no information of the band dispersion, and hence d-band center-based models, which all rely on surface localized features, lack the ability to fully account for asymmetries and distortions in the electronic structure introduced by alloying22,23,24,25,26,27,28. One model based entirely on changes in coordination-number (CN) has been shown to effectively describe the geometric (structural) impact on the chemisorption strength on pure metals22,23. However, since CN-based models ignore the inherent electronic structure of the metal and adsorbate and are independent of the chemical environment of the site, they are not transferable to complex alloys, which is essential when screening for promising multi-metallic transition metal-based catalysts. Xin et al.24 introduced an approach based on an orbitalwise coordination number, where the coordination-number of the active site and the orbital overlap with neighboring atoms were considered. Regardless of this added feature, the physical link between composition, CN, and chemisorption remains unclear. Roling et al.25, proposed a model based on metal-metal coordination (known as the α-parameter scheme) that addresses the impact of both local coordination as well as changes in the local composition. The model accurately describes the site stability and the chemisorption energy of various adsorbates on transition metals surfaces, alloys and alloy nanoparticles. However, linking the adsorption energy entirely to variations in the surface stability fails to account for effects induced explicitly by the adsorbate. Most recent developments in predicting accurate chemisorption energies are based on data-driven approaches29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57, signifying the immense challenges associated with the combinatorial complexity encountered once one extends the materials space beyond simple mono-metallic or dilute bimetallic systems. The use of multi-component catalysts, intermetallic alloys, and high entropy alloys offers a promising path towards designing formulations based on abundant low-cost materials with ideal catalytic properties58,59,60,61,62,63 but it also highlights the need for efficient methods to guide us.

In this paper, we derive a model based entirely on physical arguments that accurately describes the adsorption energies of O, N, CH, and Li (MAE of 0.13 eV) on bi- and tri-metallic surface and subsurface alloys consisting of Au, Ag, Cu, and Pt host metals alloyed with 3d, 4d, and 5d metals. This physical model is based on the 1st and 2nd moments of the d-band, as well as the d-band filling of the atoms in the alloy: variables that can readily be computed and tabulated. Model parameters for a specific adsorbate and adsorption site are subsequently optimized based on a limited set of DFT training data. As opposed to previous methods, our approach takes into consideration adsorbate-induced changes to the adsorption site and how these changes interact with variation in the chemical environment around the site. These variations explicitly give rise to deviations from the typical linear behavior of the adsorption energy with electronic structure descriptors such as the d-band center. The discrepancies in the total energies are small, of the order ± 1 eV, however still significant for chemical processes. The model is derived in the limit of small perturbations and, therefore, valid only when minor changes to a given site motif are considered. Nevertheless, by proper calibration of a representative site motif, the method can be fully generalized for the entire series of transition metals as well as for multi-component systems. The introduced method offers a promising avenue towards the discovery of improved catalysts consisting of complex multi-metal materials.

Results and discussion

Relation between chemisorption strength and surface electronic structure

In order to further deepen our understanding of surface interactions, we first critically examine existing models. In particular, we revisit the seminal work of Nørskov and co-workers on the d-band model which has widespread use in surface science and catalysis3,5,21. Regardless of its success, the d-band model also has well-known limitations in describing, e.g., noble metals, bimetallic alloys and multi-component intermetallics such as high entropy alloys17,20,64. We believe that the short-comings of the d-band model originate from the constraints under which the model was derived. We intend to lessen some of these constraints to arrive at an expression that more accurately describes a broader range of material surfaces than the d-band model while still using characteristic properties of the d-band as features or descriptors. In particular, we shall show how to account for the influence of the adsorbate on the substrate, which results in a crucial fine-tuning of the interaction model.

In general, the adsorption process occurs when the electronic states of an atom, ion, or molecule encounter the electronic states on the material surface. If the repulsive and attractive interaction contributions are overall favorable a chemical bond is formed. At the electronic level, the energy levels of both the adsorbate and surface will be perturbed (rehybridized) by the mutual overlap of electronic states and the exchange of valence electrons. In most models used to approximate chemical bonding, only the perturbation of the electronic states of the adsorbate due to interaction with the surface is considered, whereas the perturbation of the surface states due to interaction with the adsorbate is treated as negligible3,6. Given that many surfaces tend to reconstruct upon interaction with certain intermediates, metallic spin-states gets quenched, and surface segregation is induced, one can conclude that a significant perturbation of surface electronic states is induced by the adsorbate interaction65,66,67,68,69,70,71.

Exactly how the interaction with a continuum of surface electronic states affects a single adsorbate energy level can be described by applying the Newns-Andersen model. Transition metals and transition metal alloys share certain commonalities, all having delocalized sp-electronic states with a broad structureless energy distribution (see the sp-states in Fig. 1a) and localized d-states with a much narrower distribution (see the d-states in Fig. 1a). It can be shown that any adsorbate electronic energy level fully integrated in the band representing the sp-electronic states when interacting leads to a broadened single resonance state (renormalized state, εa), with an energy lower than the original gas phase adsorbate level and located below the Fermi-level primarily defined by the transition metal. On the other hand, in the Newns-Anderson model, the interaction between a single adsorbate energy level and the localized d-states, centered at εd, generally leads to two distinct solutions above and below the d-band. The lowest of the two energy states is the bonding state, which for transition metal interaction lies below the Fermi-level, whereas, the higher lying state called the antibonding state in many cases has density both above and below the Fermi-level. In models based on one-electron energies like the Newns-Anderson and the d-band model, the total energy is approximated by the sum of interactions with the sp- and d-states separately, allowing us to write the adsorption energy of adsorbate A as1:

$${{\Delta }}{E}^{A}={{\Delta }}{E}_{sp}^{A}+{{\Delta }}{E}_{d}^{A}$$
(1)

The \({{\Delta }}{E}_{sp}^{A}\) contribution from the sp-electrons is usually large and attractive, whereas the \({{\Delta }}{E}_{d}^{A}\) contribution from the interaction with the transition metal d-electrons is smaller and weakens as we move from left to right in the 3d, 4d, and 5d transition metal series. Because of the broad and structureless distribution of the sp-electronic states for transition metals and transition metal alloys, the interaction can be presumed approximately constant and, consequently, the difference in adsorbate-surface interaction between two slightly dissimilar surfaces will be completely governed by the difference in interaction between the d-states and the renormalized absorbate state. In the following, we shall use ε+ and ε− to define the bonding and antibonding states, respectively (Fig. 1b, c). These states result from the interaction of the renormalized state centered at εa with the metal d-states. The narrow distribution of both the renormalized as well as the metal d-states allows for a discrete treatment of the energy levels as illustrated in Fig. 1c and, in this approximation, we can write the expression for ε+ and ε− as1,6:

$${\varepsilon }_{\pm }=\frac{{\varepsilon }_{a}+{\varepsilon }_{d}}{2}\mp \sqrt{{V}^{2}+{\left(\frac{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}{2}\right)}^{2}}-SV$$
(2)

which can be rearranged to:

$${\varepsilon }_{\pm }=\frac{{\varepsilon }_{a}+{\varepsilon }_{d}}{2}\mp \frac{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}{2}\sqrt{1+{\left(\frac{2V}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}\right)}^{2}}-SV$$
(3)

In these equations, S and V define the intersection and interaction between the states representing the adsorbate and the metal, called the overlap and coupling matrix elements, respectively. In the limit of weak coupling (\(| V| \ll \left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|\)), we can derive an approximate expression for the energy of the bonding and antibonding states by expanding the square root in Eq. (3) to second order in \(\frac{V}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}\) :

$${\varepsilon }_{+}\simeq {\varepsilon }_{a}-\frac{{V}^{2}}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}-SV$$
(4)
$${\varepsilon }_{-}\simeq {\varepsilon }_{d}+\frac{{V}^{2}}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}-SV$$
(5)

These approximate solutions to Eq. (2) diverge when εa is equal to εd. However, it is clear that there exists a finite solution to the posed problem for all values of εa and εd as Eq. (2) is fully solvable also when εa = εd. Embracing the perturbative nature of the derivation in Eqs. (4) and (5), we can write the reaction energy associated with the interaction between the states εa and εd as a sum of the occupied electronic energies. The state occupations are clearly linked to the second moment of the distributions and their position relative to the Fermi-level. However, assuming evenly distributed energies around the mean, we can approximate the cumulative energy of each state, i, as fiεi, where fi is the filling the state centered at εi. The electronic energy associated with the bond between a state centered at εa and states centered at εd then becomes:

$${{\Delta }}{E}_{d}^{A}={f}_{+}{\varepsilon }_{+}+{f}_{-}{\varepsilon }_{-}-{f}_{a}{\varepsilon }_{a}-{f}_{d}^{{{{\rm{site}}}}}{\varepsilon }_{d}$$
(6)

In the schematics in Fig. 1b, the renormalized state (fa) and the bonding states (f+) are both completely filled, i.e., fa = f+ = 1 such that:

$${{\Delta }}{E}_{d}^{A}={\varepsilon }_{+}+{f}_{-}{\varepsilon }_{-}-{\varepsilon }_{a}-{f}_{d}^{{{{\rm{site}}}}}{\varepsilon }_{d}$$
(7)

where f− and \({f}_{d}^{{{{\rm{site}}}}}\) are the filling of the antibonding and the metal d-states, respectively. Inserting the approximations from Eqs. (4) and (5) and rearranging terms, we arrive at the following:

$$\begin{array}{rlr}{{\Delta }}{E}_{d}^{A}&\simeq {\varepsilon }_{a}-\frac{{V}^{2}}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}-SV+{f}_{-}\left({\varepsilon }_{d}+\frac{{V}^{2}}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}-SV\right)-{\varepsilon }_{a}-{f}_{d}^{{{{\rm{site}}}}}{\varepsilon }_{d}&\\ &\simeq -(1-{f}_{-})\frac{{V}^{2}}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}-(1+{f}_{-})SV+({f}_{-}-{f}_{d}^{{{{\rm{site}}}}}){\varepsilon }_{d}\end{array}$$
(8)

Expanding \(\frac{1}{\left|{\varepsilon }_{d}-{\varepsilon }_{a}\right|}\) to first order, we finally obtain an approximate expression for the electronic energy contained in the adsorbate-surface bond in terms of electronic-structure parameters associated with the adsorbate and the model surface:

$$\begin{array}{ll}{{\Delta }}{E}_{d}^{A}\simeq -\left(\frac{(1-{f}_{-}){V}^{2}-{\varepsilon }_{a}^{2}({f}_{-}-{f}_{d}^{{{{\rm{site}}}}})}{{\varepsilon }_{a}^{2}}\right){\varepsilon }_{d}\\\qquad\qquad-(1+{f}_{-})SV+\frac{(1-{f}_{-}){V}^{2}}{{\varepsilon }_{a}}\end{array}$$
(9)

Clearly, the model expression in Eq. (9) contains information that reflects the adsorbate-surface interaction and, therefore, cannot readily be extracted from the separate systems. To completely eliminate any dependency on the composite system, it has been assumed in earlier derivations, that the electron density of the adsorbate post interaction follows the electron density of the clean surface, i.e., \({f}_{-}\approx {f}_{d}^{{{{\rm{site}}}}}=f\). Furthermore, if we use that S ∝ − V and introduce a factor of 2 to account for both electronic spin-channels, Eq. (8) reduces to the expression derived in the original d-band model1:

$${{\Delta }}{E}_{d}^{A}\simeq -2(1-f)\frac{{V}^{2}}{| {\varepsilon }_{d}-{\varepsilon }_{a}| }+2(1+f){\alpha }_{s}{V}^{2}$$
(10)

where αs is introduced as an adsorbate-specific proportionality constant with the unit energy−1.

Fig. 1: Interaction between an adsorbate and transition metal surface states.
figure 1

a Schematic diagram showing the change in the electronic states of the adsorbate upon interaction with the delocalized sp-states and localized d-states on the metal surface6. Upon interaction with the metal sp-states, the sharp atomic states of the gas phase are broadened into resonances and shifted down. b, c Subsequent interaction of these renormalized states with the localized d-states gives rise to a splitting into bonding (ε+) and antibonding (ε−) states below and above the initial adsorbate and surface d-states.

Full size image

To demonstrate the range of the two different models given by Eqs. (9) and (10), respectively, it is instructive to look at the coinage metals where the metal d-states are fully occupied, i.e., \({f}_{d}^{{{{\rm{site}}}}}=f=1\). Clearly, the first term in Eq. (10) vanishes and, hence, in this picture bonding is determined entirely by the second term representing the repulsion due to overlap between states3. If we assume that αs is constant for a given adsorbate, then the adsorption strength on the coinage metals should scale with V2. In Fig. 2a, we have plotted the DFT calculated binding energy of O, N, and CH in 3-fold fcc sites on the close-packed (111) surface of Ag and Au relative to Cu(111). Only the binding of atomic oxygen shows a clear correlation with V2 and, as we shall see, the failure to describe the binding of N and CH is a consequence of the assumption that \({f}_{-}\approx {f}_{d}^{{{{\rm{site}}}}}=f\). To illustrate why the adsorbates have different chemisorption trends on the coinage metals, we use atom-projected density of states plots for N and O as shown in Fig. 2b, c (see Supplementary Fig. 5 for CH). Evidently, the antibonding states formed between the adsorbates (N, O, and CH) and the coinage metals Cu, Ag, and Au are not completely filled and, hence, we cannot neglect the term that involves hybridization between the metal d-states and the renormalized state on the adsorbates.

Fig. 2: Chemisorption on pure transition metal surfaces.
figure 2

a DFT calculated adsorption energies of O, N and CH on Cu(111), Ag(111), and Au(111) as a function of the coupling matrix element (V2). Adsorption energies are shown relative to Cu. The coupling matrix elements are taken from Hammer et al.1,4. b, c p-projected density of states for O (red line) and N (blue line), respectively. The d-projected density of states for surface atoms on Cu(111), Ag(111), and Au(111) are shown in turquoise, brown, and yellow, respectively. The fractional filling (f−) of the resulting antibonding states are shown in each case (for details on how to determine f−, see Supplementary Note 1 and see Supplementary Figs. 1–4).

Full size image

The above results can be understood by taking a closer look at Eq. (9). Considering coinage metals, we have \({f}_{d}^{{{{\rm{site}}}}}=1\) and furthermore for a given adsorbate, f− stays approximately constant for the different metals (Fig. 2b, c). It is also reasonable to assume that the change in the numerator of the first term in, \((1-{f}_{-}){V}^{2}+{\varepsilon }_{a}^{2}({f}_{-}-{f}_{d}^{{{{\rm{site}}}}})\), is approximately constant for all the coinage metals. In addition, the variation in V2 will be small in comparison to changes in \({\varepsilon }_{a}^{2}\) of the denominator. Hence, the repulsive interaction between the adsorbate and the coinage metals Ag and Au relative to Cu, as given by the first term in the approximate expression in Eq. (9), is proportional to \(-\frac{\delta {\varepsilon }_{d}}{{\varepsilon }_{a}^{2}}\), where δεd is the d-band center on the metal relative to the d-band center on Cu. Sweeping from right to left in the periodic table (here from O via N to C), one effectively empties the atomic orbitals on the adsorbate, thereby pushing up the center of the energy distribution for that orbital relative to the Fermi-level. Therefore, as the orbital population decreases, εa approaches the Fermi-level. Since Ag has much lower lying d-states than Au (Fig. 2b, c), the repulsive interaction contribution from the first term of Eq. (9) relative to Cu will always be larger for Ag than for Au. Eventually, when εa gets close enough to the Fermi-level, the first term of Eq. (9) will dominate the interaction behavior on the metal as seen in Fig. 2a. Oxygen, having the lowest lying εa, has a close to linear adsorption behavior on the coinage metals with respect to V2, a result of the first term in Eq. (9) being small. In contrast, N and CH have εa closer to the Fermi level and thus non-negligible contributions from the first term in Eq. (9), hence, explaining the non-linear adsorption behavior with V2. From the above, it is clear that effects due to the adsorbate cannot be neglected in determining chemisorption trends. Therefore, we explore more generally under which circumstances adsorbate-induced perturbations of the adsorption site and its surroundings can affect surface interactions.

Determining adsorbate-induced local effects on the surface

In the following, we shall discard the assumption that the interaction between an adsorbate and a metal surface only leads to changes in the adsorbate states. Instead, we shall also account for adsorbate-induced perturbations of the electronic states on the metal adsorption site. Although it should be emphasized that the effect of the metal on the adsorbate is still considered the dominating term in bonding, we argue that the effect of the adsorbate on the metal site plays a non-negligible role in chemisorption. To support this, we have performed an electronic-structure analysis (shown in Supplementary Note 2 and Supplementary Fig. 6) demonstrating a substantial adsorbate-induced effect on the adsorption site as the electronic-structure parameters εd and Wd (the width of the d-band) change upon adsorption.

We will employ an analysis of the stability of the adsorption site to help quantify the induced effects on the site upon adsorption. It has been suggested experimentally, and shown theoretically, that changes in the adsorbate chemisorption energy to first approximation scale linearly with variations in the stability of the metal site (BEM)26,27,72,73,74,75. To establish a link between site stability and electronic structure, we shall compare the variation of BEM with variations in the d-band center, εd, of the adsorption site for different coordination numbers (CN) on several mono-metallic fcc transition metal surfaces (Cu, Au, Rh, Ir, and Pd). The local site morphology corresponding to different coordination numbers, ranging from three to nine, is accomplished by arranging adatoms on the (100), (111) and (211) surfaces (see Supplementary Note 3 and Roling et al.25).

In Fig. 3a, the site stability is plotted and, clearly, it correlates with changes in εd as well as with changes in the CN, i.e., δεd ∝ δCN ∝ δBEM. The site with the highest CN has the lowest lying εd, which leads to increased site stability and hence weaker interaction with the adsorbate. In addition, we have shown that application of strain (compression and expansion) in the surface leads to continuous changes in the CN, thus resulting in variations in BEM along the correlations introduced in Fig. 3a76. However, if we change the local chemical composition around a specific site, a relationship appears that deviates from the original. In Fig. 3b, the correlation between the adsorption energy of a metal-adsorbate complex (\(B{E}_{{{{\rm{Au-OH}}}}}^{* }\)) and the site stability (\(B{E}_{{{{\rm{Au}}}}}^{* }\)) is seen to deviate from the linear correlation found when varying the CN (blue circles in Fig. 3b). This difference can only be a result of the adsorbate influencing the local electronic structure at the adsorption site. In the following, we shall explain how one can account for changes in chemisorption caused by local electronic structural changes by explicitly studying small perturbations around a fixed adsorption site.

Fig. 3: Variation in site stability as a function of d-band center, coordination number and change in site environment.
figure 3

a Metal atom stability as a function of the change in d-band center (εd) and coordination number (CN) on fcc (100), (111), and (211) surfaces with varying number of surface adatoms. Linear fits to the points for each metal are shown with dashed lines. b The binding energy of metal-adsorbate complexes (ontop site adsorption) as a function of the binding energy of the metal atom for pure Au (blue points) and Au surface alloys (red squares). The blue points reflect changes in CN, whereas red squares designate changes in the local chemical environment72.

Full size image

In Fig. 4, we illustrate this, by progressively adding more complexity to our model. In Fig. 4a, we only consider how the electronic states in the substrate influence the electronic energy levels on the adsorbate, A. This lowest order approximation towards understanding surface interaction has been applied in most chemisorption models, including Newns-Anderson, the d-band model, the generalized coordination number model, and the α-parameter scheme1,12,13,22,23,25. In all cases, this approximation leads to models that only depend on parameters strictly associated with the metal surface, i.e. the d-density of states, the d-band center, a surface weighted coordination number, and the site stability, and they all rely on Eq. (1). If we allow the adsorbate to influence the surface atoms in the adsorption site (pink arrows in Fig. 4b), then, presumably, the induced electronic state changes in the surface will lead to changes in the surface stability that will depend on the properties of A and the character of the nearest neighboring atoms (black arrows). In principle, the term (II) accounts for changes due to additional effects, whereas (I) accounts for the standard contributions.

Fig. 4: Schematic illustration of adsorbate-surface interaction.
figure 4

The top panel shows a general representation of the adsorbate-surface interaction. In the bottom panel, a shows the zeroth order approximation to bonding which is based on the assumption that the adsorbate induces no changes in the surface electronic states and hence, that interaction is fully governed by the changes induced in the adsorbate states. b, c Introduce first and second order effects, respectively. Here the double headed pink arrows indicate that the interaction perturbs the states at both ends, leading to explicit changes in the interaction with the surroundings as indicated by the black arrows. The change in color scale (dark to light gray) indicates the decreasing interaction away from the adsorption center.

Full size image

This exercise can be continued to include effects of increasing coordination shells (III), (IV), and so on (Supplementary Note 4). In the following, we shall assume that the screening length of the d-electronic states only includes nearest neighbor interactions, hence interactions due to next-nearest neighbor atoms and beyond can be neglected.

The second order term introduced in Fig. 4b is given by the difference in the total energies of the surface perturbed by the adsorbate and the clean metal surface. The α-parameter scheme introduced by Roling et al.25 provides an ideal framework to quantify this effect. The α-parameter scheme offers a way to attain site-resolved stabilities of alloy systems via analysis of the local coordination and the composition of the nearest neighbors of a site. The model is based on constructed features, i.e., parameters \({\alpha }_{{{{\rm{CN}}}}}^{{Z}_{1}({Z}_{1},{Z}_{2})}\) that carry information about both structure and stability but cannot be extracted directly from DFT simulation. Each \({\alpha }_{{{{\rm{CN}}}}}^{{Z}_{1}({Z}_{1},{Z}_{2})}\) corresponds to the average energy gained upon forming a bond leading to the coordination number, CN, and the parameters have to be trained against a sufficiently large number of DFT computed binding energies. For a fcc metal, the maximum CN (in bulk) is 12, hence there are 12 \({\alpha }_{{{{\rm{CN}}}}}^{{Z}_{1}({Z}_{1},{Z}_{2})}\) parameters for a mono-metallic fcc system. For alloys, additional parameters are needed. The superscript Z1(Z1, Z2) indicates that the alpha parameters are obtained for a metal of type Z1 in an alloying environment of metal types Z1 and Z2. Nonetheless, for simplicity, we will use \({\alpha }_{{{{\rm{CN}}}}}^{{Z}_{1}({Z}_{1},{Z}_{2})}\) = \({\alpha }_{{{{\rm{CN}}}}}\) in the following. The BEM of a site is determined by accounting for all the bonds formed by the atom(s) of the site, i, and the bonds formed by its neighbors, j (NNB neighbors in total)27,72:

$$B{E}_{M}=\mathop{\sum }\limits_{i=1}^{{{{\rm{CN}}}}}{\alpha }_{i}+\mathop{\sum }\limits_{j=1}^{{N}_{{{{\rm{NB}}}}}}{\alpha }_{{{{{\rm{CN}}}}}_{j}}$$
(11)

The alpha parameters can be readily employed to obtain stability, reactivity, and cohesive energies of the system26,27,72,77. Here the alpha parameters are used to determine the energy difference (\(\delta B{E}_{M}^{{{{\rm{II}}}}}\) = \(B{E}_{M}^{{{{\rm{perturbed}}}}}-B{E}_{M}^{{{{\rm{clean}}}}}\)) associated with the clean (\({\alpha }_{{{{{\rm{CN}}}}}_{i}}\)) and perturbed (\({\tilde{\alpha }}_{{{{{\rm{CN}}}}}_{i}}\)) metal surfaces in order to quantify the adsorbate-induced effect on the site stabilty:

$$\begin{array}{ll}\delta B{E}_{M}^{{{{\rm{II}}}}}=\left[\mathop{\sum }\limits_{i=1}^{{{{\rm{CN}}}}}{\tilde{\alpha }}_{i}+\mathop{\sum }\limits_{j=1}^{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{b=1}^{B}{\alpha }_{[{{{{\rm{CN}}}}}_{j}-b+1]}\right]-\left[\mathop{\sum }\limits_{i=1}^{{{{\rm{CN}}}}}{\alpha }_{i}+\mathop{\sum }\limits_{j=1}^{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{b=1}^{B}{\alpha }_{[{{{{\rm{CN}}}}}_{j}-b+1]}\right]\\\qquad\quad\,=\mathop{\sum }\limits_{i=1}^{{{{\rm{CN}}}}}\left({\tilde{\alpha }}_{i}-{\alpha }_{i}\right)\end{array}$$
(12)

where the sum over b helps to account for the fact that a neighbor (j) can bind to multiple atoms of the adsorption site (up to B). In Eq. (12), the two terms correspond to the total energies of the perturbed and clean surface, respectively. This energy difference is nothing but the change in the cohesive energy (i.e., \(\delta B{E}_{M}^{{{{\rm{II}}}}}=\delta {\varepsilon }_{{{{\rm{coh}}}}}\)) of the two systems. As for other chemisorption models, the α-parameter scheme assumes that interaction between delocalized sp-electronic states inside the metal are unaffected by the perturbations induced by the adsorbate. Thereby any variation in cohesive energy is related to variation in the localized d-states. From the above, we can thus understand the effect in terms of the influence the adsorbate has on the individual bonds of the metal adsorption site to its neighboring atoms.

It is known that the cohesive energy difference can be expressed in terms of electronic properties. If one assumes a simple rectangular shape of the d-band (see Supplementary Fig. 7) with a width, Wd, where Wd is proportional to the coupling matrix element V, we can write the cohesive energy as (see Hammer et al.1 for a more detailed discussion):

$${\varepsilon }_{{{{\rm{coh}}}}}=\int\nolimits_{{\varepsilon }_{d}-\frac{{W}_{d}}{2}}^{0}{\rho }_{d}(\varepsilon )({\varepsilon }_{d}-\varepsilon )d\varepsilon \simeq {W}_{d}{f}_{d}(1-{f}_{d})$$
(13)

where

$${f}_{d}=\frac{\int\nolimits_{{\varepsilon }_{d}-\frac{{W}_{d}}{2}}^{0}{\rho }_{d}(\varepsilon )d\varepsilon }{\int\nolimits_{{\varepsilon }_{d}-\frac{{W}_{d}}{2}}^{{\varepsilon }_{d}+\frac{{W}_{d}}{2}}{\rho }_{d}(\varepsilon )d\varepsilon }$$
(14)

is the filling of the d-states of the metal atoms in closest vicinity to the adsorption site and ρd(ε) the energy density function of the d-states (density of states, DOS), see Supplementary Note 5 for a complete derivation. Hence, for a metal system, the binding energy—or the cohesive energy, εcoh – depends on the number of electrons in the ρd(ε) distribution and its variation with energy (ε). The expression in Eq. (13) is derived for mono-metallic systems, but can easily be extended to bimetallic or multi-metallic alloys by summing over the (varying) contributions from all neighbors, j, where NNB is the total number of neighbors:

$${\varepsilon }_{{{{\rm{coh}}}}}={W}_{d}\frac{1}{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{j = 1}^{{N}_{{{{\rm{NB}}}}}}{f}_{d}^{j}(1-{f}_{d}^{j})$$
(15)

Using Eq. (15) to determine the difference between the clean and perturbed surface, and combining it with Eq. (12), we get:

$$\delta {\varepsilon }_{{{{\rm{coh}}}}}=\mathop{\sum }\limits_{i=1}^{{{{\rm{CN}}}}}\left({\tilde{\alpha }}_{i}-{\alpha }_{i}\right)={\widetilde{W}}_{d}(A)\frac{1}{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{j=1}^{{N}_{{{{\rm{NB}}}}}}\left({f}_{d}^{j}(1-{f}_{d}^{j})\right)$$
(16)

where \({\widetilde{W}}_{d}(A)\) = Wd(A/S) − Wd(S) indicates that we have taken the difference between the d-band width of the adsorbate perturbed (A/S) and the clean surface site (S). Note that we have considered the specific adsorbate effect; therefore, \({\widetilde{W}}_{d}(A)\) depends on properties of A. The change in cohesive energy is thus proportional to changes in the d-band structure \({\widetilde{W}}_{d}(A)\) of the adsorption site and also depends on the degree of filling of the d-states of the neighboring atoms. Finally, we can write a modified expression for the chemisorption energy \({{\Delta }}{E}_{{{{\rm{I+II}}}}}^{A}\) by combining Eqs. (1), (9), and (16):

$$\begin{array}{rlr}{{\Delta }}{E}_{{{{\rm{I+II}}}}}^{A}\simeq &{{\Delta }}{E}_{sp}^{A}-\left(\frac{(1-{f}_{-}){V}^{2}-{\varepsilon }_{a}^{2}({f}_{-}-{f}_{d}^{{{{\rm{site}}}}})}{{\varepsilon }_{a}^{2}}\right){\varepsilon }_{d}-(1+{f}_{-})SV&\\ &+\frac{(1-{f}_{-}){V}^{2}}{{\varepsilon }_{a}}+{\widetilde{W}}_{d}(A)\frac{1}{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{j=1}^{{N}_{{{{\rm{NB}}}}}}{f}_{d}^{j}(1-{f}_{d}^{j})\end{array}$$
(17)

Equation (17) shows how the variation in chemisorption energy depends on several properties of the substrate d-band, but also on the characteristics of the adsorbate. These include: the d-filling (\({f}_{d}^{{{{\rm{site}}}}}\)) of the site, the d-band center (εd), the width of the d-band (Wd) with and without adsorbate, the filling of the antibonding state (f−), the coupling matrix element (V), the energy of the renormalized adsorbate state relative to the Fermi-level (εa), and the filling of the d-states (\({f}_{d}^{j}\)) of the atoms in closest vicinity of the active site.

Chemisorption on near-surface alloys (NSAs)

To further investigate the effects induced by compositional changes we shall study a special group of materials, so-called near-surface alloys (NSAs). NSAs introduce compositional changes occurring in the subsurface layer of the material. Therefore, the adsorption site is fixed and all changes in adsorption with respect to a reference state can be assigned entirely to the small perturbations in the electronic structure induced by the surroundings relative to the reference. To investigate periodic trends in greater detail, we have examined chemisorption energies on NSAs based on the coinage metals Au, Ag, and Cu, where fd = 1 and εd lies well below the Fermi-level, and Pt NSAs where the d-band is tied to the Fermi-level. We consider NSAs where the second layer atoms span the full range of elements in the 3d, 4d, and 5d series of transition metals as seen in Fig. 5a. To understand the adsorbate dependence, we have considered chemisorption of O, N, CH, and Li in the fcc-hollow site on the close-packed (111) surface. These adsorbates have considerably different electronegativites and span a wide range of εa.

Fig. 5: Chemisorption on NSAs as function of d-band center.
figure 5

a Top and side view of adsorption on NSAs. b Variations in the adsorption energy for O on Au NSAs as a function of d-band center (εd). Only the guest atoms in the subsurface layer have been labeled.

Full size image

By keeping the adsorption site unchanged (i.e., \({\widetilde{W}}_{d}(A)\), V2, and SV are the same) and we study relative changes in the adsorption energy, δΔE, with respect to a reference state. We expect the energy changes to be small and well represented by the approximate expression in Eq. (18), such that

$$\begin{array}{ll}\delta {{\Delta }}{E}_{{{{\rm{I+II}}}}}^{A}\simeq -\left(\frac{(1-{f}_{-}){V}^{2}-{\varepsilon }_{a}^{2}({f}_{-}-{f}_{d}^{{{{\rm{site}}}}})}{{\varepsilon }_{a}^{2}}\right)\delta {\varepsilon }_{d}\\\qquad\qquad\quad+\frac{{\widetilde{W}}_{d}(A)}{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{i=1}^{{N}_{{{{\rm{NB}}}}}}\delta ({f}_{d}^{i}(1-{f}_{d}^{i}))\end{array}$$
(18)

It is important to note that \({\widetilde{W}}_{d}(A)\) refers to changes induced by an adsorbate onto the isolated adsorption site prior to interactions with its chemical environment. Therefore, the second term of Eq. (18) reflects the stabilizing or destabilizing effect of the neighboring atoms on the adsorption site (with the adsorbate adsorbed) relative to the reference state.

Initially, we focus on the chemisorption onto coinage metal-based NSAs. In Fig. 5b, we have plotted the adsorption energies of atomic oxygen, \(\delta {{\Delta }}{E}_{{{{\rm{ads}}}}}^{{{{\rm{O}}}}}\), against the d-band center εd for the d-projection onto the atoms making up the adsorption site on all Au-based NSAs (AuM). We observe that in most cases the d-band center shifts down relative to the Fermi-level, hence f = 1 is maintained and according to the d-band model no additional bonding through the first term in Eq. (10) is expected. Similar results are observed for NSAs based on Ag and Cu (see Supplementary Fig. 8). From Fig. 5b, we also conclude that the binding energy of oxygen gets stronger as the d-band center shifts down, which not only disagrees with the d-band model predictions but it also justifies a model derivation (Eq. (9)) that explicitly takes into account the filling of the antibonding state, f− and makes it clear that the observed variation in bonding cannot be explained by variations in the d-band center alone. In the previous section, we hypothesized that when there is significant compositional changes in the vicinity of the adsorption site then electronic changes induced in the adsorption site by the adsorbate could have a non-negligible impact on the adsorption strength. In fact, the term would be proportional to a weighted sum of fd(1 − fd) for all atoms in the vicinity of the adsorption site and with a strength or proportionality factor, \({\widetilde{W}}_{d}(A)\), that depends strongly on features related to the adsorbate A.

To test the hypothesis, we have calculated the adsorption energy for a number of adsorbates with varying features (O, N, CH and Li) on the coinage metal NSAs with 5d metals in the second layer. In Fig. 6a, the change in adsorption energy is plotted for these Au NSAs relative to clean Au as a function of the filling factor (fd) of the atoms in the subsurface layer. We notice a second order behavior in fd for all adsorbates in agreement with the second term introduced in Eq. (18). According to the model expression, the strength of this interaction will be additive in the total number of neighbors NNB. To confirm this, we computed the change in adsorption strength of oxygen while reducing the number of guest atoms in the subsurface layer from 4 to 0, for the AuRe NSA as shown in Supplementary Fig. 9. With decreasing number of Re atoms in the subsurface layer the relative energy approaches zero, clearly indicating the additive nature of the interaction resulting from the compositional changes around the active site. Depending on the adsorbate, we observe significant differences in the relative adsorption energies, an effect which is a consequence of the non-zero value of the first term in Eq. (18). As we recall, the renormalized energy level, εa is pushed towards the Fermi-level as we move to the left in the periodic table, i.e., when we empty an atomic orbital. Consequently, f− decreases gradually from O to Li as seen in Fig. 6b and shown graphically in Fig. 6c through the projection onto the p-orbital for all the adsorbates on the AuRe NSA. Hence, depending on the sign of δεd, this leads to a linearly increasing or decreasing repulsion as we move from O to Li. For the Au NSAs, δεd is generally pushed down relative to the Fermi-level leading to increasing repulsion as we sweep through the transition metals and furthermore, the degree of repulsion should increase with changing adsorbates O, N, CH, and finally Li. This effect is clearly observed in Fig. 6a.

Fig. 6: Chemisorption as a function of fd.
figure 6

a The DFT calculated adsorption energy of O, N, CH, and Li as a function of fd on Au NSAs having 5d metals in the second layer. All energies are relative to pure Au. b The filling factor (f−) of antibonding states of the adsorbates as function of fd. c The pDOS is shown for O, N, CH, and Li adsorbed on the AuRe NSA. The region enclosed by the dotted box identifying the filled antibonding states serves as a guide to the eye.

Full size image

Comparing chemisorption on coinage and non-coinage metal NSAs

For NSAs where Pt is the host material, the d-electronic states are tied explicitly to the Fermi-level; hence, the position of the antibonding state is pushed upwards in energy compared to any of the coinage metal NSAs. Consequently, the first term in Eq. (18) becomes more dominant, a dominance that is expected to increase as the renormalized state moves closer to the Fermi-level and the antibonding state is depleted, \({f}_{-}^{{{{{\rm{O}}}}}_{2p}} \,>\, {f}_{-}^{{{{{\rm{N}}}}}_{2p}} \,>\, {f}_{-}^{{{{{\rm{C}}}}}_{2p}}\). Here \({f}_{-}^{{A}_{2p}} >\) indicates the filling of the antibonding part of the p-projection onto atom A. In Fig. 7a–c and d–f, we have plotted the relative adsorption energies of O, N, and CH over Au-based and Pt-based NSAs, respectively. The same data for Ag and Cu-based NSAs can be found in the Supplementary Fig. 10. To discuss trends among all transition metals in the 3d, 4d and 5d series and address issues related to spin, all NSAs are shown in Fig. 7, and the cases where results due to spin deviate from the spin-paired calculations are indicated (open symbols).

Fig. 7: Variation in chemisorption over Au and Pt-based NSAs.
figure 7

a–c Variation of the DFT calculated adsorption energies of O, N, and CH (\({{\Delta }}{E}_{{{{\rm{ads}}}}}^{{{{\rm{A}}}}}\)) on Au NSAs as a function of the fractional number of d-electrons (fd) in the second layer metal, respectively. Adsorption energies are shown relative to the values on pure Au. The filled symbols indicate spin-polarized (spin) calculations, whereas the unfilled symbols are for spin-paired (non-spin) calculations. d–f Show the variation of the adsorption energies of O, N, CH on Pt NSAs relative to pure Pt, respectively.

Full size image

Except for the magnetic metals, 3d V through Ni and Ru in the 4d transition metal series, all the adsorption energies essentially follow the parabolic behavior on the Au NSAs. This again shows that for the coinage metal NSAs, the second term in Eq. (18) dominates the bonding trend. However, the shift in the parabolic behavior for adsorbates O, N, and CH indicates the effect of the increasing linear term as the renormalized energy level is pushed upwards. Further, we note that the binding increases going from 3d to 4d and finally to 5d metals, suggesting that the structural variations should correlate with the spatial extension of the d-orbitals (without magnetization), and thus the overlap between them. As seen in Fig. 7 and Supplementary Fig. 11, for the cases where spin and spin transitions due to adsorption play a role, the effect is much larger on the coinage metal NSAs. These results indicate that spin quenching is screened by the d-states on the non-coinage systems. In this study, we will not consider the spin transition cases; they will be addressed in a future study.

When we consider NSAs where the host is a metal with non-filled d-states, such as Pt, it introduces an extra degree of complexity. Results for Pt NSAs, in Fig. 7d–f, differ considerably from those on Au NSAs. It is observed that the relative binding energy on Pt NSAs also display a parabolic behavior as you sweep through the transition metals in the different d-series and the parabolic behavior is strongest for O and diminishes when looking at N and finally CH. As we found for Au and the remaining coinage metal NSAs, the linear term becomes more dominant for the N and CH cases and this is also true for Pt (Fig. 7e, f). Furthermore, since \({f}_{d}^{{{{\rm{site}}}}}\ne 1\) then f− is much smaller for the different adsorbates, therefore, for non-coinage NSAs the first term in Eq. (18) becomes the dominant term and the linear behavior is enhanced. Our model clearly captures this behavior and hence, we can apply it to address the adsorption energies for all NSAs. Effectively, this suggests that our model can be parametrized to predict chemisorption energies from electronic structure information of the bare NSAs only.

Generalizability across different adsorbates, host and guest atoms

We set up a generalized expression based on Eq. (18) and performed a parametric fit to the calculated data to develop a predictive model. The model expression reads:

$$\delta {{\Delta }}{E}_{{{{\rm{I+II}}}}}^{A}={\beta }_{d}\cdot \delta {\varepsilon }_{d}+{\beta }_{{f}_{d}}\mathop{\sum }\limits_{i=1}^{{N}_{{{{\rm{NB}}}}}}\delta \left({f}_{d}^{i}(1-{f}_{d}^{i})\right)+\gamma$$
(19)

The constant γ contains information about the clean reference system and the fitting coefficients βd and \({\beta }_{{f}_{d}}\) depend on site and adsorbate properties. Hence, for a given site and adsorbate the coefficients are constant and the variables relative to the reference system are adsorbate independent. Accordingly, once the coefficients have been fitted, we can easily obtain the bond strength for systems outside the training set. Nonetheless, it is essential to point out that Eq. (18) is an approximate expression which is most accurate within small variations around the reference system. The coefficients in Eq. (19) are optimized for combinations of adsorbate (O, N, or CH) and adsorption site (Au, Ag, Cu, or Pt), individually, ensuring small perturbations with respect to the reference system. Hence, the coefficients for a particular set (site composition and/or adsorbate) will have to be trained through a limited set of full DFT calculations. Note that the parameters reported are obtained from linear regression of 27 DFT calculations. The MAE changes by around 0.01 eV when the entire training set is reduced to half, 14 DFT calculations. In the future, a more thorough investigation on the sensitivity of the model with dataset size will be carried out. Beyond that, the only input needed is the d-band center (εd) of the site and the formal filling (\({f}_{d}^{i}\)) of the neighboring atoms. In Fig. 8a, we have made a parity plot of actual DFT calculations against predicted data for all our NSAs. Our model is able to predict adsorption energies with a total MAE of 0.16 eV.

Fig. 8: Prediction of adsorbation energies using the model.
figure 8

a, b Predicted adsorbation energies of O, N, and CH on Au, Ag, Cu, and Pt NSAs using the models in Eqs. (19) and (20), respectively. The predicted energies are plotted against those obtained with DFT. The black dotted line indicates parity. c, d Show the variation of O adsorption energies on Pt NSAs as a function of εd and εd + Wd/2, respectively. All the values for εd and εd + Wd/2 are listed in Supplementary Tables 1–4.

Full size image

Upon inspection of Fig. 8a, it is apparent that compared to the other NSA hosts, Pt data points are significantly more scattered from the parity line and, in particular, the alloys involving early transition metals. There is an underlying physical explanation for this deviation related to the use of εd as an independent variable in the model expression. According to Eq. (18), the bond strength of an adsorbate on a NSA surface based on non-coinage metals is dominated by the term linear in εd. Furthermore, as the d-shells are emptied of the elements in the second layer the energy distribution is shifted relative to the Fermi-level and the overlap with the d-states in the Pt surface layer increases leading to an overall downshift in εd. The effect is a weakening of the bond-strength between the surface and the adsorbate A which is fully captured by the model expression. Regardless of the compositional changes in the vicinity of the adsorption site, the number of electrons in the atomic orbitals stay constant, hence εd and Wd are not independent parameters for a given adsorption site. If the d-center is pushed down, the state density must broaden out to ensure that \({f}_{d}^{{{{\rm{site}}}}}\) remains constant. Since Wd is a physical feature related to the overlap of d-electronic states between neighboring metal atoms, then, to a large degree, it controls the position of the d-band center in a non-linear fashion. In Supplementary Fig. 12, we have shown that the late transition metal has a distinct linear correlation16. However, as we approach the early transition metals, it begins to bend off. Therefore, the adsorption energies cannot maintain a linear correlation with εd as shown in Fig. 8c64. The effect can be remedied by introducing a parameter that contains information about the d-band structure and therefore ties stronger to the Fermi-level. The upper d-band edge given approximately by εd + Wd/2 is a better descriptor since it effectively identifies the position of the antibonding state (Fig. 8d)17,64. The model expression in Eq. (19) is easily modified to incorporate this effect by replacing εd with (εd + Wd/2) to obtain:

$$\delta {{\Delta }}{E}_{{{{\rm{I+II}}}}}^{A}={\beta }_{d}\cdot \delta \left({\varepsilon }_{d}+\frac{{W}_{d}}{2}\right)+{\beta }_{{f}_{d}}\frac{1}{{N}_{{{{\rm{NB}}}}}}\mathop{\sum }\limits_{i=1}^{{N}_{{{{\rm{NB}}}}}}\delta \left({f}_{d}^{i}(1-{f}_{d}^{i})\right)+\gamma$$
(20)

Refitting the βd and \({\beta }_{{f}_{d}}\) coefficients with εd + Wd/2 variable significantly improves the agreement between actual DFT calculations and model results with a total MAE for all data points of 0.13 eV as seen in Fig. 8b. The fitted values of βd and \({\beta }_{{f}_{d}}\) coefficients are listed in Supplementary Table 5.

Transferability to mixed surface and subsurface alloys

For the model to be used as a guide in the design of catalytic materials, it needs to be able to treat generic multi-metallic systems of different structure and composition. In the following, we shall address the modifications needed to extend the model.

In Fig. 9a, we have plotted \(\delta {{\Delta }}{E}_{{{{\rm{ads}}}}}^{{{{\rm{O}}}}}\) for Au-based surface alloys (see insert in Fig. 9a) against the ideal d-filling, fd for the 5d transition metal surrounding the adsorption site. We expect to see a similar trend as for the NSAs since identical arguments concerning electronic modifications induced by the adsorbate holds true for the surface alloys. The trend shown in Fig. 9a is clearly second order in fd but contrary to the NSAs, the adsorption site in a surface alloy is stabilized, thus suggesting that interaction within the surface layer is significantly different from the interaction with the subsurface layer. This is not surprising and can be explained by simple surface science arguments. Cutting a bulk crystal along any direction designated by the Miller indexes [hkl] will leave a number of unterminated dangling bonds on each surface atom. These will interact with the surrounding surface and subsurface atoms to lower the overall surface energy as illustrated in Fig. 9b and d. If we imagine that the surplus of electronic states on one atom in the (111) surface was pushed to the neighboring surface layer atoms, then we would introduce differences in the charge density among the atoms on an otherwise completely symmetric surface. Consequently, the atoms with higher density would have to shift to compensate for the extra charge, thus leading to a buckled surface as seen in Fig. 9c(i). If the extra charge on the other hand was pushed towards the atoms in the subsurface the consequence would be a contraction/expansion of the atoms in the surface layer as seen in Fig. 9c(ii). We are not arguing that the former case is not possible. In fact, many alloy surfaces are seen to reconstruct quite significantly due to exactly this effect. We have just limited ourselves to surface alloys with negligible reconstruction and, hence, cases where electronic states in the surface exclusively are pushed into the second layer. Mixed alloys comprised of surface and near-surface alloys can be addressed in a similar way (Fig. 9d). We have demonstrated that surface and subsurface guest atoms will interact differently with the adsorption site and, hence, we need to treat them independently in the model such that:

$$\begin{array}{rlr}\delta {{\Delta }}{E}_{{{{\rm{I+II}}}}}^{A}=&{\beta }_{d}\cdot \delta \left({\varepsilon }_{d}+\frac{{W}_{d}}{2}\right)+{\beta }_{{f}_{d,1}}\frac{1}{{N}_{{{{\rm{NB}}}},{{{\rm{surf}}}}}}\mathop{\sum }\limits_{i = 1}^{{N}_{{{{\rm{NB}}}},{{{\rm{surf}}}}}}\delta \left({f}_{d}^{i}(1-{f}_{d}^{i})\right)&\\ &+{\beta }_{{f}_{d,2}}\frac{1}{{N}_{{{{\rm{NB}}}},{{{\rm{sub-surf}}}}}}\mathop{\sum }\limits_{i = 1}^{{N}_{{{{\rm{NB}}}},{{{\rm{sub-surf}}}}}}\delta \left({f}_{d}^{i}(1-{f}_{d}^{i})\right)+\gamma \end{array}$$
(21)

where explicit summation over all subsurface NNB,sub-surf and surface NNB,surf layer atoms has been introduced.

Fig. 9: Chemisorption on surface alloys.
figure 9

a Shows the variation of O adsorption energies on Au surface alloys (SAs) as a function of fd. The guest atoms M` are selected from the 5d metals. The black dotted line is added as a guide to the eye. b–d Illustrates the effect of charge transfer within the surface and subsurface layer atoms. Host atoms, guest atoms in the subsurface layer, and guest atoms in the surface layer are shown in yellow, gray, and teal, respectively.

Full size image

To test our model expression, we have randomly chosen a set of mixed surface and near-surface alloys, as shown in Fig. 10a. We fit the data using Eq. (21) with separate coefficients \({\beta }_{{f}_{d,1}}\) and \({\beta }_{{f}_{d,2}}\) (see Supplementary Tables 6 and 7 for listed values). All calculated and modeled data points follow the parity line (Fig. 10b) with a MAE of 0.13 eV. Hence, this corroborates the conceptual framework derived herein and can therefore be applied to predict and interpret chemisorption onto a broad range of alloy surfaces. Our model, which is based only on four parameters (βd, \({\beta }_{{f}_{d,1}}\), \({\beta }_{{f}_{d,2}}\), γ) shows a similar or better performance than models based purely on data-driven approaches (MAE: 0.10–0.25 eV) which rely on a large number of descriptors and a significant amount of training data30,31,44,45,51,53,54,55,56,78. The simple expression introduced in Eq. (21) holds promise for application in the design of catalysts with improved physical insight and interpretability.

Fig. 10: Parity plot showing the transferability of our model across mixed surface alloys.
figure 10

a Depicts the randomly selected Au-based mixed surface and subsurface alloys, AuMM`. Here, M and M` represent the guest atoms in the subsurface and surface layer, respectively. b Chemisorption energies (δΔEads,model) predicted by the model in Eq. (21) plotted against adsorption energies obtained from DFT calculations, (δΔEads,DFT). Both NSA, SA and mixed alloys are included in the plot. The black dotted line indicates parity. We have provided all DFT energies of the bare and adsorbed systems and adsorbates in Supplementary Tables 8–15.

Full size image

In this work, we have proposed a physics-based chemisorption model that incorporates adsorbate-induced perturbations of the surface-site. Our model successfully accounts for the observed variations in chemisorption strength when tuning the active site, e.g., upon changes in the chemical environment around the site. On near-surface alloys these compositional changes result in second-order variation in the chemisorption energy of O, N, CH, and Li. Additionally, the shift in the parabolic behavior and increasing linear behavior from O to Li is governed by the properties of the adsorbates. We have explained this in terms of substrate electronic factors and adsorbate characteristics, namely the fractional d-filling (fd), the d-band center, the width of the d-states, the filling of the antibonding state, and the position of the renormalized adsorbate state. For a given adsorbate, we have parameterized the model based on four parameters only, to effectively forecast the chemisorption energies in terms of readily derived surface characteristics: i.e., 1st and 2nd moment of the d-band and the idealized d-band filling of the neighboring atoms. The model parameters for a particular adsorbate and adsorption site can be optimized using a small set of DFT training data. The method has been generalized for the entire series of transition metals and multi-component systems. Our model predicts the chemisorption energy on a variety surface alloys with a MAE of 0.13 eV, reflecting a remarkable predictive performance combined with physical transparency and interpretability. We believe this model can be used to advance the field in guiding the engineering of complex alloys with desired catalytic properties.

Methods

Computational details

All the calculations have been performed using density functional theory (DFT) and the Vienna ab initio Simulation Package (VASP.5.4.4)79,80,81 within the atomic simulation environment (ASE)82. The Revised Perdew-Burke-Ernzerhof (RPBE)83 functional form of the generalized-gradient approximation (GGA) for electron exchange and correlation energies is used throughout. All calculations are performed with spin-polarization unless stated otherwise. The projector-augmented wave (PAW) method is employed to describe the electron-core interaction with a kinetic energy cutoff of 600 eV for the surface calculations, and PBE PAW potentials used to represent the core electronic states84,85. Geometry optimizations for all models are performed with a force convergence criterion of 0.01 eV Å−1. A Methfessel-Paxton smearing of 0.20 eV is used for all the calculations.

Model structures are created from fcc host metal structures of Au, Ag, Cu, and Pt with optimized lattice constants obtained from our previous work (experimental values are shown in parentheses and all values are given in Å): Ag 4.22 (4.09), Au 4.20 (4.08), Cu 3.68 (3.61), Pt 3.99 (3.92)25. Model surfaces for near-surface alloys (NSAs) are studied on 2 × 2 surface (111) slabs constructed with six layers of fcc metal atoms. The NSAs are generated by replacing all of the host atoms in the second layer by a single transition metal element from the 3d, 4d and 5d series. Surface alloys (SAs) are formed by replacing the atoms in the top layer in direct coordination with the atoms comprising the adsorption site on a 4 × 4 supercell containing five atomic layers. For all model structures, the top two layers were allowed to relax in the z-direction (lateral relaxation is constrained to avoid excessive reconstruction due to size differences between the guest and host atoms). All adsorbate atoms (O, N, CH and Li) are allowed to fully relax. In all surface calculations, a vacuum distance of at least 20Å between repeating slab images is used and dipole corrections86 are applied to avoid spurious interaction between the slabs. For the geometry optimization, 4 × 4 × 1 and 2 × 2 × 1 k-meshes were used for the sampling of the Brillouin zone using a Gamma centered Monkhorst-Pack scheme87 for the NSAs and SAs, respectively. Following geometry optimization, a 12 × 12 × 1 and 6 × 6 × 1 k-point mesh is used for single point calculations to determine the adsorption energies and electronic-structure properties for the NSAs and SAs, respectively. Gas phase energies of adsorbate species have computed using 21 Å × 22 Å × 23 Å unit cells with 1 × 1 × 1 k-points. Crystal orbital Hamilton populations (COHP) calculations are performed using the Lobster package88,89. We performed a supercell convergence test showing that the adsorbate interactions between neighboring cells vanish when we consider larger cells as shown in Supplementary Fig. 13. We have considered the 2 × 2 supercell that is sufficient to avoid adsorbate-adsorbate effects and address the local chemical interactions resulting from the coupling between the adsorbate and the surface. All obtained structures and energies from this study are uploaded to the Catalysis Hub90, and can be accessed directly via https://www.catalysis-hub.org/publications/SainiElectronic2022.