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

On-the-fly kinetic Monte Carlo simulations with neural network potentials for surface diffusion and reaction

We develop an adaptive scheme in the kinetic Monte Carlo simulations, where the adsorption and activation energies of all elementary steps, including the effects of other adsorbates, are evaluated "on-the-fly" by employing the neural network potentials. The configurations and energies evaluated during the simulations are stored for reuse when the same configurations are sampled in a later step. The present scheme is applied to hydrogen adsorption and diffusion on the Pd(111) and Pt(111) surfaces and the CO oxidation reaction on the Pt(111) surface. The effects of interactions between adsorbates, i.e., adsorbate-adsorbate lateral interactions, are examined in detail by comparing the simulations without considering lateral interactions. This study demonstrates the importance of lateral interactions in surface diffusion and reactions and the potential of our scheme for applications in a wide variety of heterogeneous catalytic reactions.

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.

Modeling Single-Atom Catalysis

Electronic structure calculations represent an essential complement of experiments to characterize single-atom catalysts (SACs), consisting of isolated metal atoms stabilized on a support, but also to predict new catalysts. However, simulating SACs with quantum chemistry approaches is not as simple as often assumed. In this work, the essential factors that characterize a reliable simulation of SACs activity are examined. The Perspective focuses on the importance of precise atomistic characterization of the active site, since even small changes in the metal atom's surroundings can result in large changes in reactivity. The dynamical behavior and stability of SACs under working conditions, as well as the importance of adopting appropriate methods to solve the Schrödinger equation for a quantitative evaluation of reaction energies are addressed. The Perspective also focuses on the relevance of the model adopted. For electrocatalysis this must include the effects of the solvent, the presence of electrolytes, the pH, and the external potential. Finally, it is discussed how the similarities between SACs and coordination compounds may result in reaction intermediates that usually are not observed on metal electrodes. When these aspects are not adequately considered, the predictive power of electronic structure calculations is quite limited.

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.

Deep Reinforcement Learning Environment Approach Based on Nanocatalyst XAS Diagnostics Graphic Formalization

The most in-demand instrumental methods for new functional nanomaterial diagnostics employ synchrotron radiation, which is used to determine a material's electronic and local atomic structure. The high time and resource costs of researching at international synchrotron radiation centers and the problems involved in developing an optimal strategy and in planning the control of the experiments are acute. One possible approach to solving these problems involves the use of deep reinforcement learning agents. However, this approach requires the creation of a special environment that provides a reliable level of response to the agent's actions. As the physical experimental environment of nanocatalyst diagnostics is potentially a complex multiscale system, there are no unified comprehensive representations that formalize the structure and states as a single digital model. This study proposes an approach based on the decomposition of the experimental system into the original physically plausible nodes, with subsequent merging and optimization as a metagraphic representation with which to model the complex multiscale physicochemical environments. The advantage of this approach is the possibility to directly use the numerical model to predict the system states and to optimize the experimental conditions and parameters. Additionally, the obtained model can form the basic planning principles and allow for the optimization of the search for the optimal strategy with which to control the experiment when it is used as a training environment to provide different abstraction levels of system state reactions.

Iterative multiscale and multi-physics computations for operando catalyst nanostructure elucidation and kinetic modeling

Modern heterogeneous catalysis has benefitted immensely from computational predictions of catalyst structure and its evolution under reaction conditions, first-principles mechanistic investigations, and detailed kinetic modeling, which are rungs on a multiscale workflow. Establishing connections across these rungs and integration with experiments have been challenging. Here, operando catalyst structure prediction techniques using density functional theory simulations and ab initio thermodynamics calculations, molecular dynamics, and machine learning techniques are presented. Surface structure characterization by computational spectroscopic and machine learning techniques is then discussed. Hierarchical approaches in kinetic parameter estimation involving semi-empirical, data-driven, and first-principles calculations and detailed kinetic modeling via mean-field microkinetic modeling and kinetic Monte Carlo simulations are discussed along with methods and the need for uncertainty quantification. With these as the background, this article proposes a bottom-up hierarchical and closed loop modeling framework incorporating consistency checks and iterative refinements at each level and across levels.

Exact distributed kinetic Monte Carlo simulations for on-lattice chemical kinetics: lessons learnt from medium- and large-scale benchmarks

Kinetic Monte Carlo (KMC) simulations have been instrumental in multiscale catalysis studies, enabling the elucidation of the complex dynamics of heterogeneous catalysts and the prediction of macroscopic performance metrics, such as activity and selectivity. However, the accessible length- and time-scales have been a limiting factor in such simulations. For instance, handling lattices containing millions of sites with 'traditional' sequential KMC implementations is prohibitive owing to large memory requirements and long simulation times. We have recently established an approach for exact, distributed, lattice-based simulations of catalytic kinetics which couples the Time-Warp algorithm with the Graph-Theoretical KMC framework, enabling the handling of complex adsorbate lateral interactions and reaction events within large lattices. In this work, we develop a lattice-based variant of the Brusselator system, a prototype chemical oscillator pioneered by Prigogine and Lefever in the late 60s, to benchmark and demonstrate our approach. This system can form spiral wave patterns, which would be computationally intractable with sequential KMC, while our distributed KMC approach can simulate such patterns 15 and 36 times faster with 625 and 1600 processors, respectively. The medium- and large-scale benchmarks thus conducted, demonstrate the robustness of the approach, and reveal computational bottlenecks that could be targeted in further development efforts. This article is part of a discussion meeting issue 'Supercomputing simulations of advanced materials'.

Site Communication in Direct Formation of H2O2 over Single-Atom Pd@Au Nanoparticles

Single atom alloy catalysts offer possibilities to obtain turnover frequencies and selectivities unattainable by their monometallic counterparts. One example is direct formation of H₂O₂ from O₂ and H₂ over Pd embedded in Au hosts. Here, a first-principles-based kinetic Monte Carlo approach is developed to investigate the catalytic performance of Pd embedded in Au nanoparticles in an aqueous solution. The simulations reveal an efficient site separation where Pd monomers act as active centers for H₂ dissociation, whereas H₂O₂ is formed over undercoordinated Au sites. After dissociation, atomic H may undergo an exothermic redox reaction, forming a hydronium ion in the solution and a negative charge on the surface. H₂O₂ is preferably formed from reactions between dissolved H⁺ and oxygen species on the Au surface. The simulations show that tuning the nanoparticle composition and reaction conditions can enhance the selectivity toward H₂O₂. The outlined approach is general and applicable for a range of different hydrogenation reactions over single atom alloy nanoparticles.

Assessing the Impact of Adlayer Description Fidelity on Theoretical Predictions of Coking on Ni(111) at Steam Reforming Conditions

Methane steam reforming is an important industrial process for hydrogen production, employing Ni as a low-cost, highly active catalyst, which, however, suffers from coking due to methane cracking. Coking is the accumulation of a stable poison over time, occurring at high temperatures; thus, to a first approximation, it can be treated as a thermodynamic problem. In this work, we developed an Ab initio kinetic Monte Carlo (KMC) model for methane cracking on Ni(111) at steam reforming conditions. The model captures C-H activation kinetics in detail, while graphene sheet formation is described at the level of thermodynamics, to obtain insights into the "terminal (poisoned) state" of graphene/coke within reasonable computational times. We used cluster expansions (CEs) of progressively higher fidelity to systematically assess the influence of effective cluster interactions between adsorbed or covalently bonded C and CH species on the "terminal state" morphology. Moreover, we compared the predictions of KMC models incorporating these CEs into mean-field microkinetic models in a consistent manner. The models show that the "terminal state" changes significantly with the level of fidelity of the CEs. Furthermore, high-fidelity simulations predict C-CH island/rings that are largely disconnected at low temperatures but completely encapsulate the Ni(111) surface at high temperatures.

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.

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.

A probabilistic microkinetic modeling framework for catalytic surface reactions

We present a probabilistic microkinetic modeling (MKM) framework that incorporates the short-ranged order (SRO) evolution for adsorbed species (adspecies) on a catalyst surface. The resulting model consists of a system of ordinary differential equations. Adsorbate-adsorbate interactions, surface diffusion, adsorption, desorption, and catalytic reaction processes are included. Assuming that the adspecies ordering/arrangement is accurately described by the SRO parameters, we employ the reverse Monte Carlo (RMC) method to extract the relevant local environment probability distributions and pass them to the MKM. The reaction kinetics is faithfully captured as accurately as the kinetic Monte Carlo (KMC) method but with a computational time requirement of few seconds on a standard desktop computer. KMC, on the other hand, can require several days for the examples discussed. The framework presented here is expected to provide the basis for wider application of the RMC-MKM approach to problems in computational catalysis, electrocatalysis, and material science.

A Kinetic Monte Carlo Approach to Model Barite Dissolution: The Role of Reactive Site Geometry

Barite (Ba[SO4]) is one of the promising candidates for sequestration of radioactive waste. Barite can incorporate radium (Ra) and form ideal solid solutions, i.e., (Ba,Ra)[SO4]. Together with isostructural celestite (Sr[SO4]), ternary solid solutions, (Ba,Sr,Ra)[SO4], may exist in natural conditions. Our fundamental understanding of the dissolution kinetics of isostructural sulfates is critically important for a better risk assessment of nuclear waste repositories utilizing this mineral for sequestration. So far, the barite-water interface has been studied with experimental methods and atomistic computer simulations. The direct connection between the molecular scale details of the interface structure and experimental observations at the microscopic scale is not yet well understood. Here, we began to investigate this connection by using a kinetic Monte Carlo approach to simulate the barite dissolution process. We constructed a microkinetic model for the dissolution process and identified the reactive sites. Identification of these sites is important for an improved understanding of the dissolution, adsorption, and crystal growth mechanisms at the barite–water interface. We parameterized the molecular detachment rates by using the experimentally observed etch pit morphologies and atomic step velocities. Our parameterization attempts demonstrated that local lattice coordination is not sufficient to differentiate between the kinetically important sites and estimate their detachment rates. We suggest that the water structure and dynamics at identified sites should substantially influence the detachment rates. However, it will require more work to improve the parameterization of the model by means of Molecular Dynamics and ab initio calculations.