Bayesian Learning of Chemisorption for Bridging Complexities of Electronic Descriptors

6 Building upon the d -band reactivity theory in surface chemistry and catalysis, we develop a Bayesian 7 learning approach to probing chemisorption processes at atomically tailored metal sites. With representative 8 species, e.g., *O and *OH, the Bayesian model trained with ab initio adsorption properties of transition 9 metals predicts site reactivity at a diverse range of intermetallics and near-surface alloys while naturally 10 providing uncertainty quantiﬁcation from posterior sampling. More importantly, this conceptual framework 11 sheds light on the orbitalwise nature of chemical bonding at adsorption sites with d -states characteristics 12 ranging from bulk-like semi-elliptic bands to free-atom-like discrete energy levels, bridging complexities 13 of electronic descriptors for the prediction of novel catalytic materials.

Adsorption of molecules or their fragments at transition-metal surfaces is a fundamental process for many technological applications, such as chemical sensing, molecular self-assembly, and heterogeneous catalysis.Because of the convoluted interplay between electron transfer and orbital coupling, chemical bonding can be formidably complex.Recent decades have brought major advances in spectroscopic tools [1,2] which reveal orbitalwise information of chemisorbed systems and concurrently in predicting chemical reactivity at sites of interest via electronic factors, e.g., number of valence d-electrons [3], density of d-states at the Fermi level [4], d-band center [5], and d-band upper edge [6,7].Compared with a full quantum-mechanics treatment of many-body systems, the simplicity of physics-inspired descriptors comes at a cost of limited generalization, particularly for high-throughput materials screening in which variations of site composition and configuration are sufficiently large to invalidate the perturbation approximation.Incorporation of multi-fidelity site features into reactivity models with machine learning (ML) algorithms has shown early promise for the prediction of adsorption energies with an accuracy comparable to the typical error (∼0.1−0.2 eV) of density functional theory (DFT) calculations [8][9][10][11][12].However, the approach is largely black-box in nature, prohibiting its physical interpretation.Developing a theory-based, generalizable model of chemisorption that bridges complexities of electronic descriptors and predicts the binding affinity of active sites to key reaction intermediates with uncertainty quantification represents one of the biggest challenges in fundamental catalysis.
Here we present a Bayesian inference approach to probing chemisorption processes at metal sites by learning from ab initio datasets.The model is built upon the basic framework of the d-band reactivity theory [5] while employing a Newns-Anderson-type Hamiltonian [13,14] to capture essential physics of adsorbate-substrate interactions.Such types of simplified Hamiltonians were originally used for describing magnetic properties of impurities in a bulk metallic host [13] and later extended with success by Newns and Grimley to chemisorption at surfaces [14,15].A basis set of orbitals consisting of the adsorbate and substrate states was used for solving the hybridization problem within a self-consistent Hartree-Fock scheme [14].Despite a remarkable success in advancing the basic understanding of adsorption phenomena at surfaces, particularly for d-block metals [6], its application in materials design remains limited due to the lack of accurate model parameters and meaningful error estimates.Bayesian inference produces the posterior probability distribution of model parameters under the influence of observations and prior knowledge [16].
With representative species, e.g., *O and *OH, we demonstrate the predictive performance and physical interpretability of the Bayesian models for chemical bonding at a diverse range of in-termetallics and near-surface alloys, bridging complexities of electronic descriptors in search of novel catalytic materials.

Results
The d-band reactivity theory.Within the basic framework of the d-band reactivity theory for transition-metal surfaces, the formation of the adsorbate-metal bond conceptually takes place in two consecutive steps [5].First, the adsorbate orbital (or orbitals) |a at ǫ 0 a couples to the delocalized, free-electron-like sp-states of the metal substrate, leading to a Lorenzian-shaped resonance state at ǫa .Second, the adsorbate resonance state interacts with the localized, narrowly-distributed metal d-states, shifting up in energies due to the orthogonalization penalty for satisfying the Pauli principle and then splitting into bonding and anti-bonding states.The first step interaction contributes a constant ∆E 0 albeit the largest part of chemical bonding.The variation in adsorption energies from one metal to another is determined by the metal d-states.This part of the interaction energy ∆E d can be further partitioned into orbital orthogonalization and hybridization contributions [17].As a first-order approximation, the orbital hybridization energy can be evaluated by the changes of integrated one-electron energies [18].The orbital orthogonalization cost is con-sidered simply as proportional to the product of interatomic coupling matrix and overlap matrix, V S, or equivalently αV 2 ad , where V 2 ad is the interatomic coupling matrix element squared when the interacting atoms are aligned along z-axis and α is the orbital overlap coefficient.The absolute value of V  [20,21], the Schmickler model of electron transfer has been developed to understand H 2 evolution/oxidation and OH − adsorption at metal-electrolyte interfaces.However, the deterministic fitting of adsorption properties from a single surface is prone to overfitting or trapping into a locally optimal region, limiting its application in catalysis.
Bayesian learning.We instead employ Bayesian learning to infer the vector of model parameters θ = (∆E 0 sp , ǫa , ∆ 0 sp , α, β) ′ from the evidence, i.e., ab initio adsorption properties, along with prior knowledge if available [16].In Bayes' view, those parameters are not deterministic point values, but rather a probabilistic distribution reflecting not necessarily the random nature of physical variables but rather the uncertainty.The use of parameter distributions as opposed to computationally-derived point values has obvious advantages for uncertainty quantification.In the chemical sciences, Bayesian learning has been used for calibration and validation of thermodynamic models for the uptake of CO 2 in mesoporous silica-supported amines [22], designing the Bayesian error estimation functional with van der Waals correlations [23], and identifying potentially active sites and mechanisms of catalytic reactions [24], just to name a few.The Bayesian approach allows one to infer the posterior probability distribution P ( θ|D) for latent variables based on the prior P ( θ) as well as the likelihood function P (D| θ) subject to the observation D. The (1.67 eV −1 ) is smaller than that of p z (1.77 eV −1 ), consistent with the symmetry analysis that the p xy orbitals that are parallel to a surface form π-bonds with the d-states, while the p z orbital can interact through a stronger σ bond.A weaker coupling manifests itself in a narrower orbital splitting of π/π * than that of σ/σ * , which has been previously observed using the angle-resolved photoemission spectroscopy on Cu and Ni [28].In Fig. S3 and S4, it shows that the model-constructed projected density of states onto symmetry-resolved orbitals closely resemble the DFT-calculated distributions and the predicted values of *O adsorption energies have a MAE ∼0.17 eV, suggesting the robustness and generalizability of the approach.
To test prediction capability of the Bayesian model for unseen systems, we took the *OH species as a case study because of its fundamental importance in understanding the nature of chemical bonding [29] and practical interests as a key reactivity descriptor in metal-catalyzed electrochemical O 2 reduction [30], CO 2 reduction [31], H 2 oxidation in alkaline electrolytes [32], etc. Three frontier molecular orbitals, i.e., 3σ, 1π, and 4σ*, are assumed to be involved in chemical bonding [29].Symmetry-resolved, molecular orbital density of states projected onto OH along with adsorption energies are used as the DFT ground truth Y in Eq. 6.With the Bayesian model developed here (see Fig. S5-7 for posterior parameter distributions, model-predicted adsorption energies and projected density of states on training samples), we predict *OH binding energies at a diverse range of intermetallics and near-surface alloys.Specifically, we included A 3 B, A ′ @A ML , A-B@A ML , A 3 B@A ML , A@A 3 B, and A@AB 3 , where A (A ′ ) represents 10 fcc/hcp metals used in the model development and B covers d-metals across the periodic table (see ref [33] for structural details and tabulated data).The A sites of above-mentioned surfaces exhibit diverse characteristics of the metal d-states ranging from bulk-like semi-elliptic bands to free-atom-like discrete energy levels [34], as illustrated in Fig. 4  We demonstrated that the Bayesian models of descriptor species, e.g., *O and *OH, optimized with pristine transition-metal data predicts adsorption energies at a diverse range of atomicallytailored metal sites with a MAE ∼0.1−0.2 eV while providing uncertainty quantification.Incorporation of physics-based models into data-driven ML algorithms, e.g., deep learning, might hold the promise toward developing highly accurate while interpretable reactivity models.Furthermore, this conceptual framework can be broadly applied to unravel orbital-specific factors governing adsorbate-substrate interactions, paving the path toward design strategies to go beyond adsorption-energy scaling limitations in catalysis.

Methods
DFT calculations Spin-polarized DFT calculations were performed through Quantum ESPRESSO [36] with ultrasoft pseudopotentials.The exchange-correlation was approximated within the generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof (PBE) [37].{111}-terminated metal surfaces were modelled using (2 × 2) supercells with 4 layers and a vacuum of 15 Å between two images.The bottom two layers were fixed while the top two layers and adsorbates were allowed to relax until a force criteria of .1 eV/ Å.A plane wave energy cutoff of 500 eV was used.A Monkhorst-Pack mesh of 6×6×1 was used to sample the Brillouin zone, while for molecules and radicals only the Gamma point was used.Gas phase species of O and OH were used as the reference for adsorption energies of *O and *OH, respectively.The projected atomic and molecular density of states were obtained by projecting the eigenvectors of the full system at a denser k-point sampling (12×12×1) with a energy spacing 0.01 eV onto the ones of the part, as determined by gas-phase calculations.Further details and tabulated data can be found in ref [9].
The d-band reactivity theory To revisit the d-band theory of chemisorption along with new developments, let's consider a metal substrate M in which electrons occupy a set of continuous states with one-electron wavefunctions |k and eigenenergies ǫ k , and an isolated adsorbate species A with a valence electron described by an atomic wavefunction |a at ǫ 0 a , see Fig. 1.When the adsorbate is brought close to the substrate, the two sets of states will overlap and hybridize with each other.The strength of such interactions is determined by the coupling integral V ak = a| Ĥ|k , where Ĥ is the system Hamiltonian.Within the Newns-Anderson-Grimley model of chemisorption [13][14][15], Ĥ is defined as, where σ denotes the electron spin, n is the orbital occupancy operator, and c † and c represent the creation and annihilation operator, respectively.The first two terms in Eq. 1 are the one-electron energies from the adsorbate and the substrate when they are infinitely separated in space.The last term captures the coupling, or intuitively electron hopping, between the adsorbate orbital |a and a continuum of substrate states |k .If the one-electron states of the whole system can be described as a linear combination of the unperturbed adsorbate and substrate states, the one-electron Schrödinger equation can be solved using the Green's function approach [14].In Fig. 1, we illustrate the chemisorption process of a simple adsorbate onto a d-block metal site characterized by delocalized sp-states and localized d-states [17].The interaction of the adsorbate state at ǫ 0 a with the structureless sp-states, typically accompanied with electron transfer from/to the Fermi sea, results in a broadened resonance (or so-called renormalized adsorbate state) at a perturbed energy level ǫa .Conceptually viewing chemical bonding as consecutive steps in Fig. 1, the renormalized adsorbate state then couples with the narrowly distributed d-states, shifting up in energies due to orbital orthogonalization that increases the kinetic energy of electrons and splitting into bonding and anti-bonding states.One important information from this framework is the evolving density of states projected onto the adsorbate orbital |a upon adsorption in which spin is neglected for simplicity.The effective adsorbate energy level, denoted by ǫ a , is determined by the image potential of a charged particle in front of conducting surfaces and the Coulomb repulsion between electrons in the same orbital [14].The chemisorption function ∆(ǫ) includes contributions from the sp-states and the d-states To simplify the matter, only the 2 nd step interaction, i.e., the coupling of the renormalized adsorbate state with the substrate d-states, is explicitly considered in Eq. 2. As a new development in our approach, we include an energy-independent constant ∆ 0 sp along with ∆ d as the chemisorption function ∆(ǫ).The inclusion of ∆ 0 sp provides a lifetime broadening of the adsorbate state, serving as a mathematical trick to avoid burdensome sampling of the resonance, i.e., the Lorentzian distribution ρa from the 1 st step interaction in Fig.  [27,38]: The constant 2 considers spin degeneracy of the orbital, ña is the occupancy of the renormalized adsorbate state by integrating the Lorentzian distribution ρa up to the Fermi level ǫ F (taken as 0), and f is the idealized d-band filling of the metal atom.The tan −1 is defined to lie between −π to 0 since ∆ 0 sp is a non-zero constant across the energy scale [-15, 15] eV.Thus there is no need to explicitly include localized states even if present below or above the d-band.As a good approximation, the overlap integral S ad is linearly proportional to the coupling integral for a given adsorbate, i.e., S ad ≈ α|V ad |, in which α is termed the orbital overlap coefficient.Similarly, the effective coupling integral squared V 2 ad can be written as β V 2 ad , where β denotes the orbital coupling coefficient and V 2 ad characterizes the interorbital coupling strength when the bonding atoms are aligned along the z-axis at a given distance [39].Its values of d-block metals relative to that of Cu are readily available on the Solid State Table [19].∼ N (4, 1).We assume that the predicted adsorption properties from Eqs. 2 and 4 are subject to independent normal errors.Specifically, for the property Y and the surface i we have where ǫ i is an independent and standard normal random variable and σ is the standard deviation, allowing for a mismatch between the model prediction Ŷi ( θ) and the DFT ground truth Y i .In this approach, we define the likelihood function of the property Y from n observations [41] P where the sum runs over n training samples for the property Y , which is either the projected density of states onto an adsorbate orbital or adsorption energies.For adsorption energies, Y i and Ŷi are scalar values with no ambiguity.For projected density of states, it is a vector of paired values, i.e., the one-electron energy of a state and its probability density, thus deserving a clarification.The mean squared residuals of model prediction from Eq. 2 for the surface i is used as {Y i − Ŷi ( θ)} 2 in Eq. 6.To compute the transition probability of each MCMC step, we define the sum of the (negative) logarithm of the likelihood functions corresponding to projected density of states onto each adsorbate orbital and binding energies with a hyperparameter λ adjusting the weight of two contributing metrics, i.e., − ln(P ρa ) − λ ln(P ∆E ).To optimize this parameter, we varied it on a grid of 1, 10, 100, and 1000, and found that 100 is the optimal value to

Supplementary Files
This is a list of supplementary les associated with this preprint.Click to download. supp.pdf

Figure 1 .
Figure 1.Illustration of chemical bonding at transition-metal surfaces within the d-band reactivity theory.An adsorbate A with a valence electron at a discrete energy level ǫ 0 a first interacts with the free-electronlike sp-states of the substrate M, forming a broadened resonance at ǫa accompanied with electron transfer.Conceptually, it further overlaps and hybridizes with the narrowly distributed d-states, which leads to a splitting into bonding and anti-bonding states.The work function φ and Fermi level ǫ F of M are marked.

Figure 2 .Figure 3 .
Figure 2. Bayesian parameterization.a The co-variance of the joint posterior distribution for each parameter pair and the 1D histogram of model parameters (∆E 0 sp , ǫa , ∆ 0 sp , α, and β) from MCMC simulations for *O adsorption at the fcc-hollow site of the {111}-terminated transition-metal surfaces.A top view of the model structure is shown in inset.b Schematic illustration of the MCMC sampling in a multi-dimensional parameter space.
Fig. 4(c) shows the partition of *OH adsorption energies resulting from the 2 nd step interaction (∆E d ) into orbital orthogonalization and hybridization.As we can see, for 3d, 4d, and 5d series

Figure 4 .
Figure 4. Model test and interpretation.a The d-states of a transition-metal site exhibit diverse characteristics ranging from bulk-like semi-elliptic bands to free-atom-like discrete energy levels (Pt and Ag 3 Pt as examples).b DFT-calculated vs. model-predicted adsorption energies of *OH at the atop site of {111}terminated intermetallics and near-surface alloys.c Partition of *OH adsorption energies at the M site of Ag 3 M into orbital hybridization and orthogonalization of 3σ, 1π, and 4σ* orbitals with the metal d-states.
Bayesian learning Due to the computationally intensive nature of the MCMC algorithm, there is a need for a more efficient implementation of the Newns-Anderson-Grimley model than what is obtained by Python and standard libraries like SciPy and NumPy.We make extensive use of Cython, a C++ extension to the standard Python, to speed up the performance (10−1000 times) of some CPU-intensive functions in the model, e.g., Hilbert transform.To perform MCMC sampling, we use PyMC, a flexible and extensible Python package which includes a wide selection of built-in statistical distributions and sampling algorithms[40], e.g., Metropolis-Hastings.A "burn-in" of the first half of the samplings and then thinning (1 out of 5 samplings) was performed to ensure that subsequent ones are representative of the posterior distribution.Convergence of our MCMC-based sampling was verified using parallel chains[25].The MCMC sampling results can be directly visualized using corner, a open-source Python module.We took Normal for floatingpoint variables unrestricted in sign, LogNormal for non-negative parameters, and Uniform for others.∆E 0 sp and ǫa can be estimated from DFT calculations of the adsorbate on a simple metal, e.g., sodium (Na) at the face-centered cubic (fcc) phase.Specifically, for *O, we used ∆E 0 sp ∼ N (−5.0, 1), ǫa ∼ N (−5, 1), ∆ 0 sp ∼ LN (1, 0.25), β ∼ LN (2, 1), and α ∼ U (0, 1).For *OH, we used ∆E 0 sp ∼ N (−3.0, 1), ǫ3σ a ∼ N (−6, 1), ǫ1π a ∼ N (−2, 1), and ǫ4σ * a

Figures Figure 1
Figures

Figure 3 Model
Figure 3

Figure 4 Model
Figure 4 [19]can be written as β V 2 ad , in which the standard values of V 2 ad relative to Cu are readily available on the Solid State Table[19].The overall adsorption energy ∆E can then be written as the sum of the energy contributions from the sp-states ∆E 0 and the d-states ∆E d , with the latter 1. Accordingly, ǫ a will be replaced by the renormalized adsorbate state at ǫa .Attributed to the narrowness of a typical metal d-band, ∆ d can be simplified as the projected density of d-states onto the metal site ρ d (ǫ) modulated by an effective coupling integral squared V 2 ad , i.e., ∆ d ≃ πV 2 ad ρ d (ǫ).Λ(ǫ) is the Hilbert transform of ∆(ǫ).In this framework, the interaction energy between the adsorbate and the substrate can be partitioned into two contributions, i.e., ∆E sp and ∆E d .∆E sp is the energy change due to the interaction of the unperturbed adsorbate orbital(s) with the delocalized sp-states, while ∆E d is the energy contribution from further interactions with the localized d-states of the substrate.Since all d-block metals have a similar, free-electron-like sp-band, ∆E sp can be approximated as a surface-independent constant ∆E 0 sp albeit the largest contribution to bonding [17].To calculate ∆E d , we include both the attractive orbital hybridization