Kinetic Monte Carlo (KMC) simulations in combination with first-principles (1p)-based calculations are rapidly becoming the gold-standard computational framework for bridging the gap between the wide range of length scales and time scales over which heterogeneous catalysis unfolds. 1p-KMC simulations provide accurate insights into reactions over surfaces, a vital step toward the rational design of novel catalysts. In this Perspective, we briefly outline basic principles, computational challenges, successful applications, as well as future directions and opportunities of this promising and ever more popular kinetic modeling approach.

## I. INTRODUCTION

Heterogeneous catalysts are capable of accelerating chemical reactions, thereby achieving high conversion and product selectivity at low cost.^{1} 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.^{8} 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).^{9} 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.^{33} On the other hand, KMC simulations are very accurate but computationally demanding.^{34} 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.

## II. THE KMC FRAMEWORK: WHY IS IT USEFUL AND WHY DOES IT WORK?

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.^{43} 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.^{44} 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.^{53} 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)^{69} 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.

### A. The temporal coarse-graining underlying KMC simulations

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.^{77} This hopping probability is independent of the previous history of state *ω* and is calculated via TST, taking into account the properties of the PES.^{78}

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.^{79}

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}

*ω*to a new possible state. Thus, under Poisson statistics, the probability that the system has escaped from state

*ω*at some time less than or equal to

*τ*(i.e., the cumulative probability of the time to escape from state

*ω*to any other state) is given by

*ω*due to all possible transitions (elementary events). Then, the probability density function that a transition will occur at time

*t*+

*τ*, given that the system is at state

*ω*at time

*t*, is given by the following exponential distribution:

*ω*to another new state

*ω*′ are uncorrelated, the escape time of each separately follows an exponential distribution, given as

^{76,77}one obtains

*u*∈ (0, 1) is a random number from the uniform distribution and

*τ*

_{ωω′}is the escape time sought. In KMC parlance, these random times are also referred to as occurrence, waiting, or inter-arrival times. Although at a given state each possible elementary event has its own individual escape time, the first event to happen is, of course, the one with the shortest time.

^{55,76,77}

*u*

_{1}∈ (0, 1) is also a random number from a uniform distribution. In this case, however, the escape time only depends on the total rate constant and is independent of the elementary event that brings the system out of state

*ω*. Such an elementary event, which propagates the system to a new state

*ω*′, is picked with probability

*u*

_{2}into an integer that indicates the next event.

^{54,55,76,77,80}

^{55,81}In this case, the escape time for each elementary event is obtained by solving the following non-linear equation:

*k*

_{ωω′}is now changing over time. The time advancement and the next elementary event to occur are then given by the minimum escape time among all events. Alternatively, one can obtain the time advancement from

*k*

_{tot}is also changing over time. Then, the integer random variable indicating the next event to happen would follow a probability distribution similar to that of Eq. (7) but with the rate constants integrated over time. If the rate constants do not change over time, Eqs. (8) and (9) reduce to Eqs. (5) and (6), respectively.

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.

### B. Detailed balance, rate constants, and lattice energetics

^{55,84}Given that these elementary events are treated as independent Poisson processes, the evolution of the catalytic system is described by the so-called Markovian master equation (MME),

^{75}

*k*

_{ωω′}represents, as mentioned above, the rate constant for a transition from state

*ω*to state

*ω*′ due to an elementary event and

*k*

_{ω′ω}represents the rate constant of the reverse transition. The sums in the right-hand side of this equation are taken over all possible states of the state space Ω, keeping always in mind that if there is no elementary event connecting two states, the corresponding rate constant is zero.

^{55}

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.

*st*indicates stationarity. This equality expresses that, under stationary conditions, the total probability flux leaving a state must be balanced by the total probability flux toward that state. Moreover, if this stationary solution corresponds to a state of thermodynamic equilibrium, microscopic reversibility imposes the even stronger constraint of detailed balance. The detailed balance condition demands that each microscopic process has a corresponding reverse process and that the average rate of every elementary event is equal to the average rate of its reverse event. This, together with the fact that the system at thermodynamic equilibrium in a particular state follows a Boltzmann distribution, allows us to conclude that

*k*

_{B}is the Boltzmann factor,

*T*is the temperature, and

*ɛ*

_{ω}and

*ɛ*

_{ω′}are the free energies associated with states (or PES basins)

*ω*and

*ω*′, respectively.

^{24,28}Because the free energy contains contributions for the electronic and vibrational energies of the lattice species and the rotational, vibrational, translational, and electronic energies from the gas species, Eq. (12) can be rewritten as

*Q*

_{ω}and

*Q*

_{ω′}are quasi-partition functions for the initial and final states, respectively.

^{24}Each of these partition functions incorporates the vibrational contributions of adsorbates, as well as the translational, rotational, and vibrational contributions of gas species. Regarding the vibrational contributions, we are assuming that the vibrational energy levels are with respect to the bottom of the potential energy well; if the ground state is chosen as a reference, zero-point energy corrections need to be accounted for in the activation and reaction energies (but not in the vibrational partition functions anymore, to avoid double-counting). Δ

*E*

_{rxn}is the change of electronic energy of the system due to a transition from

*ω*to

*ω*′ brought about by removing all the reactant molecules from the lattice and the gas phase (depending on the nature of the elementary event), as well as adding the product molecules onto the lattice and gas phase. In other words, Δ

*E*

_{rxn}, is the reaction energy,

*E*and

*E*

^{gas}denote the electronic energy of the given lattice configuration and gas species, respectively.

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.^{24} 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.

*ω*to the final state

*ω*′ is calculated from TST,

^{24,28,69}

*κ*is known as the dynamic transmission coefficient. It is a correction term, typically calculated from short MD simulations, that captures the possibility of trajectories recrossing the transition state region of the PES back to the initial state.

^{1}

^{,}

*h*is Planck’s constant,

*Q*

_{‡}is the quasi-partition function of the transition state, and $E\omega \omega \u2032\u2021$ is the activation energy barrier of the elementary event (i.e., the difference in the potential energy between the transition and initial states). The latter is typically computed from the PES using DFT calculations combined with efficient methods for locating the transition state. Examples of these methods are the dimer method and the nudged elastic band method.

^{43,96}In most studies, the partition functions are obtained from first-principles calculations within the harmonic TST (hTST) approximation.

^{1,30}However, recent efforts to account for anharmonic effects, either directly or via the hindered translator and rotor models, have been put forward in the literature.

^{53,97–100}Invoking the detailed balance condition together with the TST rate constants for the forward (

*ω*→

*ω*′) and reverse (

*ω*′ →

*ω*) elementary events gives

^{24}

^{,}

^{28,87}However, since there may be several different configurations of the spectators, it is impractical to calculate distinct DFT activation energy barriers for every possible configuration. Instead, a convenient way to model the impact of lateral interactions on the activation energy is by using the so-called Brønsted–Evans–Polanyi (BEP) relations, which represent a linear correlation between the activation barrier and the reaction energy of an elementary step in a reaction mechanism.

^{87,88}Such BEP relations are obtained by linear fitting to appropriate DFT calculations.

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.

## III. FROM KMC ALGORITHMS TO CODES

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.

### A. The general flow chart of KMC algorithms

KMC is a computational framework that generates stochastic trajectories fulfilling the fundamental premises of the MME (see Sec. II B).^{75} 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).^{24}

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).^{101} 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.*^{28} and Prats *et al.*^{27}

### B. Computational codes

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)^{106} and the catalytic CO oxidation on Pt surfaces.^{107}

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*,^{108} 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.^{114} 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.^{119} 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.^{120} 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.^{123} Finally, the software package SuSmoST (SUrface Science MOdeling and Simulation Toolkit)^{124} 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.

## IV. COMPUTATIONAL CHALLENGES AND METHODOLOGICAL DEVELOPMENTS

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.

### A. Scheduling and executing elementary events

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.*^{130} 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.^{131}

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.^{102} 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^{101} reviewed several KMC simulation methods and explored the implementation and efficiency of several search and update algorithms.^{128} 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].^{132} 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),^{135} a simplified model of the water–gas shift (WGS) reaction on Pt(111),^{136} and a TPD model of CO on Cu(111).^{137} 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.

### B. Treating complicated energetic models of non-ideal adlayers

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.^{28} 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.*^{87} developed an implementation of this idea, which uses open multi-processing (OpenMP)^{141} 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.*^{109} 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.^{95} 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,^{142} as opposed to the traditional linear search implementation of the DM,^{101} 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_{2}(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.

### C. Treating large surface domains with distributed simulations

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^{2}. 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.^{101} 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,^{148} the synchronous relaxation,^{149} the optimistic synchronous relaxation,^{150,151} the optimistic Time-Warp,^{159} 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),^{160} 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*.^{82} 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.^{161} 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).^{83} Crucially, Ravipati *et al.*^{82} 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,^{146} and lead to novel applications in the area of spatiotemporal control of reactive systems.^{162}

### D. Treating event frequency disparity

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),^{40} 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.^{173} 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.^{174} 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.*^{172} 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.*,^{172} 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.*^{13} 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.*^{113} 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.^{175} 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.^{176} 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.^{177} 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.^{190}

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.

### E. Steady-state detection, sensitivity analysis, and uncertainty quantification

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\omega (t1)$, $P\omega (t1),\omega (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.*^{192} 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.^{193} 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.^{194} 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.*^{195} 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.*^{196} 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.*^{197} 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.^{198}

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.

### F. Coupling with larger scales

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.^{204}

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.^{213}

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.

## V. CHEMICAL MECHANISMS AND PHYSICAL PHENOMENA UNRAVELED BY KMC SIMULATIONS

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.

### A. Adlayer non-ideality and effects of lateral interactions

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.^{216} 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_{2}, another model reaction and a key step in NO_{x} elimination from exhaust gases, has also been recently explored by the CEH approach.^{93} 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_{2} on the HCL oxidation reaction over RuO_{2} (i.e., the so-called Sumitomo–Deacon process).^{217,218} This reaction is a green chemistry route to recover high purity Cl_{2} from HCL waste. Extensive KMC simulations revealed that, in contrast to previous suggestions, neither the adsorption of O_{2} 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.*^{185} 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_{2}) 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.*^{185} 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.^{186}

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.

### B. The complexity of chemical pathways

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_{2}. 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_{2} 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_{2}O_{3} catalyst on the WGSR can be closely correlated with the Cu surface area.^{226} In such a catalyst, large Cu particles were present, predominantly exhibiting (111) facets. This finding motivated Prats *et al.*^{227} 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_{2}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_{2}.

Recently, Chutia *et al.*^{228} 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_{2}O and OH decomposition are the most common events. The split of H_{2}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_{2}, or it may break down to oxygen, which then reacts with CO to form CO_{2} 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_{2} 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_{2} from this reaction is a challenge because formic acid can either (i) dehydrate into CO + H_{2}O or (ii) dehydrogenate to CO_{2} + H_{2}. In this regard, some recent experimental and theoretical studies have found that Au supported nanoclusters decompose formic acid with complete selectivity toward H_{2}.^{231} However, the reasons behind this selectivity and activity remain unclear. Intending to clarify these Au nanoclusters properties, Chen *et al.*^{229} recently performed KMC simulations of formic acid decomposition over Au_{18}. The objective was to determine the nature of the active sites and the reaction mechanism. Their simulation results showed that indeed Au_{18} 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_{18}, 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_{2} by the HCOO pathway, in agreement with aforementioned theoretical and experimental studies.

1p-KMC simulations have also delivered insights into the heterogeneously catalyzed CO_{2} hydrogenation reaction.^{214,232,233} This reaction has gained attention as an approach to mitigate the greenhouse effect of CO_{2} and as an economical source of single carbon (C_{1}) products (e.g., CO, methane, methanol, and formic acid).^{234} 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_{2} partial reduction using Pd-, Au-, or Cu-based catalysts.^{232} 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_{2} 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.*^{232} recently carried out 1p-KMC simulations to investigate the molecular mechanism of CO_{2} 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,^{236} 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_{2} 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_{2} reduction to methanol over Cu-based catalysts has also been investigated by using the 1p-KMC approach.^{233} 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_{2}) has also been studied by the 1p-KMC framework.^{237} The purpose of this work was to understand the role of metal support, which affects the electronic structure and geometry of defect-rich MoS_{2}. The combination of first-principles thermodynamics with 1p-KMC simulations showed that the support boosts the reactivity and product selectivity of defect-rich MoS_{2}, 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.*^{241} 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_{3}) on O pre-covered Au(111), with CO_{2} 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_{2} and HCOOCH_{3} 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.*^{242} predicted below-room-temperature release of formaldehyde from methanol over CeO_{2}. 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_{2} is an inexpensive and abundant natural oxide. The mechanism behind the catalytic conversion of ethanol over the reducible La_{0.7}Sr._{0.3}MnO_{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.^{243} 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.^{244} 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_{2} as soft oxidant on a bimetallic PtNi(111) model surface.^{251} The simulations were used to explore the selectivity of PtNi(111) toward syngas (CO + H_{2}) and ethylene (CH_{2}CH_{2}) 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.*^{252} 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_{2}O_{3}(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).^{255} 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_{2}O_{3}(0001),^{254} 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_{2}O_{3}(0001)^{253} 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.*^{261} 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.*^{260} implemented 1p-KMC simulations to show that Pt atoms doped into a CeO_{2} 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_{2}.

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^{264} 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.^{266}

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.

### C. Structural changes of the catalyst surface and structure sensitivity effects

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.^{272} 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.*^{273} 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_{2} hydrogenation to methanol,^{274} under experimental conditions, was also recently investigated using 1p-KMC simulations by Kopač *et al.*^{275} 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_{2}COH hydrogenation, resulting in a higher CH_{3}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.*^{278} 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.*,^{278} 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.^{282} 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_{2} 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^{283} 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).^{284} 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.^{285} 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_{2} via the associative (carboxyl) mechanism.

The structure sensitivity of the RWGSR in SrTiO_{3}-based perovskite-supported copper catalysts was also recently investigated by Kopač *et al.*^{286} 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,^{287} the same authors explored the synergistic effect of bifunctional Cu/perovskite catalysts on the catalytic hydrogenation of CO_{2} to methanol. In this investigation, KMC simulations were performed on a Cu phase with four perovskite substrate materials (i.e., Cu/CaTiO_{3}, Cu/SrTiO_{3}, Cu/BaTiO_{3}, and Cu/PdTiO_{3}). It was found that all systems outperformed the pure Cu, with Cu/PbTiO_{3} and Cu/SrTi_{3} 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.*^{288} 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_{8}, Au_{18}, and Au_{20}.^{289,290} Thus motivated, Nikbin *et al.* investigated the catalytic behavior of $Aun\u2212$ clusters with *n* = 6, 8, and 10 atoms. They used negatively charged clusters because Au clusters up to the Au_{13} 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\u2212$ is inert, $Au8\u2212$ is active, and $Au10\u2212$ is less active than $Au8\u2212$ (Fig. 15). More specifically, the KMC simulations predicted that $Au6\u2212$ gets poisoned by carbonate, while $Au8\u2212$ and $Au10\u2212$ exhibit sustained activity, via pathways involving CO–O_{2} 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^{297} and oxygen reduction on oxide-supported PtNi nanoalloys.^{298}

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.*^{143} 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_{2}(111) in a CO atmosphere as a case study, since CO is the most common probe molecule in infrared (IR) spectroscopy.^{299} 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.*^{300} 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.*^{301} investigated how Cr_{2}O_{3} 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.

## VI. FUTURE 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.^{87} 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.^{82} 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,^{300} 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.^{170} 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.

## VII. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

## REFERENCES

*Fundamental Concepts in Heterogeneous Catalysis*

*Physics of Surface, Interface and Cluster Catalysis*

*Introduction to Computational Chemistry*

*Theoretical Surface Science: A Microscopic Perspective*

*Modeling and Simulation of Heterogeneous Catalytic Reactions*

*Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory*

*Electronic Structure: Basic Theory and Practical Methods*

*Introduction to the Kinetic Monte Carlo Method*

*An Introduction to Kinetic Monte Carlo Simulations of Surface Reactions*

*Monte Carlo Methods in Statistical Physics Chapter*

*Multiscale Simulation Methods in Molecular Sciences*

*Markov Processes: An Introduction for Physical Scientists*

*Data Structures and Algorithms*

*Introduction to Algorithms*