Kinetic modelling of heterogeneous catalytic systems

The importance of heterogeneous catalysis in modern life is evidenced by the fact that numerous products and technologies routinely used nowadays involve catalysts in their synthesis or function. The discovery of catalytic materials is, however, a non-trivial procedure, requiring tedious trial-and-error experimentation. First-principles-based kinetic modelling methods have recently emerged as a promising way to understand catalytic function and aid in materials discovery. In particular, kinetic Monte Carlo (KMC) simulation is increasingly becoming more popular, as it can integrate several sources of complexity encountered in catalytic systems, and has already been used to successfully unravel the underlying physics of several systems of interest. After a short discussion of the different scales involved in catalysis, we summarize the theory behind KMC simulation, and present the latest KMC computational implementations in the field. Early achievements that transformed the way we think about catalysts are subsequently reviewed in connection to latest studies of realistic systems, in an attempt to highlight how the field has evolved over the last few decades. Present challenges and future directions and opportunities in computational catalysis are finally discussed.


Introduction
The importance of heterogeneous catalysis in modern applications that enhance the quality of life cannot be overstated. A well-known example is the catalytic converter used in vehicles to oxidize CO and reduce NO x , thereby controlling emissions and improving air quality [1]. Iron-based catalysts for the production of ammonia by the Haber process is another notable example, the importance of which is evidenced by the fact that fertilizers, generated from ammonia thus produced, are responsible for sustaining one-third of the Earth's population [2]. Heterogeneous catalysts are commonly used in industrial practice, being Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. involved in the production of several key chemical compounds, such as hydrogen, sulfuric acid, nitric acid, hydrogen cyanides, and polymers [3].
However, the discovery of catalysts for chemistries of industrial interest is a non-trivial task. Traditional practice involves trial-and-error experimentation which is timeconsuming, costly, and does not guarantee success. On the other hand, theoretical studies at the molecular-and catalystlevel can provide a detailed understanding of how catalysts work and what affects their performance [4]. Moreover, modelling can be used to screen for potentially active catalytic materials before synthesizing the most promising candidates in the lab, thereby improving the efficiency of catalyst discovery. These ideas are central to the 'rational catalyst design' approach, which has been steadily evolving from a proof-of-concept to a revolutionary paradigm in the catalysis field [4][5][6][7][8][9][10].
To simulate the various aspects of catalytic behaviour, there are nowadays several theoretical and computational approaches at one's disposal. For instance, first-principles methods, such as the widely used density functional theory (DFT), can be employed at the molecular scale to understand the elementary events and the reaction mechanisms giving rise to catalytic activity. First-principles calculations thus provide important information about the chemistry of interest; yet, it may often be difficult to make statements about the performance of the catalytic material under operating conditions, especially in the presence of competing pathways with similar barriers. To this end, kinetic modelling studies are required to estimate activity and selectivity at the catalyst scale, which may employ Langmuir-Hinshelwood-type models, Sabatier analysis [11], mean-field micro-kinetic models [12], or more sophisticated statistical mechanical treatments, in particular kinetic Monte Carlo (KMC) simulation [13,14]. The latter is motivated by the structural complexity exhibited by real-life heterogeneous catalysts, which expose several different types of sites that may have distinct functionalities. For instance, sites at the interface between metal and support may behave differently than these two phases, planar versus low coordinated sites or defects may exhibit disparate catalytic activity, and alloys may expose a variety of sites at the interface of the two (or more) components. First-principles based KMC models have indeed been successfully used to unravel such complexities and aid in catalyst discovery; for instance, recent studies on ethanol synthesis from syngas identified effective promoters for Rh-based catalysts and made connections with experiments [15,16]. From an industrial/practical perspective, one is further interested in simulating the behaviour of chemical reactors, a task for which computational fluid dynamics (CFD) can be employed. More importantly, multiscale modelling frameworks can provide unified descriptions ranging from the molecular to the reactor scale [17,18], and are therefore playing an increasingly prominent role in theoretical catalysis and reaction engineering.
Thus, in the following section, we give an overview of the most popular methods at each scale, briefly discussing the main concepts and ideas that enable coarse-graining. We then focus on molecular-and catalyst-scale modelling using firstprinciples KMC methods, which are regarded as exceptionally promising in the detailed understanding of catalytic systems. We highlight the timeliness of this approach by showcasing a number of recent studies that have elucidated the complexity encountered in heterogeneous catalysts, and discuss current challenges and opportunities in the field.

Catalysis: a multiscale problem
Catalytic activity arises as a 'collective' phenomenon, in which myriads of elementary events (adsorption, desorption, diffusion, surface or Eley-Rideal reactions) happening on different sites of the catalytic surface, result in an increased rate for a net chemical reaction. This increase is with respect to the reaction rate in the absence of the catalyst, a material that is neither consumed nor generated during the process. Catalysis is inescapably a problem involving multiple spatial and temporal scales. The former encompass the electronic and molecular scales (Å), the catalytic nanoparticle scale (nm−µm), and the reactor scale (mm−m). The temporal scales encompass the timescales for the occurrence of an elementary event (ns−µs), as well as for the evolution of the dynamics on single nanoparticles (µs−s) and reactors (s-h). The theoretical understanding of catalytic phenomena at the fundamental level has thus motivated the development of bottom-up multiscale approaches (table 1). At the lowest level of these approaches lie detailed descriptions of the elementary events of adsorption, desorption, diffusion and reactions on the catalytic surface. Each such elementary event is the result of the rearrangement of atoms from a reactant configuration to a product configuration, which entails bond breaking and/or formation. From a physics perspective, such an event can be fully described by a timedependent Schrödinger equation for the system of interacting nuclei and electrons. In practice, this is a formidable task, and one resorts to simpler approaches: the Born-Oppenheimer approximation [19,20] can thus be used to parametrize the energy of the system with respect to the positions of the nuclei E (R 1 , R  , . . . , R N ), thereby yielding the so-called potential energy surface (PES) [19]. To calculate the latter, electronic structure methods [21,22] can be used, out of which DFT has emerged as highly popular in the catalysis community due to its computational efficiency. The PES provides a wealth of information about the catalytic system under investigation: local minima in this function correspond to stable species or binding configurations, and a reaction event is the transition from such a minimum to a neighbouring one by crossing the intermediate barrier ( figure 1). In our discussion, we will assume that electronic excitations do not take place and we will focus entirely on the ground state PES.
Conceptually, one could simulate a catalytic system, consisting of numerous reacting molecules, by employing a first-principles molecular dynamics (MD) approach, in which the force-field is prescribed by the aforementioned energy function, E (R 1 , R  , . . . , R N ).
To obtain this function, one could for instance use ab initio calculations to parametrize a reactive force-field (ReaxFF [23,24]), and there has been a significant amount of work that makes use of such methods on catalytic systems [24][25][26][27][28][29][30][31]. Yet, this approach becomes extremely demanding computationally when simulating catalytic systems for long timescales. The main cause of this challenge is that 'interesting' barrier crossings are rare events: the simulation algorithm would thus spend a lot of computational time wandering around a local minimum before a reaction happens. Even if one was able to simulate a few reaction events, this would hardly be enough to provide an overview of the catalytic function of the material, since the reaction mechanisms encountered therein often contain several elementary steps and parallel or competing pathways. Thus, a proper sampling of several such events would be needed, which would necessitate rather long simulations. KMC algorithms were designed to overcome this challenge by focusing on the timescale of the barrier crossings, while treating the wandering motion around a PES minimum in a coarse-grained sense. To this end, one only needs to locate the minima of the PES (stable molecular configurations) as well as the transition states; for the latter, several methods have been developed, a comparison of which has been performed by Klimeš and Michaelides [32]. Subsequently, transition state theory (TST) arguments are employed to describe the statistics of the transitions rather than the exact times of occurrence, durations and trajectories thereof. Central to TST is the assumption that the reactant and transition states in a reactive system are at quasi-equilibrium [33,34]. Thus, by making use of fundamental laws in statistical mechanics, one can derive the rate constant, k TST , for the crossing (see section 2.2.3 below). This constant provides the statistical description of the reaction occurrence: it does not reveal exactly when a reaction will take place in the system, it rather implies that within a time interval [t ini , t fin ], on average k TST · (t fin − t ini ) reaction events will be observed (Poisson statistics) per reacting molecule (or per pair of molecules for a bimolecular reaction). The average is taken over an ensemble of systems that have the same number of reactant molecules but different positions and momenta which satisfy the quasi-equilibrium condition of TST.
Rather than talking about the average number of events within a time interval, the theory of Poisson processes enables one to also make statements about the statistics of the waiting times between successive reaction events [35]. Under constant reaction conditions (T , P , compositions) these waiting times are exponentially distributed with rate parameter equal to k TST . This result lies in the heart of a KMC simulation: given the reaction mechanism of the chemistry under investigation, namely the elementary events and the rate constants thereof, KMC proceeds with the simulation by generating exponentially distributed random times for the occurrence of these events and executing them in the order they arise. This randomness gives rise to sampling error which can be reduced to arbitrarily small values by performing multiple simulations, or just long runs for ergodic systems under stationary conditions. A KMC simulation would in principle yield the same macroscopic observables (e.g. catalytic activity) with a corresponding MD simulation if both methods use the same PES, to the extent that TST faithfully captures the statistics of barrier crossings, providing accurate values for the rate constants. The predictive power of TST is of course a non-trivial subject: there exist several variants of TST approximations [36][37][38], and the development of novel, highly accurate methodologies has long been the subject of intense research.
It is possible to coarse-grain KMC simulations even further to derive deterministic models incorporating ratebased descriptions of reactions. System-expansion type of arguments [39][40][41] and mean-field approximations [42][43][44][45] can be employed to develop 'well-mixed' micro-kinetic (MK) models [12] widely used in the catalysis community, or spatially distributed models [46][47][48]. Such descriptions are derived on the basis of the following premises: (i) for large systems, reaction events occur extremely frequently on the macroscopic timescale, because of the large populations of molecules and consequently the large number of reactive collisions among them; (ii) the contribution of each individual reaction event in changing the number of molecules in the system is rather small due to the large populations of molecules therein. One can subsequently use scaling arguments to introduce a set of intensive variables (for instance, coverages or concentrations) and derive the macroscopic rate-based model as a limiting case of the underlying KMC model. In practice, the rate-based models incorporate further simplifying approximations. MK models, for instance, make use of the mean-field approximation which neglects spatial correlations, treating the adsorbate overlayer as homogeneous. Such an assumption limits the accuracy of MK models, since adsorbate ordering due to adsorbate lateral interactions is wellestablished experimentally and strongly supported by detailed theoretical models [45,[49][50][51]. Moreover, finite diffusion rates also result in spatial correlations which have been shown to exert a strong influence on the catalytic rate [52][53][54]. Improved MK models that take into account such effects have thus been proposed [51][52][53]; yet, a general-purpose framework is still under development.
Deterministic macroscopic models do have, however, significant advantages over the previous approaches discussed: their computational efficiency makes it possible to use them in reactor-scale simulations employing CFD frameworks. Moreover, coupling CFD with deterministic MK models (versus stochastic KMC models) avoids problems with accuracy [55,56]. Since CFD is considered as a mature field, employing tested and thoroughly validated modelling approaches and computational algorithms, recent research efforts have focused on the development of frameworks that integrate chemistry modules with CFD codes, culminating in the implementation thereof into computational packages such as CatalyticFOAM [57]. These advances pave the road towards the design of chemical reactors and processes from first-principles-based multiscale methods [58].
Still though, there are several challenges that need to be overcome, primarily related to the accuracy of the models at the individual scales, their ability to treat the several sources of complexity in catalytic systems, as well as the accurate coarsegraining and seamless coupling of models across different scales, so that no loss of important information is incurred. In our discussion, we will focus on the recent advances in KMC simulation methods, which couple the molecular with the catalyst scale, and appear to be rather promising in dealing with the complex features of heterogeneous catalysts.

Fundamental concepts behind KMC algorithms
Spatial KMC algorithms have steadily gained popularity in the catalysis community over the past two decades, being able to simulate spatial scales ranging from nm to µm, and temporal scales from ms to hours). Coupled with first-principles calculations in a multiscale type of approach [13,14,[59][60][61], KMC has matured into an indispensable tool for the kinetic modelling of heterogeneous catalysts. The most important advantage of the KMC method, lies in that it integrates several sources of complexity into a comprehensive framework: it thus models the spatial inhomogeneity of the catalytic surface due to the different facets or defects encountered in catalytic nanoparticles, faithfully represents the adlayer structure arising from spatial correlations, and captures the effect of neighbouring environment on reaction rate. Below, we examine the technical aspects of KMC that make it possible to model these phenomena.

The Master equation.
A KMC simulation captures stochastically the changes in the numbers and locations of molecules bound on a catalytic surface. The latter is represented as a lattice, namely a collection of discrete sites, onto which molecules can adsorb, diffuse, and react, and from which adsorbates can desorb. These elementary events are treated as Poisson processes under constant reaction conditions (see discussion in section 2.1). Let W denote the state space of the system, encompassing all the possible adsorbate configurations on the lattice and numbers of gas phase molecules specified up to additive constants (this will become clearer shortly), and let w ∈ W denote a state thereof. Moreover, let k w →w w represent the rate for the transition from state w to w due to an elementary event, and k w→w w represent the rate of the reverse transition. Then the stochastic equation that describes the dynamics of the system, known as the Master equation, can be written as follows [39,62,63]: The above equation expresses a balance for the probability of finding the system at state w at time tP (w, t). The sum in the right-hand side is taken over all the possible states of the system, with the understanding that if there is no elementary event that brings the system to state w from state w ∈ W, then k w →w w = 0. For instance, consider an empty lattice of N L sites in contact with a gas phase containing a single species that adsorbs on a single site (without dissociating). There will thus be N L possible lattice processes, each of which brings the lattice to a configuration where one specific site is occupied by a molecule of that species. There is, however, no elementary event that would result in the simultaneous adsorption of two molecules on the surface. KMC simulates a sequence of elementary events, which is nothing more than a random walk on W that obeys equation (1) and is also referred to as a stochastic realization of the corresponding random process. At stationary conditions, the probability at each state must stay constant (dP (w, t) /dt = 0) and thus the total probability flux leaving a state must be balanced by the total probability flux towards that state. Microscopic reversibility and the principle of detailed balance impose an even stronger condition [64,65]: each microscopic process has a corresponding reverse process, and the average individual rates of the two processes must balance at thermodynamic equilibrium: k w →w w · P eq w − k w→w (w) · P eq (w) = 0. (2) Imposing the detailed balance condition by individually specifying the rate constants of the transitions can be rather cumbersome. It is however possible to simplify things by introducing the microscopic energy of the system with respect to configuration E (w). This energy has electronic and vibrational contributions from the lattice species [66], as well as translational, rotational and vibrational contributions from the gas species that have not yet adsorbed on, or have already desorbed from the lattice: where E, E vib are the electronic energy and vibrational free energy of the given lattice configuration, whereas E gas , E gas transl , E gas rot , E gas vibr are the electronic energy, as well as the translational, rotational and vibrational free energies, respectively, of the gas species. In the above, we have assumed that the vibrational energy levels are with respect to the bottom of the potential energy well; if the ground state is used as the reference, zeropoint energy correction terms need to be added in equation (3).
The gas species reservoirs are implicitly taken into account in the KMC, since elementary events may consume or produce a given number of gas molecules. Thus, one needs to track the differences in the total energy of the system as a result of such changes. Since we are interested only in differences, the numbers of gas phase species are specified only up to additive constants (for instance they can all be set to zero initially).
At thermodynamic equilibrium, according to the Boltzmann distribution: and therefore, the transition rate constants must satisfy: The above equation can also be written as: where q w and q w are quasi-partition functions for the products and reactants, respectively, incorporating the vibrational contributions of surface species, as well as the translational, rotational and vibrational contributions of gas species, and E w →w rxn is the change in the electronic energy of the system, brought about by the reaction/transition from state w to w: Equations (6) and (7) can be used as guidelines in specifying the functional relations of the rate constants of the forward and reverse transitions. In the next two sections, we will discuss how this is achieved in the context of KMC simulations of catalytic surfaces.

Modelling lattice energetics.
For a reversible event that entails transitions between states w and w, the expression E w →w rxn in equation (7) accounts for the energy changes brought about by removing all the reactant molecules from the lattice and the gas phase (depending on the nature of the event), as well as adding the product molecules onto the lattice and gas phase. Accounting for the contributions of gas species or non-interacting lattice species is straightforward: given the electronic and vibrational energies of the species of interest, calculated for instance via DFT, one merely needs to count the number of molecules consumed or produced by the reaction for each species. However, depending on the catalytic system, lattice species may exert strong lateral interactions which lead to the formation of structured adlayers [45,49,50] and strongly influence the rates of elementary events [50,67]. These adsorbate-adsorbate lateral interactions thus have to be taken into consideration when building a KMC model, since they can have a significant effect on the quality of the predictions thereof. It is typically assumed that the major contribution of these interactions is on the electronic energy, since the vibrational modes of adsorbates are to a good approximation independent of coverage [68]. Therefore, the objective is to find E(w) in equation (3), which may be a complicated function of the state variable w (the latter of course incorporates all information about the coverage).  [13,59,69] developed first-principles-based KMC models that took into account two types of interactions: (i) through-space electrostatic contributions (attractive or repulsive) modelled by an explicit function of the distance between the interacting adsorbates; (ii) through-surface interactions that were modelled with either the use of pairwise additive contributions or within the bond-order conservation (BOC) framework [70,71]. For the pairwise additive interaction model, in the case of repulsive interactions, the number of contributions due to the neighbouring adsorbates are subtracted from the intrinsic binding energy found from ab initio calculations. In the case of attractive interactions, such contributions are added to the intrinsic binding energy. The total energy of the lattice E(w) is then calculated as the sum of the binding energies of all adsorbates.
The BOC framework, on the other hand, employs a different parametrization of the adsorbates' binding energy. BOC, which was later reformulated and renamed as unity bond index-quadratic exponential potential (UBI-QEP) [71], assumes that in a many-body system: (i) all forces are functions of distance only (spherical symmetry), (ii) interaction between two bodies (for instance an adsorbate and a metal atom) can be captured by a Morse potential, and (iii) the bond order between two bodies is an exponential function of bond distance, and the total bond order (sum of bond orders for all interacting pairs of bodies) is conserved at unity. Thus, the total bond energy in this framework is a sum of two-body contributions, and the binding energy of an adsorbate can be expressed as where BE A,0 is the binding energy at the top site in the absence of neighbouring adsorbates, n is the coordination number of the adsorbate (1 for an atop site, 2 for bridge, 3 or 4 for (1 1 1) or (1 0 0) hollow site, respectively), x A,i is the order of the ith bond: where r is the distance between the two bodies, r 0 the bond length at equilibrium, and b a scaling constant. It is then possible to estimate the interactions between an adsorbate and its neighbours by quantifying the metal atom sharing: this amounts to using an expression of the form (8) in which the bond order is expressed by considering the bonds of the metal atom [59,69]: where BE A is the adsorbate binding energy in the presence of neighbours, BE metal the binding energy of the metal, and m i the coordination of metal atom i. Given the binding energies of all adsorbates, the total lattice energy is the sum thereof, similarly to the case of the pairwise additive interaction model. Hansen and Neurock [59] observed that for several of the systems they considered, BOC overestimated the energetic interactions compared to ab initio calculations. To overcome this issue they used a linear rescaling of the BOC approach [59]. A generally applicable modified parametrization scheme that markedly improves the accuracy of the UBI-QEP method has also recently been developed by Maestri and Reuter [72].

2.2.2.2.
Cluster expansion-based models. The cluster expansion Hamiltonian [73] is another class of models that can be used to describe lattice energetics. Such models have long been used to study the thermodynamics of adsorption, adlayer ordering, phase change, and reconstruction on surfaces [74]. It is only recently that they have been incorporated into kinetic models [75][76][77], and have already shown remarkable potential in elucidating subtle features in catalytic systems, such as the inhomogeneity in the chemical reactivity of even structurally simple single-crystal surfaces [75].
The cluster expansion formulation is the most general method for capturing lattice energetics. This mathematical construct first devised by Sanchez et al [73] is based on the definition of a complete basis set using discrete Chebyshev polynomials n (σ ) for the description of a single lattice site with state variable σ = 0, ±1, ±2, . . . , and products thereof for constructing the basis functions for sets of lattice sites, referred to as clusters in this approach. A function αs (σ ) for cluster α containing lattice sites {p, p , . . . , p 2 } is formed as the following product: αs (σ ) = n σ p · n σ p · . . . · n σ p (11) where s = {n, n , . . . , n 2 }, and the configurational variable σ = (σ 1 , σ 2 , . . . , σ N L ) expressing the occupation of all sites. The variable σ is a component of the KMC state variable w, along with the numbers of gas species molecules N gas . αs (σ ) is thus a basis function on the configurational space of the lattice, and the set of all such functions can be used to represent any function of the lattice configuration, for instance the lattice energy: where α ranges over all clusters, namely combination of sites on the lattice, and s ranges over all indexes specifying Chebyshev polynomials. This approach has been recently generalized to lattices with multiple site types [78]. In a more intuitive formulation of the above principle, the total energy of the lattice can be decomposed into single-, pair-and many-body contributions, encompassing short-and long-range interactions. The total energy of the lattice is then given as the sum of all of these contributions [66,68,75,76,79]: where k ranges over all 'figures', which are sets of lattice sites with specific coverage patterns representing elementary interactions, n c k is the number of instances of such figures found on the lattice for the configuration w, m k is a multiplicity constant that prevents over-counting, and ECI k is the effective cluster interaction of figure k. In practice, the basis set is truncated to include interactions of up to a specific length (e.g. 8th nearest-neighbour [68]) and between a limited number of adsorbates (e.g. three bodies [68]).
Fitting the ECI parameters to ab initio data is conceptually straightforward: the problem is linear, and thus a least-squares formulation is sufficient. However, complications arise due to the large number of parameters that would have to be fitted if the procedure was done in a brute-force way [80]. This large number of parameters emanates from the fact that adsorbates can bind to several different site types, and adsorbate overlayers are multicomponent systems in catalysis applications; these factors result in a combinatorial explosion of the possible interaction patterns. Thus, instead of a one-off fitting, an iterative training procedure can be employed [81][82][83] (this can in principle also be done on the fly during the KMC simulation). Higher order patterns (for instance three-body terms) are not considered until the smaller and same-range lower order patterns have been exhausted, and the accuracy of the expansion is continuously evaluated by computing the cross-validation score [83]. The latter is a statistical measure of the error, obtained by leaving one configuration out of the training set, predicting that configuration's energy, repeating the 'leave-one-out and predict' procedure for all configurations, and summing the squared errors of these predictions. Moreover, Seko et al [84] have stressed the necessity of including states away from the ground state in the training set, a practice particularly important in systems with high entropy, which sample an extensive region of the configurational space.
Since cluster expansion-based models have only recently been employed in kinetic simulations in catalysis, there is still a lack of studies comparing them with earlier BOCbased models. Yet, the following remarks could guide the choice of a method for a particular application. With regard to representing lattice energetics, the ability of the cluster expansion to mathematically represent any function on a lattice is certainly a desirable feature, even though in practice one may not need to include terms beyond three-body interactions [68]. BOC-based models on the other hand employ a specific parametrization of the binding energy, and have been shown to deviate from the DFT data necessitating corrections [59]. Yet, an advantage of the BOC approach is that it can be used to also estimate activation energies, making it possible to build models for both energetics and kinetics from DFT-computed binding energy data. The refined scheme by Maestri and Reuter [72] has actually been shown to be rather successful in capturing the activation energy of dissociation reactions encountered in the water-gas shift and steam reforming chemistries: the predictions show a 10% variation around the DFT data for Rh(1 1 1) and Pt(1 1 1) at various oxygen coverages (within a 2 × 2 cell). It would further be of interest to investigate how the scheme generalizes for more complex surfaces, long-range interactions and multicomponent adlayers.  (1) are the result of microscopic processes such as adsorption, desorption, diffusion and elementary reaction events. The rate constants, or propensities, of these events can be calculated from TST as first proposed by Eyring in his seminal paper 80 years ago [85]. Central to this theory is the idea of treating the barrier crossing as a one-dimensional problem along the reaction coordinate, and averaging over all other degrees of freedom. TST further assumes that the reactant and transition states are in quasi-equilibrium: this enables the calculation of the probability flux through the transition state by considering the (equilibrium) probability of finding the system at the transition state, as well as the (equilibrium) distribution of momenta. These ideas lead to the following expression for the rate constant of an elementary event [33,86].
where κ is the so-called 'transmission coefficient', a correction factor accounting for trajectories that do not lead to product formation due to re-crossing the barrier; k B denotes Boltzmann's constant; h Planck's constant; T the temperature; q ‡ and q reac the quasi-partition functions of transition state and reactants, containing translational, rotational and vibrational contributions, E ‡ the activation energy of the event (the electronic energy of the transition state minus that of the reactants). For examples of how this equation is used to calculate the rates of specific events (for instance non-activated adsorption or Eley-Rideal reactions) the reader is referred to [61] and the supplementary material of [80]. As discussed in section 2.2.1, microscopic reversibility [65] dictates that for every forward event a corresponding reverse event is also possible. The rate of a reverse event must be such that the detailed balance condition holds (equation (6)), which can be written as follows by making use of equation (14): assuming that the transmission coefficients for the forward and reverse events are equal. Equation (15) has an important consequence: since the energy change for an elementary event E rxn is related to forward and reverse activation energies, it follows that any change in E rxn brought about by interaction of reactants or products with other adsorbates in the vicinity of the reaction event, must also affect the kinetics of the elementary event (figure 2(a)). A convenient way to model this effect, is by assuming a linear relationship between the activation energy and the reaction energy (figure 2(b)), as originally proposed by Brønsted, Evans and Polanyi [87]: where the constants α and β can be found by fitting firstprinciples data into the above relations. A parametrization that relates these constants with the activation and reaction energies in the zero-coverage limit (E ‡ f wd,0 or E ‡ rev,0 , and E rxn,0 respectively) as well as the 'proximity factor' ω of the reaction (a parameter quantifying how reactant-or product-like the transition state is [88]) has been recently implemented in KMC by Nielsen et al [76]: (16) or (17) ensure thermodynamic consistency in a KMC simulation [89] and allow for a quantitative (or semi-quantitative) treatment of the effect of lateral interactions on the reaction rate. The predictions of the barriers from these equations are of course limited by how well the linear functions fit the first-principles data, and recent studies that have implemented this methodology show errors of about 0.2 eV at most [75,77] (which is within the acceptable DFT error range of a few tens of an eV [90,91]). Alternatively, one can opt for more detailed approaches, for instance the energy of the transition state may be parametrized by a cluster expansion [92], or even explicitly defined for specific configurations of the neighbouring species [93,94]. The latter approach can quickly become intractable due to the combinatorial complexity of the possible configurations, especially for multicomponent adsorbate layers on multi-site lattices. The cluster expansion approach is more promising; however, it necessitates a large amount of calculations in order to fit each transition state's cluster expansion.
Given a set of elementary events whose kinetic constants are known from equation (14) (subject to (15)), one can generate random waiting times for the occurrence of these events on the lattice by noting that these times follow the distribution (see chapter 4.3: 'The next jump density function' in [62]): The above equation expresses the probability density that elementary event j will occur at time t +τ , given that the system is at state w at time t. Note that we have explicitly indicated the dependence of the rate constant on the state of the system, since for instance, a desorption event will exhibit different activation energies for the various possible local environments in the presence of lateral interactions. Moreover, we have denoted a temporal dependence of the rate constant, mainly to take into account cases in which the temperature changes over time, for instance in temperature programmed desorption (TPD) or annealing experiments. In the case that the rate constant is time independent, equation (18) simplifies to the exponential distribution: One can generate samples from the above distributions by appropriately mapping uniformly distributed variables (see chapter 1.8: 'Random number generating procedures' in [62]). For instance, using the inversion generating method for equation (18) necessitates solving the following non-linear equation:

Reaction Coordinate
where u is a uniform random number and τ j the waiting time sought. For time independent rate constant the above equation simplifies to which provides sample values of distribution (19).

Implementing KMC simulation computationally
2.3.1. General structure of KMC algorithms. The energetic and kinetic components discussed in the two previous subsections are essentially the foundation of KMC. The role of a KMC algorithm is to orchestrate the event simulation procedures such that the resulting trajectory follows the statistics prescribed by the master equation (1). The general structure of a KMC algorithm is shown in figure 3. After defining the simulation setup and parameters, with the specification of the lattice structure, the energetics model, and the reaction mechanism, the lattice is initialized with a desired configuration and the simulation starts. At each KMC step of the simulation, a lattice process is executed, the KMC time is advanced by an appropriately calculated random time, all data structures are updated and the simulation proceeds with the next KMC step. The result is a stochastic trajectory which contains information about the sequence of events executed, the resulting lattice configurations, as well as the instantaneous loss or appearance of reactant or product molecules, respectively, in the gas phase. This information can be post-processed to yield a wealth of information about the catalytic system considered: one can thus estimate catalytic performance metrics such as activity or selectivity, perform pathway analysis to identify the dominant pathway and the rate-determining step, and calculate species coverages to detect the most abundant reaction intermediate.
Most of the computational effort in a KMC simulation is directed towards the selection of the events simulated, as well as the update of the data structures holding information about the lattice state and the elementary lattice processes. For carrying out these operations, various techniques have been proposed, which exhibit different levels of efficiency. We will briefly outline the general principles and considerations behind these techniques; for a more detailed presentation the reader can consult [95], and chapter 3: 'Kinetic Monte Carlo algorithms' in [96].

Null-event method.
For determining which lattice process will take place, early studies employed the random selection of a lattice site and subsequently an event to be executed, as well as a neighbouring site for a two-site event (for instance the Ziff-Gulari-Barshad model for CO oxidation [97]). This approach is referred to as random selection [93] or null-event method [95], since there is always the possibility that the event selected cannot be realized, for example adsorption in an already occupied site, thereby resulting in a null-event.
With regard to when the executed event occurs, early studies used the number of Monte Carlo steps as a rough estimate of time. On the other hand, one can rigorously consider the time for the occurrence of any lattice process, which is generated by an equation similar to (20), but with the sum of rate constants entering the integrant: In the above equation, index j ranges over a set of N ∅ proc lattice processes, including the ones that are not active and would result in null-events. Also, we have used the fact that if u 1 is a uniformly distributed random number, 1 − u 1 will also be as such.
The computational effort for the null-event method, normalized by the number of sites, is lattice-size independent. However, the efficiency of this category of implementations strongly depends on the probability that a non-realizable event (null-event) will be selected. They perform rather well for simple cases, for instance, if there is a type of process that may occur practically anywhere on the lattice. The existence of several types of elementary events that occur infrequently as well as the effect of lateral interactions on the rate constant can often undermine the performance of the random selection methods (see section 3.9: 'Practical considerations concerning algorithms' in [98]).

Direct method.
Approaches that avoid the selection of non-realizable events belong to the class of rejectionfree methods [95]. Two such fundamental techniques that give statistically equivalent results were introduced by D Gillespie back in the 70s. The first one is referred to as the direct method [99], or variable step size method [100], and necessitates the generation of two random numbers: the first number represents the time for the occurrence of any lattice process, and is generated by an equation similar to (22), with the difference that only realizable (active) processes are now counted (the number of these processes is denoted by N proc ). The second random number indicates which process j exe will be executed, and is chosen such that the following inequality holds: Instead of considering each lattice process individually, one can group events of the same type into classes. Consider, for instance the case of adsorption of a single species on a lattice with a single site type and coordination number 4, and assume there are only pairwise nearest-neighbour interactions.
One can thus have adsorption on a site where the number of neighbouring adsorbates ranges from 0 to 4. This creates five classes of events; the multiplicity n q of each class q is equal to the number of empty sites with the corresponding number of occupied neighbours. This is referred to as the N -fold way [101] and allows for a more efficient search through the smaller list of classes. This approach, however, becomes impractical for systems with several species and complicated energetics models (long-range interactions or many-body terms), due to the high number of classes that need to be specified.
To reduce the computational expense of the direct method, one can alternatively perform the search for j exe by utilizing divide-and-conquer types of methods [94,[102][103][104]]. An efficient implementation thereof employs a binary tree storing the partial sums of the rate constants appearing in equation (23). Then the summation is done in constant time, since the head node of the tree will always contain the updated value of the sum. Searching through the tree for the process j exe in equation (23) can be done in computational time of O(log 2 (N proc )) at worst [103], as opposed to linear time for a brute-force implementation. Updating the partial sums after the execution of an event will also incur a O(log 2 (N proc )) cost.

First reaction method.
A conceptually different family of methods relies on the generation of random times for each of the possible lattice processes using equation (20) or (21), followed by the execution of the most imminent process, namely the one with the shortest inter-arrival time: This method is referred to as the first reaction method [99].
As originally proposed, one would have to generate as many random times at each KMC step as the possible processes; however, this is not necessary for independent processes [102], such as those involving species further from the maximum lateral interaction length. One can thus compute the absolute time t + τ j TST in which process j will be executed, and keep it intact until either of following occurs: (i) process j is simulated; (ii) another process takes place that consumes any one of the reactants of process j (in which case process j can no longer happen); (iii) the rate constant k j TST changes due to lateral interactions emerging as a result of another lattice process being executed within the vicinity of j (in which case τ j TST needs to be updated).
Efficient implementations of the first reaction method employ inverted lists for an efficient book-keeping of the dependencies between different events [13,59,76,105]. For instance, they keep track of the processes in which an adsorbate participates, thereby enabling an efficient deletion of the processes that can no longer happen in scenario (ii) just noted. Moreover, they incorporate binary heaps [76,102,105] in which the search for the most imminent lattice process can be done in constant time [103]. However, updating, inserting or deleting elements in the heap are now O(log 2 (N proc )) operations, and can be expensive if numerous lattice processes are to be stored therein. On the other hand, a major advantage of the first reaction method is that it is the most efficient approach for time-dependent rate constants [98].
The aforementioned algorithms simulate one lattice process per KMC step: during each such step, the state of only a small number of sites is manipulated, and any updates (of rates or possible processes) remain within the local neighbourhood. In the cases where massive updates are needed, for instance because of the existence of long-range interactions, one can trivially parallelize the computation of new rate constants. In a recent study, an OpenMP algorithm implementing such a scheme exhibited a speedup of about 10× for 16 processors in the case of a KMC model of NO oxidation with a cluster expansion Hamiltonian containing up to 8th nearest-neighbour interactions [76]. The efficiency exhibits a plateau for a high number of processors, since the overhead of scheduling and communication becomes significant compared to the workload of each processor.
Obviously, one can trivially parallelize KMC runs in the case of parameter sweeps, or by invoking multiple copies of the same simulation with different seeds for the random number generator. The latter approach would be useful in order to decrease the statistical error by increasing the sample size. Such a scheme would work well for a spatially homogeneous system; however, in certain physical problems this may not be the case. Thus, the underlying catalytic surface may be spatially inhomogeneous, such as in the case of nanoparticles which expose multiple facets and low coordinated sites. Moreover, the adsorbate layer may exhibit pattern formation [106]. In such cases, one may be interested in simulating large domains.
Thus, a different approach to parallelizing KMC simulations is via spatially decomposing the simulated domain and assigning each domain to a processor. This is a much more challenging approach because of the need for communication between neighbouring domains to avoid causality errors. The origin of the latter can be understood by a simple example (figure 4): suppose that neighbouring domains A and B are simulated by processors P1 and P2. If no attempt is made to synchronize the efforts of the two processors, it can happen that P2 finds itself at simulation time t B > t A , the simulation time of P1. If P1 executes a diffusion event that sends an adsorbate to domain B, any KMC events simulated by processor P2 for times after that diffusion event took place (t conflict = t A up to t B ) are now incorrect.
Tackling this type of problem in a rigorous way has motivated the development of algorithms that either prevent such issues, or detect and correct them by using 'rollback' procedures. Such methods are yet to be applied in the simulation of catalytic systems; however, we will briefly highlight them, as this is a potentially promising area of research. Thus, a conservative approach has been suggested by Lubachevsky [107,108] and later improved by Korniss and co-workers [109,110], for the simulation of Ising spin systems. Lubachevsky's algorithm allows the simulation of an event in a subdomain only under the condition that the waiting time of the first lattice process to occur therein is smaller than the corresponding waiting times of all neighbouring subdomains. This condition makes sure that no causality errors are incurred and therefore requires no rollbacks. Relaxation algorithms, on the other hand assume that causality errors are not too frequent, and incorporate rollback methods for correcting them. A widely known such algorithm is Jefferson's 'Time Warp' [111] which allows the multiple processors to propagate the simulated system asynchronously, and when a causality violation arises, the conflicted threads are rolled-back and propagated again in revised trajectories. Synchronous relaxation algorithms have also been developed [112,113] which have been used in the simulation of film growth [114,115].
The aforementioned algorithms can be complicated to implement and provide limited speedup. To tackle these problems, semi-rigorous algorithms have been proposed: the algorithm by Martínez et al [116], is a generalization of the N -fold way [101] and relies on the partitioning of a lattice into subdomains and computing the total rate constant for each of them. The largest such aggregate rate constant is used to compute rate constants of null-events for all the other subdomains (this step is necessary to ensure synchrony). A global waiting time is then generated and event selection proceeds in parallel for each subdomain. Appropriate rules for resolving boundary conflicts would result in the simulation of the correct statistics; however, such conflicts are treated in an approximate way in favour of computational efficiency. To this end, the synchronous sublattice method developed by Shim and Amar [117] can be used [118] thereby incurring a quantifiable error in the KMC simulation [119]. A comparison of rigorous synchronous algorithms and the semi-rigorous synchronous sublattice method for the simulation of island coarse-graining on Ag(1 1 1) revealed that the sublattice method provides the highest efficiency [120].
With the need for accurate simulations of realistic catalytic systems it is expected that the aforementioned methods will attract significant interest within the catalysis community. P2 has to roll-back in time (to point t conflict ) and re-simulate.

Computational codes.
Despite the power of the KMC method, there is currently a lack of widely used publicly available computational codes and software packages, in contrast to the abundance of several prevalent packages that perform ab initio simulations. As the focus of investigations of catalytic systems continuously expands to larger scales, it is expected that existing KMC packages will gain popularity and result in the creation of established user communities. We thus review some of the main KMC codes that have been developed in recent years. SPPARKS (Stochastic Parallel PARticle Kinetic Simulator) is a general C++ software package for on-and offlattice KMC simulations that has been developed by Steve Plimpton, Aidan Thompson, Alex Slepoy and co-workers at Sandia National Laboratories [121][122][123][124][125]. SPPARKS is open source publicly available, incorporates an approximate parallelization technique similar to the sublattice method discussed in section 2.3.5 (for more information see section 3.2 'Parallelism' of [125]), and employs a range of algorithms for KMC, in particular, rejection-free (direct) methods as well as null-event KMC (see sections 2.3.1-2.3.4). In addition, Metropolis (equilibrium) MC simulations can be performed. A Python interface for SPPARKS has also been developed, allowing the user to run the program interactively, and visualize the results. To simulate a specific system, the user has to write a SPPARKS script that involves commands which are subsequently executed by the program. SPPARKS employs a modular architecture such that a developer can expand and extend the code with new solvers, sweeping rules for event detection, as well as new applications (for the simulation of specific systems of interest). In fact, a SPPARKS module for heterogeneous catalysis simulations has recently been developed by Van den Bossche and applied to the hydrogenation of benzene on Pt(1 1 1) [126]. Yet, SPPARKS has mostly been used for materials problems, such as grain or thin film growth, diffusion, and sintering phenomena [127][128][129].
CARLOS is a general-purpose program developed by Johan Lukkien and Tonek Jansen at Technische Universiteit Eindhoven, for the simulation of chemical reactions on catalytic crystal surfaces [93,98,130,131]. CARLOS is distributed by the software company Dassault Systèmes BIOVIA (formerly Accelrys), as a Materials Studio package called Kinetix. CARLOS can simulate a variety of catalytic surfaces (with periodic boundaries) including stepped planes and defects, and also tackle time-dependent conditions (for instance while simulating TPD spectra). The effect of lateral interactions on reaction rate is treated by explicitly defining the different rate constants for the various configurations of the spectator species. This 'proliferation' of reaction events can, however, result in a large list of possible events, especially for multicomponent adlayers and long-range interactions. In addition, rotational symmetry is not taken into account by default, for example, a two-site event has to be defined explicitly for each and every neighbour on the lattice. A simulation is defined in a 'system file' that is subsequently parsed by the program, which can visualize the results on the fly. CARLOS has been optimized for speed and memory utilization, and employs a variety of solvers with different computational efficiencies [93] (first reaction method, random selection method, and variable step size method; see sections 2.3.2-2.3.4).
The kmos package [94] developed in Fortran90 and Python by Max J Hoffmann, Sebastian Matera and Karsten Reuter at Technische Universität München, is able to simulate periodic or aperiodic 2D and 3D systems, reaction events that involve multiple active sites, and multi-dentate species that bind to more than one sites. As in the case of the CARLOS code, the symmetries of the lattice are not taken into account, and thus one has to explicitly define elementary events involving equivalent collections of sites in different orientations. However, the kmos application programming interface (API) can tackle the definition of such events in a straightforward way. Lateral interactions are taken into account in this package by explicitly defining the rate constants for specific occupancies of distant sites (similarly to the CARLOS code). The solvers implemented are the ones described by Lukkien et al [93], and model definition, runtime and post-processing is handled by the API offered with kmos. The code is available from the CPC program library [132] and also from Github [133]. A number of projects are already using kmos, including the investigation of quantum effects in hydrogen activation [134], and the simulation of catalytic hydrogenation of benzene on Pt(1 1 1) [126].
Zacros is a Fortran2003 software package developed by Michail Stamatakis at University College London (UCL), based on the Graph-Theoretical KMC framework introduced by Stamatakis and Vlachos [105]. This package can tackle periodic or aperiodic 2D lattices of arbitrary complexity, elementary events involving multi-dentate species and several site types, as well as long-range lateral interactions and manybody effects. The program allows the user to define a cluster expansion Hamiltonian (see section 2.2.2.2) as the energetics model, independently of the reaction mechanism. Then, the influence of lateral interactions on reaction rate is treated within the BEP framework as described in section 2.2.3 [76]. If desired, the user can still explicitly define the kinetic parameters of a reaction for a specific arrangement of spectator species on the distant sites. Moreover, the symmetries of the lattice are inherently accounted for in Zacros, while also giving the ability to restrict the occurrence of reaction events to specific directions, for instance on a 4-fold lattice, dissociative adsorption of O 2 can be specified on only vertically arranged pairs of sites. Keyword-based text files are used to set up a simulation, whereas MATLAB scripts can be used for post-processing the results. OpenMP parallelization has been implemented to tackle the detection and update of the processes affected by lateral interactions after each reaction event. Significant speedup can be achieved for systems involving long-range interactions between adsorbates. Zacros is distributed free of charge for academics via e-Lucid, the online licensing portal of UCL [135]. The graph-theoretical KMC and the Zacros implementation have already been used in a variety of applications, such as the simulation of structure sensitivity effects for the water-gas shift reaction on Pt [136]; CO oxidation on Au nanoclusters [137,138]; dissociation, spillover and associative desorption dynamics on bimetallics [8,139,140]; NO oxidation and NO 2 reduction on Pt [76], as well as CO oxidation on Pd [77].
Finally, a very recent addition to the list of available KMC codes comes from Mikael Leetmaa and Natalia V Skorodumova, in KTH-Royal Institute of Technology and Uppsala University, Sweden [141]: KMCLib is a Python module with a highly optimized C++ library that incorporates the N -fold way for event selection and execution, thereby being able to simulate the reaction and diffusion of millions of particles in up to 3D geometries. MPI parallelization is employed in the detection of lattice processes (by matching a process with a local geometry) and the rate calculations thereof, and has been shown to be efficient whenever these procedures are computationally expensive, for instance in the presence of long-range interactions. Analysis modules that can perform on-the-fly calculations of diffusion rates for instance, are implemented as plugins via an API. KMCLib is well-suited for applications involving diffusion on surfaces or in solids, and its capabilities and performance have been demonstrated in two systems: (i) a conceptual one-dimensional model involving three types of particles, two of which can flip locations reversibly while particles of the third type are inert; (ii) oxygen-vacancy diffusion in a 3D fluorite structure of a metal oxide [141].

Physics unravelled by heterogeneous catalyst models
Kinetic models of heterogeneously catalysed reactions have shed light into the underlying chemical mechanisms and physical phenomena that give rise to catalysis, or affect the performance of a catalyst. A review of chemical systems simulated in detail with first-principles KMC has been presented before [80], whereas here we focus on the physical phenomena unravelled with KMC. Our discussion presents conceptual studies on model systems but also first-principlesbased KMC models of realistic catalysts, in order to make the link between the two and hopefully stimulate further studies on the topics highlighted.

Complexity and emergent behaviour
Complexity and emergent behaviour are ubiquitous in nature, exhibited by a plethora of physical systems including catalytic reactions. Thus, kinetic phase transitions and bistability [142], as well as oscillations and wave formation [143] have been observed in a variety of catalytic systems. Such phenomena have long been the focus of both experimental and theoretical studies.
A kinetic phase transition is a qualitative change in the dynamic behaviour of the system as a result of altering a parameter. A simple KMC simulation scheme of kinetic phase transitions in a conceptual model of the CO oxidation reaction was first proposed in the seminal work by Ziff, Gulari and Barshad (ZGB) in 1986 [97]. The ZGB model, assumes irreversible adsorption of CO and O 2 (the latter adsorbs dissociatively), and rapid surface reaction between CO and O adatoms toward CO 2 gas. For CO-lean feeds, the surface is poisoned by O, and as the molar fraction of CO (y CO ) increases, a 2nd order transition occurs into a reactive region ( figure 5). This region disappears at a 1st order transition, after which the catalytic surface gets poisoned by CO ( figure 5(a)). In the ZGB model, the adsorbed species are immobile, a limitation that was overcome in a subsequent model by Jensen and Fogedby [144]. This work showed that accounting for diffusion extends the range of the reactive region by shifting the 2nd and 1st order transition points to lower and higher values of y CO , respectively, and increasing the rate of CO 2 production. Moreover, fast diffusion makes the O-poisoned phase eventually disappear. These effects were attributed to the homogenization achieved by diffusion, which reduces the density of empty site pairs and breaks the islands of O and CO on the surface. Since the original publication by Ziff et al [97], several studies employing KMC and mean-field models have investigated kinetic phase transitions in variants of the ZGB and other conceptual models (reviewed in [142,[145][146][147]). In order to tackle the more realistic case of automotive catalysts, in which the exhaust mixture contains nitric oxides in addition to CO, Basit and Ahmad considered adding NO to the gas phase [148]. Thus, to study the effect of this species on the kinetic phase transitions of the CO oxidation reaction, Basit and Ahmad [148] combined the ZGB model for CO oxidation by O 2 and a model by Yaldram and Khan for the NO-CO reaction (YK model) [149]. It was thus shown that the 1st order transition of the ZGB model gets converted to a continuous 2nd order transition, even for small NO concentrations. On the other hand, Buendía and Rikvold studied the effect of gas phase impurities (non-reactive) on the ZGB model with O 2 adsorption on next-nearest-neighbour sites (to eliminate the Opoisoned state) [150]. If impurities are not allowed to desorb from the surface, Buendía and Rikvold's model shows that the CO and O coverages vary continuously with respect to y CO , no CO 2 is produced at steady state, and the 1st order transition is no longer observed. In the case that impurities can desorb, the reactive region reappears, and its range depends on the concentration and desorption rate of the impurities.
Zhdanov and co-workers studied the effect of surface oxide formation on the bistable behaviour and the hysteresis of CO oxidation on Pt [151,152]. The model incorporates reversible formation of surface metal oxide by adsorbed oxygen adatoms, and subsurface oxide by diffusion into the subsurface layer. A number of kinetic parameters of this KMC model were fixed to physically reasonable values and an exploration of the remaining parameter space was performed [151]. Thus, the oxide formation propensity and its effect on reaction rate were investigated with respect to the energetics of surface and subsurface oxide. The model predicts that surface oxide is formed in islands on the nanometre scale and suppresses CO adsorption and reaction, thereby making CO impingement much faster than the CO oxidation process, in agreement with experimental observations [152]. The range where CO coverages exhibit bistability is greatly influenced by oxide formation, for instance the bistable regime diminishes when surface (but not subsurface) oxide is stable. In the case where surface and subsurface oxide formation is possible, the shape of the hysteresis loop is markedly affected, making the transition from high to low CO coverages smoother [151].
The phase behaviour of the CO oxidation has also been investigated by first-principles KMC models; in particular, Reuter et al [14,60,61] studied this reaction extensively on the RuO 2 (1 1 0) surface. In these studies, DFT was employed to investigate the binding of CO and O species on the bridge and cus sites of the oxide, as well as all the possible diffusion and elementary reaction events therein. Subsequently, KMC simulation was used to map the catalytically active regions as well as the poisoned states on a 2D phase diagram with the partial pressures of CO and O 2 as the control variables. The transitions between the different regions were found to occur continuously over a range of pressures, rather than abruptly (as in the case of 1st order transitions).
Kinetic oscillations have also been observed in catalytic systems, and thus, studies of kinetically driven oscillations have been performed. Noussiou and Provata have developed a surface oxide model for CO oxidation, showing that a large difference in the reactivity of the metallic versus the oxidized surface can result in oscillatory behaviour [153]. Zhdanov has studied the reduction of NO by CO on Pt(1 0 0) at low temperatures, under which no reconstruction occurs [154]. The model accounts for CO and NO adsorption, NO decomposition to N and O, CO oxidation by the O adatom, as well as associative desorption of N. Parametrization was done in an empirical way, in line with experimental observations. Sustained oscillations were reproduced by this model for approximately equal CO and NO impingement rates, slow desorption thereof, and fast NO decomposition and CO oxidation. The qualitative features of the oscillations were further studied for various sizes of the four-fold (rectangular) lattice: it was found that lattices of size 15×15 or larger  [154]. Such dynamical behaviours are reminiscent of the oscillations exhibited by the Belousov-Zhabotinsky reaction [155] and similar systems. Still though, oscillatory dynamics in catalytic systems are to a large extent attributed to catalyst restructuring and will be discussed in the following section.

Catalyst reconstruction
Even when a catalytic reaction takes place on a welldefined crystal plane, the catalytic surface may not be static. Under certain conditions (temperature, composition of the gas phase), restructuring (also referred to as reconstruction) can occur, which changes the coordination numbers of the surface atoms and can profoundly affect the reaction kinetics. One of the most interesting effects of catalyst reconstruction is perhaps the emergence of oscillations and chemical waves [106]. A number of KMC models that incorporate catalyst reconstruction have been proposed to explain such phenomena. A review of such modelling studies has been authored by Zhdanov [156]; below, we highlight seminal and recent studies.
Zhdanov has thus proposed a conceptual KMC model to study the effect of catalyst reconstruction on the CO-NO reaction on Pt(1 0 0) [157]. This model assumes a square lattice of the Pt atoms and captures reconstruction by considering two different energy states for each Pt atom, referred to as 'stable' and 'metastable'. The actual lattice remains static during the simulation, and therefore, changes in the density of sites are not modelled. The energetic interactions between the stable and unstable states are, however, accounted for, to enable the propagation of the reconstructed phase along the lattice (hence, same states exhibit attractive interactions, different states repulsive). For the switching between the two states, a Metropolis criterion is applied, which takes into account the adsorbate-catalyst energetic interactions. The reaction mechanism is the same as in the later work by Zhdanov [154], discussed in the previous section, but with NO dissociation occurring only on the (1×1) reconstructed phase, modelled by the metastable atoms. The main findings of this model pertain to the qualitative characteristics of the oscillations and the islands of reconstructed catalyst. Thus, sustained oscillations (figure 6(a)) and island formations of reconstructed catalyst (figure 6(b)) were observed for a relatively narrow, yet not unrealistically so, range of CO pressures. The oscillation features and the size of the islands were found to exhibit a weak dependence on the relative rates of NO and CO diffusion with respect to the rest of the elementary steps: faster diffusion results in oscillations with larger amplitude and longer periods, and also promotes the formation of slightly larger islands [157]. To assess the stability of the observed patterns in the presence of defects on the catalytic surface, Pt atoms were forced to stay in the metastable state; this however, did not affect the stability of the oscillations and islands [157]. In addition to the CO-NO reaction, Zhdanov has also applied this KMC simulation scheme to study the NO-H 2 reaction on Pt(1 0 0) [158][159][160], and the CO oxidation on the same catalytic surface [161], also investigating the effect of system size [162].
More recent studies of the CO-NO reaction on Pt(1 0 0) have been presented by Zgrablich and co-workers [163,164]. In these models, the change in the density of sites is also neglected, and the emergence and propagation of the reconstructed phase is modelled by rule-based mechanisms: the nucleation mechanism enables the hex → square (1 × 1) reconstruction when five adsorbed molecules (CO or NO) cluster together. The trapping mechanism on the other hand, allows for the reconstruction to propagate to occupied sites that are nearest-neighbours to a reconstructed cluster. Moreover, Zgrablich and co-workers [163,164] took into account the different reactivities of the two catalytic phases towards NO dissociation, and incorporated a novel elementary step in the reaction mechanism according to which adsorbed NO and N react giving rise to an activated complex (N-NO) • that subsequently decays to N 2 gas and an O adatom. The existence of sustained oscillations was explained on the basis of the sequential emergence of reactive and recovery stages. In particular, during an oscillatory cycle, CO and NO adsorb on the hex phase triggering reconstruction. NO dissociates on the reconstructed phase only, and the products, N and O, react with CO and NO (or desorb associatively in the case of N). These reactions are fast, thereby resulting in an increase of the density of empty sites, which triggers the recovery to the hex phase, in preparation for the onset of the next cycle. As a result of this mechanism, reaction fronts can be observed on large catalytic surfaces, when the mobility of species is appreciable. The NO dissociation rate was found to modulate the period of the oscillations, with higher (lower) rates resulting in shorter (longer) periods, or even the cessation of the oscillations for rates beyond a specific range. More importantly, Zgrablich and co-workers [163,164] discovered that the emergence of oscillations and reactive fronts are rather sensitive to the rates with which the surface reconstructs: too fast hex → (1 × 1) reconstruction rates or too slow rates for the opposite transition, both have a detrimental effect to the regularity of the oscillations. Finally, the effect of temperature was investigated, and it was shown that higher temperatures result in higher oscillation frequencies, in agreement with experiments.
Rafti and Vicente [165] used Zhdanov's [158] approach of modelling reconstruction in order to study the reduction of NO by NH 3 on the Pt(1 0 0) surface. The NO induced reconstruction, hex → (1 × 1), was found to result in oscillations in this system as well: for temperatures in the range of 390 to 410 K, aperiodic oscillations were obtained, whereas for higher temperatures, namely 410-480 K, periodicity was observed. Outside of this temperature range, the system exhibits negligible or no reactivity: for T < 390 K, the surface exists predominantly in the hex phase, onto which NH 3 does not adsorb, and NO coverages are low. For T > 480 K, the system exhibits a NO-covered (1 × 1) phase and negligible turnover frequency [165]. It is important to note the qualitative differences with the case where no reconstruction takes place: by repeating the simulations after forcing the surface to always be at the (1 × 1) state, Rafti and Vicente [165] observed a NH 3 -rich surface for T < 420 K, and a NO-rich surface for T > 540 K, both non-reactive. In the intermediate range of temperatures (420-540 K) the system exhibits a finite turnover frequency, but no oscillatory features in its temporal behaviour. These findings are in agreement with experiments [166] and underscore the importance of reconstruction on pattern formation in catalytic systems.
Reconstruction can also occur in catalytic structures other than surfaces: in a recent study of the CO oxidation reaction on the Au − 6 nanocluster supported on MgO, DFT calculations revealed a 'breathing' mechanism for catalyst reconstruction that results in carbonate (CO 3 ) formation [137]. The latter process was not detected when the Au − 6 cluster atoms were frozen in the calculation setup. The energetic and kinetic data obtained from DFT for all CO oxidation and CO 3 formation pathways, were subsequently incorporated into KMC simulations, which showed that carbonate is rapidly formed on the Au − 6 cluster and, due to its high stability, poisons the catalyst. This study thus highlights the effect of reconstruction on catalytic activity at the nanoscale.
In all previous studies, reconstruction processes were taken into account in the energetics and kinetics of the elementary steps of the reaction mechanisms considered; the lattices remained static. While such approaches have provided useful insight, they are unable to account for changes in the density of active sites, or capture the emergence of new sites upon roughening or re-faceting. Conceptually different approaches which overcome these challenges by taking into account the dynamics of the lattice have thus been developed.
In particular, Noussiou and Provata [167] studied reconstruction effects on the ZGB model [97], by switching between 4-and 6-fold lattices thereby mimicking the behaviour of the (1 × 1) to hex reconstruction. This switching was done globally in the simulation domain, a scheme that was justified on the basis of the sizes of the reconstructed islands, which are found experimentally to be on the order of the simulation domain size [167]. As a result of this reconstruction mechanism, the model exhibits sustained oscillations.
Moreover, Cherezov et al [168] have presented a conceptual model of reconstruction on a bcc metal lattice, employing Metropolis MC for metal atom hops, and KMC to capture the chemistry. The latter involves dissociative adsorption of a species B 2 , single-site adsorption/desorption of a species A, and the bimolecular reaction between A and B similar to the ZGB model [97]. It is found that a rough surface is smoothened out by stabilization of flat areas due to adsorption of molecules. Moreover, the mobility of the catalyst atoms shifts the 1st and 2nd order transition points to lower partial pressures of component A (thus, higher B fractions). These shifts were attributed to the hindering of adsorption of the bimolecular species by the surface roughness, as this type of adsorption requires a flat site of six catalyst atoms.
Finally, Monine et al [169,170] have introduced hybrid models for simulating CO oxidation on reconstructing Pt(1 1 0) surfaces. In these models, the dynamics of the Pt atoms are treated within a KMC framework (parametrized using DFT and experimental data), whereas the surface coverages of CO and O species are treated as continuous mean-field variables averaged over the surface. This setting is justified by the higher adsorbate diffusion rates compared to the Pt atom hopping rates. Thus, the KMC model for reconstruction is coupled with a system of deterministic mean-field rate equations capturing the chemistry. The models successfully reproduce the conditions (CO pressure) under which roughening and faceting occurs, as well as the emergence of a periodic hilland-valley structure.
This and the aforementioned studies are particularly promising, and pave the road towards accurate, unified descriptions of the structural dynamics and the chemistry of catalytic systems.

Effects of adsorbate lateral interactions
When adsorbates bind to a catalytic surface, the resulting overlayers are not uniform; rather, ordered structures are formed due to the repulsive and/or attractive lateral interactions between the adsorbed entities [50]. These interactions are also known to affect the rates of elementary events [67,[171][172][173][174], by stabilizing or destabilizing the initial, transition or final states thereof, resulting in different activation barriers (figure 2).
For instance, repulsive interactions would destabilize the bound state and therefore promote desorption of an adsorbate from the surface. In addition, ordering in the adlayer can affect the density of 'reactive sites' [75]; for example, pairwise repulsions would result in a lower density of nearest-neighbour occupied sites on the catalytic surface.
Early studies employing analytical and computational methods were able to unravel the richness of such phenomena in catalysis, by focusing on specific elementary events [173]. For instance, analytical approximations to lattice-gas model solutions explained the drop in the heat of adsorption at high coverages as well as qualitative features of TPD spectra, such as multiple peaks [175]. In an early computational study employing KMC, Weinberg and co-workers [176] modelled N 2 dissociative adsorption and associative desorption via a precursor on Ru(0 0 1). Nearest-neighbour pairwise additive interactions between precursor and chemisorbed species, as well as between two chemisorbed molecules, were taken as repulsive, whereas next-nearest-neighbour interactions were taken as attractive. This model reproduced the experimentally observed initial increase of the adsorption coefficient with respect to coverage, due to the formation of energetically favourable ( √ 3 × √ 3)R30 • islands, as well as the subsequent decrease of this coefficient at high coverages due to crowding effects, which on average intensify the repulsions. TPD spectra were also qualitatively captured: at low coverages, the spectra exhibit one peak attributed to desorption from the peripheries of the islands. At high coverages, a low temperature peak emerges, because of desorption from grain boundaries between antisense domains, whereby lie loosely bound species due to repulsive interactions.
Several further studies have been carried out, for instance, Lombardo and Bell [177] investigated the effect of metal-adsorbate and adsorbate-adsorbate interactions using BOC-based models for the energetics, incorporated into KMC simulation. These models successfully captured the experimentally observed numbers and locations of TPD peaks for CO desorbing from Pd(1 1 1), as well as H 2 from Mo(1 0 0) and Ni(1 1 1). Meng and Weinberg [178] employed KMC simulation with a lattice-gas model and analysed the TPD features for attractive, repulsive or absent next-nearestneighbour interactions, in the presence of repulsive nearestneighbour interactions. Models of co-adsorbed species have also been developed by Gupta and Hirtzel [179] using latticegas energetics, as well as Lombardo and Bell using BOC [180], and were shown to exhibit complex features, not present in the case where pure species desorb from the catalytic surface. Thus, peaks can broaden or shift to lower temperatures due to repulsions; multiple peaks can merge, and new peaks can form as well. Attractions can lead to sharper peaks, shift the location thereof to higher desorption temperatures and markedly affect the overlayer structure, for instance, they can lead to segregation.
In several of the aforementioned studies adsorption/ desorption and diffusion events were incorporated in the simulations. It is, however, worth mentioning that the effect of lateral interactions on diffusion alone has been characterized in the early studies by Bowker and King [171,172]. These studies focused on the effect of lateral interactions on the diffusivity of adsorbed molecules on a square 2D lattice (with four-fold symmetry), but also in 1D geometries, relevant for example in the case of stepped surfaces. KMC simulations and comparisons with an analytical model, based on order-disorder theory in the quasi-chemical approximation framework, were performed for the ideal case of no interactions, as well as in the presence of repulsive or attractive interactions. Pairwise additive nearest-neighbour [171] and next-nearestneighbour interactions [172] were considered. In the case of no interactions, the model reproduces Fickian diffusion, as expected. For repulsive interactions, the coverage profile of the diffusing adsorbate becomes asymmetric, with a low slope in the high-coverage regime and a sharp front in the low-coverage region, giving the impression of a moving boundary. Attractive interactions result in sharper, slightly asymmetric coverage profiles. The coverage dependence of the diffusion coefficient was also investigated: for no interactions this coefficient remains constant. However, in the presence of repulsive interactions, higher coverages result in increased mobility (higher diffusion coefficients). The opposite effect is found for attractive interactions, where higher coverages slow down diffusion. By combining nextnearest-neighbour attractive with nearest-neighbour repulsive interactions [172], Bowker and King were able to provide a qualitative explanation for the maximum in the diffusion coefficient with respect to coverage, observed experimentally for O adatoms on W(1 1 0) [181]. Subsequent studies on the diffusion of interacting adsorbates have elaborated on the effect of different types of interactions (attractive, repulsive, nearest-and next-nearest-neighbour) [182][183][184][185], correlated non-monotonic trends of the diffusion coefficient with the phase transitions of the overlayer [186,187], focused on lattices with multiple site types [188][189][190], analysed memory effects in tracer and collective diffusion [191], and investigated the correlation function of fluctuations in the ordered structure of adlayers [192,193].
Conceptual models of catalytic reactions between interacting adsorbates have been developed as well. Starting with the ZGB model for CO oxidation [97], Wolf and co-workers [194] incorporated nearest-neighbour interactions by including normalization factors in the rate constants that depend on the coverage of neighbouring sites. Luque et al [195] incorporated interactions into the ZGB model by expressing the transition probabilities of CO/O 2 adsorption and CO+O oxidation via Arrhenius equations, with barriers that depend on nearest-neighbour coverages. Satulovsky and Albano [196] employed a Maxwell-Boltzmann type dependence of the rate constants on the interaction energies. In the studies just noted, no diffusion was taken into account. Kaukonen and Nieminen [197] did consider diffusion, and accounted for energetic interactions through Boltzmann terms that captured the changes in the activation energies of elementary events. These studies revealed that for repulsive interactions that facilitate desorption (or hinder adsorption), the 'interacting ZGB' model exhibits the same kinetic phase transitions as the original one, namely a 2nd order transition from the O-poisoned phase to the reactive phase, and a 1st order transition from the reactive to the CO-poisoned phase. Yet, the reactive window is wider in the presence of interactions. Attractive interactions that facilitate adsorption (or hinder desorption) result in a narrowing of the reactive region and eventually the loss thereof, in which case the system exhibits only a 1st order transition from the O-poisoned to the COpoisoned regime.
More recent studies have aimed at developing realistic KMC models of reaction kinetics on catalytic surfaces. Völkening and Wintterlin have presented an elegant model of CO oxidation on Pt (1 1 1), motivated by scanning tunnelling microscopy experiments [198]. A titration protocol was employed in which pre-adsorbed O is exposed to gas phase CO and allowed to react. The simulations employed a lattice with a single site type of coordination number 6, and a variety of energetic models, yielding results that were compared to experiments. It was thus found that the adsorbate domain morphologies can be adequately captured by a KMC model that includes attractions between O adatoms and their effect on the activation energy of the CO oxidation reaction. The diffusion of O and CO (the latter being fast) were also necessary to include in the model.
An ab initio KMC study of a similar system, namely a titration experiment on O-precovered Pd(1 1 1), has recently been conducted by Piccinin and Stamatakis [77]. This study was also motivated by experiments, which suggested that different overlayer phases exhibit different reactivity, with the p(2×2) being inert, and the ( √ 3 × √ 3)R30 • and p(2×1) exhibiting distinct kinetic behaviours [199]. To capture the lateral interactions, a detailed cluster expansion Hamiltonian was developed, incorporating pairwise additive terms up to the 3rd nearest-neighbour, as well as three-body terms, for CO and O species bound to face-centred cubic (fcc) and hexagonal close packed sites. The activation energies in the reaction mechanism were expressed via Brønsted-Evans-Polanyi (BEP) relations (section 2.2.3), in order to explicitly account for the influence of lateral interactions on the rate constants of the elementary events (adsorption/desorption, diffusion and reaction). These simulations were able to reproduce the experimental coverage profiles (figure 7(a)), as well as the quantitative behaviour of the system, in particular, the reaction orders and apparent activation energies. A sensitivity analysis with respect to the next-nearest-neighbour O-O interaction parameter was performed and suggests that, while there may be a correlation between the emergence of different phases in the adlayer and the variable reactivities, the presence or absence of ordering ( figure 7(b)) leaves the catalytic properties of the system practically unaffected. Hence, the variable reactivity was attributed to the different local environments around a CO+O pair that exist at different adsorbate coverages.
It is thus interesting to note that among the several reactive configurations present at any given time on the catalytic surface, only those with the lowest activation energies contribute to the overall activity (figure 7(c)). A similar effect has also been observed in an earlier first-principlesbased MC model of NO oxidation on Pt(1 1 1) by Schneider and co-workers [75], in which O adsorbed on fcc sites was taken into account. A very detailed cluster expansion model of lateral interactions was utilized, incorporating up to 8th nearest-neighbour interactions and two triplets (three-body terms). By analysing the statistics of O 2 dissociation events, which are rate limiting for this reaction under the conditions studied, Schneider and co-workers showed that the dominant contribution to the turnover frequency comes from events occurring on pairs of empty sites that have a low number of occupied neighbours. These site pairs are precisely the ones with lower activation energies, given that the dominant interactions between O adatoms are strongly repulsive.
Schneider and co-workers have very recently extended this methodology to a wide range of metals, in particular Ru, Os, Co, Rh, Ir, Ni, Pd, Pt, Cu, Ag and Au [200]. By comparing models that incorporate lateral interactions with models that neglect such interactions, it was shown that the traditional 'volcano plot' of activity [11], in line with the Sabatier principle, has to be refined: most metals investigated actually exhibit rates that are uncorrelated with those of their non-interacting (conceptual) counterparts, and their activities depend on the deviations from mean-fieldlike behaviour [200]. The authors further argue that strong repulsive lateral interactions are critical for the catalyst retaining its functionality over a wide range of operating conditions, by stabilizing the oxygen coverage, but should not be too strong such that they negatively affect reaction rates.
First-principles kinetic models incorporating lateral interactions can also shed light into puzzling experimental observations regarding kinetic phase transitions of certain catalytic systems. In particular, Zhdanov and Matsushima [201], presented a DFT-parametrized KMC model of the reaction between N 2 O and CO on Pd(1 1 0), to explain a counterintuitive feature of this system, pertinent to the existence of a first-order kinetic phase transition without hysteresis being observed [202]. The proposed reaction mechanism takes into account decomposition of gaseous N 2 O into adsorbed O and N 2 gas, CO adsorption, and reaction between bound CO and O towards gaseous CO 2 . A mean-field model of this system clearly exhibits hysteresis; however, MC simulations revealed that in the presence of attractive interactions the hysteresis regime shrinks, providing a possible explanation why such a feature was not detected experimentally. It is interesting to note that in this model, the lateral interactions between adsorbed O-O pairs were taken to be anisotropic, to implicitly model the missing row reconstruction on the Pd(1 1 0) surface. Thus, the interactions along the horizontal and vertical directions were taken as repulsive and attractive, respectively [201].
Even though lateral interactions have been found to have a marked effect on the observed catalytic behaviour, not all observables may exhibit a strong dependence thereon. For instance, Mei et al [203] studied the selective hydrogenation of acetylene on Pd(1 1 1) using the first-principles KMC method by Neurock and co-workers [13,59] of the full model, it was found that the overall apparent activation energy decreases only slightly in the absence of lateral interactions. Still though, the surface coverages differ considerably between these two cases, and taking into account the interactions is important for a correct prediction of reaction orders and selectivity, the latter being severely underestimated in the absence of interactions.

Structure sensitivity
One of the most intriguing problems in catalysis is how the geometric structure of a material affects its catalytic performance, if at all. It was indeed observed early on, that for some reactions, the size and/or shape of catalytic nanoparticles markedly affects activity [204]. This effect was attributed to nanoparticles of various morphologies exposing different Miller planes (with distinct surface atom coordination numbers), and exhibiting different densities of low-coordinated sites. Defect sites may also be present, depending on the catalyst preparation. A reaction may proceed on these sites via different pathways, characterized by diverse energetic landscapes, and may therefore occur preferentially on a specific subset of catalytic sites. Such reactions are referred to as 'structure sensitive' and require detailed investigations to unravel the dependence of activity (or possibly selectivity) on structure, with the ultimate aim of optimizing catalytic performance.
Theoretical investigations can facilitate this task and aid in obtaining a fundamental understanding of the underlying physics of structure sensitivity.
Early studies have thus investigated the behaviour of models of catalytic reactions as a function of lattice geometry. Focusing on the kinetic phase transitions of the ZGB model, Meakin and Scalapino [205] pointed out the importance of geometric factors pertinent to the dissociation of O 2 on a pair of empty nearest-neighbour sites. In particular, while the original ZGB simulations were conducted on a rectangular lattice (coordination number 4), Meakin and Scalapino [205] extended this work to a hexagonal lattice (coordination number 6), as well as 'confined geometries' of 1D lattices and narrow strips. The qualitative behaviour of the hexagonal lattice is similar to the rectangular, with a 2nd and a 1st order transition for increasing CO molar fraction, y CO . Yet, the location of these transitions shifts to lower and higher reactive y CO values; as a result, the reactive window becomes wider for this lattice. More importantly, on the 1D lattice there is no reactive window, and the transition between an O-poisoned state to a CO-poisoned state occurs abruptly. Meakin and Scalapino [205] further carried out simulations on strips, by introducing 'inert' sites along the long sides of rectangular lattices. It was found that for strips of three lattice units and above, the reactive window is restored. Changes in the reactivity of the surface, depending on the combinations of the neighbouring sites onto which the O adatoms can be located upon dissociation, have also been reported by Tammaro and Evans [206] and Cortés and Valencia [207]. These investigations are perhaps more pertinent to the adsorption mechanism than the geometric structure of the surface.
A conceptually similar study to that of [205] was performed by Yaldram and Khan [149] on the NO-CO reaction. The mechanism simulated incorporates adsorption of CO, dissociative adsorption of NO into N and O, CO oxidation by O adatoms towards CO 2 gas, adsorbed NO reduction by a N adatom (leaving O on the surface and producing N 2 gas), and finally N 2 associative desorption. The simulations revealed no kinetic phase transitions (and no reactive window) for a lattice with coordination number 4; whereas for coordination number 6, a reactive state exists. The latter emerges at a critical CO concentration via a 2nd order transition and disappears via a 1st order transition, for a fixed NO concentration. Yaldram and Khan [149] also showed that the reactive window becomes narrower for lower values of the NO dissociation rate, until it disappears for a critical value thereof. More recently, Luque et al [208] investigated the effect of lattice connectivity in the presence of diffusion of N adatoms and lateral interactions (attractive or repulsive between CO, O and N). In agreement with the studies by Yaldram and Khan [149], a reactive window (with respect to CO concentration) was observed only in the six-fold lattice in the absence of interactions and diffusion (note that Luque et al [208] refer to the lattice with coordination number 6 as triangular). Diffusion of N enables the four-fold lattice to exhibit a reactive window, whereas it leaves the reactive window of the six-fold lattice practically unaffected. Repulsive lateral interactions destroy reactivity on both lattices, whereas attractive interactions produce reactive windows on both lattices for an extended range of CO concentrations, still though with the six-fold lattice having a broader window [208].
KMC simulations of the NO-CO catalytic reaction on a series of percolation clusters, ranging from a uniform surface to the incipient percolation cluster fractal, have also been performed by Cortés and Valencia [209][210][211][212]. Interestingly enough, these studies revealed that the rate of CO 2 production correlates linearly with the mean coordination number of the lattice for high gas CO concentrations, whereas for low concentrations linearity is lost and the dependence is less pronounced (figure 8). In the former case, the NO gas concentration is low, making the NO dissociative adsorption and the N and O adatom coverages critical for the catalytic activity. In the latter case, the high availability of N and O adatoms causes the observed insensitivity of CO 2 production rate with respect to catalyst morphology [211]. These results reveal that structure sensitivity (or lack thereof) has to be considered with reference to the operating conditions, a Average Number of Neighbours Activity towards CO conclusion also reached by a more recent KMC study on the water-gas shift reaction [136], discussed later in this section. The aforementioned studies considered relatively simple chemistries and investigated the effect of the geometric structure of the lattice, assuming a single type of active sites. The latter assumption is true only in highly ideal situations concerning model catalysts, whereas real catalysts expose low coordinated as well as defect sites, and they involve synergy with neighbouring support sites. Such sources of complexity attracted the interest early on: variants of the ZGB system have thus been developed, in which special types of sites, mimicking the function of a reducible support, act as an infinite reservoir of O adatoms in the periphery of a metal island [213] or interspersed with the metal phase [214]. In both cases, the existence of these sites was shown to enhance the reactivity of the surface, due to the synergy between O-rich sites and the remaining metal sites. Later, Savchenko et al [215] considered surface morphologies involving patches of two site types that exhibit different sticking coefficients towards CO and O 2 . KMC simulations on this system showed that CO spillover (diffusion between the different patches) enhances the reactivity even higher, in addition to the interfacial synergistic effect shown by the previous studies. Subsequent studies on the same type of reaction mechanism, involving dissociative adsorption of one of the two reactants, generalized these observations for nanometre-sized catalytic particles, emphasizing that the observed kinetics may not necessarily be captured by just superimposing the individual facet kinetics [216,217].
In addition, Vlachos et al [218] considered the case of defect sites onto which adsorbates can bind strongly, for a different reaction mechanism. The latter incorporated reversible adsorption of a species A, and decomposition thereof in the presence of one or two neighbouring vacant sites towards a gas phase species B. Even though in that study no reaction or desorption was allowed on defect sites, the latter constitute nucleation centres in the presence of attractive lateral interactions. As a result, the kinetic phase transition from (almost) clean surfaces to A-rich configurations shifts to lower partial pressures of the reactant species A. Moreover, the range of the observed bistable region decreases for higher defect densities, indicating a sharper transition to the A-rich state [218].
The studies discussed above revealed several interesting effects of the morphology of the catalytic surface on the observed kinetic behaviour thereof. With the advent of ab initio KMC methodologies, new avenues are opening for in-depth comprehensive investigations of the origin and consequences of structure sensitivity effects on catalytic systems.
Thus, Stamatakis et al [136] recently investigated structure sensitivity effects for the water-gas shift reaction (CO+H 2 O ↔ CO 2 +H 2 ) on Pt catalysts using first-principlesbased KMC. The simulated reaction mechanism incorporated reversible adsorption of CO and H 2 O, the decomposition of the latter into OH, O and H, and different pathways for the water-gas shift: direct oxidation of CO via O adatoms, or alternatively, formation and decomposition of carboxyl (COOH) or formate (HCOO), all of these processes followed by H 2 associative desorption. Different operating conditions were also investigated: CO-rich or CO-lean feeds for a range of temperatures on three different surface morphologies. The latter were Pt(1 1 1), a flat surface, Pt(3 2 2) with a step : terrace site ratio of 4 : 1, and Pt(2 1 1) with a step:terrace site ratio of 2 : 1 (figure 9(a)). It was found that CO-rich feeds result in negligible structure sensitivity with the Pt(1 1 1) being slightly more reactive than stepped surfaces ( figure 9(b)). Under these conditions, the model predicts that step sites remain almost poisoned by CO, whereas the water-gas shift proceeds via the carboxyl pathway on terraces or at the step-terrace interface ( figure 9(d)). On the other hand, the stepped surfaces exhibit much higher activities for CO-lean feeds and high temperatures (figure 9(c)), an effect that was attributed to the lower coverages and higher activities on steps, and the dominant effect of the direct CO oxidation pathway thereon (figure 9(d)).
It is worth mentioning that different rates on different exposed surfaces may not necessarily result in a measurable structure sensitivity effect. In particular, Mei et al [219] studied NO oxidation and reduction on Pt(1 0 0) and Pt(1 1 1), both planes exposed by cuboctahedral nanoparticles. Spillover of adsorbates was explicitly taken into account in the KMC simulations, and the edge and corner sites were treated as (1 1 1) sites for simplicity. By comparing the activities for a range of particle sizes, it was shown that NO oxidation appears as structure sensitive while NO reduction does not. In fact, a close inspection of the specific rates of these two chemistries on the two facets reveals that NO reduction does not occur on (1 1 1) facets, whereas NO oxidation rates can take place on both facets with a slightly higher rate on Pt (1 1 1). The apparent insensitivity of NO reduction was attributed to the fraction of (1 0 0) sites staying almost constant as the particle size increased for the sizes/shapes investigated. This study underlines the importance of nanoparticle morphology (shape) in addition to size, in characterizing the structure sensitivity of a catalytic reaction.

Summary and perspective
Computational catalysis is an active and exciting field of research that has already yielded results important in the fundamental understanding and future development of catalytic materials. Motivated by these developments, we briefly discussed the different scales encountered in catalytic systems and explored recent methodological advances in the KMC simulation thereof. In this topical review, we further highlighted some of the interesting physical phenomena that have thus been unravelled, whereas in a previous paper, we had summarized recent findings from the perspective of a chemist [80]. A striking observation is that, in less than two decades, a marked transition was achieved, from simulations of conceptual systems and simple chemistries, to comprehensive studies of realistic kinetics, able to provide semi-quantitative or even quantitative predictions for model catalysts. Still, though, there are major outstanding challenges that need to be overcome in order to capture and analyse the complexity of industrial catalysts [58].
The large numbers of reaction intermediates for realistic chemistries, in conjunction with the structural complexity of supported catalysts, necessitates detailed investigations of a multitude of reaction pathways on the different site types: support, metal, as well as metal-support interface. For catalysts with multicomponent metal or support phases (for instance bimetallics) proposing a full reaction mechanism and parametrizing the model via DFT calculations, can be even more challenging, as the number of elementary events may quickly become intractable. Recent research on selflearning KMC, also referred to as adaptive KMC, which detects possible reaction events 'on the fly' as the simulation proceeds, may prove invaluable for such applications [220][221][222]. Alternatively, this challenge could be overcome by developing software packages that would facilitate (or automate) the workflow of detecting the elementary events a priori, and subsequently incorporate the latter into the KMC simulation. Equally important would be methods for parameterizing the energetics and kinetics models on the basis of descriptors [10], or by using group additivity approaches [223][224][225][226][227][228].
Capturing the dynamic changes of the catalytic surface is also a major challenge, as we have already discussed in section 3.2 on catalyst reconstruction. Changes in the composition of the catalytic entities may take place and markedly affect catalytic activity, for instance in the methanol synthesis chemistry, Zn atoms migrating from the ZnO support onto the edges of Cu nanoparticles have been proposed as the active sites for this chemistry [229]. Moreover, sintering of nanoparticles may lead to the deactivation of the catalyst [230]. To model such changes, one would need to explicitly capture the temporal evolution of the catalytic surface, a task that could be achieved by algorithms combining a submodel for the catalytic surface dynamics and another one for the kinetics thereon (similar, for instance, to the work by Monine et al [169,170]).
A perhaps more fundamental challenge has to do with the quality of the predictions, given that the vast majority of published first-principles KMC studies use DFT calculations for the parametrization of the models. The error of such calculations on estimating the energy of transition states is expected to be on the order of a few tens of an eV [90,91], which is rather high, considering that an error of 0. the observed behaviour and characterize such influences, with the ultimate goal of optimizing the performance of the system. Such a task is achieved through sensitivity analysis methods. Both uncertainty quantification and sensitivity analysis, in the context of lattice KMC simulations, are attracting significant interest in the catalysis and applied mathematics communities [237][238][239], and they may hold the key to obtaining an unprecedented understanding of the underlying catalytic mechanisms in complex systems. Finally, bridging the gaps between the supported metal nanoparticle, the catalytic pellet and ultimately the reactor would be a major achievement in the development of bottomup multiscale modelling frameworks with predictive power. CatalyticFOAM, a CFD software package able to simulate complex catalytic chemistries based on MK models has recently been developed by Maestri and Cuoci [57]. It is desirable to also couple CFD with KMC models, in order to benefit from the higher fidelity and accuracy of KMC. Recent studies have focused on this challenge, developing hybrid algorithms that incorporate mass and heat transfer effects into KMC simulations [55,56,[240][241][242][243]; yet, deviations from the predictions of corresponding deterministic models have been observed [55], as well as violations of mass conservation, which have to be corrected during the course of the simulation [56]. Future development efforts for treating mass transfer could make use of the reaction-diffusion master equation [244] (also known as multivariate master equation [245]), essentially a stochastic finite volume formulation. One could use this approach in combination with surface KMC, or by coupling to both CFD and KMC, if computational expense becomes an issue. Note that this master equation is in principle able to treat diffusion under low pressures or in confined geometries, for which the continuum approximation breaks down. For instance, one would be able to model selective access of desired reactants to active sites within a zeolite pore. The aforementioned advancements are extremely promising and, with the necessary future developments, they are envisioned to eventually enable the purely first-principlesbased catalyst discovery, as well as reactor design and optimization.