Kinetic Monte Carlo simulations for heterogeneous catalysis: Fundamentals, current status, and challenges

Heterogeneous catalysts are capable of accelerating chemical reactions, thereby achieving high conversion and product selectivity at low cost.¹ For this reason, this type of catalysis is at the heart of many processes of critical importance in the chemical industry.^2,3 It not only has played a pivotal role in shaping our current society but also will play an essential role in ensuring its sustainable future. This role involves improving the efficiency of catalytic materials and processes leading to significant energy and raw material savings, as well as pollution reduction. While the development of novel improved catalysts has typically been an empirical trial-and-error process, computational modeling approaches are playing an increasingly important role in the quest for the “rational catalyst design.”^4–7 

However, due to the wide range of length scales and time scales involved in heterogeneous catalysis (Fig. 1), developing pertinent computational models is not an easy task.⁸ The elementary processes happening on the catalytic surface are determined by the structure and composition of sites onto which the reactants are adsorbed and where the breaking and making of chemical bonds occur (i.e., the active sites).⁹ At the same time, the structure and composition of active sites evolve in response to local variations in temperature and concentrations, which are themselves dictated by the overall mass and heat transport at the reactor scale. In other words, there is a dynamic interplay between the reactive elementary processes and the large-scale environment in which they evolve.

First-principles multiscale modeling approaches have been implemented to tackle this challenge.^8,10 One of the central pillars of these approaches is to exploit the disparities in the time scales at which processes at different levels unfold, from the electronic to the reactor level. In doing so, it is possible to hierarchically couple the different levels of theoretical descriptions. In its bottom-up version, a first-principles multiscale modeling approach starts from quantum mechanical electronic structure calculations toward building the atomic-scale model for the catalyst and exploring the energetics of the elementary steps in the reaction mechanism. Then, the elementary steps are used as inputs into either the mean-field microkinetics modeling (MF-MKM) approach^11–21 or the so-called Kinetic Monte Carlo (KMC) framework.^22–32 The MF-MKM approach is computationally efficient but of questionable accuracy.³³ On the other hand, KMC simulations are very accurate but computationally demanding.³⁴ Finally, the multiscale approach can integrate either of these two microkinetic modeling approaches with mass and heat transport models to account for inhomogeneities in the reactant concentrations and temperature in the reactor.

It would not be an overstatement to claim that KMC is an essential computational tool to bridge the gap between the atomistic and the catalyst scale toward high-fidelity (predictive) models of catalytic kinetics. Unlike the more widely adopted MF-MKM, KMC allows for incorporating spatial inhomogeneities and correlations in the distribution of reactants on the catalytic surface, as well as detailed information about their configurations on different types of active sites. In addition, it takes into account the discrete nature of the coverage for nanoparticles or clusters, which expose a small number of sites. Thus, it is capable of explaining phenomena that are beyond the limitations of the MF-MKM approach.

This Perspective aims to give a general perspective on the 1p-KMC modeling approach for heterogeneous catalysis. We will focus on general concepts, main computational challenges, and recent progress in implementing this versatile computational approach. The rest of this Perspective is structured as follows: First, we briefly present motivations and concepts behind the 1p-KMC framework for heterogeneous catalysis. Second, we introduce the general structure of its computational implementation and give an overview of prominent software packages and codes that build upon it. Third, we discuss outstanding computational challenges faced by the 1p-KMC framework. We continue by reviewing applications to catalytic systems of practical relevance. Finally, we discuss future challenges and directions for further research and development and end with our conclusions.

The accurate computational modeling of heterogeneous catalysis is certainly a challenging and vital endeavor.^{4,6,8,10,35–42} The acceleration of chemical reactions by the catalyst and its ability to direct the paths to yield particular products arise as a consequence of a highly complex interplay between a large number of elementary events (i.e., adsorption, desorption, diffusion, and reaction events) occurring at the active sites of the catalytic surface. A significant challenge in modeling such events is to find efficient ways to treat the time scale disparity arising from the fact that such elementary events are commonly thermally activated and are thus rare on the typical femtosecond scale of smallest atomistic vibrations.⁴³ The effective treatment of such a problem is essential to capture the wide range of temporal and spatial scales involved in heterogeneous catalysis, as outlined in Fig. 1.⁴⁴ In this respect, the KMC framework has been demonstrated to be a powerful and versatile computational modeling approach that, by taking advantage of the aforementioned time scale disparity, allows us to develop models that couple the electronic, atomistic, and catalytic nanoparticle scales in a relatively straightforward way.^{22–25,27–29} This coupling opens the window for the accurate prediction of catalytic performance metrics (i.e., selectivity, activity, and stability) critical for the successful design of new catalysts and the improvement of chemical processes at the reactor scale.

The elementary events happening on the catalytic surface are transformations in which the atoms of a system, including those of the active sites and their surroundings, rearrange from one configuration to another by breaking and making chemical bonds. The characterization of such systems requires, at the most detailed level, a fully quantum mechanical description. However, obtaining solutions of Schrödinger’s time-dependent many-body equation remains a challenging task.

Over the years, this practical problem has motivated the development of several approximations, of which the Born–Oppenheimer (BO) one is the most widely used in the catalysis research community.^45–50 The BO approximation, along with the solution of the time-independent quantum mechanical problem for electrons, given different configurations of the nuclei,^51,52 leads to the so-called potential energy surface (PES). The PES provides the energy of the system of electrons and nuclei as a function of nuclei positions and thus contains valuable information about reaction paths, adsorption energies, vibrational frequencies, and the existence of barriers for the elementary events.⁵³ A minimum in the PES represents a stable species, and the basin of attraction around such a minimum contains all atomistic configurations that will perform a vibrational relaxation into that minimum. Moreover, a transition from a basin to a neighboring one by crossing the intermediate energetic barrier constitutes an elementary event of adsorption, desorption, diffusion, or reaction.^54,55 As in many other branches of computational chemistry, the PES is a fundamental concept in KMC simulations (Fig. 2).

Although this static quantum mechanical description provides valuable information about individual chemical transformations in terms of the PES, to understand how a catalytic reaction mechanism proceeds and to calculate the reaction rate of the entire process, dynamical simulations are required. However, quantum-dynamical studies are still computationally very costly. To alleviate this issue, it is often a good approximation to neglect the quantum mechanical effects in the motion of the nuclei. This approximation entails treating the system classically by solving Newton’s equations of motion with suitable boundary conditions for all the nuclei with the forces calculated from the PES, which, in turn, is obtained from electronic calculation methods.^51,52 Such an approach is the well-known first-principles molecular dynamics (1p-MD),^47–49 and despite its strong predictive power, it becomes computationally inefficient when simulating catalytic systems for practically relevant time scales and system sizes.

Indeed, a demanding step in this method is the computation of the PES for every time step of a MD run. Several efficient methods to obtain good and efficient representations of the PES during MD simulations have been introduced,^56–62 each having its advantages and drawbacks. However, the main obstacle in this approach is that, for catalytic systems, the energetic barriers of the PES are typically much higher than the system’s thermal energy. Therefore, if only thermal energy drives the dynamics, the system performs many unsuccessful attempts before the corresponding barrier is crossed due to a large enough energy fluctuation (see Fig. 2 for an example of a typical MD trajectory). In many cases, the time between two successful elementary events can reach several nanoseconds or even longer time scales. Under these circumstances, the MD algorithm, with its time steps of the order of femtoseconds, spends most of the time resolving irrelevant thermal vibrations inside the basin of the PES before a catalytically relevant barrier crossing occurs. Although the barrier crossings per se happen quickly and could be, in principle, followed by 1p-MD simulations,^63,64 the infrequent occurrence thereof means that the full functionality of a catalytic system cannot be currently studied by this computational method. Developing frameworks to extend the accessible time scale of MD simulations has been a long-standing challenge for the atomistic simulation research community.^65–67 

The KMC framework is specifically tailored to address the time scale problem just discussed, and in doing so, it allows us to access the long time dynamics necessary to explore heterogeneously catalyzed reactions.^{24,28,54,55,68} As its name indicates, KMC combines the Monte Carlo (MC) method with a kinetic approach that focuses on the time scale of the barrier crossings and uses transition state theory (TST)⁶⁹ arguments to describe the statistics of the transitions between basins of the PES. Therefore, in contrast to traditional equilibrium MC simulations where the time variable is absent,^70,71 in KMC, the temporal evolution of the catalytic system is taken into account but just in a coarse-grained sense, namely, the simulations pertain to the time scale of the barrier crossings rather than that of atomic vibrations, as in MD.

KMC simulations in combination with density functional theory (DFT)-based calculations have grown in popularity in the catalysis research community over the past few decades.^{10,22,24,28,72–74} This 1p-KMC approach can cover temporal scales ranging from ms to hours and spatial scales from nm to μm. Let us provide the fundamental basics of such a 1p-KMC simulation approach in the rest of this section.

In KMC simulations, the long-term time evolution of the catalytic system is said to be governed by successive state-to-state transitions, with a state corresponding to a single basin of the PES. A key observation is that because the system stays inside each basin (or KMC state) for a very long time (relative to the time scale of the fastest atomistic vibration), it “forgets” how it got there. Then, it is possible to assume that such state-to-state dynamics is Markovian.^{54,55,75–77} This assumption allows assigning to each elementary event leading from a state ω to a different state ω′ a rate constant k_ωω′ that represents the probability, per unit time, that the transition to state ω′ occurs.⁷⁷ This hopping probability is independent of the previous history of state ω and is calculated via TST, taking into account the properties of the PES.⁷⁸ 

Therefore, during a KMC simulation, the system passes through a large number of states. The collection of all possible such states define the state space, Ω, to which ω and ω′ belong, and thus, a KMC trajectory is simply a random walk on Ω. In Fig. 3, we present an example of a KMC state space together with two possible ways to mathematically represent it. If all elementary events are known for every state the KMC trajectory goes through, the rate constant expressions are exact, and the Markovian approximation is valid, a state-to-state trajectory generated by KMC simulation will be statistically identical to the trajectory generated by the MD method projected onto the states—provided, of course, both methods use the same PES.⁷⁹ 

The Markovian approximation underpinning KMC is typically valid because once the catalytic system reaches a particular PES basin, it equilibrates very rapidly within that basin, which represents a KMC state. This equilibration breaks any connection between the state visited before the current state and any state the system will visit next. Moreover, as the system vibrates within a basin, it also repeatedly loses memory about where it was vibrating before. As a result, the probability of waiting for a certain amount of time until a transition occurs to another new state is the same at every moment spent in the current state (assuming, of course, no prior information about how much time has passed since the last transition). This behavior gives rise to a Poisson process with exponential decay statistics.^{54,55,75–77}

It is now evident that both the random selection of an escape time and the identification of the associated elementary event are the foundation of any KMC algorithm. This information is almost all that we need to discuss the computational implementation of the KMC framework. However, before doing that, let us discuss in Sec. II B important aspects required for implementing KMC in combination with first-principles methods.

The state space Ω encompasses all the possible adsorbate configurations on the lattice and the associated numbers of gas-phase molecules (produced or consumed during elementary events), which are specified up to an additive constant.^24,30 The MME is a set of coupled differential equations whose solution gives the probability of being in state ω at time t. However, obtaining analytical or even numerical solutions of this set of equations is a challenging task due to the high-dimensional state space of typical catalytic systems. Instead, it is more convenient to perform a sampling of trajectories that adhere to the MME via KMC simulation. One can then obtain the correct probability P_ω(t) by averaging over many KMC trajectories or, if the system is ergodic, by time-averaging over a long enough KMC trajectory.

Because the reaction energy, ΔE_rxn, is obtained from Eq. (14), the task is to find the electronic energies of the initial and final states. To achieve this, we need an appropriate energetic model of the system. Accounting for the contributions of the gas species is straightforward because the gas reservoir is typically assumed to be an ideal gas.²⁴ However, modeling the lattice energetics (e.g., electronic energy of the lattice configuration) is a more delicate issue.^72,85,86 In this respect, the so-called cluster expansion Hamiltonian (CEH) approach has gained much popularity in recent years.^30,87–93 CEHs are typically fitted to a limited dataset derived from DFT calculations, which, however, can become costly and tedious. However, in recent years, efforts have been dedicated to reducing the computational cost of getting accurate CEHs.^94,95 Hence, Eqs. (12)–(14) enable the implementation of thermodynamically consistent models in 1p-KMC, and we now turn our attention to the calculation of the rate constants.

Having described the critical aspects behind the 1p-KMC approach, we now proceed to briefly discuss how the KMC framework works and how one can implement it.

Let us now briefly outline the computational implementation of the KMC framework, as well as software packages and codes particularly devoted to simulate catalytic surface reactions.

KMC is a computational framework that generates stochastic trajectories fulfilling the fundamental premises of the MME (see Sec. II B).⁷⁵ The general flow chart of a KMC algorithm is presented in Fig. 6. To begin with a KMC simulation, it is necessary to specify in advance the lattice structure of the catalytic surface, the chemical processes that can happen on the lattice (i.e., reaction mechanism), the energetics of the system (i.e., lateral interactions), and the simulation parameters (i.e., pressure, temperature, and composition of the gas phase). Once the simulation setup and parameters are established, one can start the KMC simulation from a desired initial lattice configuration and the number of gas molecules specified up to additive constants. Then, at each time step of the simulation, the KMC algorithm scans the lattice and detects all the elementary events that can happen on it. Based on this detection, a lattice queue of lattice processes/events is created, and the next elementary event is selected. The execution of such an event is accompanied by the appropriate update of the lattice configuration and time clock. The lattice process queue is then updated so that the elementary events that cannot happen anymore are deleted and newly appeared events are added. The simulation ends once a predetermined termination criterion is fulfilled (e.g., the KMC simulation time exceeds a target value). The information contained in a KMC trajectory can be post-processed to obtain, for instance, catalytic performance metrics, calculate surface coverages, or perform reaction pathway analysis (see Fig. 4 for examples of some typical KMC outputs).²⁴ 

A large part of the computational cost in a KMC simulation comes from selecting the elementary event to be executed and updating the data structures carrying information about the possible lattice events and lattice configuration. The way to decide which event will happen next and when gives rise to two broad algorithms, namely, the null-event and rejection-free algorithms.^24,55,101 In the former, not all selected elementary events are realizable, and if a non-realizable event is selected, the KMC time is advanced but no event is executed. On the contrary, in the latter, all selected elementary events are executed. Although these algorithms are, in principle, exact and yield statistically equivalent results, they typically do so at different computational costs, as will be discussed in more detail later (see Sec. IV A).¹⁰¹ The null-event algorithm has minimum bookkeeping and was one of the early approaches employed to simulate catalytic surface reactions. However, it becomes quite inefficient when the number of null events is significant, which is quite likely when dealing with realistic systems having multiple site types, many-site events, and mechanisms with several elementary steps. Additionally, KMC clock updates in a null-event algorithm are complicated in the presence of lateral interactions, since one would have to consider all possible spectator arrangements for each reaction event. Due to these issues, the null-event method has given way to the rejection-free method, which has become, over the years, the gold standard in modeling the kinetics of heterogeneous catalysts.

The rejection-free algorithm is implemented either by the direct method (DM) or the first reaction method (FRM).^24,55,102 The DM generates the next elementary event and its occurrence time directly from Eqs. (6) and (7). On the other hand, the FRM generates a putative time for each possible elementary event using Eq. (5) and selects, as the next event, the one having the shortest time of all. Both methods are equivalent and reproduce the statistics of the master equation exactly (up to numerical accuracy, of course).^55,101,102 However, their efficient implementation depends heavily on the data structure employed to handle the selection of the elementary event, the execution of this event, and the post-execution updates. The latter can also be computationally intensive, especially for systems with long-range lateral interactions or, more generally, complicated energetics.

Detailed discussions on KMC algorithms have been carried out by several authors.^24,55,101 Therefore, rather than extending our discussion here, let us move on to comment on some of the KMC software packages and codes that build on these algorithms in Sec. III B. Readers interested in practical guides for constructing and evaluating KMC models for the simulation of catalytic surface reactions are referred to, e.g., the recent tutorials by Andersen et al.²⁸ and Prats et al.²⁷ 

KMC simulations of surface reactions have been steadily gaining popularity since KMC was coupled with first-principles calculations in the late 1990s.^72,85 Although the computational packages that perform first-principles calculations are more numerous and mature, at the present time, several general-purpose KMC software packages are at the disposal of researchers working in the areas of heterogeneous catalysis and surface science.

To the best of our knowledge, the first two general-purpose KMC software packages that were developed are CARLOS^55,103 and SPPARKS^104,105 (Stochastic Parallel PARticle Kinetic Simulator). CARLOS was probably the first user-friendly KMC code primarily devoted to simulating molecular phenomena on catalytic surfaces.^55,103 CARLOS employs null-event and rejection-free algorithms and is optimized for memory and speed utilization. The influence of lateral interactions on the system’s catalytic performance is treated by explicitly defining the different rate constants for the various configurations of the lattice. On the other hand, SPPARKS has been mostly implemented to investigate material science problems and biochemical reaction networks.^104,105 This KMC code was the first to implement parallelization techniques and spatial decomposition of the simulation domain. SPPARKS uses null-event and rejection-free algorithms as well and is designed to facilitate its modification and extension. The code has been sporadically used to study catalytic systems. Examples include the hydrogenation of benzene on Pt(111)¹⁰⁶ and the catalytic CO oxidation on Pt surfaces.¹⁰⁷ 

The last decade has witnessed the appearance of new powerful KMC software packages and codes. We briefly describe some representative such tools, noting that our list is not intended to be exhaustive. The first in this list is the versatile software package Zacros,¹⁰⁸ which was developed at University College London (UCL) and made available in 2013. Zacros is written in Fortran 2003, implements the FRM, and couples the so-called graph-theoretical KMC (GT-KMC) framework with CEHs for the adlayer energetics.^{22,24,26,30,84} Moreover, to deal with the environment-dependent activation energies of elementary events, the code uses BEP relations that enable the calculation of rate constants on-the-fly (OTF). With Zacros, one can perform simulations that consider species that bind to more than one active site, complex reaction patterns, spatial correlations and ordering arising from lateral interactions involving many-body contributions, as well as changes in the activation energies of elementary events due to interactions of reactants and spectators. Zacros applies OpenMP parallelization along with a sophisticated caching scheme to tackle the high computational cost arising from the presence of lateral interactions, especially when the CEH involves long-range interactions or a large number of contribution terms.^87,109 In addition, the code implements domain decomposition and lookahead-rollback parallel algorithms, built on the message passing interface (MPI) framework, to enable the distributed simulation of very large domains.^82, Zacros can be used for KMC simulation but is also capable of performing equilibrium MC simulations.^110,111 An interactive post-processing and visualization tool is also available (Zacros-post graphical user interface), and a Python wrapper has been developed and is made available independently by Núñez et al.^112,113 This wrapper can run multiple simulations with parallel processing, rescale rate constants of fast equilibrated reactions to accelerate the simulation, and perform parametric sensitivity analyses.

About a year after the release of Zacros, the software packages kmos and KMCLib were introduced. kmos uses the DM to generate KMC trajectories.^28,114,115 The code is an application programming interface (API) that simulates reactions involving multiple active sites and species that bind to more than one site. Although the standard backend of kmos handles lateral interactions in a way similar to CARLOS, to run complex catalytic systems containing lateral interactions, kmos also implements the OTF backend, which uses BEP relations to modify rate constants during runtime.¹¹⁴ On the other hand, KMCLib is a general KMC framework that can simulate the time evolution of systems of up to millions of particles in one, two, or three dimensions.^116–118 The code uses MPI parallelization to detect elementary events on the lattice and calculate the associated rate constants. KMCLib uses the DM as well.

The software package MoCKA (Monte Carlo Karlsruhe) also appeared around the same time as Zacros, kmos, and KMCLib.¹¹⁹ MoCKA provides a general KMC framework combining graph-theoretical and lattice-based approaches to model molecular phenomena on nanoparticles. The code can efficiently handle the simulation of nanoparticles exposing several nanofacets, which may exhibit different chemical environments. To this end, the code uses multiple lattices to represent the nanofacets and handles the communication effects between catalytic particle nanofacets or between support and particle facets. MoCKA uses an efficient implementation of the DM and makes the FRM also available. MoCKa addresses only pairwise lateral interactions but, in principle, could be extended to three-body and higher-order interactions as they can be derived from CEHs.

A recent addition to the list of user-friendly software packages is MonteCoffee.¹²⁰ This code is an open-source object-oriented programmable (OOP) application, which exploits similar ideas to the graph-theoretical (or GT) approach and has the advantage of quick extension. Several challenging tasks are possible with this code, such as evaluation of lateral interactions, sensitivity analysis (SA), and descriptor-based energy landscapes. The code uses neighbor lists to represent the lattice connectivity. The FRM is the primary KMC simulation driver, but the code can be extended in a straightforward way to incorporate other algorithms.

Moreover, Excimontec^121,122 is a recent KMC software package developed in modern C++ and optimized for efficient execution on high-performance computing clusters using MPI. It allows users to create simulation models on a cubic lattice and combines the DM and several variations of the FRM. This software package uses object-oriented design and extends the KMC_Lattice framework.¹²³ Finally, the software package SuSmoST (SUrface Science MOdeling and Simulation Toolkit)¹²⁴ consists of computer programs and libraries intended to support surface investigation, focusing on adsorption systems and phase transitions. The code implements mainly transfer-matrix and tensor-network methods, as well as Metropolis Monte Carlo simulations, with KMC support as a recently added feature.

The computational codes and software package just described are broadening the adoption of 1p-KMC simulations for heterogeneous catalysis. Although this could indicate certain maturity reached by this simulation technique, several computational challenges still have to be overcome to fully exploit its predictive power. Key such challenges and recent attempts to address them are the subjects of Sec. IV.

KMC is a numerical simulation framework, and as such, it faces computational challenges when dealing with realistic catalytic systems. These challenges, but also the increasing popularity of 1p-KMC as a tool for understanding and predicting catalytic performance metrics, have motivated the development of novel methods, as discussed in this section.

As mentioned in Sec. III A, at every step of a KMC simulation (i.e., KMC iteration), we have to create a list of all of the possible elementary events that may occur on the lattice. Then, we need to randomly select one of these events and its time of occurrence and execute it. Finally, we conduct the necessary updates on both the list of possible events and the lattice configuration. These procedures are at the heart of any KMC simulation and can become a significant computational bottleneck for systems with many elementary events. The bookkeeping of these procedures depends on the KMC method implemented and is usually performed employing appropriate data structures.

In a series of early studies aiming at finding efficient KMC methods for the simulation of catalytic surface reactions, Jansen and co-workers implemented binary trees and binary search to store the elementary events of the KMC simulation.^{55,81,125–129} In a similar effort, Reese et al.¹³⁰ explored the implementation and efficiency of rejection-free and null-event approaches with local and global search and update algorithms. In a global such algorithm, the entire library containing the information of interest is sequentially searched and updated, whereas in a local one, only a portion of the library is searched and updated. Reese et al. concluded that when the reaction mechanism involves events occurring at disparate frequencies, e.g., fast diffusion, slow reactions rejection-free algorithms can be much faster than null-event ones (this disparity of frequencies/time scales is referred to as “stiffness;” see Sec. IV D for more details). However, Reese et al. also found that null-event algorithms could, under some conditions, outperform rejection-free approaches due to the fewer operations per KMC iteration that the former entails. Along similar lines, Dooling and Broadbelt explored a KMC tool that maintains detailed information about the current state of the catalytic surface and only updates the information locally, thereby greatly increasing the efficiency of the simulations.¹³¹ 

In a general context beyond heterogeneous catalysis, Gibson and Bruck employed a dependency graph for locally updating only the rate constants that are actually affected when a particular event is executed.¹⁰² Furthermore, they implemented an indexed priority queue to identify the elementary events with the least waiting time when using the FRM. Such an indexed priority queue enabled the retrieval of an arbitrary elementary event in constant time (though, of course, updates are more costly). Moreover, for both the DM and FRM, a complete tree data structure was proposed in which the amount of time per KMC iteration is proportional to the logarithm of the number of elementary events on the lattice, not to the number of events itself.

More recently, Chatterjee and Vlachos¹⁰¹ reviewed several KMC simulation methods and explored the implementation and efficiency of several search and update algorithms.¹²⁸ They concluded that rejection-free approaches, which use efficient search and update algorithms, are very efficient but require much more complicated coding than null-event approaches. They reported that implementations involving linear search on the data structure scale as O(N_queue), where N_queue is the total number of elements in the queuing data structure. They also found that extensions of the linear search, such as the two-level and n-level linear searches, scale as $O(Nqueue)$ and $nO(Nqueuen)$⁠, respectively. The so-called binary and hash-table searches^80,128 were also discussed in this review. This work also discussed techniques for updating the lattice configuration, lattice energetics, and rate constants.^102,130

In a latest effort, Savva and Stamatakis implemented and compared four data structures (i.e., the unsorted list, the binary heap, the pairing heap, and the one-way skip list),^128,133,134 as alternative queuing systems to handle the elementary event queue during the implementation of the FRM, and further developed the two-way skip list-based queuing data structure [see Fig. 7(a) for the schematic representation of three of the data structures implemented by Savva and Stamatakis].¹³² These five approaches were implemented in the Zacros software package and benchmarked against one another using a CO oxidation model adapted from the seminal work by Ziff, Gulari, and Barshad (the ZGB system),¹³⁵ a simplified model of the water–gas shift (WGS) reaction on Pt(111),¹³⁶ and a TPD model of CO on Cu(111).¹³⁷ They also investigated the effect of compiler optimizations on the performance of these data structures. This work found the unsorted list to be, as expected, impractical. However, as presented in Fig. 7(b), they observed a 3× speedup of the binary or pairing heaps compared to the one-way skip list. They also found that compiler optimization delivers a speedup of up to 1.8×.

Since a KMC simulation involves the iterative scheduling and executing elementary events, efficient implementations of these operations are of paramount importance. A fundamental part of this endeavor consists in adopting the latest approaches or developing novel ones. In this respect and to the best of our knowledge, data structures, such as the Brodal queue or the Fibonacci heap, have not yet been implemented for kinetic simulations of catalysts. Both approaches are known to have better time complexities than the binary or pairing heaps, though large memory requirements and high constant factors might make them competitive for only specific applications. Therefore, these non-trivial considerations might represent an opportunity for further algorithm development and optimization.

In reaction mechanisms for which adsorbate–adsorbate lateral interactions are significant, the non-ideal mixing in the adlayer could result in the formation of islands and ordered structures on the catalytic surface. Such phenomena can dramatically affect the observed reaction kinetics and should be properly taken into account when performing 1p-KMC simulations (see Sec. II B).

The incorporation of lateral interactions into a comprehensive kinetic modeling framework has been a long-standing challenge. Even at the present time, for simplicity and computational efficiency, such interactions are often disregarded altogether with the use of a “hard-sphere” (lattice) model, also referred to as the “site blocking” model.²⁸ The minimum-effort approach beyond this approximation considers simple pairwise additive nearest neighbor (1NN) interactions within a lattice-gas model.^138,139 In a more elaborate approach, early 1p-KMC simulations implemented various versions of the so-called bond-order conservation (BOC) method to account for 1NN and longer-range interactions in a realistic way.^72,85 Currently, a popular methodology to accurately treat lateral interaction within the 1p-KMC framework is the CEH approach,^{30,87–95,140} which is based on a general formalism for representing a real-valued function of a lattice configuration variable.

In this approach, the electronic energy of the lattice of a state is represented by a sum of interaction clusters (patterns). Each of them represents an energetic interaction pattern or figure in which the adsorbates are arranged, for instance, single body, two-body, three-body, and larger contributions. Interaction patterns typically included in the different CEH models are presented in Fig. 8(a). While this approach enables, in principle, the most accurate modeling of lateral interactions, it dramatically increases the computational cost (or runtime) of KMC simulations. For this reason, a number of strategies to reduce such computational burden have been implemented in recent years.

Parallelizing the detection of the interaction clusters has been a successful approach to this end. For instance, Nielsen et al.⁸⁷ developed an implementation of this idea, which uses open multi-processing (OpenMP)¹⁴¹ in the update of kinetic constants after the execution of a reaction event. By simulating a kinetic model for NO_x oxidation/reduction on Pt(111) as a benchmark, this approach was found to significantly reduce computational time by about 10× for 16 threads with a CEH containing up to eighth-NN (8NN) interactions. However, the results of these benchmarks also revealed that the performance improvement reaches a plateau as the number of threads increases. This plateau was attributed to a synchronization overhand when collecting the updated rate constants from the different threads. Motivated by the need for further improvement, Ravipati et al.¹⁰⁹ recently developed a novel and exact scheme that implements a sophisticated caching data structure along with OpenMP for faster updating of the rate constants. The approach is based on caching information about the energetic interaction patterns associated with the products of each possible elementary event. The scheme was also benchmarked for the same NO_x oxidation/reduction system and yielded acceleration factors of up to 20× when comparing single-thread runs without caching to runs on 16 threads with caching for simulations with up to 8NN interactions [Fig. 8(b)]. The two approaches just mentioned have been implemented in the software package Zacros. It is worth noting that OpenMP parallelization alone merely distributes the workload to several threads (with the unavoidable overhead), while only the second approach (caching) actually reduces it. However, due to the complicated nature of the caching algorithm, no formal estimates of computational complexity have yet been derived.

Aiming at reducing the computational complexity of 1p-KMC simulations with lateral interactions, Hess recently proposed the supercluster approach, subtraction schemes, and the supersite algorithm.⁹⁵ The supercluster approach improves the computational performance of evaluating the CEH by essentially reducing the number of terms in the CEH, and when tested on several compilers and central processing unit (CPU) architectures, it showed a significant speedup compared to the traditional CEH approach. The subtraction scheme was introduced to optimize the calculation of the sum of rate constants when implementing the DM. This scheme exploits the well-established strategy of performing local updates, i.e., recomputing rate constants only in the vicinity of an elementary event and within lateral interaction range. This way, the recalculation of the total rate constant scales as O(1) with the number of sites of the lattice. Finally, the supersite search implements Maksym’s two-level algorithm,¹⁴² as opposed to the traditional linear search implementation of the DM,¹⁰¹ in which the sites of the lattice are visited successively. Overall, it was demonstrated that one must combine the three algorithms to reach the best performance. In such a case, it was found that the computational cost when including lateral interactions increases by less than a factor of 3 as compared with the case without lateral interactions. The performance studies were mainly carried out with simple models, However, simulations of a full model of HCl oxidation over RuO₂(110) confirmed the improvements in computational savings.

At the present time, lateral interactions are an integral part of most of the existing KMC software packages. However, the computational burden associated with its implementation remains quite high. Although a number of acceleration schemes have appeared in recent years,^87,95,109 further developments are necessary.^{38,86,90,91,94,143,144} Along these lines, further improvements in the efficiency of lateral interaction pattern detection, e.g., via fingerprinting and precomputing, could be a way forward.

Conventional KMC algorithms are sequential in the sense that they handle one elementary event at a time. For this reason, KMC does not scale well with the lattice size. In particular, even if the execution of a KMC event happens at constant time (i.e., is independent of the lattice size), thermodynamic scaling laws imply that the number of events needed to reach a target/final KMC time scales linearly with system size. This poor scalability has limited the application of 1p-KMC simulations to relatively small lattices, typically on orders of up to 10 × 10 nm². Although KMC simulations of such lattices can yield good estimates of catalytic performance metrics, understanding long-range spatiotemporal phenomena, such as pattern formation in the context of surface reconstruction, necessitates the simulation of larger lattices. The size of these lattices is dictated by the wavelength of the pattern, which could be in the order of μm to mm.^145–147 These large-sized lattice simulations pose challenges in terms of computational time and memory footprint.¹⁰¹ With the hardware available at present, using the sequential KMC framework to simulate complex chemistries on catalytic surfaces that span μm in size is impractical or infeasible.

A way to address this challenge is by implementing distributed parallelization techniques for discrete event simulation.^148–158 Applications of these techniques are possible because of the typically localized character of the elementary events on the catalytic surface. Among these approaches, we note as representative the synchronous,¹⁴⁸ the synchronous relaxation,¹⁴⁹ the optimistic synchronous relaxation,^150,151 the optimistic Time-Warp,¹⁵⁹ and the semi-rigorous synchronous algorithms.^152–157 In general terms, they rely on domain decomposition and event execution protocols that either keep the different processors synchronized as the simulation proceeds, thereby avoiding boundary conflicts, or carry on simulating in an “optimistic spirit” until such conflicts arise, at which point they employ rollback and re-simulation procedures to rectify the pertinent causality violations. Thus, in the latter approach, provisional KMC trajectories are progressively amended and finally incorporated into the global simulation history, provided they are mutually consistent.

Although an approach similar to the semi-rigorous synchronous algorithm is already implemented in SPPARKS^{104–106,154} and features of the Time-Warp algorithm have also been incorporated into the software package SPOCK (Scalable Parallel Optimistic Crystal Kinetics),¹⁶⁰ it was only recently that a general and validated approach was made available to the computational physics community by the coupling of the optimistic Time-Warp algorithm with the GT-KMC framework, as implemented in Zacros.⁸² In the Time-Warp-GT-KMC framework, the lattice is decomposed into domains that are assigned to different “processing elements” (each of which may be handled by a single CPU core or involve several threads). Each processing element carries out the KMC algorithm for the assigned domain and communicates with the neighbor domains, if necessary. This way, elementary events that occur far from domain boundaries are handled privately and asynchronously, while events close to those boundaries must be appropriately communicated. The communication among processes is handled by using the message passing interface (MPI) framework.¹⁶¹ Because each processing element follows its own simulation time, when boundary events happen, causality violations may arise. Such violations are resolved using the protocols of the Time-Warp algorithm (i.e., state-saving, rollback, and message/anti-message sending protocols). Further protocols for the computation of the collective KMC time of all processing elements (or global virtual time) and the corresponding termination of the distributed run are implemented as well. Scaling benchmarks of the Time-Warp-GT-KMC framework revealed that the overhead associated with the Time-Warp-related procedures can be significant. However, the approach was found to scale well with the lattice size and outperformed the sequential KMC for large lattices. The simulation runs were performed for two simple models and a more realistic one that captured CO oxidation dynamics on Pt(111) and made use of an elaborate CEH energetic model (Fig. 9).⁸³ Crucially, Ravipati et al.⁸² devised a serial simulation protocol that delivers, by construction, the same results as the Time-Warp-GT-KMC, thereby validating the implementation in a general way.

Approaches such as the Time-Warp-GT-KMC are very promising and will allow the investigation of heterogeneous catalysis at spatial and temporal scales of unprecedented magnitude. The research topic of distributed parallel KMC for heterogeneous catalysis and surface science applications is still in its infancy, and we envision further exciting developments, both from a method development and from an applications standpoint, in the years to come. Thus, the development and implementation of more efficient exact simulation protocols, as well approximate schemes whose accuracy can be assessed by comparisons with exact ones, could pave the road for unraveling the molecular origins of emergent behaviors at the macro-scale, such as pattern formation on catalytic surfaces,¹⁴⁶ and lead to novel applications in the area of spatiotemporal control of reactive systems.¹⁶² 

By performing simulations on the long time scale of the elementary events, KMC achieves an enormous speedup over MD simulations. However, a time scale disparity problem can emerge even at that level of description because, in many cases, these events can occur at vastly different time scales (i.e., their rate constants can span multiple orders of magnitude). For instance, it is common to encounter surface reaction mechanisms that involve very fast surface diffusion events but also slow reactions. Then, because in KMC simulations, elementary events are selected based on their relative rate constants (see Sec. II A), it is far more likely that a diffusion event will occur at each KMC step. Typically, these slow reactions are the events of interest (i.e., the ones directly related to the catalytic performance metrics). In contrast, the fast diffusion events lead to a quasi-equilibrated adlayer but do not contribute to the net evolution of reactants and products. Furthermore, because these fast (frequent) diffusion events are the limiting factor in the advancement of the simulation, numerous KMC steps will be needed in order to sample over a sufficient number of reactions so as to obtain accurate estimates of catalytic activity or other metrics. This problem, which is also known as KMC stiffness (similar to stiffness in differential equations),⁴⁰ can severely undermine KMC simulation performance and has motivated the development of acceleration schemes that address it.

An early such scheme is the so-called absorbing Markov chain KMC (AMC-KMC) framework.^163–165 Fundamental to this scheme is the concept of a superbasin, i.e., a set of basins (KMC states) connected by fast (low energy barrier) quasi-equilibrated elementary events. The AMC-KMC framework formulates the escape from a superbasin as an absorbing Markov chain, with the repeatedly visited states within the superbasin referred to as “transient states” and the states outside the superbasin, to which the system exits, called “absorbing states.” This framework is exact and relies on creating the so-called Markov matrix to describe the transitions among transient and absorbing states. The determination of the absorbing state that the system enters (thereby escaping the superbasin), as well as the time when this happens, requires the diagonalization or inversion of the Markov matrix. AMC-KMC works well with and without time scale separation and has been implemented in the software package EON for off-lattice (adaptive or OTF) KMC simulations.^43,166–171 However, a downside of this approach is that the matrix operations can become computationally intensive when dealing with large Markov matrices (i.e., large superbasins). Moreover, the method does not explicitly specify how to group states to form a superbasin.

An approximate method that does not require expensive matrix manipulations is the accelerated superbasin KMC (AS-KMC) introduced by Chatterjee and Voter.¹⁷³ AS-KMC relies on the construction of a database of elementary events and counting the number of times each event has occurred. A larger number of occurrences for an event indicates that a superbasin may be present. If the system is considered to be in a superbasin, the rate constants of all quasi-equilibrated elementary events that connect the basins of the superbasin are lowered (by, e.g., raising the activation barriers). This scaling increases the probability for a non-equilibrated elementary event to be selected, leading to a transition to a new superbasin (see Fig. 10). In doing so, AS-KMC offers a significant computational speedup over “standard” KMC simulations, incurring an error that has been shown to increase asymptotically linearly with the factor used to downscale the fast events.¹⁷⁴ This factor can be taken such that the error of the approximation is negligible compared to the error of the KMC sampling. However, a potential drawback of this method is that the identification of even a single superbasin can be computationally too expensive for complex catalytic systems, where the state space is typically so large that many states may only be visited very infrequently during the entire simulation.

More recent KMC acceleration methods are essentially based on the fundamental concepts and approximations just discussed and largely focus on overcoming the computational shortcomings of the AS-KMC framework. For instance, Dybeck et al.¹⁷² introduced a method that automatically partitions the elementary events into equilibrated and non-equilibrated and scales down the rate constants of the quasi-equilibrated ones (i.e., the elementary events inside the current superbasin). The partitioning and scaling are applied collectively to all lattice events in a given reaction channel rather than to the individual lattice events (a reaction channel includes all the events of, e.g., adsorption of some species at a given site type that have been detected on the lattice, independently of where they occur). The classification of events according to reaction channels allows identifying those events efficiently. Once a transition to a new superbasin is executed, the rate constants are unscaled again to their original values to allow for sufficient sampling of the new superbasin. This approach was used to model the complex Fisher–Tropsch synthesis reaction over ruthenium nanoparticles and exhibited computational saving spanning several orders of magnitude. In this work, for simplicity, surface diffusion and lateral interactions were not taken into account.

Following the development of the AS-KMC framework by Dybeck et al.,¹⁷² Andersen et al.^13,28 implemented it on the software package kmos and applied it for the simulation of CO methanation over stepped transition metal surfaces using scaling-relation-based rate constants. Although the scheme was found to perform quite well, some challenging cases were encountered, for which the scheme exhibited artifacts. For instance, problems may arise when simulating reactions between two low-coverage species, which are both produced by independent quasi-equilibrated elementary events. In this case, the method may scale the rate constants of the quasi-equilibrated events too aggressively for an adequate sampling of states where the two low-coverage species are found at neighboring lattice sites for their reaction to proceed. The problem was linked to the fact that the method does not track system states and therefore cannot verify if all states within a superbasin have been sufficiently sampled. To address the breakdown, Andersen et al.¹³ proposed a correction to the original version of the scheme that takes into account lattice configurations of the nearest neighbor sites in the definition of the reaction channels. However, such a corrected version was shown to work well only for the simple case of a reaction between two low-coverage species formed directly at neighboring sites. The correction does not apply if the low-coverage species are created at distant lattice sites and rely on diffusion events in order to “meet” and react.

Other methods, similar to the ones described above, have appeared in the literature. For instance, Núñez et al.¹¹³ addressed the challenge by employing rate constant rescaling techniques. Motivated by the problem of establishing a direct correspondence between the MF-MKM and KMC, Hoffmann and Bligaard presented a KMC acceleration scheme in which, before rescaling the rate constants of the fast elementary events, the algorithm automatically makes sure the system is at steady-state.¹⁷⁵ A common feature between these two schemes is that the rate constants of fast elementary events are not restored (“unscaled”) once the system enters a new superbasin.

The “staggered quasi-equilibrium rank-based throttling for steady-state” (SQERTSS) algorithm is another recent example in which the idea of rate constant rescaling is implemented to accelerate KMC simulations.¹⁷⁶ The SQERTSS algorithm was designed to decrease the occurrence of the rapid elementary events that do not significantly progress the system toward a steady-state and increase the occurrence of slower but relevant ones. A variant of SQERTSS is the “staggered quasi-equilibrium rank-based throttling geared toward transient kinetics” (SQERT-T) algorithm, which deals with ameliorating the time disparity problem in the context of transient kinetics simulations.¹⁷⁷ The SQERT-T algorithm is helpful in the context of simulating temperature programming desorption and temperature programming reaction.

The accelerated KMC algorithms described above automatically treat the time scale disparity problem (to the extent that this is feasible). However, a manual adjustment and verification of the fast elementary events is also possible, provided the catalytic surface reaction is simple enough.^83,178,179 There exist also algorithms that require the user to supply in advance information about fast and slow elementary steps or identify the quasi-equilibrated ones.^180,181 For surface reaction mechanisms with a small number of fast diffusion elementary events coexisting with slow reactions, such a piece of information is, in principle, easy to establish. However, having such information in advance does not always help because the quasi-equilibrated elementary events may change as the simulation progress. In other cases, the system can be too complex for the user to correctly classify the elementary events into fast vs slow.

A method that adopts the strategy of separating in advance fast diffusion events and slow reaction events is the recently introduced “extended phenomenological kinetics” (XPK) approach.^182–187 This approach is a hybrid between a diffusion-only KMC explicitly taking place on the lattice, which enables the evaluation of the reaction propensities for elementary events, and a subsequent implicit-lattice KMC that evolves the number of species on the surface. Another method that addresses the problem of fast diffusion is the recently introduced fast species redistribution (FSR) method.^188,189 The FSR-KMC method builds on the idea that fast diffusion steps hardly affect the system’s time evolution but result in the redistribution of the fast diffusive species. A prerequisite for applying the FSR-KMC method is the identification of the fast and slow diffusion elementary events, which is performed automatically. The FSR-KMC method takes into account detailed spatial and energetic information to ensure a proper species redistribution. Finally, accelerated KMC algorithms to address KMC stiffness in the presence of catalyst deactivation by site blocking have been also introduced in the literature.¹⁹⁰ 

Although we have focused our discussion on accelerated KMC frameworks within the context of heterogeneous catalysis, the time scale disparity problem is a generic computational challenge of the KMC framework. It is also important to reiterate that, with the exception of the AMC-KMC, the schemes discussed are approximate, and due to the heuristics involved, it is quite difficult to estimate the magnitude of the error incurred in the simulation observables. Thus, even though progress has been made and several numerical approaches and algorithms have been developed to tackle the issue of event frequency disparity, important issues remain unresolved. The challenge remains to develop an acceleration approach that provides approximate solutions within pre-specified error tolerances and that is also robust and generic or transferable, i.e., applicable to, in principle, any (type of) KMC simulation.

In KMC simulations of catalytic surface reactions, we are typically interested in estimating observables at steady-state. Typical examples of such observables are the average surface coverage of certain species and the average rate of production of a certain molecule per lattice site, referred to as the turnover frequency (TOF). We are also often required to identify the most influential rate constants of the system and understand how the uncertainty in the rate constants influences the estimates of the observables of interest; these aims are achieved via sensitivity analysis (SA) and uncertainty quantification (UQ), respectively. These two aspects are important in assessing the accuracy of 1p-KMC simulations, and since the pertinent analyses have to be done at “steady-state” conditions, identifying the latter is also important in this context.

It is important to stress that the term “steady-state” is used somewhat loosely in our discussions (and in the cited literature), and such a “state” should not be confused with any of the KMC states visited during the simulation. The mathematically acceptable term referring to time-invariant conditions in the context of random processes is “stationarity,” and several definitions thereof exist, capturing different levels of rigor. For instance, if the mean and autocovariance are constant and the second moment is finite at all times, then the stochastic process is referred to as “wide-sense stationary.” The much more restrictive definition of a “strictly stationary” processes requires all joint probability distributions, $Pω(t1)$⁠, $Pω(t1),ω(t2)$⁠, … to be time-invariant. Of interest, in practice, is the probability of the state variable at stationary conditions, $Psts$⁠, which is known to be the long-time limit of P_ω(t), at least for “well-behaved” random processes.

Thus, for an ensemble of KMC runs starting from the same state (e.g., empty lattice), the sharp probability corresponding to the initial condition evolves toward the stationary probability. Once “steady-state” (stationarity) has been reached, it is possible to invoke the assumption of ergodicity to calculate an average quantity of interest as a time average from a single KMC run instead of an ensemble average. This is the strategy typically followed in practice; however, when performing multiple simulations, e.g., for parametric analyses, verifying that the system has reached steady-state can be tedious. Indeed, the “steady-state” behavior of the system depends on the pressure, the temperature, and the composition of the gas phase in contact with the catalyst. Given the large number of possible such conditions, the number of KMC steps required to reach steady-state can extend over a wide range, and manual verification of steady-state behavior for each simulation would be impractical. What would be needed is to implement parallel processing for simultaneously conducting large numbers of simulations, coupled with automated algorithms able to decide with enough confidence when a KMC simulation has reached steady-state and exclude the transient period from sampling. As an example, Fig. 4(c) shows a fluctuating KMC trajectory exhibiting a transient up to about 10 s, before approaching what appears to be a “statistical steady-state” behavior. Multiple types of different fluctuating profiles like that can be produced by KMC simulations.

It is only very recently that attempts to develop robust and computationally efficient steady-state detection (SSD) algorithms for KMC simulations have been made, based on the implementation of appropriate statistical tests. For instance, Núñez et al.^113,191 implemented a criterion for having achieved steady-state, which uses a batch means test complemented with a t-test. In another work, Nellis et al.¹⁹² developed the F-t-Pj-RG method, which relies on the subsequent application of an F-test, a t-test, and a projection test on adjacent windows of the time series while rolling and growing the windows when any of the tests fail. Passing the F-t-Pj-RG test requires that all three individual tests pass. Upon passing the test, a statistical steady-state is considered to have been detected within a confidence level and tolerances provided by the user. The F-t-Pj-RG method determines the appropriate window size on the fly and does not require advance empirical knowledge of the system under consideration. The computational cost of the method was proved to be on average less than one percent of the computational cost of the KMC simulations, and it was also shown to be suitable for various types of data trends that can occur in real KMC applications. Despite all that, further testing will be needed to verify whether the F-t-Pj-RG method can be applied in more general settings.

In recent years, the development of comprehensive and accurate SA and UQ approaches to 1p-KMC simulations of catalyzed surface reactions has also gained attention. This development is motivated by the high computational cost of calculating rate constants from DFT calculations. Ideally, one would like to know in advance how those rate constants influence the simulated behavior of the catalytic system. Then, one could dedicate computational resources to the accurate calculation of the most influential kinetic constants using expensive electronic structure methods, while the remaining kinetic constants could be estimated by low-cost methods. In addition, such analyses can help us quantify the propagation of errors from the electronic structure calculations of choice to the quantities of interest.¹⁹³ In this respect, both local and global SA and UQ approaches have been advocated.

One of the most popular methods of local SA in heterogeneous catalysis is the Campbell’s degree of rate control (DRC) that evaluates the partial derivative of the average reaction rate with respect to a rate constant parameter while keeping other parameters constant.¹⁹⁴ Despite its successful application to first-principles MF-MKM (1p-MF-MKM), the implementation of such a local SA method is computationally demanding in KMC simulations. In this respect, Meskine et al.¹⁹⁵ were among the first to explore the usefulness of implementing various definitions of the DRC to understand the propagation of errors from electronic structure calculations to 1p-KMC simulations. Although the results of this seminal work were encouraging, the application of DRC to perform SA has been hampered by the high computational effort required to accurately sample numerical derivatives of quantities obtained from stiff KMC simulations. To address these issues, Hoffmann et al.¹⁹⁶ proposed an efficient and robust approach that was proven capable of performing a reliable evaluation of the sensitivity measures for stiff KMC simulations. Moreover, local SA of KMC simulations based on log-likelihood estimators have been recently carried out. For instance, Nuñez and Vlachos introduced the singularly perturbed steady-state likelihood ratio (SP-SSLR) method to address the problem of SA in stiff KMC simulations.^113,191 Techniques based on the likelihood ratio reduce the computational cost of sensitivity analysis by obtaining all gradient information in a single KMC run. A novel element of this method is also that it enables the study of how parameter influences propagate from fast time scales to slow dynamics.

In the local SA approach, the results of sensitivity measures are valid locally in the input parameter space (i.e., an assumption of linearity is made). However, kinetic models of catalytic surface reactions often exhibit highly non-linear behavior and strong parameter correlations. For such cases, global approaches that focus on the entire parameter space of the system are needed. These approaches have been applied mostly to 1p-MF-MKM models but also to KMC simulations. For instance, by performing a correlative derivative-based global sensitivity analysis, Sutton et al.¹⁹⁷ showed that neglecting correlations in reaction and species’ energies can lead to incorrect identification of key reaction intermediates and influential parameters when performing KMC simulations. However, the drawback of such an approach is again that derivatives are difficult to estimate from 1p-KMC simulations. In this respect, Döpking and Matera proposed a probability distribution-based approach that allows quantifying the error propagation in 1p-KMC simulation with reasonable cost.¹⁹⁸ 

The development of ever refined approaches that assess how errors from electronic structure calculations propagate to essential outcomes of 1p-KMC simulations, such as catalytic activity and selectivity, is an active and vital topic of research. The challenge is to develop dynamical schemes that tune the rate constants’ accuracy depending on the kinetic importance of the associated chemical pathways in the reaction mechanism.

Chemical kinetics has a central role in the modeling of reactors and chemical processes; yet, in practice, simplistic approaches based on empirical kinetic models are typically employed. Over the last decade, the vision of the bottom-up reactor design has been steadily gaining attention; to this end, the equations describing mass transport need to be coupled with high-fidelity kinetic models, and it is attractive to use 1p-KMC for the latter. The coupling of 1p-KMC simulations to computational fluid dynamics (CDF) or continuous models described by partial differential equations (SDEs) is a relatively young area of research that opens up interesting new avenues for method development. The ultimate goal of such a coupling is to properly integrate the descriptions of the momentum, heat, and mass transport phenomena occurring in either small-scale laboratory reactive chambers or industrial reactors, with the molecular-scale processes occurring on the catalyst surface. To this end, continuous models are used to simulate the gas-phase reaction mixture, while the KMC-derived reaction rates (or TOF values) are used in the boundary conditions of the problem. However, the stochastic character and computational cost of KMC simulations poses challenges. For instance, the noise in the KMC outputs could negatively affect the stability of the numerical methods used to model the transport phenomena, and the computational costs for obtaining those outputs may render the coupling challenging. The situation is complicated further by the inherent computational burden of CFD simulations.

This coupling can be carried out in a direct or an indirect way.^{8,10,138,199–214} In the direct coupling, at each time step of the simulation, the KMC outputs are used as the in situ boundary conditions for the continuous model. The temperature and pressure obtained from the continuous model simulations are used as the operating conditions for the KMC simulations in the next time step. Thus, this approach establishes a concurrent coupling with a direct exchange of information between the two scales. The direct coupling has been achieved for simple continuous models and/or surface reaction systems.^199–204 Most of them follow a domain decomposition approach in which the entire heterogeneous catalytic reaction system is decomposed into the catalyst surface domain, which is simulated with the lattice KMC framework, and a gas phase domain, which is captured by a continuous model. If the descriptions of these two domains are simple enough, the direct coupling is relatively straightforward.^199,200,202 However, KMC simulations of the whole surface domain may become computationally prohibitive for complicated spatially heterogeneous systems with strong concentration gradients. For such scenarios, patch dynamics approaches have been proposed, such as the so-called gap-tooth scheme.^201,203,204 In this framework, the entire surface domain is subdivided into periodically spaced patches, referred to as “teeth,” separated by spaces referred to as “gaps.” Each “tooth” is simulated as an independent KMC lattice that gives the surface coverage at its location within the domain. Then, the surface concentration across the gaps is approximated using interpolation. Such attempts for direct coupling have been carried out for steady and non-steady-state conditions in the fluid phase.^199–204 Several methods to reduce the impact of the noise of KMC simulations on the numerical stability of the schemes employed have also been proposed for this type of coupling.²⁰⁴ 

In recent years, interesting indirect coupling attempts have also been carried out, most based on the so-called steady-state approximation.^{138,205–214} The latter relies on the fact that upon a change of the gas-phase conditions (pressure, temperature, and gas species molar fractions), the surface kinetics typically relaxes rapidly to the steady-state corresponding to those conditions. Thus, under this approximation, the KMC reaction rates required at each time step for the boundary conditions of the continuous model can be replaced by the steady-state KMC reaction rates for the current gas-phase conditions. The drawback of this approach is that the constant generation of steady-state KMC simulations is a computationally demanding task for most systems of interest. A strategy to address this issue is to construct surrogate models from several pre-calculated data points obtained via independent KMC simulations. To this end, KMC simulations are first carried out to determine the steady-state reaction rate for a wide range of surface conditions. Different interpolation schemes, such as splines or Shepard interpolation variants, are then applied to the resulting data toward the development of a continuous representation, which, in turn, provides the required necessary boundary conditions for the continuous model.^209,212 This methodology has been successfully used to compare simulations with experiments (Fig. 11)^138,210 and has allowed the integration of KMC into several CDF software packages [e.g., CatalyticFOAM and Multiphase Flow with Interphase eXchanges (MFIX)].^138,209 However, traditional interpolation techniques typically suffer from accuracy and efficiency issues as the number of variables increases. In this respect, machine learning techniques, such as random forests, are already making a game-changing contribution.²¹³ 

The coupling of 1p-KMC simulations to transport models, in order to capture the local environment experienced by the catalytic surface, is still at a relatively early stage of development. Although progress has been accomplished, direct coupling approaches are still only applicable to simple systems, while indirect coupling ones necessitate advanced surrogate modeling techniques able to handle noisy datasets and high-dimensional spaces. The generation of the training dataset for the surrogate model by pre-computing can also be challenging in such cases, and it may happen that a large portion of the dataset corresponds to conditions that are never actually realized in the reactor-level CFD simulation. As a potential remedy, OTF training of the surrogate model (while the CFD simulation is running) could reduce the number of samples needed for training and focus on the relevant conditions, thereby improving efficiency. However, this could introduce computational challenges if the KMC simulations needed to generate new data points are quite slow. The acceleration methods discussed in Secs. IV A and IV D could address such issues. Overall, the efforts toward the efficient coupling between 1p-KMC simulations and full CFD descriptions are expected to draw from multiple fields, including machine learning, software engineering, and high-performance computing.

The 1p-KMC framework has played a key role in helping to elucidate the underlying chemical mechanisms and physical phenomena of many heterogeneous catalytic reactions of practical relevance. We will briefly discuss some of these contributions in this section, with special attention to the last five to six years. Readers interested on older contributions can consult Refs. 22 and 24.

As discussed in Secs. II B and IV B, lateral interactions typically lead to the formation of non-ideal adlayers that strongly influence the kinetic behavior of catalytic surface reactions. Such interactions are nowadays mostly modeled by the CEH approach, and CEH-based 1p-KMC (1p-CEH-KMC) models have already demonstrated the accurate simulation of surface reactions of increasingly complexity. In this section, we will comment on some representative works that have implemented such models to rationalize phenomena induced by non-ideal adlayers and lateral interactions.

Recent studies by Stamatakis and Piccinin explored the impact of lateral interactions on the CO oxidation on Pd(111),^83,178,215 a model catalytic reaction with various practical applications, including the treatment of automotive gas exhausts. These studies rationalized the different reaction orders with respect to oxygen coverage at different temperatures that were observed in the context of titration experiments.²¹⁶ In the latter, oxygen-precovered Pd(111) surfaces were exposed to gas CO, and the CO oxidation rates as well as the adlayer ordering were monitored as the reaction progressed (see Fig. 12 for a schematic of the process). The reaction was shown to exhibit half order kinetics with respect to O coverage at high temperatures and first order kinetics at low temperatures, a behavior that was explained on the basis of island formation at high temperatures. This explanation was, however, somewhat counter-intuitive, since if island formation were possible, it should be favored at the low temperature range. The 1p-CEH-KMC simulations were shown to reproduce the experimental orders with remarkable accuracy, and analysis of the results demonstrated that subtle coverage effects due to lateral interactions are responsible for the observed reaction orders, contrary to the previous explanation, based on island formation.

The role of lateral interactions on the catalytic oxidation of NO to NO₂, another model reaction and a key step in NO_x elimination from exhaust gases, has also been recently explored by the CEH approach.⁹³ In this work, adlayer structure and lattice size effects on catalytic rates were explored in detail and, among others, it was found that highly ordered overlayer structures due to short-range lateral interactions were responsible for lattice size effects on the catalytic rates. As shown in Fig. 13, it was also found that if the lattice was non-commensurate to the lowest energy adlayer structure, defective regions and anti-phase boundaries emerge, which give rise to short-lived, highly reactive configurations that promote catalytic activity.

1p-CEH-KMC simulations have also shed light on the role of lateral interactions on the experimentally observed promoting effect of O₂ on the HCL oxidation reaction over RuO₂ (i.e., the so-called Sumitomo–Deacon process).^217,218 This reaction is a green chemistry route to recover high purity Cl₂ from HCL waste. Extensive KMC simulations revealed that, in contrast to previous suggestions, neither the adsorption of O₂ nor the associative desorption of chlorine are rate-determining during typical reaction conditions. Instead, hydrogen transfer in the water formation step is likely to determine the rate of the overall reaction. It was found that such hydrogen transfer processes are not highly activated, but they are strongly configuration controlled. Such configurations were found to be strongly determined by lateral interactions.

Recently, Chen et al.¹⁸⁵ implemented a protocol that combines the XPK method for accelerating KMC simulations (see Sec. IV D) with the CEH approach to investigate how the microscopic surface nonuniformity affect the syngas (CO + H₂) conversion on Rh(111) under operando conditions. Experimental observations of this reaction show that the selectivity toward acetaldehyde can increase with pressure, while the selectivity of methane exhibits the opposite trend (decreases with pressure).^219–222 Increasing the pressure increases the surface coverages and the impact of lateral interactions on the reaction mechanism. The operando theoretical analysis of Chen et al.¹⁸⁵ was found to agree with the experimental findings and demonstrated the power of 1p-CEH-KMC simulations in elucidating how the dynamic and intermediate-specific local coverage controls the selectivity. In a separate article, the same authors successfully implemented the methodology to elucidate the role of lateral interactions on the bistable region of the catalytic CO oxidation on platinum group metals.¹⁸⁶ 

Catalytic surface reactions typically occur under conditions in which lateral interactions affect activity and selectivity. Therefore, the accurate modeling of these interactions is of paramount importance in obtaining a detailed mechanistic understanding of the catalytic reaction mechanism. Such a relevance appears to be progressively more appreciated, as demonstrated by the fact that an increasing number of studies incorporate lateral interactions into their KMC simulations.

1p-KMC simulations are excellent for validating reaction mechanisms, identifying dominant pathways and rate-determining steps, and calculating species coverages to detect the most abundant reaction intermediates. These analyses are of paramount importance because surface reactions can be complex, and mechanisms with parallel or competing pathways are common. Moreover, understanding the surface reaction mechanisms provides insight into how changes in the operating conditions affect the overall reaction outcome and, hence, the catalyst’s performance. In this section, we discuss how 1p-KMC simulations have helped elucidate the complexity of chemical pathways in reactions of practical interest.

Let us start our discussion with the industrially important water–gas shift reaction (WGSR).^223–225 This reaction involves carbon monoxide and water vapor as reactants toward hydrogen and carbon dioxide and, thus, provides a route for the production of high-purity hydrogen. Hence, it is involved in several critical industrial and technological processes related, for instance, to the synthesis of ammonia and methanol. The interest in this reaction is also linked to the tight requirements of high purity hydrogen needed in fuel cells. The mechanism of the WGSR has been a subject of intense study during recent years; it is, however, still widely debated. Examples of reaction pathways include the redox or regenerative mechanism and the associative or carboxyl mechanism.^223,224 In the redox mechanism, the adsorbed water molecule is dissociated into H and O atoms; the resulting adsorbed oxygen atom then interacts with the adsorbed CO to form CO₂. In the associative mechanism, the adsorbed water molecule partially dissociates into OH and H adspecies. The resulting OH reacts with adsorbed CO to give a carboxyl (COOH) intermediate, which then decomposes to CO₂ and an H adatom.

Although several theoretical and experimental studies have been carried out to identify the most efficient WGSR catalysts, Cu-based catalysts are still the most widely used industrially. For this reason, a lot of effort has been devoted to understanding the mechanism of the WGSR catalyzed by Cu-based catalysts. In one of these studies, experimental evidence was presented that the activity of the CuO/ZnO/Al₂O₃ catalyst on the WGSR can be closely correlated with the Cu surface area.²²⁶ In such a catalyst, large Cu particles were present, predominantly exhibiting (111) facets. This finding motivated Prats et al.²²⁷ to perform a comprehensive 1p-KMC study of the WGSR on the Cu(111) surface. This study revealed that the reaction proceeds predominantly through the associative mechanism via a carboxyl intermediate and that the rate-limiting steps change at higher temperatures. Furthermore, a surface coverage analysis indicated that H₂O and H are the main adsorbed species at low temperatures, whereas OH and O are dominant at high temperatures. The simulations also revealed that reactant mixture compositions with high CO proportion enhance the production of H₂.

Recently, Chutia et al.²²⁸ performed a 1p-KMC study of the WGSR on the Pd(100) surface. Their focus on the Pd surface was motivated by the fact that Pd-based membranes are known to be able to isolate hydrogen in large quantities while maintaining stability during the WGSR. On the other hand, the choice of the Pd(100) surface was based on the fact that it is well characterized, relatively active, and sufficiently stable. 1p-KMC simulations allowed Chutia et al. to probe a reaction mechanism in which the redox and associative pathways operated simultaneously. These simulations indicated that H₂O and OH decomposition are the most common events. The split of H₂O is followed by the production of an H adatom and an OH species on the Pd(100) surface. Then, the OH molecule either reacts with the CO molecule to form carboxyl, which subsequently generates CO₂, or it may break down to oxygen, which then reacts with CO to form CO₂ via direct CO oxidation. Thus, the study indicated that the proposed redox-associative mechanism progress via both direct oxidation and carboxyl pathways that occur in parallel.

Decomposition of formic acid (HCOOH) to H₂ is another catalytic reaction relevant to hydrogen production, storage, and transportation, which has been recently explored with the 1p-KMC approach.^229,230 The selective release of H₂ from this reaction is a challenge because formic acid can either (i) dehydrate into CO + H₂O or (ii) dehydrogenate to CO₂ + H₂. In this regard, some recent experimental and theoretical studies have found that Au supported nanoclusters decompose formic acid with complete selectivity toward H₂.²³¹ However, the reasons behind this selectivity and activity remain unclear. Intending to clarify these Au nanoclusters properties, Chen et al.²²⁹ recently performed KMC simulations of formic acid decomposition over Au₁₈. The objective was to determine the nature of the active sites and the reaction mechanism. Their simulation results showed that indeed Au₁₈ is highly active and selective for formic acid decomposition, with triangular ensembles of atoms with a coordination number (CN) of five being the likely active sites (see Fig. 14). Interestingly, it was found that even though there are two of these active sites on Au₁₈, only one HCOOH molecule can be dehydrogenated at a time. This is because the strong stabilizing interactions between the adsorbates when they occupy both active sites lead to poisoning of the active sites by a pair of H and HCOO. Dissociation of additional HCOOH molecules is prohibited because of this transient poisoning of the cluster. These results constituted one of the first examples of heterogeneous catalysis by clusters one molecule at a time. The KMC simulations were also consistent with the selectivity of H₂ by the HCOO pathway, in agreement with aforementioned theoretical and experimental studies.

1p-KMC simulations have also delivered insights into the heterogeneously catalyzed CO₂ hydrogenation reaction.^214,232,233 This reaction has gained attention as an approach to mitigate the greenhouse effect of CO₂ and as an economical source of single carbon (C₁) products (e.g., CO, methane, methanol, and formic acid).²³⁴ By choosing the metal and type of support, it is possible to direct this reaction toward a particular product. For example, CO is mainly obtained through the reverse WGSR (RWGSR) using Rh-, Pt-, or Ni-based catalysts, methane can be produced using Ru-, Ni-, Pd-, or Pt-based catalysts via the so-called Sabatier reaction, and methanol is generated via CO₂ partial reduction using Pd-, Au-, or Cu-based catalysts.²³² Among all these catalysts, Ni is one of the most commonly used because of its low price compared to noble metals and its relatively high activity. Much investigation has been directed toward the hydrogenation of CO₂ on the ideal Ni(111) surface because it constitutes the most stable extended Ni surface. However, despite many efforts, the overall mechanism of this reaction is not fully understood, with opposing theoretical and experimental results regarding selectivity toward the RWGSR or the Sabatier reaction.^232,235

Thus motivated, Lozano-Reis et al.²³² recently carried out 1p-KMC simulations to investigate the molecular mechanism of CO₂ hydrogenation on a Ni(111) surface. This investigation aimed to determine the leading products and give insights regarding the dominant pathways governing the reaction under relevant operating conditions. The study considered several mechanisms for both the RWGSR and the Sabatier pathways. In contrast to the suggestion of early DFT calculations,²³⁶ the KMC simulations showed no methane formed on Ni(111) for any operating conditions investigated. The same simulations also showed that the RWGSR dominates mainly through the redox mechanism but also through the carboxyl mechanism to a lesser extent. Furthermore, the simulations identified the CO₂ dissociation step as the only rate-determining step. These results led Lozano-Reis et al. to speculate that the methane production typically observed experimentally on Ni-based catalysts is not due to the presence of Ni(111) facets of the Ni nanoparticles but the result of other contributions, for instance, the interplay between the nanoparticle and the support or the presence of other active sites.

The mechanism behind the CO₂ reduction to methanol over Cu-based catalysts has also been investigated by using the 1p-KMC approach.²³³ This work focused on pure Cu(111) catalysts, and KMC simulations were performed at various pressures and temperatures to study the selectivity, conversion, and TOF dependence at multiple conditions. The simulations demonstrated that methanol production is favored at low temperatures and high pressures, selectivity is highly dependent on pressure, and conversion and TOF are low. These results showed qualitative trends as obtained from experiments.

Similarly, the possibility of alcohol synthesis from syngas on Cu(111) supported defect-rich molybdenum disulfide (MoS₂) has also been studied by the 1p-KMC framework.²³⁷ The purpose of this work was to understand the role of metal support, which affects the electronic structure and geometry of defect-rich MoS₂. The combination of first-principles thermodynamics with 1p-KMC simulations showed that the support boosts the reactivity and product selectivity of defect-rich MoS₂, making it promising for ethanol synthesis. In particular, while thermodynamics appeared to favor reaction pathways whereby the Cu(111) support promotes both methanol and ethanol production, KMC simulations actually suggested a high selectivity toward the formation of ethanol.

The catalytic conversion of alcohols, an essential step in the valorization of biomass,^238–240 has also been investigated through 1p-KMC simulations. For instance, Réocreux et al.²⁴¹ recently performed KMC simulations of temperature programmed oxidation (TPO) spectra to investigate the role of oxygenated species in the mechanism of methanol oxidative coupling toward methyl formate (HCOOCH₃) on O pre-covered Au(111), with CO₂ being a key by-product due to overoxidation. In this study, a detailed comparison of the simulations with experimental TPO spectra enabled the validation of the proposed mechanism and the identification of rate-determining steps. Moreover, the simulations reproduced well the desorption temperatures of CO₂ and HCOOCH₃ and demonstrated the importance of considering van der Waals forces and adsorbate–adsorbate lateral interactions for the accurate modeling of the system.

The conversion of alcohols over oxide catalysts has also been the subject of attention. For instance, Sutton et al.²⁴² predicted below-room-temperature release of formaldehyde from methanol over CeO₂. KMC simulations showed that C–H bond breaking occurs via the disproportionation of adjacent methoxy species at such low temperatures. This finding is important because the conversion of alcohols by C–H activation is typically difficult to achieve at low temperatures, and CeO₂ is an inexpensive and abundant natural oxide. The mechanism behind the catalytic conversion of ethanol over the reducible La_0.7Sr._0.3MnO_3−x(100) surface to acetaldehyde and ethene was also recently investigated by 1p-KMC simulations in conjunction with pre-exposure temperature-programming reaction (PE-TPR) experiments.²⁴³ A branched mechanism was revealed by which ethene is produced by a β-dehydrogenation reaction and acetaldehyde is produced by a previously unknown disproportionation reaction.

Due to the increasing demand for light alkenes, catalytic dehydrogenation of light alkanes has received a lot of focus in recent years.^244–248 It constitutes an environmentally friendly alternative route to the energetically costly and unselective production of alkenes by the traditional cracking of oil products. The pertinent mechanisms can be classified as direct (non-oxidative) or oxidative, depending on whether an oxidant is used or not.^244,249 Direct dehydrogenation (DH) is currently used in industry, and conventional catalysts are platinum and chromium oxide. A problem with DH is the deactivation of the catalysts by the coke formed during the reaction.^244,250 On the other hand, the oxidant used in oxidative dehydrogenation can help prevent coke formation.²⁴⁴ However, oxidative dehydrogenation is not commonly implemented due to its low selectivity toward desired alkenes. In this context, 1p-KMC simulations were recently implemented to investigate the oxidative C–C and C–H bond cleavage of ethane with CO₂ as soft oxidant on a bimetallic PtNi(111) model surface.²⁵¹ The simulations were used to explore the selectivity of PtNi(111) toward syngas (CO + H₂) and ethylene (CH₂CH₂) under typical experimental reaction conditions. It was shown that oxidative ethane dehydrogenation to ethylene primarily occurs by two successive C–H bond scissions. The propane DH on Pt(111) was also recently studied by Lian et al.²⁵² In this work, KMC simulations were used to reveal the dominant reaction pathway to propylene formation and the origin of coke formation. It was found that the availability of active sites crucially affects both propylene and coke formation and that the quick deactivation of the catalysts occurs because most of the active sites are occupied by cracking products that are difficult to remove from the surface.

The DH of propane and butane on Cr₂O₃(0001) has also been investigated by 1p-KMC simulations.^253–255 Chromiun oxide catalysts are used in the so-called CATOFIN process for the production of olefins, such as propylene (from propane) and iso-butylene (from iso-butane).²⁵⁵ Such studies have delivered insights into the most abundant products of the reaction pathways investigated, unwanted pathways leading to side-products or poisons, and the nature of catalyst deactivation by coking. In the case of DH of propane on Cr₂O₃(0001),²⁵⁴ it was found that the accumulation of propylene and propyne in the reaction mixture adversely affects the reaction rate and selectivity. It was also reported that higher pressures increase the reaction rate and coke formation rate. The simulations also revealed that the deactivation of the catalyst has a strong temperature dependence and is caused by the accumulation of coke and coke-like intermediates. On the other hand, 1p-KMC simulations of butane DH on Cr₂O₃(0001)²⁵³ demonstrated that 2-butene is the most abundant product of dehydrogenation. KMC simulations were also used to investigate the rate at which the modeled catalyst gets deactivated. In particular, it was shown that coking is negligible at low temperatures and becomes relevant at higher temperatures and that deactivation of the catalyst was caused by the formation of coke deposits.

Finally, let us briefly comment on two new “frontier topics” in catalysis research to which the 1p-KMC framework has also been successfully applied. The first is single-atom catalysts (SACs), whereby one reduces the active site of the heterogeneous catalyst to a single isolated precious metal, typically coordinated to an oxide support.^{189,256–261} In this respect, Alexopoulos et al.²⁶¹ recently performed 1p-KMC simulations of CO oxidation over Pd atoms on α-alumina as a test case to provide insights into this area of research. The results of this investigation demonstrated that KMC analysis can help discriminate between different mechanisms as well as between different active sites. It was also shown that KMC can be a complementary tool to current spectroscopic methods typically implemented to investigate the active site(s) on reactions over SACs. In another interesting study, Su et al.²⁶⁰ implemented 1p-KMC simulations to show that Pt atoms doped into a CeO₂ surface exhibit a very high CO oxidation activity and thermodynamic stability in comparison to models involving Pt single atoms on terrace and steps of CeO₂.

The second “frontier topic” is single-atom-alloy (SAA) catalysts. SAAs are another type of SACs in which catalytically active components, such as Pt, Ni, Pd, Ru, and Rh, are atomically dispersed in more inert, but more selective, host metals, such as Cu, Au, and Ag.^{7,256,262–271} In this regard, 1p-KMC simulations in conjunction with surface science experiments have shed light into the mechanisms behind the efficient and selective C–H activation on coke-resistant PtCu SAAs²⁶⁴ and C–C coupling on coke-resistant PdAu and NiAu SAAs.^265,267 A similar approach was recently implemented to explore how to control hydrocarbon (de)hydrogenation pathways with bi-functional PtCu SAAs.²⁶⁶ 

The works described in this subsection are just a representative sample of a larger number of studies that have successfully implemented the 1p-KMC approach to validate and elucidate complex chemical pathways in heterogeneous catalysis. Let us now move onto Sec. V C, where we will discuss the first-principles multiscale modeling of structure-sensitive catalytic reactions.

It is widely recognized that the structural complexity of real catalysts can often affect the catalytic performance. Real catalysts are highly dynamic materials that adapt their morphology (geometric structure) to the constantly changing chemical environment and operating conditions. In doing so, they may exhibit a variable number of active sites (i.e., low-coordination sites and defects) at which the elementary steps of the reaction considered might proceed. Such spatial and/or temporal variability may affect the active pathways of the so-called “structure-sensitive” reactions, whose investigation is a persistence challenge.²⁷² In this section, we briefly comment on recent 1p-KMC investigations devoted to obtaining a fundamental understanding of the role of structural sensitivity on heterogeneous catalytic systems.

The structure sensitivity in catalytic activity and selectivity is typically approached from a reductionist perspective in which the catalyst is divided into isolated facets. Such an approach has been recently implemented, in conjunction with 1p-KMC simulations, to investigate several relevant structure-sensitive catalytic reactions. For instance, Prats et al.²⁷³ compared the catalytic activity of the WGSR on the stepped Cu(321) surfaces with the one on a flat Cu(111) surface. Their simulations revealed that in contrast to the prevalence of the associative mechanism for Cu(111), both the redox and associative mechanism are possible for Cu(321). It was also found that despite exhibiting lower activation energies, stepped surfaces (i.e., low coordinated sites) do not necessarily have an overall high catalytic activity. This unexpected observation was rationalized based on a detailed investigation of the coverage effects due to lateral interactions and the relative contribution of several elementary steps to the overall TOF.

The effect of the structure of Cu catalyst surfaces on the catalytic CO₂ hydrogenation to methanol,²⁷⁴ under experimental conditions, was also recently investigated using 1p-KMC simulations by Kopač et al.²⁷⁵ Flat Cu(111) and stepped Cu(533) surfaces were considered in this work. As expected, the stepped Cu(533) surface enhances the activity and selectivity toward methanol compared to the Cu(111) surface. Surface coverage studies revealed information about the most abundant intermediates. At steady-state, HCOOH is the most abundant on Cu(111), while hydrogen is the most abundant on Cu(533). Furthermore, this study also gave a clear perspective on which reactions on the methanol synthesis pathway are favored on both surfaces. In this respect, the event frequency per active site indicated that Cu(533) enables the otherwise suppressed H₂COH hydrogenation, resulting in a higher CH₃OH yield.

The industrially important ethylene epoxidation reaction on silver catalysts has also been investigated by extensive 1p-KMC simulations by Huš et al.^276,277 The reaction was modeled on three pristine silver surfaces: Ag(100), Ag(110), and Ag(111), as well as on the missing-row reconstructed Ag(110) surface. One of the objectives was to understand how oxygen coverage affects selectivity and activity. The simulation results revealed that Ag(111) maintains very low oxygen coverage while being the least active surface with a moderate selectivity, Ag(100) exhibits the highest selectivity at high oxygen coverage, and both pristine and reconstructed Ag(110) surfaces lack any appreciable selectivity but are the most active. These observations were compared with previous experimental reports and were found to be in good agreement with the data for Ag(111) and Ag(100). The agreement between predicted results and experimental data for Ag(110) was less satisfactory due to the fact that the model did not take into account reconstruction as a dynamic phenomenon. It turns out that, on this surface, reconstruction can be brought about by the adsorption of oxygen.

The reductionist approach has also been used to investigate the structure sensitivity of several prototypical catalytic reactions.^278–281 For instance, Fajín et al.²⁷⁸ conducted 1p-KMC simulations to clarify the role of the silver facets in the catalytic CO oxidation on nanoporous gold (NPG) catalysts obtained by dealloying an AuAg alloy. This catalytic system is important for cleaning hydrogen before feeding it into fuel cells. Fajín et al.,²⁷⁸ thus, performed computer simulations for the reaction on both Au(110) and Ag(110) surfaces, and based on their results, it was proposed that the small silver microfacets can be responsible for the CO activation on such NPG catalysts.

The reductionist approach of dividing the catalysts into separate components, which are studied individually, delivers significant insights into the role of the structural complexity of real catalysts but certainly has limitations. To properly unveil synergistic effects due to such complexity, a holistic approach is needed. In this respect, the 1p-KMC framework is also proving to be a valuable computational tool, for instance, the work of Guo and Vlachos is an interesting example of such a holistic 1p-KMC simulation approach.²⁸² This work entailed 1p-KMC simulations of ammonia decomposition on patched bimetallic Ni/Pt surfaces, a prototypical system for structure sensitivity. By varying the size and/or shape of Ni clusters on Pt, the bifunctional behavior of such patched bimetallic surfaces was elucidated. Among other things, it was shown that the Ni terrace sites catalyze N–H bond scission and the (110) edges of the Ni patches catalyze N₂ association. Furthermore, the computational analysis revealed that such a dual-site behavior is responsible for the higher activity of patched bimetallic surfaces compared to full Ni monolayers on Pt or the pure metal surfaces of Ni and Pt. Interestingly, it was also found that the structure sensitivity of the reaction was rather weak on these patched surface bimetallics under the chosen reaction conditions. In an interesting follow-up work, Núñez and Vlachos²⁸³ combined KMC simulations with active learning to optimize the Ni/Pt catalysts surface microstructure to enhance reaction rates of the ammonia decomposition reaction.

1p-KMC simulations have also been recently implemented to understand the mechanism behind the enhancement of the catalytic activity and/or selectivity due to the so-called strong metal support interactions (SMSIs).²⁸⁴ In several cases, SMSIs lead to catalyst deactivation. However, it has been found experimentally that SMSIs make Au nanoparticles dispersed on molybdenum carbide (MoC), a highly active catalyst for the low-temperature WGSR.²⁸⁵ In this respect, KMC simulations have unraveled the origin of the experimentally observed high activity and have provided strong evidence for a cooperative effect between the different regions of the catalysts. In particular, it was found that the clean MoC regions are responsible for adsorbing and dissociating water molecules and that the interface regions between the nanoparticles and support act as attractors for CO molecules. The latter subsequently react with OH molecules produced in the clean MoC region to form COOH, which then produces CO₂ via the associative (carboxyl) mechanism.

The structure sensitivity of the RWGSR in SrTiO₃-based perovskite-supported copper catalysts was also recently investigated by Kopač et al.²⁸⁶ Their KMC simulations revealed that considering copper and the support sites in addition to the interface sites results in higher predicted rates for the reaction compared to the case in which the interface alone is modeled. In a separate work,²⁸⁷ the same authors explored the synergistic effect of bifunctional Cu/perovskite catalysts on the catalytic hydrogenation of CO₂ to methanol. In this investigation, KMC simulations were performed on a Cu phase with four perovskite substrate materials (i.e., Cu/CaTiO₃, Cu/SrTiO₃, Cu/BaTiO₃, and Cu/PdTiO₃). It was found that all systems outperformed the pure Cu, with Cu/PbTiO₃ and Cu/SrTi₃ being the most promising copper/perovskite catalysts.

Nanometer-sized particles dispersed on oxide supports are generally ill-defined with respect to size and shape and expose a range of different interconnected active sites from which complex kinetic behavior could arise. 1p-KMC simulations demonstrating genuine nanoscale effects on catalytic activity have also been performed. For instance, Nikbin et al.²⁸⁸ carried out a 1p-KMC analysis to elucidate the experimentally observed “magic number” behavior of sub-nanoscale Au clusters toward CO oxidation. Traditionally, low coordinated sites are thought to be highly active, and thus, it would be expected that all types of Au nanoparticles would be effective catalysts. However, experimental observation on small Au_n clusters (n = 2–20 atoms) on a O-defective MgO support revealed that the most active clusters were Au₈, Au₁₈, and Au₂₀.^289,290 Thus motivated, Nikbin et al. investigated the catalytic behavior of $Aun−$ clusters with n = 6, 8, and 10 atoms. They used negatively charged clusters because Au clusters up to the Au₁₃ size, when supported on MgO with O vacancies, are known to be negatively charged. Apart from unraveling a high degree of complexity in the catalytic behavior of these Au clusters, the simulation results were in good agreement with the experimental observations that $Au6−$ is inert, $Au8−$ is active, and $Au10−$ is less active than $Au8−$ (Fig. 15). More specifically, the KMC simulations predicted that $Au6−$ gets poisoned by carbonate, while $Au8−$ and $Au10−$ exhibit sustained activity, via pathways involving CO–O₂ intermediates.

To further enrich our understanding of structure-sensitive catalytic reactions over nanoparticles, it is of paramount importance to “map out” complex reaction energy landscapes by considering all catalytically relevant pathways on a potentially large number of inequivalent active sites. A full treatment of such landscapes is as of yet impractical or even infeasible. Instead, the problem has recently been tackled by constructing KMC models incorporating DFT-based structure sensitivity scaling relations.^120,291 In such scaling relations, generalized coordination numbers are used as descriptors for adsorption energies and reaction barriers. In doing so, one can efficiently and accurately address the structural complexity of nanoparticles and the synergistic effects emerging from assemblies of active sites. This first-principles scaling-relation KMC (1p-SR-KMC) approach has been successfully applied to model structure sensitivity features of the archetypal CO oxidation reaction on nanoparticles.^{120,292–296} It has also been implemented to study more complex surface reaction systems, such as the selective acetylene hydrogenation over SAA nanoparticles²⁹⁷ and oxygen reduction on oxide-supported PtNi nanoalloys.²⁹⁸ 

Motivated by the need to understand the performance of subnanometer catalysts and explain how catalyst treatment and exposure to spectroscopic probe molecules change the structure, Wang et al.¹⁴³ developed a 1p-KMC framework that incorporates machine learning-based Hamiltonians. This framework was used to follow the evolution of subnanometer clusters at experimentally relevant time scales. Wang et al. choose Pd (n = 1–40) on CeO₂(111) in a CO atmosphere as a case study, since CO is the most common probe molecule in infrared (IR) spectroscopy.²⁹⁹ The approach gave important insights into the effect of temperature, CO partial pressure, metal loading, and initial catalyst state on cluster formation at that scale.

A catalytic phenomenon also addressed to some extent by the KMC framework is the so-called surface morphological rearrangement (or surface reconstruction). Such a dynamical behavior can have a significant impact on the observed catalytic activity and selectivity. Thus, recently, Hoffmann et al.³⁰⁰ introduced a multi-lattice 1p-KMC approach to describe the dynamics of morphological transitions on solid surfaces in the special case whereby the reconstructed and the pristine structures can be captured by commensurate lattices. As a case study, they modeled the reduction of surface oxide on Pd(100). The simulation results reproduced the observed experimental trends in the reduction rates and revealed the crucial role of elementary steps at the boundary between oxide and metal domains.

In a recent work, Huš et al.³⁰¹ investigated how Cr₂O₃ catalyzes the propane dehydrogenation reaction in oxidative and reducing environments. An interesting contribution of this work is that the catalytic system was studied on a mixed surface, consisting of equal parts of oxidized and reduced areas. It was shown that this mixed surface exhibits considerably better performance than each individual surface. Furthermore, investigation of the surface with varying degrees of oxidation showed that there exists an optimal degree of surface oxidation with respect to propene yield.

In this section, we have discussed several representative studies showcasing the versatility of the 1p-KMC approach in exploring structural changes of the catalyst surface and structure sensitivity effects. Even though much progress is being made in this area, challenges remain, which motivate exciting developments. Thus, the computational expense and the amount of effort required toward developing comprehensive models of complex chemistries (involving numerous pathways) on realistic catalysts (exposing various types of active sites) remain prohibitive. In addition, a general framework for modeling complex catalytic reconstruction effects is still lacking. Descriptor-based approaches facilitated by the implementation of machine learning could be critical in overcoming these challenges.

KMC simulations in combination with first-principles-based calculations are becoming essential in modeling heterogeneous catalysis. As highlighted in this Perspective, such a computational approach enables us to explore the wide range of length scales and time scales over which structure–function relationships unfold. Here, we briefly introduced this versatile modeling framework, discussed the main outstanding computational challenges, and commented on successful applications. The latter clearly demonstrate the power of 1p-KMC in delivering insights into reaction mechanisms and rate-determining steps, the role and interplay of different active sites, as well as the impact of operating conditions and the reaction micro-environment on catalytic performance metrics (observables).

1p-KMC is a bottom-up simulation framework, and as such, its capacity to generate reliable results depends on the accuracy and efficiency of the first-principles methods employed to obtain the rate constants fed into it. At present, the first-principles method widely used for this purpose is DFT. However, using DFT in building detailed KMC models can become computationally expensive because the number of rate constants to be obtained can be very large, particularly when lateral interactions are taken into account. To this end, a common approach is to implement CEH-based BEP relations to calculate activation energies and hTST to get prefactors.⁸⁷ The quest for efficient and accurate ways to obtain rate constants from first-principles is thus of paramount importance to eventually reach the full potential of this powerful computational tool.

Another challenge in obtaining rate constants via DFT calculations is that the errors arising from such calculations may lead to rate constants that are potentially inaccurate by several orders of magnitude. This issue underlines the need for more reliable DFT (or, in general, first principles) methods but also efficient approaches for estimating the sensitivity of the model predictions on the rate constants and quantifying the uncertainties on these predictions. In this Perspective, we have discussed recent progress in developing such methods for KMC simulations. The ultimate goal of such SA and UQ approaches is to direct computational efforts into improving the truly relevant first-principles calculations.

A major problem in modeling heterogeneous catalysis by 1p-KMC is that elementary events often happen at vastly different time scales. The development of algorithms to overcome this issue is an active area of research. While progress has been (and is being) made, the current algorithms are still not robust or reliable enough for “out-of-the-box” usage. Having generic and easy-to-use algorithms to tackle this issue is paramount for performing direct comparisons of KMC predictions with experimental data and for expanding the applicability of 1p-KMC simulations to more complex catalytic systems.

The poor scalability of KMC with respect to the lattice size has prevented its implementation in simulating large catalytic surface domains, relevant to catalytic reaction systems that exhibit pattern formation. Previously proposed approximate methods have seen limited adoption, potentially due to the fact that the errors incurred by these approximations have not been thoroughly studied and quantified. Ideally, for simulations of spatiotemporal pattern formation, exact KMC schemes would be preferable, ensuring simulations free of artifacts due to numerical error. Latest research efforts have started to address this challenge by implementing distributed parallelization techniques that properly handle causality errors arising from boundary conflicts due to the domain decomposition of the lattice.⁸² The development and application of such distributed KMC approaches is quite an exciting research area with lots of scope for further efforts in algorithm development. Equally important is developing and improving algorithms for event search and execution, as well as KMC state update, in order to evolve the state-to-state dynamics in KMC simulations more efficiently.^109,132

The 1p-KMC framework discussed throughout this Perspective is based on a rigid lattice representation of the catalytic surface. This representation poses challenges when modeling situations in which dynamical reconstruction or other morphological changes of the catalyst take place. Although some authors have addressed the kinetics of such transformations within the rigid on-lattice KMC approach,³⁰⁰ there is a growing demand for new tools and methodological developments to investigate the effects of such dynamical changes on catalytic kinetics. The so-called off-lattice or adaptive KMC approach could, in principle, be useful for this purpose.¹⁷⁰ However, its implementation is currently computationally too expensive. A boost to solve this challenge could also be provided by integrating machine learning approaches into the 1p-KMC-based multiscale modeling framework.^8,38,302

The bottom-up modeling approach for heterogeneous catalysis aims to bridge the gap between the supported metal nanoparticle, the catalytic pellet, and the reactor. From a computational standpoint, the goal is thus to reach a direct coupling between CFD models and 1p-KMC simulations. This topic is a rather new area of research full of many exciting challenges.

In the past decade, the KMC approach has steadily matured in the computational catalysis field. An indicator of this maturity is the emergence of several user-friendly software applications/packages devoted to 1p-KMC simulations for heterogeneous catalysis. These software packages have contributed to the impressive growth in the number of studies employing this computational framework. This growth will continue, driven by the need for kinetic simulations in catalysis, enabling detailed comparisons of theoretical predictions with experimental data. Moreover, this growth will potentially accelerate, facilitated by advances in algorithm development for the efficient computational implementation of the KMC framework itself, but also advances in the generation, processing, sharing, and re-use of large amounts of data on catalyst–adsorbate interactions and reaction mechanisms. This impetus is expected to spread to other branches of heterogeneous catalysis not covered in this Perspective, for instance, electrocatalysis,^303–308 and photocatalysis.^309–311 The KMC approach will thus (continue to) drive exciting breakthroughs in the quest toward understanding and predicting catalytic materials for niche applications.

The authors acknowledge funding from the European Union’s Horizon 2020 Research and Innovation Program under Grant Agreement No. 814416.

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.