Modeling the impact of structure and coverage on the reactivity of realistic heterogeneous catalysts

Main

In the typical image of a transition-metal heterogeneous catalyst, molecules adsorb, react and desorb independently and on a static catalyst surface. Yet, this is far from reality. Catalysts undergo a wide variety of structural transformations during reactions, such as surface segregation, dispersion, agglomeration and even phase transitions. Many of these transformations can be induced or accelerated by adsorbates, which can interact strongly with and partially cover the catalyst surface under realistic reaction conditions. Although simulations¹ and surface science experiments² have been invaluable for gaining insights into how the nature of a catalyst affects the binding properties of adsorbates, the reverse—how adsorbates affect the structure and especially the reactivity of their host catalyst—is proving much more challenging to understand^3,4.

These challenges largely stem from the wide variety of ways in which adsorbates can interact with and alter the catalyst. At high coverages, adsorbates can strongly modify the interactions between catalyst atoms. This may lead to the formation of dynamic, reaction-driven active species and sites^5,6,7 as well as changes in the morphologies of nanoparticles, clusters and even single-atom catalysts⁸. Predicting the structure of catalysts in the presence of adsorbates is therefore a formidable task, not to mention how these changes can affect catalytic reactivity—even on static catalyst surfaces, the stability and reactivity of catalytic intermediates can differ substantially at low and high coverages as they interact via direct through-space and surface-mediated interactions^9,10. Indeed, high coverage is crucial for enabling certain processes, such as Fischer–Tropsch synthesis¹¹. In reality, the situation is often more complex, as the catalyst structure and reactivity can be tightly intertwined. In particular, Nobel laureate Gerhard Ertl’s pioneering work was pivotal in uncovering the critical importance of accounting for a dynamic surface structure to understand catalytic reactivity. His studies on CO oxidation over Pt surfaces revealed striking spatiotemporal oscillatory behavior, which was traced back to a reversible switching between reconstructed and non-reconstructed surfaces—each with their own unique reactivity—as CO* coverage rises and falls¹².

To unravel the nature of such catalyst–adsorbate interactions and their impact on catalyst structure and reactivity, researchers have increasingly employed in situ studies, which probe catalytic behavior under reactive conditions, as well as operando studies, which combine in situ characterization with the simultaneous measurement of catalytic performance^13,14. However, they face various challenges due to the limitations of spectroscopic instruments at the high pressures and temperatures typical of experimental conditions. The primary obstacle lies in the inability of spectroscopic probes such as photons and electrons to penetrate the dense reaction environment and access the catalyst surface.

Compared with experiments, computations do not face such limitations and are theoretically well-suited to probe such atomic-scale phenomena¹⁵. However, the cost of modeling these complex phenomena is currently prohibitive. Typically, idealized surfaces corresponding to 0 K/ultrahigh vacuum (UHV) conditions are used as models, giving rise to a large pressure, temperature and materials gap between simulations and catalytic experiments. There is therefore a tremendous opportunity to develop more realistic and representative models with improved predictive power, enabling the construction of more accurate structure–property relationships and driving a more rational, data-driven approach to designing catalysts¹⁶. Critically, these higher-fidelity models would also enable us to systematically understand and harness the complex interactions between catalysts and adsorbates, which are a rich and unexplored dimension for tuning and enhancing catalytic properties.

In this Review we share our vision of how cutting-edge computations, combined with in situ and operando experiments, can help unravel the dynamic interplay between adsorbates and catalysts, as well as how we can leverage these interactions in rational catalyst design. Recent strides in computational power, artificial intelligence and spectroscopic methods are beginning to shed light on these complex molecular-level processes, making this a timely moment for exploration. We first examine the chemistry of adsorbate–surface interactions and how they can affect catalytic structure and reactivity. We then provide our perspectives on how recent advances in simulations, machine learning (ML) and artificial intelligence can help overcome challenges in catalyst modeling to make sense of these phenomena and design improved catalysts.

Although our Review focuses on and draws most examples from thermal catalysis over transition-metal nanoparticles, the concepts discussed are broadly relevant to electro- and photocatalysts^4,11,17, and can be generalized to various types of catalytic system, such as zeolite-supported metal catalysts¹⁸, metal–organic frameworks (MOFs)¹⁹ and single- and dual-atom catalysts¹⁷. We finally note that it is important to distinguish between a true ‘catalyst’—the active material responsible for chemical reactivity—and a ‘precatalyst’ that transforms in operando to generate the active catalyst. As it is challenging to definitively determine whether a material functions as a catalyst or a precatalyst, we will use ‘catalyst’ throughout this Review as an inclusive term for both cases, except when discussing materials widely accepted as precatalysts in the literature.

Adsorbate-induced structural changes

Understanding the gamut of adsorbate-induced structural changes, from the facet to the nanoparticle scale (Fig. 1a), is the first step toward devising new modeling strategies for tackling them. Fundamentally, adsorbates tend to weaken inter-catalyst interactions, increasing the mobility of catalyst atoms by lowering the barriers for their diffusion—this is also known as the skyhook effect (Fig. 2). For instance, the formation of Pt hydride (Pt–H) complexes enhances the mobility of Pt adatoms on Pt(110) by more than 500 times at room temperature, as has been directly visualized via scanning tunneling microscopy (STM)²⁰. Similarly, CO adsorption enhances the mobility of Pd clusters on Fe₃O₄ via the formation of Pt carbonyl (Pt–CO) species²¹, and analogous CO-induced mobility has also been observed for Au clusters supported on CeO₂ (ref. ²²).

**Fig. 1: Structural transformations catalysts can undergo.**

**Fig. 2: Explaining the skyhook effect.**

This increase in the atomic mobility of the catalyst is a key effect of adsorption and is the foundation for further structural transformations. For example, the resulting adatoms may combine to form metal clusters once a critical nucleation density is crossed^21,23,24, as seen in the formation of nanoclusters upon exposure of Pd²⁵ and Cu^23,26 nanoparticles to CO (Fig. 1b). Such nanoclusters can be highly active for various chemistries, such as Pt nanoclusters for the (de)hydrogenation of hydrocarbons²⁷. In another example, Ag clusters formed when Ag nanodisks were exposed to laser irradiation in a hydrogen atmosphere²⁸. In such photocatalytic reactions, the mobility of the catalyst atoms is also intrinsically higher as the illuminated catalysts are heated to temperatures that can nearly melt them, leading to substantial reconstruction, even under an inert N₂ atmosphere²⁸.

The skyhook effect also mediates strong metal–support interactions (SMSIs)²⁹, a phenomenon whereby the support migrates onto the surface of the catalyst and which can be used to tune and enhance the reactivity of the catalyst¹¹^,30 (Fig. 1c). One key mechanism for this enhancement is the increased area of the active support–metal interface³¹. Due to the increase in the atomic mobility of the support, adsorbates can also reduce the temperatures needed for SMSIs^32,33. This adsorbate-enhanced mobility has been leveraged to reduce the temperature required to activate, by SMSI, Cu/ZnO/Al₂O₃ catalysts for methanol steam reforming³⁴. Interestingly, however, not all adsorbates increase the mobility of adatoms. OH* for example, instead serves to anchor Pd–CO* species on Fe₃O₄(001) (ref. ²¹).

On alloy catalysts, adsorbates may have a stronger affinity for one metal over the other. This drives adsorbate-induced surface segregation, where typically the more reactive metal species are drawn from the subsurface to the surface^35,36. However, in the case of CuPt near surface alloys with a Pt skin, introduction of CO surprisingly induced the segregation of Cu to the surface, forming a CuPt surface alloy³⁷. Although Cu by itself binds CO more weakly than Pt, the CuPt surface binds CO more strongly than pure Pt, thus illustrating the complex and sometimes counterintuitive nature of adsorbate interactions with alloy surfaces.

Another fundamental effect of adsorption to consider is how it alters surface energies. For example, adsorption-induced strain can lead to local transformations in the catalyst structure to alleviate this strain³⁸. Changes in the surface energy of a facet can also induce substantial effects on a larger scale, such as precipitating facet reconstruction³⁹ and changes in the morphology of entire nanoparticles to minimize the system’s total energy⁴⁰. For example, just 0.15 monolayers of thallium (Tl) was able to entirely lift the conjugate honeycomb-chained-trimer reconstruction of a Au/Si(111)-√3 × √3 surface⁴¹ (Fig. 1d). This was observed for Cs, In and Na as well, and was proposed to be due to charge donation from the adatoms to the surface, which increases the Coulomb repulsion between surface atoms and thus the surface stress⁴². In contrast to the conventional view of a static catalyst, many catalysts may therefore exist in a state of fluctuating shape and morphology as adsorbate coverages rise and fall^12,43.

In liquid-phase reactions, such as electrocatalytic reactions, additional avenues for catalyst reconstruction open up due to environmental factors such as the applied electric potential⁴⁴ and solvents⁴⁵. For example, positive potentials can oxidize the catalyst during the oxygen evolution reaction (OER), resulting in the in situ formation of highly oxidized metal species, such as metal (oxy)hydroxides, which are widely recognized as the active species⁴⁶. To take advantage of this, MnN nanocuboids have been designed to be oxidized in situ to form a strained, hydroxylated Mn₃O₄ shell that is highly active for the oxygen reduction reaction (ORR)⁴⁷. During CO₂ electroreduction (CO2RR), negative potentials were found to trigger the formation of Cu nanograins via the reduction of surface Cu oxides^17,48, as well as the formation of Cu clusters⁴⁴ via increased H* and CO* coverages⁴⁴. Although commonly thought of as one-way processes, such potential-induced phenomena may be reversible and exist in a dynamic equilibrium. This is exemplified by the case of mixed NiFe hydroxides for the OER, which deactivates by forming a FeOOH secondary phase during OER conditions but can be reactivated by lowering the potential to reduce this secondary phase⁴⁹.

Beyond these physical changes, chemical phase transformations can also occur if adsorbates are able to alloy with the catalyst. For instance, Fe oxide, commonly used as a precatalyst for Fischer–Tropsch synthesis, carburizes to Fe carbide, the in situ generated active phase, during the reaction itself⁵⁰. Cu nanoparticles can also oxidize into Cu₂O during various reactions, such as propylene oxidation⁵¹ and the aqueous phase reforming of methanol to H₂ (ref. ⁵²), which leads to a decrease in their activity. Interestingly, however, these changes can be also reversed by illumination with light, which excites charge carriers to reduce the Cu oxides back into metallic Cu (refs. ^51,52).

Modeling and prediction of catalyst structure

The wide range of effects adsorbates have on catalysts poses formidable challenges for modeling and simulation, which we are beginning to be able to tackle, in part due to the rapid advancements in ML techniques and computational power. In this section, we first briefly explore first-principles methods for modeling catalysts, and then discuss approaches to model and predict how catalysts reconstruct in the presence of adsorbates, which we broadly categorize into two main approaches: equilibrium-based and dynamic methods (Fig. 3).

**Fig. 3: Equilibrium-based and dynamic approaches for the prediction of catalyst morphology in the presence of adsorbates.**

Overview of first-principles methods

Choosing a first-principles method is a complex choice that depends on the catalytic system under study and the properties of interest. Density functional theory calculations with the generalized gradient approximation (DFT-GGA) serves as the main workhorse of computational heterogeneous catalysis due to its good balance of efficiency and accuracy, especially for metallic systems. However, more sophisticated approaches may be needed for specific materials. For example, strongly correlated materials with localized d or f electrons, such as NiO and CeO_x, exhibit substantial self-interaction errors that require advanced treatment with Hubbard U corrections (DFT + U), meta-GGAs or hybrid functionals⁵³. For modeling charge transfer and excitation processes in photo- and electrocatalysis, such higher-level methods are also required to obtain accurate band structures and excited states^54,55. We refer interested readers to recent reviews^56,57 for more comprehensive discussions of these first-principles methods.

Equilibrium-based methods

On selecting a first-principles method, predicting what structural reconstructions adsorbates can induce can then be framed as a high-dimensional global optimization problem of finding energetically low-lying catalyst states over the vast number of atomic degrees of freedom contributed by the catalyst, adsorbates and even the support¹¹. Equilibrium-based methods are one way to do this. Typically, such methods involve compiling a database of candidate structures and searching across this curated set of structures to locate the most stable one.

To simplify the problem, much research has focused on finding stable adlayer structures on a fixed catalyst facet. For example, the myriad of different phases formed by co-adsorption of O* and CO* on Pd(100) can be determined by means of comprehensive configurational searches to build up a library of candidate adlayer and surface oxide phases⁵⁸. Phase diagrams can then be constructed to ascertain the most stable phases, which change as a function of temperature and CO and O₂ partial pressures as these changes affect their surface energies (Fig. 4). In electrocatalysis, analogous Pourbaix diagrams can help reveal how catalyst phases change as a function of electrode potential and pH⁵⁹. The computational hydrogen electrode provides a convenient framework to estimate electrochemical potentials, such as those of protons and halides for example, by assuming equilibrium between gas-phase species and their corresponding solvated ions⁶⁰. These theoretical predictions can be validated via experimental techniques such as STM and low-energy electron diffraction (LEED) spectra.

**Fig. 4: Phase diagram of Pd(100) in the presence of CO(g) and O₂(g).**

To model reconstructions, the constraint of a fixed catalyst facet can be removed, extending the set of candidate structures to include different catalyst structures such as reconstructed or surface-segregated phases⁶¹. Taking Fig. 4 for example, surface oxides were included in the construction of the phase diagram to capture the propensity of Pd(100) to form thin films of surface oxides. Finally, by considering combinations of both different adlayer and catalyst structures, we can arrive at a comprehensive set of structures that allows us to directly predict adsorbate-induced facet reconstruction, as showcased in studies of the O*-induced segregation of PdRh surfaces⁶².

Accurate free-energy calculations are essential for bridging temperature and pressure gaps and enabling robust phase predictions. Free energies are typically calculated using entropies derived from the harmonic oscillator approximation, but this approach may breakdown under several conditions, such as high temperatures, dense adsorbate coverages, and weakly bound adsorbates. More sophisticated methods have been employed to tackle such cases. For instance, an exact solution of the quantum mechanical translational states of monoatomic adsorbates can be derived by mapping out the potential energy surface (PES) of the system⁶³. Similarly, the hindered translator and rotor model uses information from the PES to better describe soft vibrational modes associated with adsorbate movement parallel to surfaces⁶⁴.

Building on these phase predictions from individual facets, we can then model complete nanoparticles using the Wulff construction⁶⁵, which optimizes the proportions of individual catalyst facets to minimize the total surface energy of the nanoparticle. This approach allows us to predict the structure of adsorbate-decorated nanoparticles^66,67, which can be validated and compared with spectroscopic measurements, either directly via transmission electron microscopy (TEM) imaging to capture the shapes of the nanoparticle, or indirectly via infrared spectroscopy. In particular, infrared spectroscopy uses the spectroscopic signatures of adsorbed intermediates such as CO* and H* to probe and quantify the different types of surface site, from which the structure of the nanoparticle may be deduced³. These probe molecules, however, should be an existing reactant or intermediate in the considered reaction, or they might inadvertently introduce unrelated structural changes in the catalysts. An additional limitation is that coverage effects may alter their electronic structure and potentially complicate the interpretation of the spectra³.

One possible modeling workflow to predict adsorbate-dependent nanoparticle morphology from the bottom up would be to first gather candidate structures of a standard set of low-index facets at varying surface coverages. At any given reaction condition, we can then determine the most stable phases for each facet. Although the surface energies of clean facets are typically used in the Wulff construction, we can also minimize the energy of the system using the surface energies of these different adlayer-covered facets, which will then lead us to obtain the most stable predicted nanoparticle morphologies, considering the influence of adsorbates^68,69. To simulate finite temperatures, Boltzmann statistics can also be employed to predict an ensemble of possible nanoparticle structures⁷⁰.

Although these equilibrium-based methods offer a systematic approach to predict adsorbate-influenced morphologies, they have a major defect that limits their predictive strength—their reliance on the extensivity of the predetermined set of structures used to construct these phase diagrams and Wulff constructions. The combinatorial complexity of possible placements of adsorbates on surfaces renders an exhaustive search of the configurational space impractical in most cases. Developing more efficient global minima searching tools to traverse the large unified configurational space of supports, nanoparticles and adsorbates is therefore imperative for accelerating the search for stable catalyst structures.

Examples of efforts in this area include enhanced sampling methods, such as metadynamics simulations⁷¹ and stochastic surface walking (SSW)⁷², which introduce bias potentials to accelerate sampling of rare events, as well as basin hopping⁷³, particle swarm optimization (PSO)⁷⁴ and grand canonical Monte Carlo (GCMC) simulations⁷⁵. Among these approaches, the replica-exchange grand canonical (REGC) scheme stands out for its ability to perform unbiased sampling of the configurational and compositional space by maintaining a series of molecular dynamics (MD) simulation replicas spanning different temperatures and chemical potentials⁷⁶. These replicas periodically attempt to swap their configurations to help overcome kinetic barriers for exploring the phase space. Additionally, this MD-based approach naturally captures configurational and anharmonic effects, providing reliable free energies for the different phases.

However, the cost of applying these techniques with first-principles accuracy remains prohibitive for all but the simplest systems. Approximations are commonly used, such as relying on human intuition to sample only well-known phases, or assuming that the most stable adsorption sites at the limit of clean coverage are the first to be occupied at higher coverages. In most cases, this provides a good balance between configurational space sampling and computational cost. ML-based approaches have also been developed, such as AutoSurfRecon, which uses machine learning interatomic potentials (MLIPs) for rapid energy evaluation together with the virtual surface site relaxation–Monte Carlo (VSSR-MC) algorithm to accelerate configurational sampling⁷⁷. Generative diffusion models also show substantial promise at finding new phases, as evidenced by a recent study where a novel structure for a large Ag₂₉O₂₂ boundary domain was discovered with good agreement with STM images⁷⁸.

Dynamic methods

As models of catalytic systems become more realistic, they naturally exhibit increased complexity. This is exemplified in electrocatalysis, where accurate modeling must simultaneously account for multiple phenomena: surface adsorption and reconstruction, diverse configurations of solvent molecules with similar energies, charge rearrangement at interfaces, and the presence of interfacial structures such as the electrical double layer⁷⁹. Similarly, supported nanoparticles present their own set of challenges, including symmetry mismatches between the support, and complex and dynamic interfaces between the nanoparticle and support.

For these complex systems, dynamic simulations that simultaneously sample all components of the system—including adsorbates, catalyst structure and the reaction environment—and treat them on the same level, may be a more efficient way to explore configurational space. For instance, on-lattice kinetic Monte Carlo (KMC) simulations that model surface atoms as individual lattice sites have revealed how adsorbates can induce surface segregation on alloy surfaces, showing the formation of surface ensembles (dimers/trimers) upon CO adsorption⁸⁰. In a similar fashion, KMC simulations have also been used to predict CO-induced cluster formation on Cu(111), providing a mechanistic explanation for a phenomenon commonly observed in experiments²³.

Compared with on-lattice methods such as KMC, methods that are more continuous, such as MD simulations, will be more efficient for sampling systems with high fluxionality or low symmetry, including complex adlayer configurations⁸¹. One key example is subnanometric nanoparticles⁸², which can access ensembles of low-energy metastable configurations at reaction conditions⁸³. Having a larger surface area and therefore a higher proportion of adsorbates to catalyst atoms also means adsorbates can potentially have a large impact on the structure of these nanoclusters. This is especially so due to their high curvature, which can support higher coverages than slab models⁸⁴. Ab initio MD (AIMD) simulations are particularly valuable in this context, successfully tracking the structural evolution of Pd/Au clusters with CO* (ref. ⁸⁵), as well as determining stable geometries of supported Pt cluster at predefined H* coverages⁸⁶. Beyond nanoclusters, MD simulations are also critical for modeling electrocatalytic environments, where they enable rigorous sampling of the diverse landscape of energetically similar water configurations⁸⁷.

Another key advantage of dynamic simulations as compared with equilibrium-based methods is that they can provide valuable kinetic insights into how reconstructions form and evolve during reactions, enabling us to understand how catalysts evolve as reactions proceed^85,88. For instance, reactive MD simulations driven by Bayesian MLIPs have been synergistically integrated with extended X-ray absorption fine structure (EXAFS) measurements to reveal a shared quasi-icosahedral during structural transformation of Pt nanoclusters of different morphologies when exposed to hydrogen⁸⁹. Such insights could be augmented with static DFT calculations to shed light onto the specific atomic-scale interactions driving the reconstruction⁹⁰, which will be crucial for guiding the design of custom metastable surface structures for catalysis⁹¹.

Although such time-based dynamic methods can usually cover a broader configurational space than equilibrium-based methods, much remains to be improved in terms of their sampling ability. MD simulations are prone to being trapped in local states, requiring long simulations to overcome the barriers between different basins. Furthermore, the space explored by KMC methods is typically limited by the user-defined lattice. We describe two paths to mitigate these challenges, which can also be applied to alleviate the cost of global optimization searches that involve dynamic simulations in their workflows.

The first involves the development of surrogate models^92,93,94, which can be trained to reproduce first-principles data and therefore provide first-principles-like accuracy at a fraction of the cost. Physics-based cluster expansion models^95,96, which parameterize the energetics of the system in terms of pair, triplet and higher-body atomic interactions, have been used to great effect in lattice-based simulations such as KMC and GCMC simulations. More recently, MLIPs show promise by serving as drop-in substitutes for first-principles calculations, providing quick evaluations of energetics and forces to extend the time and length scales of atomistic simulations. For instance, MLIP-accelerated basin-hopping Monte Carlo simulations have successfully simulated CO-induced nanoisland formation from Pt steps⁹⁷, and MLIP-accelerated MD simulations have been able to locate low-energy metastable structures of Pt–H clusters⁹⁸. The Au herringbone reconstruction, which previously eluded explanation by DFT simulations using smaller unit cells, was also successfully explained by MLIPs with an Au slab model consisting ~40,000 atoms⁹⁹. Such reconstructions are typically realized in the order of nanoseconds (~10 ns for Au reconstruction⁹⁹ and 50 ns for Pt nanocluster reconstruction⁸⁹), well within the reach of current MLIPs¹⁰⁰.

The second method involves developing more efficient global optimization methods. The advent of powerful generative AI tools—such as Boltzmann generators¹⁰¹, generative flow networks (GFlowNets)¹⁰², generative adversarial networks (GANs)¹⁰³ and diffusion-based generative processes¹⁰⁴—holds great promise due to their ability to pinpoint novel structures via the inverse design of desired properties within the discrete space of catalytic materials. This is a substantial advantage compared to methods such as MD and MC simulations, which are trajectory-based and may require long equilibration due to the presence of kinetic barriers. As new generative AI tools proliferate, it will be important to develop metrics to determine which tools can be more useful for problems in catalysis. For example, the ability to generate diverse, metastable structures that exist at finite reaction temperatures, as opposed to searching only near the global minima, as is the tendency for methods such as GANs, will be particularly relevant for fluxional materials such as nanoclusters.

Effects of structural transformations on catalytic reactivity

Catalyst reconstructions and the presence of spectator adsorbates can substantially influence the local chemical environment, impacting adsorption and activation energies, and thus overall catalytic reactivity. Typically, these transformations are detrimental to catalytic activity. For instance, coking, sintering or agglomeration of catalyst particles generally leads to catalyst deactivation. In another example, pure metallic Ir and Ru are known to be highly active catalysts for the acidic OER¹⁰⁵. However, they oxidize into more stable but less active IrO₂ and RuO₂ phases during the OER process¹⁰⁵.

Much research has focused on how to prevent detrimental deactivation, such as through encapsulation techniques¹⁰⁶, but it is increasingly being recognized that some types of reconstruction may actually be beneficial in nature. For instance, Cu nanograin boundaries are highly active for CO2RR^107,108. Similarly, strained Pt overlayers¹⁰⁹ and undercoordinated Pt sites left behind after the leaching of lanthanide Pt alloys¹¹⁰ and PtNi nanowires¹¹¹, respectively, also possess enhanced ORR activity. In the alkaline OER, transition-metal oxides and perovskites function as precatalysts that undergo surface or even overall structural transformation in situ as they are oxidized, forming active sites with high metal oxidation states such as metal oxyhydroxides^112,113. Leveraging this phenomenon, there have been substantial efforts in rationally designing reconstruction-prone catalysts with potentially enhanced OER activity^114,115. An instructive example is the doping of Fe into the inactive precatalyst CoAl₂O₄ spinel¹¹⁶, which facilitates its restructuring into active Co oxyhydroxides during the alkaline OER by promoting the pre-oxidation of Co and increasing the structural flexibility of the catalyst¹¹⁶.

Catalytic activity may also be enhanced by the broad range of binding environments in reconstructed catalysts, which may assume cooperative bi- and multifunctional roles. This allows different site types to catalyze successive steps in cascade reactions: for example, intermediate products may desorb from one type of site and subsequently adsorb onto another for further reaction. This approach, also known as tandem catalysis¹¹⁷, leverages different sites to perform different tasks in the overall reaction. Interestingly, this concept was also proposed to be responsible for the enhanced activity of high-entropy alloys¹¹⁸, for instance.

Another increasingly important concept in understanding how structural changes affect reactivity is the idea of self-adaptive catalysts, which possess the ability to generate their own active sites via coverage-induced transformations^117,119. By intentionally designing for reconstructions, new catalytic sites can be generated dynamically, offering novel reaction paths on the same catalytic material. Take, for instance, the design of dual-metal-site pairs (Cu–Ni DMSPs) in MOFs¹²⁰. The flexibility of the ethylenediaminetetraacetic acid (EDTA) skeleton allows the distance between the Cu–Ni DMSPs to continuously adapt to different reaction intermediates in the reduction of CO₂ to CH₄ as the reaction progresses, allowing the reaction to progress with higher selectivity¹²⁰.

Modeling reactivity of low-symmetry catalysts

The diverse local environments on low-symmetry reconstructed catalysts, compounded by the presence of spectator adsorbates, imposes significant complexity on the modeling of catalytic reaction mechanisms. These diverse environments can come from both the low intrinsic symmetry of the catalyst surface or from high adsorbate coverages. Therefore, instead of the traditional most stable site approach commonly used for idealized slab surfaces, it becomes crucial to consider the distribution of available binding sites and how each of them contributes to overall reactivity.

The first hurdle to investigating the reactivity of low-symmetry catalysts lies in identifying the numerous unique sites of interest and generating the required simulation inputs, which can be tedious and impractical to perform by hand. Automated frameworks for sampling such structures will therefore be required. For example, surface scanning methods¹²¹ can be used to map out the PESs for adsorption of various species with minimal human intervention. Alternatively, heuristics-based determination of adsorption sites—such as those based on Voronoi tessellation algorithms in CatKit¹²², for example—may be more efficient. However, such heuristics are currently designed for slab systems and have yet to be adapted to handle nanoparticles and rugged surfaces reliably.

ML-based techniques are also increasingly employed to determine possible adsorption sites. Notably, AdsorbDiff uses diffusion models to predict the optimal adsorption site and adsorbate orientation, followed by refinement via MLIP-based optimization¹²³. AdsorbML uses a stochastic sampling approach to sample up to 100 random site placements for subsequent MLIP optimization¹²⁴. With an added DFT single-point calculation, such an approach can achieve up to >85% accuracy in determining binding energies within 0.1 eV of the DFT global minimum. This can be especially efficient for low-symmetry slabs where heuristic approaches may generate many possible binding sites.

An additional challenge is the high cost of first-principles calculations, which may pose constraints when characterizing the energetics of these numerous binding sites. As such, surrogate models are increasingly important, either to aid directly with property prediction or to sieve out important candidates for DFT refinement¹²⁵. Such models include descriptor-based approaches¹²⁶ such as scaling relations¹²⁷, generalized¹²⁸ and orbital-wise¹²⁹ coordination number, d-band centers¹³⁰ and electron counting¹³¹. However, despite considerable work in this area, whether these descriptors can generalize from simple ideal models to complex reconstructed surfaces in the presence of spectator adsorbates remains an open question.

As an alternative to these hand-crafted descriptors, graph neural networks (GNNs), a natural and general way to model discrete systems where atoms and bonds are represented as vertices and edges of a graph, respectively, have emerged as a promising method for predicting a wide variety of properties, including energetics. Due to their high inherent model complexity, they can generalize better, especially when trained on large datasets (>10,000 data points). For instance, the adsorbate chemical environment-based graph convolution neural network (ACE-GCN) accurately ranked the energetics of various high-coverage configurations of adsorbates on mono- and bimetallic transition-metal surfaces⁹³. MLIPs based on GNNs or other architectures could also help to perform pre-optimization or preliminary screening and ranking of high-coverage configurations¹²⁴.

Besides adsorption calculations, transition state (TS) searches, already a costly and complex task for model surfaces, will be especially challenging when considering the additional degrees of freedom contributed by the catalyst atoms, spectator adsorbates and potential solvent molecules in the reaction coordinate. The use of low-coverage scaling and Brønsted–Evans–Polanyi relations¹²⁷ obtained from idealized surfaces have been proposed as possible solutions, but they may not generalize well to high-coverage scenarios due to changes in the possible reaction pathways. For example, co-adsorbates capable of hydrogen bonding can shift the preferred reaction pathway for formate decomposition by ‘hooking’ the formate adsorbate from above¹³². Such a system exhibits characteristics of a three-dimensional rather than a two-dimensional system, which can potentially break scaling relations¹²⁷. These co-catalytic intermediates might also be metastable, making them difficult to elucidate via experimental techniques and traditional DFT approaches.

The large number of possible binding sites also poses difficulties in sieving out those that should be prioritized for TS searches, especially as metastable configurations tend to be more reactive than their more stable counterparts. One way to overcome this challenge is to use enhanced sampling methods, such as metadynamics or umbrella sampling, which can smartly sample the PES to locate low-energy pathways in complex energetic landscapes¹³³ and at solvated interfaces¹³⁴. Yet, this requires extensive prior knowledge of the reaction, such as the possible products and the collective variables leading to them. This prevents them from being employed in a high-throughput manner.

The above challenges underscore the need to develop tools for fully autonomous TS calculations for heterogeneous catalysis, which will be especially useful for reactions with large reaction networks, such as for CO2RR and biomass reforming. Such tools currently lag behind those available for organic chemistry¹³⁵ and homogeneous catalysis¹³⁶, perhaps due to the unique complexities imposed by the surface. Early tools developed include the adaptive KMC (aKMC) method, which uses the dimer min-mode following method to populate a rate table with possible reaction processes¹³⁷. These processes are ultimately applied to coarse grain state transitions in MD simulations. Another tool is the stochastic surface walking reaction sampling (SSW-RS) method⁷², which automatically identifies reactant–product pairs for a given reactant and connects them with the double-ended surface walking TS search¹³⁸. Combined with a NN-based MLIP for rapidly evaluating reaction energetics, the SSW-RS method has been used to map out the reaction network for glucose pyrolysis, totaling 6,407 elementary reactions¹³⁹ (Fig. 5). The scale of the reaction network for this commonplace molecule highlights the crucial need for such methods to tackle complex reactions. To enable the full potential of these methods, enhancements such as integrating them with state-of-the-art MLIPs to alleviate their high computational demands, and automated processing of their raw information to determine unique reaction pathways¹⁴⁰, would be welcome.

**Fig. 5: Mapping of the reaction network of glucose pyrolysis via NN potential-accelerated SSW-based global optimization.**

Even with tools for automated TS calculations, an exhaustive probe of the reactivity of all possible metastable structures and the active sites within remains impractical. One way to tackle this is via importance sampling techniques¹⁴¹, which can be used to reduce the number of calculations needed by using inexpensive surrogate models to quickly identify promising active sites. These sites will then be validated with more costly, but accurate, ab initio calculations. This leverages the fact that only a small number of active sites (<20%) dominate the total rate.

Alternately, a paradigm shift in how reaction studies are conducted could be advantageous. In traditional approaches, one typically first determines the most stable and near most stable adsorption intermediates, using these as initial starting points for TS searches. However, the global activity search (GAS) technique, designed to probe the reactivity of nanocluster catalysts, inverts this process¹⁴². It directly samples different cluster configurations at the TS to locate the one that minimizes the energy of the TS. Once this most stable TS is located, it is then connected to its corresponding minima, completing the reaction pathway. When applied to model alkane dissociation on Pt nanocatalysts, GAS revealed that concomitant reconstruction of the cluster was vital for accessing the most stable TS for C–H bond dissociation. By bypassing the need to probe multiple intermediates on different cluster configurations, many of which may be unreactive, this paradigm shift substantially simplifies the modeling of nanocatalysts.

Once the relevant minima and TSs are identified, kinetic modeling can be used to obtain macroscopic kinetic data, such as reaction rates. However, the mean-field microkinetic models (MKMs) typically used for idealized surfaces will not be suitable for more realistic, inhomogeneous surfaces with high adsorbate coverages due to a breakdown in the mean-field approximation. If the catalyst structure can be approximated as fixed at steady state, the use of site-averaging techniques¹⁴¹ may be one way to estimate the reaction rate, assuming the active sites contribute independently. On-lattice KMC simulations, which are higher fidelity and can model interactions between the active sites, are another option. However, they require meticulous specification that will likely involve the development of automated tools and scripting. For cases where catalyst degrees of freedom are significant, off-lattice KMC simulations, such as aKMC¹³⁷, may be employed. These kinetic models can be coupled with the advanced TS searching methods described above to build up a catalog of reactive events on-the-fly, allowing them to autonomously explore the reaction landscape to build up the reaction network.

Future perspectives in modeling

To bridge the gap between idealized slab models and realistic catalysts at reaction conditions, innovative and efficient simulation approaches are imperative. We envision that a data-driven approach, embracing both theory and experiments, will be an effective method for dealing with the complexity of these systems. Here we suggest four key themes where concentrated efforts are essential to advance the modeling of these types of system and ultimately enable the overarching goal of rational catalyst design.

Integration of AI and ML

By reducing the reliance on human intuition and teasing out previously unseen patterns in catalytic datasets, the adoption of AI and ML has already begun to yield promising results in the development of new catalysts^{94,143,144,145}, and even in discovering entirely new catalytic reactions¹⁴⁶. Given the high complexity of realistic catalysts, we believe it is inevitable that AI and ML will play key roles in their modeling. However, we should avoid ‘black box’ models that are non-interpretable and whose outputs cannot be meaningfully traced back to underlying chemical or structural features¹⁴⁷. Even if such models may provide accurate predictions, they would have limited physical insights, preventing us from advancing our understanding of the underlying phenomena.

To avoid this, we advocate for two complementary modeling approaches. First, AI and ML models can serve as surrogates for computationally intensive first-principles calculations, which are used as inputs for subsequent multiscale modeling. Such ML-accelerated physics simulations retain the ability to trace the catalyst’s reactivity back to fundamental energetic and mechanistic insights, striking a balance between computational efficiency and scientific interpretability. Second, simulations can help provide catalysis-relevant electronic and structural descriptors that can augment experimental datasets, which can then be used to train interpretable ML models that directly predict reactivity^148,149. Such models can provide chemical insights via feature importance analyses such as SHAP (SHapley Additive exPlanations). However, care must be taken to validate such models against experiments and in-depth mechanistic studies, as well as to avoid drawing spurious correlations that lack rigorous mechanistic grounding.

Predictive ML models, especially MLIPs, show substantial potential to advance the modeling of dynamic catalysts. They have already demonstrated success in elucidating the ground-state structures of supported Cu clusters (Fig. 6a), modeling of complex systems such as multiply-promoted catalysts that present numerous possible adlayer arrangements (Fig. 6b), as well as understanding surface segregation phenomena of large alloy nanoparticles with over 1,000 atoms (Fig. 6c). Their growing contributions are attributable to rapid advances in the field, from early Behler–Parrinello neural networks (BPNNs)¹⁵⁰ to kernel approaches such as the Gaussian approximation potential (GAP)¹⁵¹, and now to sophisticated message-passing GNNs (MP-GNNs), such as MACE¹⁵² and NEquIP¹⁵³. These latest MP-GNNs inherently possess a larger receptive field due to their message-passing architecture, while their graph-based representation naturally accommodates rotational, translational and permutational symmetries.

**Fig. 6: MLIPs for accelerating the exploration of vast configurational spaces.**

Yet there is still much room for improvement, especially in the treatment of long-range magnetic and electrostatic interactions, which are crucial, especially in electrocatalytic systems¹⁵⁴. Notable efforts in this direction include charge equilibration frameworks, which solve for charge distributions and explicitly account for electrostatic interactions via classical Coulomb terms, as implemented in the 4G-NNP architecture¹⁵⁵. Alternate models, such as SpookyNet¹⁰⁰, directly predict charges, and incorporate charge and spin embeddings to account for global charge and spin. However, the computational cost of solving for the charge distribution and building the Coulomb matrix can be highly significant, especially for large systems¹⁵⁶. Improving on this, So3krates employs a spherical self-attention mechanism to capture long-range nonlocal effects, avoiding the need for explicit Coulomb terms¹⁵⁷. The success of attention mechanisms in treating long-range interactions in other domains, such as AlphaFold’s modeling of protein folding¹⁵⁸, suggests that such architectures may be a crucial piece of the puzzle for future work in capturing long-range interactions in catalytic systems.

Another exciting frontier for MLIPs lies in their training on data from higher-level theoretical methods. Potential improvements extend even to metallic systems: for instance, the adsorption site of CO on Pt(111) is not correctly described by DFT-GGA but can be reproduced with the random phase approximation (RPA)¹⁵⁹. However, the computational expense of such calculations makes developing such potentials prohibitive. To address this challenge, multi-fidelity modeling approaches such as Δ-ML¹⁶⁰, transfer learning¹⁶¹ and explicit fidelity encoding¹⁶² can leverage abundant low-level theoretical data to minimize the required high-level calculations. Alternately, training could focus only on the most crucial rate-determining steps and intermediates of a given reaction¹⁶³. This is particularly important given the steep computational scaling of higher-level methods like RPA, which formally scales as O(N⁶) with the number of electrons N, as compared to DFT-GGA, which scales as O(N³).

Although increasing accuracy is crucial, another emerging bottleneck for MLIPs is their computational speed. Some structural changes in catalytic systems, such as deactivation processes, may occur over seconds or even days¹⁶⁴, revealing the challenges of relying on MLIPs and molecular simulations alone for navigating the vast configurational landscapes of catalyst structures. Even with the fastest force fields and powerful computational architectures, MD simulations struggle to cross the millisecond timescale barrier¹⁶⁵. This temporal gap necessitates new methods for sampling and generating relevant catalytic structures.

In this context, the recent emergence of powerful generative AI methods presents exciting opportunities for accelerating catalysis simulations by enabling direct inverse design, which can be much more efficient at sampling. For instance, the diffusion-based model MatterGen¹⁰⁴ can generate stable, unique and novel inorganic materials with specified chemical and mechanical properties. We note that generative AI models should respect fundamental physical laws and constraints to avoid producing unphysical structures. Early models such as MatGen¹⁶⁶ and CDVAE¹⁶⁷ struggled to generate structures with high space-group numbers due to a lack of symmetry constraints in their architectures. To tackle this, recent approaches explicitly embed physical and geometric constraints ranging from simple periodic boundary conditions (PBC) to more sophisticated implementations of force and energy equivariance and invariance with respect to permutation, and rotation, and translation. For example, MatterGen uses GNNs that are equivariant in the special Euclidean group in three dimensions (SE(3)-equivariant) to guide the generation of lattices, atom positions and elemental compositions in its denoising process¹⁶⁶. Also, the physics-guided crystal generative model (PGCGM) incorporates physics-guided explicit loss functions that enforce realistic bond distances and structural symmetry¹⁶⁸. Although these models were primarily developed for materials discovery, their adaptation to catalysis could radically transform the way we create catalyst models for subsequent simulations.

Large language models (LLMs) also have potential to be powerful tools in catalysis modeling due to their flexibility and intuitive human–computer interface¹⁶⁹. CatBERTa—a bidirectional encoder representations from transformers (BERT) model for catalysis—demonstrates an innovative approach for adsorption energy prediction by training on structures transformed into textual representations¹⁷⁰. Beyond direct predictions, LLMs can also serve as research co-pilots—given their strong multimodal capabilities, they may help seamlessly integrate computations and experiments, and, in the near future, could serve at the core of goal-oriented, agentic AI platforms to orchestrate complex workflows. For instance, we envision that such platforms could be tasked to autonomously elucidate surface structures by iteratively analyzing STM images, generating candidate surface structures, and validating with first-principles DFT software.

Training robust and transferable ML and generative AI models that can accurately span the wide range of possible catalytic structures, is, however, a major challenge. Not only does it require the curation of large amounts of expensive first-principles simulations data, but these data must also be diversified enough for the models to generalize well. For instance, MLIPs trained on local minima will struggle with simulating reactive events. To meet this challenge, various strategies for selecting diverse initial training data have been developed, including dimensionality reduction with stratified sampling¹⁷¹, entropy optimization¹⁷² and diversity-maximizing metrics¹⁷³.

Such strategies cannot, however, prevent generalization errors when MLIPs encounter extrapolated configurations outside their training domain, which can be especially dangerous if this goes undetected. Active learning provides a possible solution via its iterative approach cycling between training and uncertainty-based data selection and labeling¹⁷⁴. Not only does this safeguard against extrapolation errors if carried out on-the-fly, but it is also highly data-efficient as it judiciously selects only the most informative data points for labeling (Fig. 6d). Recent innovations in this direction include the use of adversarial attacks with a specially crafted loss function to sample uncertain structures that are also physically relevant¹⁷⁵.

The success of such an active learning framework hinges crucially on efficient uncertainty quantification (UQ), which guides the selection of the next structure for labeling¹⁷⁶. An effective UQ metric should correlate well with prediction errors, allowing new configurations to be ranked and selected by the magnitudes of their uncertainties. Current approaches include variance estimates from model ensembles, as well as mean-variance estimation from single models, which output predictions together with an uncertainty estimate. However, these methods are not universally reliable for all use cases^176,177. Developing more predictive UQ measures remains essential to efficiently curate datasets and enhance the adoption of active learning in the community.

Catalysis informatics and databases

Apart from model and workflow development, building up robust informatics infrastructures will immensely support the AI and ML efforts in catalysis. Sharing of data related to commonly employed models, such as (111) slabs, in open databases and repositories can help streamline research efforts by reducing the number of redundant computations. These databases could therefore serve as a vital foundation for jumpstarting investigations into more complex phenomena. Although a great deal of useful data have been tabulated in journal publications, enabling for example, the prediction of alloy segregation¹⁷⁸ and adatom formation²³, efforts to digitalize the wealth of information from these previous studies would improve their accessibility and applications in future research. Adoption of common standards, such as the Open Databases Integration for Materials Design (OPTIMADE) application programming interface (API)¹⁷⁹, will also be useful as we begin to develop such databases.

Initiatives such as the Open Catalyst Project¹⁸⁰, nomad-lab.eu¹⁸¹ and the Catalysis-Hub platform¹⁸² have made great strides in making catalysis data more accessible. However, they focus on elementary quantities such as low-coverage binding energies, and more comprehensive efforts to characterize catalyst materials remain lacking. One fruitful area for possible exploration would be databases focusing on high coverage surfaces, which would complement the existing low-coverage data. Such databases could offer pre-computed adsorbate configurations at various coverages, which would serve as basic building blocks to generate temperature- and pressure-dependent phase diagrams, and ultimately coverage-dependent Wulff constructions.

To generate such databases, we propose a two-pronged approach. First, pretrained ML models can be used to perform zero-shot inference on the adsorption energies of high-coverage configurations and rank them by their stability, which can be used as a metric of importance to prioritize subsequent DFT evaluation. Second, we identify a set of maximally diverse configurations that span the relevant configurational space using featurization techniques such as cluster expansion models. This carefully selected initial dataset can then be augmented by active learning, enabling efficient exploration of the vast configuration space while ensuring comprehensive coverage of physically relevant structures. Similar approaches can be taken for the building up of other complex databases.

Synergistic experimental and computational characterization

Together with the curation of simulation data, concomitant curation of high-fidelity experimental databases will be vital to calibrate, validate and refine the increasingly complex computational models used. For example, single-crystal adsorption microcalorimetry (SCAC) data have proven valuable in benchmarking different DFT functionals¹⁸³. Such databases should comprehensively document synthesis methods, characterization protocols and data, as well as reaction kinetics performance, to enable robust and meaningful comparisons with simulations⁴⁰. Initiatives such as digcat.org and CatTestHub¹⁸⁴ are important steps, but much work remains to be done to develop consistent standards for generating and managing diverse, high-quality data.

A crucial part of this effort will be to integrate and corroborate in situ experimental information from a wide range of complementary surface-sensitive techniques, such as high-pressure STM¹⁸⁵, liquid-phase TEM¹⁸⁶, X-ray absorption spectroscopy¹³, infrared³ and Raman spectroscopy¹⁸⁷. To ensure data consistency, measurements should be performed simultaneously on the same catalyst sample, as opposed to the use of different facilities to obtain the data at different times¹⁸⁸. Moreover, it will be essential to advance the capabilities of in situ methods by pushing the boundaries of their operating pressures and temperatures to encompass realistic reaction conditions, as well as to develop higher spatial and temporal resolution, some efforts in this space being the development of time-resolved infrared reflection absorption spectroscopy and quick-EXAFS¹⁶⁴.

Structure–activity relationships can be constructed via operando studies, where the reaction kinetics are tracked by simultaneously quantifying the product distribution via, for example, a gas analyzer. Such relationships will be instrumental for rationalizing catalytic activity and developing a fundamental understanding of how dynamic changes in catalyst structure affect their activity¹⁸⁹. However, the transient nature of many of these catalytic states poses challenges for in situ/operando studies, underscoring the need to develop more accurate techniques with finer time resolution, as well as methods leveraging atomistic simulations to help interpret experimental data. Fundamental studies with shape- and size-selected nanoparticles can also provide useful activity trends for more direct comparison with simulations.

Advances in experimental data curation will open new opportunities in areas such as computational structure elucidation. One promising approach involves inverting experimental spectra via ML-based approaches. For example, reinforcement learning has been used in combination with image-recognition techniques to determine the atomistic structure of reconstructed rutile SnO₂(110)-(4 × 1) from STM images¹⁹⁰. Diffusion models have also been used to determine the structure of amorphous carbon by matching with X-ray absorption near-edge structure (XANES) spectra¹⁹¹. ML algorithms that combine data from different spectroscopic sources, as well as from DFT, would be a particularly promising direction for future research.

Adaptive multiscale modeling frameworks

To accurately model catalytic performance—including activity, selectivity and stability—end-to-end simulation frameworks that are able to self-consistently adapt to possible changes in adlayer and catalyst structure will be essential. In traditional catalysis modeling, the choice of model is fixed throughout: first-principles thermodynamic and kinetic data are collated with an assumed base model consisting of a certain surface and adlayer coverage, and then fed into multiscale models to obtain macroscopic properties such as reaction rates and selectivities. In coverage and structure-cognizant modeling, however, the base model is constantly challenged and refined by incorporating the interrelated effects of composition and coverage on local reactivity to obtain a more accurate and dynamic picture of a catalyst’s stability, structure and reactivity¹⁹² (Fig. 7a). This makes sure that our predicted catalysts are not only theoretically active, but also practically viable.

**Fig. 7: Structure- and coverage-consistent framework capable of accounting for adsorbate-induced reconstruction when modeling catalytic reactivity.**

Specifically, coverage-cognizant modeling involves a self-consistent loop where, if a divergence appears between the coverage predicted by the kinetic model and the assumed base model, the catalyst should be modeled anew with the freshly predicted coverages and then re-evaluated by kinetic modeling. This continues iteratively until the modeled coverage matches the kinetic model outputs; that is, we have reached coverage self-consistency. Taking this concept a step further, structure self-consistency can also be achieved with the help of a parallel model that informs us of the most stable catalyst structure at a given coverage, which may be different from the initial structure.

Such a combined structure and coverage self-consistency approach has been used successfully in the prediction of adsorbate-induced nanoparticle reconstruction via microkinetic modeling⁶². Incorporating such approaches into volcano plot-based discovery campaigns will help improve their predictive power and hit rate for novel catalysts by fully mapping out the predicted catalyst structure and coverage in the descriptor space (Fig. 7b), which can be used to more accurately guide subsequent experimental synthesis.

Finally, the integration of nanoscale atomistic models with mesoscale transport models, such as first-principles calculations with computational fluid dynamics (CFD) simulations, can provide critical insights into the local reaction environments throughout the reactor and a better understanding of intrinsic reactivity with heat and mass-transfer limitations^193,194. Operando experiments on single crystals, for instance, commonly face severe heat and mass-transport limitations at high conversion, leading to substantially altered local concentrations, such as the formation of boundary layers of product molecules above the catalyst surface¹⁹⁵. The close coupling between reaction kinetics, heat dissipation and mass transport has been exemplified by a combined first-principles KMC and CFD simulation for CO oxidation over RuO₂(110)¹⁹⁶. Under adiabatic conditions—that is, with no heat dissipation into solid parts of the reactor—the catalyst’s temperature rose to 35 K above the inlet temperature. Yet, due to intensified mass-transfer limitations, adiabatic operation resulted in lower turnover frequencies (TOFs) than isothermal operation at the same temperature as the inlet¹⁹⁶.

Such integrated multiscale studies remain rare, as important technical challenges remain in seamlessly integrating picosecond nanoscale elementary reactions with millisecond mesoscale diffusion processes^194,197. This is especially true for electrocatalysis, where reactions can occur simultaneously in both the liquid phase and on the catalyst surface, thus dynamically coupling the two media together. Additionally, traditional models commonly approximate the catalyst as a homogeneous surface rather than using more realistic probability distributions that would better represent real catalysts¹⁹⁷. Overcoming these challenges will be critical for revealing the true nature of the active site while on stream.

Conclusions

Modeling can play a key role in understanding the nature of the active site under industrially relevant reaction conditions. The materials, temperature and pressure gaps between models and real-life catalysts are currently large, and there is an enormous configurational space of the combined catalytic system where the adsorbates can affect catalyst structure, and catalyst structures in turn can affect adsorbate stability. Simulating such behavior will require the treatment of adsorbates, the catalyst and its reaction environment on an even footing. Overcoming the much-increased complexity of such models necessitates the development of new data-driven approaches for intelligently navigating the catalytic space and elucidating possible active sites that might be formed in situ during the reaction, as well as adaptive modeling frameworks that iteratively refine simulation models to ensure coverage and structural self-consistency. Crucially, these approaches will benefit from the integration of ML, generative AI and catalysis informatics infrastructure to accelerate data curation and foster collaborative sharing. To validate and finetune these increasingly intricate simulations, it will be imperative to corroborate them with multiple in situ experimental characterization and operando techniques, which also have much room for development in terms of their temporal and spatial resolution, as well as their range of operating conditions.

Ultimately, obtaining a deeper understanding of the structural transformations catalysts can undergo, and their effects on catalytic performance, will be vital for enabling us to design better catalysts. We envision that the intertwined and complex dependence between catalyst structure and its adlayer can eventually be harnessed as a powerful lever for rational catalyst design.