Unraveling the Complexity of Catalytic Reactions via Kinetic Monte Carlo Simulation: Current Status and Frontiers

Modified Activation-Relaxation Technique (ARTn) Method Tuned for Efficient Identification of Transition States in Surface Reactions

Exploring potential energy surfaces (PES) is essential for unraveling the underlying mechanisms of chemical reactions and material properties. While the activation-relaxation technique (ARTn) is a state-of-the-art method for identifying saddle points on PES, it often faces challenges in complex energy landscapes, especially on surfaces. In this study, we introduce iso-ARTn, an enhanced ARTn method that incorporates constraints on an orthogonal hyperplane and employs an adaptive active volume. By leveraging a neural network potential (NNP) to conduct an exhaustive saddle point search on the Pt(111) surface with 0.3 monolayers of surface oxygen coverage, iso-ARTn achieves a success rate that is 8.2% higher than the original ARTn, with 40% fewer force calls. Moreover, this method effectively finds various saddle points without compromising the success rate. Combined with kinetic Monte Carlo simulations for event table construction, iso-ARTn with NNP demonstrates the capability to reveal structures consistent with experimental observations. This work signifies a substantial advancement in the investigation of PES, enhancing both the efficiency and breadth of saddle point searches.

Basics of the reaction kinetics on metallic alloys

Kinetics of heterogeneous catalytic reactions are often complicated by various factors, and from the perspective of statistical physics the development of the corresponding models is frequently challenging especially in the case of practically important alloy catalysts. To extend the basics in this area, I use a generic kinetic model of N_{2} formation in the NO reaction with such species as CO or H_{2}. The focus is on the reaction on the surface of a bimetallic alloy with low integral fraction of one of the metals so that the alloy is formed only at the surface. The main goal is to illustrate the specifics of the reaction kinetics and to clarify the accuracy of various approximations which are often inevitable for the analysis of such systems. The key results are as follows. (i) In the baseline one-metal case with the Langmuir-type equations, the model predicts the existence of the optimal adsorbate binding energy corresponding to the maximal reaction rate. (ii) With suitable substitution of the symbols, the equations employed for a uniform surface [item (i)] can be used for the mean-field description of reaction on the surface of a random alloy. In particular, the maximum in the dependence of the reaction rate on the binding energy is converted to the maximum with respect to the alloy composition. (iii) The reaction kinetics on the random-alloy surface have been described exactly. The corresponding maximum in the reaction rate is lower and somewhat smeared compared to that predicted in the mean-field approximation for an alloy or for the optimal monometallic catalyst. (iv) The reaction kinetics have been scrutinized also in the case of two-dimensional (2D) segregation of metal atoms at the surface in the limits of slow and rapid adsorbate diffusion between the spots. The kinetics are shown to be qualitatively different in these limits and also compared to the random-alloy case. (v) The thermodynamic criteria for 2D segregation of metal atoms at the surface have been derived as well.

Adaptive tau-leaping methods for microscopic-lattice kinetic Monte Carlo simulations

Traditional Kinetic Monte Carlo (KMC) approaches, rooted in Gillespie's stochastic simulation algorithm, become computationally demanding in systems with a large range of timescales. The goal of this work is to propose and study new adaptive lattice-KMC time integration strategies for spatially non-uniform systems. To that end, two novel adaptive tau-leaping methods and their corresponding time integration strategies are developed based on the idea of the "n-fold" direct KMC method. These strategies allow for the simultaneous execution of multiple reactions, advancing time by adaptively selected coarse increments. We present numerical experiments comparing the proposed methods with existing approaches in a catalytic surface kinetics application involving ammonia decomposition.

Kinetic Monte Carlo Analysis Reveals Non-mean-field Active Site Dynamics in Cu-Zeolite-Catalyzed NO x Reduction

Copper-exchanged chabazite (Cu-CHA) zeolites are the preferred catalysts for the selective catalytic reduction of NO x with NH3. The low temperature (473 K) SCR mechanism proceeds through a redox cycle between mobile and ammonia-solvated Cu(I) and Cu(II) complexes, as demonstrated by multiple experimental and computational investigations. The oxidation step requires two Cu(I) to migrate into the same cha cage to activate O2 and form a binuclear Cu(II)-di-oxo complex. Prior steady state and transient kinetic experiments find that the apparent rate constants for oxidation (per Cu ion) are sensitive to catalyst composition and follow nonmean-field kinetics. We develop a nonmean-field kinetic model for NO x SCR that incorporates a composition-dependent Cu(I) volumetric footprint centered at anionic [AlO4]- tetrahedral sites on the CHA lattice. We use Bayesian optimization to parameterize a kinetic Monte Carlo model against available experimental composition-dependent SCR rates and in situ Cu(II) fractions. We find that both rates and Cu(II) fractions of a majority of catalyst compositions can be captured by single oxidation and reduction rate constants combined with a composition-dependent Cu(I) cation footprint, highlighting the contributions of both Cu and Al densities to steady-state SCR performance of Cu-CHA. The work illustrates a pathway for extracting robust molecular insights from the kinetics of a dynamic catalytic system.

On the Capabilities of Transition Metal Carbides for Carbon Capture and Utilization Technologies

The search for cheap and active materials for the capture and activation of CO2 has led to many efforts aimed at developing new catalysts. In this context, earth-abundant transition metal carbides (TMCs) have emerged as promising candidates, garnering increased attention in recent decades due to their exceptional refractory properties and resistance to sintering, coking, and sulfur poisoning. In this work, we assess the use of Group 5 TMCs (VC, NbC, and TaC) as potential materials for carbon capture and sequestration/utilization technologies by combining experimental characterization techniques, first-principles-based multiscale modeling, vibrational analysis, and catalytic experiments. Our findings reveal that the stoichiometric phase of VC exhibits weak interactions with CO2, displaying an inability to adsorb or dissociate it. However, VC often exhibits the presence of surface carbon vacancies, leading to significant activation of CO2 at room temperature and facilitating its catalytic hydrogenation. In contrast, stoichiometric NbC and TaC phases exhibit stronger interactions with CO2, capable of adsorbing and even breaking of CO2 at low temperatures, particularly notable in the case of TaC. Nevertheless, NbC and TaC demonstrate poor catalytic performance for CO2 hydrogenation. This work suggests Group 5 TMCs as potential materials for CO2 abatement, emphasizes the importance of surface vacancies in enhancing catalytic activity and adsorption capability, and provides a reference for using the infrared spectra as a unique identifier to detect oxy-carbide phases or surface C vacancies within Group 5 TMCs.

Dynamic Active Sites In Situ Formed in Metal Nanoparticle Reshaping under Reaction Conditions

Understanding the nonequilibrium transformation of nanocatalysts under reaction conditions is important because metastable atomic structures may be created during the process, which offers unique activities in reactions. Although reshaping of metal nanoparticles (NPs) under reaction conditions has been widely recognized, the dynamic reshaping process has been less studied at the atomic scale. Here, we develop an atomistic kinetic Monte Carlo model to simulate the complete reshaping process of Pt nanoparticles in a CO environment and reveal the in situ formation of atomic clusters on the NP surface, a new type of active site beyond conventional understanding, boosting the reactivities in the CO oxidation reaction. Interestingly, highly active peninsula and inactive island clusters both form on the (111) facets and interchange in varying states of dynamic equilibrium, which influences the catalytic activities significantly. This study provides new fundamental knowledge of nanocatalysis and new guidance for the rational design of nanocatalysts.

Robust Gaussian Process Regression Method for Efficient Tunneling Pathway Optimization: Application to Surface Processes

Simulation of surface processes is a key part of computational chemistry that offers atomic-scale insights into mechanisms of heterogeneous catalysis, diffusion dynamics, and quantum tunneling phenomena. The most common theoretical approaches involve optimization of reaction pathways, including semiclassical tunneling pathways (called instantons). The computational effort can be demanding, especially for instanton optimizations with an ab initio electronic structure. Recently, machine learning has been applied to accelerate reaction-pathway optimization, showing great potential for a wide range of applications. However, previous methods still suffer from numerical and efficiency issues and were not designed for condensed-phase reactions. We propose an improved framework based on Gaussian process regression for general transformed coordinates, which has improved efficiency and numerical stability, and we propose a descriptor that combines internal and Cartesian coordinates suitable for modeling surface processes. We demonstrate with 11 instanton optimizations in three representative systems that the improved approach makes ab initio instanton optimization significantly cheaper, such that it becomes not much more expensive than a classical transition-state theory rate calculation.

Comprehensive Density Functional and Kinetic Monte Carlo Study of CO2 Hydrogenation on a Well-Defined Ni/CeO2 Model Catalyst: Role of Eley-Rideal Reactions

A detailed multiscale study of the mechanism of CO2 hydrogenation on a well-defined Ni/CeO2 model catalyst is reported that couples periodic density functional theory (DFT) calculations with kinetic Monte Carlo (kMC) simulations. The study includes an analysis of the role of Eley-Rideal elementary steps for the water formation step, which are usually neglected on the overall picture of the mechanism, catalytic activity, and selectivity. The DFT calculations for the chosen model consisting of a Ni4 cluster supported on CeO2 (111) show large enough adsorption energies along with low energy barriers that suggest this catalyst to be a good option for high selective CO2 methanation. The kMC simulations results show a synergic effect between the two 3-fold hollow sites of the supported Ni4 cluster with some elementary reactions dominant in one site, while other reactions prefer the another, nearly equivalent site. This effect is even more evident for the simulations explicitly including Eley-Rideal steps. The kMC simulations reveal that CO is formed via the dissociative pathway of the reverse water-gas shift reaction, while methane is formed via a CO2 → CO → HCO → CH → CH2 → CH3 → CH4 mechanism. Overall, our results show the importance of including the Eley-Rideal reactions and point to small Ni clusters supported on the CeO2 (111) surface as potential good catalysts for high selective CO2 methanation under mild conditions, while very active and selective toward CO formation at higher temperatures.

Combined DFT and Kinetic Monte Carlo Study of UiO-66 Catalysts for γ-Valerolactone Production

Zr-based metal-organic frameworks (MOFs) are excellent heterogeneous porous catalysts due to their thermal stability. Their tunability via node and linker modifications makes them amenable for theoretical studies on catalyst design. However, detailed benchmarks on MOF-based reaction mechanisms combined with kinetics analysis are still scarce. Thus, we here evaluate different computational models and density functional theory (DFT) methods followed by kinetic Monte Carlo studies for a case reaction relevant in biomass upgrading, i.e., the conversion of methyl levulinate to γ-valerolactone catalyzed by UiO-66. We show the impact of cluster versus periodic models, the importance of the DF of choice, and the direct comparison to experimental data via simulated kinetics data. Overall, we found that Perdew-Burke-Ernzerhof (PBE), a widely employed method in plane-wave periodic calculations, greatly overestimates reaction rates, while M06 with cluster models better fits the available experimental data and is recommended whenever possible.

Speeding up the Detection of Adsorbate Lateral Interactions in Graph-Theoretical Kinetic Monte Carlo Simulations

Kinetic Monte Carlo (KMC) has become an indispensable tool in heterogeneous catalyst discovery, but realistic simulations remain computationally demanding on account of the need to capture complex and long-range lateral interactions between adsorbates. The Zacros software package (https://zacros.org) adopts a graph-theoretical cluster expansion (CE) framework that allows such interactions to be computed with a high degree of generality and fidelity. This involves solving a series of subgraph isomorphism problems in order to identify relevant interaction patterns in the lattice. In an effort to reduce the computational burden, we have adapted two well-known subgraph isomorphism algorithms, namely, VF2 and RI, for use in KMC simulations and implemented them in Zacros. To benchmark their performance, we simulate a previously established model of catalytic NO oxidation, treating the O* lateral interactions with a series of progressively larger CEs. For CEs with long-range interactions, VF2 and RI are found to provide impressive speedups relative to simpler algorithms. RI performs best, giving speedups reaching more than 150× when combined with OpenMP parallelization. We also simulate a recently developed methane cracking model, showing that RI offers significant improvements in performance at high surface coverages.

Surface coverage dynamics for reversible dissociative adsorption on finite linear lattices

Dissociative adsorption onto a surface introduces dynamic correlations between neighboring sites not found in non-dissociative absorption. We study surface coverage dynamics where reversible dissociative adsorption of dimers occurs on a finite linear lattice. We derive analytic expressions for the equilibrium surface coverage as a function of the number of reactive sites, N, and the ratio of the adsorption and desorption rates. Using these results, we characterize the finite size effect on the equilibrium surface coverage. For comparable N's, the finite size effect is significantly larger when N is even than when N is odd. Moreover, as N increases, the size effect decays more slowly in the even case than in the odd case. The finite-size effect becomes significant when adsorption and desorption rates are considerably different. These finite-size effects are related to the number of accessible configurations in a finite system where the odd-even dependence arises from the limited number of accessible configurations in the even case. We confirm our analytical results with kinetic Monte Carlo simulations. We also analyze the surface-diffusion case where adsorbed atoms can hop into neighboring sites. As expected, the odd-even dependence disappears because more configurations are accessible in the even case due to surface diffusion.

Operando Modeling of Zeolite-Catalyzed Reactions Using First-Principles Molecular Dynamics Simulations

Within this Perspective, we critically reflect on the role of first-principles molecular dynamics (MD) simulations in unraveling the catalytic function within zeolites under operating conditions. First-principles MD simulations refer to methods where the dynamics of the nuclei is followed in time by integrating the Newtonian equations of motion on a potential energy surface that is determined by solving the quantum-mechanical many-body problem for the electrons. Catalytic solids used in industrial applications show an intriguing high degree of complexity, with phenomena taking place at a broad range of length and time scales. Additionally, the state and function of a catalyst critically depend on the operating conditions, such as temperature, moisture, presence of water, etc. Herein we show by means of a series of exemplary cases how first-principles MD simulations are instrumental to unravel the catalyst complexity at the molecular scale. Examples show how the nature of reactive species at higher catalytic temperatures may drastically change compared to species at lower temperatures and how the nature of active sites may dynamically change upon exposure to water. To simulate rare events, first-principles MD simulations need to be used in combination with enhanced sampling techniques to efficiently sample low-probability regions of phase space. Using these techniques, it is shown how competitive pathways at operating conditions can be discovered and how broad transition state regions can be explored. Interestingly, such simulations can also be used to study hindered diffusion under operating conditions. The cases shown clearly illustrate how first-principles MD simulations reveal insights into the catalytic function at operating conditions, which could not be discovered using static or local approaches where only a few points are considered on the potential energy surface (PES). Despite these advantages, some major hurdles still exist to fully integrate first-principles MD methods in a standard computational catalytic workflow or to use the output of MD simulations as input for multiple length/time scale methods that aim to bridge to the reactor scale. First of all, methods are needed that allow us to evaluate the interatomic forces with quantum-mechanical accuracy, albeit at a much lower computational cost compared to currently used density functional theory (DFT) methods. The use of DFT limits the currently attainable length/time scales to hundreds of picoseconds and a few nanometers, which are much smaller than realistic catalyst particle dimensions and time scales encountered in the catalysis process. One solution could be to construct machine learning potentials (MLPs), where a numerical potential is derived from underlying quantum-mechanical data, which could be used in subsequent MD simulations. As such, much longer length and time scales could be reached; however, quite some research is still necessary to construct MLPs for the complex systems encountered in industrially used catalysts. Second, most currently used enhanced sampling techniques in catalysis make use of collective variables (CVs), which are mostly determined based on chemical intuition. To explore complex reactive networks with MD simulations, methods are needed that allow the automatic discovery of CVs or methods that do not rely on a priori definition of CVs. Recently, various data-driven methods have been proposed, which could be explored for complex catalytic systems. Lastly, first-principles MD methods are currently mostly used to investigate local reactive events. We hope that with the rise of data-driven methods and more efficient methods to describe the PES, first-principles MD methods will in the future also be able to describe longer length/time scale processes in catalysis. This might lead to a consistent dynamic description of all steps-diffusion, adsorption, and reaction-as they take place at the catalyst particle level.

Formation of active sites on transition metals through reaction-driven migration of surface atoms

Adopting low-index single-crystal surfaces as models for metal nanoparticle catalysts has been questioned by the experimental findings of adsorbate-induced formation of subnanometer clusters on several single-crystal surfaces. We used density functional theory calculations to elucidate the conditions that lead to cluster formation and show how adatom formation energies enable efficient screening of the conditions required for adsorbate-induced cluster formation. We studied a combination of eight face-centered cubic transition metals and 18 common surface intermediates and identified systems relevant to catalytic reactions, such as carbon monoxide (CO) oxidation and ammonia (NH₃) oxidation. We used kinetic Monte Carlo simulations to elucidate the CO-induced cluster formation process on a copper surface. Scanning tunneling microscopy of CO on a nickel (111) surface that contains steps and dislocations points to the structure sensitivity of this phenomenon. Metal-metal bond breaking that leads to the evolution of catalyst structures under realistic reaction conditions occurs much more broadly than previously thought.

Ethane Dehydrogenation over the Core-Shell Pt-Based Alloy Catalysts: Driven by Engineering the Shell Composition and Thickness

Pt-based catalysts as the commercial catalysts in ethane dehydrogenation (EDH) face one of the main challenges of realizing the balance between coke formation and catalytic activity. In this work, a strategy to drive the catalytic performance of EDH on Pt-Sn alloy catalysts is proposed by rationally engineering the shell surface structure and thickness of core-shell Pt@Pt₃Sn and Pt₃Sn@Pt catalysts from a theoretical perspective. Eight types of Pt@Pt₃Sn and Pt₃Sn@Pt catalysts with different Pt and Pt₃Sn shell thicknesses are considered and compared with the industrially used Pt and Pt₃Sn catalysts. Density functional theory (DFT) calculations completely describe the reaction network of EDH, including the side reactions of deep dehydrogenation and C-C bond cracking. Kinetic Monte Carlo (kMC) simulations reveal the influences of the catalyst surface structure, experimentally related temperatures, and reactant partial pressures. The results show that CHCH* is the main precursor for coke formation, and Pt@Pt₃Sn catalysts generally have higher C₂H₄(g) activity and lower selectivity compared to those of Pt₃Sn@Pt catalysts, which is attributed to the unique surface geometrical and electronic properties. 1Pt₃Sn@4Pt and 1Pt@4Pt₃Sn are screened out as catalysts exhibiting excellent performance; especially, the 1Pt₃Sn@4Pt catalyst has much higher C₂H₄(g) activity and 100% C₂H₄(g) selectivity compared to those of 1Pt@4Pt₃Sn and the widely used Pt and Pt₃Sn catalysts. The two descriptors C₂H₅* adsorption energy and reaction energy of its dehydrogenation to C₂H₄* are proposed to qualitatively evaluate the C₂H₄(g) selectivity and activity, respectively. This work facilitates a valuable exploration for optimizing the catalytic performance of core-shell Pt-based catalysts in EDH and reveals the great importance of the fine control of the catalyst shell surface structure and thickness.

Nature of Catalytic Behavior of Cobalt Oxides for CO2 Hydrogenation

Cobalt oxide (CoO x ) catalysts are widely applied in CO2 hydrogenation but suffer from structural evolution during the reaction. This paper describes the complicated structure-performance relationship under reaction conditions. An iterative approach was employed to simulate the reduction process with the help of neural network potential-accelerated molecular dynamics. Based on the reduced models of catalysts, a combined theoretical and experimental study has discovered that CoO(111) provides active sites to break C-O bonds for CH4 production. The analysis of the reaction mechanism indicated that the C-O bond scission of *CH2O species plays a key role in producing CH4. The nature of dissociating C-O bonds is attributed to the stabilization of *O atoms after C-O bond cleavage and the weakening of C-O bond strength by surface-transferred electrons. This work may offer a paradigm to explore the origin of performance over metal oxides in heterogeneous catalysis.

Large-scale benchmarks of the time-warp/graph-theoretical kinetic Monte Carlo approach for distributed on-lattice simulations of catalytic kinetics

Motivated by the need to perform large-scale kinetic Monte Carlo (KMC) simulations, in the context of unravelling complex phenomena such as catalyst reconstruction and pattern formation, we extend the work of Ravipati et al. [S. Ravipati, G. D. Savva, I.-A. Christidi, R. Guichard, J. Nielsen, R. Réocreux and M. Stamatakis, Comput. Phys. Commun., 2022, 270, 108148] in benchmarking the performance of a distributed-computing, on-lattice KMC approach. The latter, implemented in our software package Zacros, combines the graph-theoretical KMC framework with the Time-Warp algorithm for parallel discrete event simulations, and entails dividing the lattice into subdomains, each assigned to a processor. The cornerstone of the Time-Warp algorithm is the state queue, to which snapshots of the simulation state are saved regularly, enabling historical KMC information to be corrected when conflicts occur at subdomain boundaries. Focusing on three model systems, we highlight the key Time-Warp parameters that can be tuned to optimise performance. The frequency of state saving, controlled by the state saving interval, δ_snap, is shown to have the largest effect on performance, which favours balancing the overhead of re-simulating KMC history with that of writing state snapshots to memory. Also important is the global virtual time (GVT) computation interval, Δτ_GVT, which has little direct effect on the progress of the simulation but controls how often the state queue memory can be freed up. We also find that pre-allocating memory for the state queue data structure favours performance. These findings will guide users in maximising the efficiency of Zacros or other distributed KMC software, which is a vital step towards realising accurate, meso-scale simulations of heterogeneous catalysis.

Enhancing Catalytic Performance through Subsurface Chemistry: The Case of C2H2 Semihydrogenation over Pd Catalysts

Subsurface chemistry in heterogeneous catalysis plays an important role in tuning catalytic performance. Aiming to unravel the role of subsurface heteroatoms, C₂H₂ semihydrogenation on a series of Pd catalysts doped with subsurface heteroatom H, B, C, N, P, or S was fully investigated by density functional theory (DFT) calculations together with microkinetic modeling. The obtained results showed that catalytic performance toward C₂H₂ semihydrogenation was affected significantly by the type and coverage of subsurface heteroatoms. The Pd-B_0.5 and Pd-C_0.5 catalysts with 1/2 monolayer (ML) heteroatom coverage, as well as Pd-N, Pd-P, and Pd-S catalysts with 1/16 ML heteroatom coverage, were screened to not only obviously improve C₂H₄ selectivity and activity but also effectively suppress green oil. The essential reason for subsurface heteroatoms in tuning catalytic performance is attributed to the distinctive surface Pd electronic and geometric structures caused by subsurface heteroatoms. In the Pd-B_0.5 and Pd-C_0.5 catalysts, the Pd surface electronic and geometric effects play the dominant role, while the geometric effect plays a key role in the Pd-N, Pd-P, and Pd-S catalysts. The findings provide theoretically valuable information for designing high-performance metal catalysts in alkyne semihydrogenation through subsurface chemistry.

Nonlinear Reactor Design Optimization With Embedded Microkinetic Model Information

Despite the success of multiscale modeling in science and engineering, embedding molecular-level information into nonlinear reactor design and control optimization problems remains challenging. In this work, we propose a computationally tractable scale-bridging approach that incorporates information from multi-product microkinetic (MK) models with thousands of rates and chemical species into nonlinear reactor design optimization problems. We demonstrate reduced-order kinetic (ROK) modeling approaches for catalytic oligomerization in shale gas processing. We assemble a library of six candidate ROK models based on literature and MK model structure. We find that three metrics—quality of fit (e.g., mean squared logarithmic error), thermodynamic consistency (e.g., low conversion of exothermic reactions at high temperatures), and model identifiability—are all necessary to train and select ROK models. The ROK models that closely mimic the structure of the MK model offer the best compromise to emulate the product distribution. Using the four best ROK models, we optimize the temperature profiles in staged reactors to maximize conversions to heavier oligomerization products. The optimal temperature starts at 630–900K and monotonically decreases to approximately 560 K in the final stage, depending on the choice of ROK model. For all models, staging increases heavier olefin production by 2.5% and there is minimal benefit to more than four stages. The choice of ROK model, i.e., model-form uncertainty, results in a 22% difference in the objective function, which is twice the impact of parametric uncertainty; we demonstrate sequential eigendecomposition of the Fisher information matrix to identify and fix sloppy model parameters, which allows for more reliable estimation of the covariance of the identifiable calibrated model parameters. First-order uncertainty propagation determines this parametric uncertainty induces less than a 10% variability in the reactor optimization objective function. This result highlights the importance of quantifying model-form uncertainty, in addition to parametric uncertainty, in multi-scale reactor and process design and optimization. Moreover, the fast dynamic optimization solution times suggest the ROK strategy is suitable for incorporating molecular information in sequential modular or equation-oriented process simulation and optimization frameworks.

Challenges of modeling nanostructured materials for photocatalytic water splitting

Understanding the water splitting mechanism in photocatalysis is a rewarding goal as it will allow producing clean fuel for a sustainable life in the future. However, identifying the photocatalytic mechanisms by modeling photoactive nanoparticles requires sophisticated computational techniques based on multiscale modeling. In this review, we will survey the strengths and drawbacks of currently available theoretical methods at different length and accuracy scales. Understanding the surface-active site through Density Functional Theory (DFT) using new, more accurate exchange-correlation functionals plays a key role for surface engineering. Larger scale dynamics of the catalyst/electrolyte interface can be treated with Molecular Dynamics albeit there is a need for more generalizations of force fields. Monte Carlo and Continuum Modeling techniques are so far not the prominent path for modeling water splitting but interest is growing due to the lower computational cost and the feasibility to compare the modeling outcome directly to experimental data. The future challenges in modeling complex nano-photocatalysts involve combining different methods in a hierarchical way so that resources are spent wisely at each length scale, as well as accounting for excited states chemistry that is important for photocatalysis, a path that will bring devices closer to the theoretical limit of photocatalytic efficiency.

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

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

Autonomous Reaction Network Exploration in Homogeneous and Heterogeneous Catalysis

Autonomous computations that rely on automated reaction network elucidation algorithms may pave the way to make computational catalysis on a par with experimental research in the field. Several advantages of this approach are key to catalysis: (i) automation allows one to consider orders of magnitude more structures in a systematic and open-ended fashion than what would be accessible by manual inspection. Eventually, full resolution in terms of structural varieties and conformations as well as with respect to the type and number of potentially important elementary reaction steps (including decomposition reactions that determine turnover numbers) may be achieved. (ii) Fast electronic structure methods with uncertainty quantification warrant high efficiency and reliability in order to not only deliver results quickly, but also to allow for predictive work. (iii) A high degree of autonomy reduces the amount of manual human work, processing errors, and human bias. Although being inherently unbiased, it is still steerable with respect to specific regions of an emerging network and with respect to the addition of new reactant species. This allows for a high fidelity of the formalization of some catalytic process and for surprising in silico discoveries. In this work, we first review the state of the art in computational catalysis to embed autonomous explorations into the general field from which it draws its ingredients. We then elaborate on the specific conceptual issues that arise in the context of autonomous computational procedures, some of which we discuss at an example catalytic system.

GRAPHICAL ABSTRACT

SUPPLEMENTARY INFORMATION

The online version contains supplementary material available at 10.1007/s11244-021-01543-9.

Ab initio molecular dynamics with enhanced sampling in heterogeneous catalysis

Enhanced sampling ab initio simulations enable to study chemical phenomena in catalytic systems including thermal effects & anharmonicity, & collective dynamics describing enthalpic & entropic contributions, which can significantly impact on reaction free energy landscapes.

Kinetic Monte Carlo simulations of the dry reforming of methane catalyzed by the Ru (0001) surface based on density functional theory calculations

This study shows the main pathways for the DRM reaction and the competitive RWGS reaction upon changing reaction conditions, displaying the importance of including lateral–lateral interactions to describe the reaction in agreement with the experiment.

Structure-Reactivity Relationship for Nano-Catalysts in the Hydrogenation/Dehydrogenation Controlled Reaction Systems

For the activity of a nano-catalyst, a general and quantitative solution to building direct structure-reactivity relationship has not yet been established. On top of the first-principle-based kinetic Monte Carlo (KMC) simulations, we developed a model to build the adsorption site dependence of the activity. We applied this model to study the nano effects of Cu catalysts in the water-gas shift reaction. By accumulating the activities of different adsorption sites, our model satisfactorily reproduced the experimental apparent activation energies for catalysts with sizes over hundreds of nanometers, which were out of reach for conventional KMC simulations. Our results disclose that, even for a cubic catalyst with size of 877 nm, its activity can still be closely related to the activity of edge sites, instead of only the exposed Cu(100) facets as might be expected. The present model is expected to be useful for systems that are controlled by the hydrogenation/dehydrogenation processes.

Effects of surface diffusion in electrocatalytic CO2 reduction on Cu revealed by kinetic Monte Carlo simulations

Kinetic Monte Carlo (KMC) methods are frequently used for mechanistic studies of thermally driven heterogeneous catalysis systems but are underused for electrocatalysis. Here, we develop a lattice KMC approach for electrocatalytic CO₂ reduction. The work is motivated by a prior experimental report that performed electroreduction of a mixed feed of ¹²CO₂ and ¹³CO on Cu; differences in the ¹³C content of C2 products ethylene and ethanol (Δ¹³C) were interpreted as evidence of site selectivity. The lattice KMC model considers the effect of surface diffusion on this system. In the limit of infinitely fast diffusion (mean-field approximation), the key intermediates ¹²CO* and ¹³CO* would be well mixed on the surface and no evidence of site selectivity could have been observed. Using a simple two-site model and adapting a previously reported microkinetic model, we assess the effects of diffusion on the relative isotope fractions in the products using the estimated surface diffusion rate of CO* from literature reports. We find that the size of the active sites and the total surface adsorbate coverage can have a large influence on the values of Δ¹³C that can be observed. Δ¹³C is less sensitive to the CO* diffusion rate as long as it is within the estimated range. We further offer possible methods to estimate surface distribution of intermediates and to predict intrinsic selectivity of active sites based on experimental observations. This work illustrates the importance of considering surface diffusion in the study of electrochemical CO₂ reduction to multi-carbon products. Our approach is entirely based on a freely available open-source code, so will be readily adaptable to other electrocatalytic systems.

Real-time dynamics and structures of supported subnanometer catalysts via multiscale simulations

Understanding the performance of subnanometer catalysts and how catalyst treatment and exposure to spectroscopic probe molecules change the structure requires accurate structure determination under working conditions. Experiments lack simultaneous temporal and spatial resolution and could alter the structure, and similar challenges hinder first-principles calculations from answering these questions. Here, we introduce a multiscale modeling framework to follow the evolution of subnanometer clusters at experimentally relevant time scales. We demonstrate its feasibility on Pd adsorbed on CeO2(111) at various catalyst loadings, temperatures, and exposures to CO. We show that sintering occurs in seconds even at room temperature and is mainly driven by free energy reduction. It leads to a kinetically (far from equilibrium) frozen ensemble of quasi-two-dimensional structures that CO chemisorption and infrared experiments probe. CO adsorption makes structures flatter and smaller. High temperatures drive very rapid sintering toward larger, stable/metastable equilibrium structures, where CO induces secondary structure changes only.

The role of oxygenated species in the catalytic self-coupling of MeOH on O pre-covered Au(111)

The oxidation of alcohols plays a central role in the valorisation of biomass, in particular when performed with a non-toxic oxidant such as O₂. Aerobic oxidation of methanol on gold has attracted attention lately and the main steps of its mechanism have been described experimentally. However, the exact role of O and OH on each elementary step and the effect of the interactions between adsorbates are still not completely understood. Here we investigate the mechanism of methanol oxidation to HCOOCH₃ and CO₂. We use Density Functional Theory (DFT) to assess the energetics of the underlying pathways, and subsequently build lattice kinetic Monte Carlo (kMC) models of increasing complexity, to elucidate the role of different oxygenates. Detailed comparisons of our simulation results with experimental temperature programmed desorption (TPD) spectra enable us to validate the mechanism and identify rate determining steps. Crucially, taking into account dispersion (van der Waals forces) and adsorbate-adsorbate lateral interactions are both important for reproducing the experimental data.

Bragg Coherent Diffraction Imaging for In Situ Studies in Electrocatalysis

Electrocatalysis is at the heart of a broad range of physicochemical applications that play an important role in the present and future of a sustainable economy. Among the myriad of different electrocatalysts used in this field, nanomaterials are of ubiquitous importance. An increased surface area/volume ratio compared to bulk makes nanoscale catalysts the preferred choice to perform electrocatalytic reactions. Bragg coherent diffraction imaging (BCDI) was introduced in 2006 and since has been applied to obtain 3D images of crystalline nanomaterials. BCDI provides information about the displacement field, which is directly related to strain. Lattice strain in the catalysts impacts their electronic configuration and, consequently, their binding energy with reaction intermediates. Even though there have been significant improvements since its birth, the fact that the experiments can only be performed at synchrotron facilities and its relatively low resolution to date (∼10 nm spatial resolution) have prevented the popularization of this technique. Herein, we will briefly describe the fundamentals of the technique, including the electrocatalysis relevant information that we can extract from it. Subsequently, we review some of the computational experiments that complement the BCDI data for enhanced information extraction and improved understanding of the underlying nanoscale electrocatalytic processes. We next highlight success stories of BCDI applied to different electrochemical systems and in heterogeneous catalysis to show how the technique can contribute to future studies in electrocatalysis. Finally, we outline current challenges in spatiotemporal resolution limits of BCDI and provide our perspectives on recent developments in synchrotron facilities as well as the role of machine learning and artificial intelligence in addressing them.

Dynamic Interconversion of Metal Active Site Ensembles in Zeolite Catalysis

Catalysis science is founded on understanding the structure, number, and reactivity of active sites. Kinetic models that consider active sites to be static and noninteracting entities are routinely successful in describing the behavior of heterogeneous catalysts. Yet, active site ensembles often restructure in response to their external environment and even during steady-state catalytic turnover, sometimes requiring non-mean-field kinetic treatments to describe distance-dependent interactions among sites. Such behavior is being recognized more frequently in modern catalysis research, with the advent of experimental methods to quantify turnover rates with increasing precision, an expanding arsenal of operando characterization tools, and computational descriptions of atomic structure and motion at chemical potentials and timescales increasingly relevant to reaction conditions. This review focuses on dynamic changes to metal active site ensembles on zeolite supports, which are silica-based crystalline materials substituted with Al that generate binding sites for isolated and low-nuclearity metal site ensembles. Metal sites can become solvated and mobilized during reaction, facilitating interactions among sites that change their nuclearity and function. Such intersite communication can be regulated by the zeolite support, resulting in non-single-site and potentially non-mean-field kinetic behavior arising from mechanisms of catalytic action that combine elements of those canonically associated with homogeneous and heterogeneous catalysis.We discuss recent literature examples that document dynamic active site behavior in metal-zeolites and outline methodologies to identify and interpret such behavior. We conclude with our outlook on future research directions to develop this evolving branch of catalysis science and harness it for practical applications.

Directing reaction pathways via in situ control of active site geometries in PdAu single-atom alloy catalysts

The atomic scale structure of the active sites in heterogeneous catalysts is central to their reactivity and selectivity. Therefore, understanding active site stability and evolution under different reaction conditions is key to the design of efficient and robust catalysts. Herein we describe theoretical calculations which predict that carbon monoxide can be used to stabilize different active site geometries in bimetallic alloys and then demonstrate experimentally that the same PdAu bimetallic catalyst can be transitioned between a single-atom alloy and a Pd cluster phase. Each state of the catalyst exhibits distinct selectivity for the dehydrogenation of ethanol reaction with the single-atom alloy phase exhibiting high selectivity to acetaldehyde and hydrogen versus a range of products from Pd clusters. First-principles based Monte Carlo calculations explain the origin of this active site ensemble size tuning effect, and this work serves as a demonstration of what should be a general phenomenon that enables in situ control over catalyst selectivity.

Single-Atom Alloys for the Electrochemical Oxygen Reduction Reaction

Single-atom alloys (SAAs) consisting of isolated transition-metal atoms doped in the surface of coinage metal hosts exhibit unique catalytic properties, harnessing the high activity of the dopant metals with the selectivity of the coinage metal hosts. Here we use density functional theory (DFT) to study SAAs comprised of Ni, Pd, Pt, Co and Rh doped into Ag and Au hosts, as candidate electrocatalysts for the oxygen reduction reaction (ORR) in proton-exchange membrane (PEM) fuel-cells. Our calculations reveal that the PdAu SAA exhibits a slightly lower theoretical overpotential, enhanced selectivity for 4-e- ORR, and tolerance to CO-poisoning compared to Pt(111). While the number of active sites of PdAu SAA is lower than that of Pt(111), the aforementioned desirable properties could bring the overall catalytic performance thereof close to that of Pt/C, indicating that the PdAu SAA could be a viable material for electrocatalytic ORR in PEM fuel-cells.