PaperlayerView exact revision on arXiv

Source-linked AI summary

A Practical Guide to Surface Kinetic Monte Carlo Simulations

Mie Andersen, Chiara Panosetti, Karsten Reuter

arXiv:1904.02561v1physics.comp-phcond-mat.mtrl-sci

TL;DR

KMC needs a practical framework for simulating complex surface and interface dynamics without explicitly resolving every vibrational motion or configuration. This review presents lattice-KMC modeling practices, worked examples, rate-constant derivations, and acceleration methods, showing both the method’s broad utility and its sensitivity to model choices and timescale disparity. Its guidance emphasizes careful process selection, uncertainty assessment, and treatment of lateral interactions.

  • Problem

    Surface KMC models must handle enormous configuration spaces, widely separated event timescales, and rate constants that can be uncertain by orders of magnitude.

  • Method

    The review provides practical lattice-KMC guidance through kmos examples covering lattice mapping, rate-constant construction, initial states, lateral interactions, and acceleration algorithms.

  • Results

    The examples demonstrate that included processes and rate constants can substantially alter predictions, while acceleration algorithms address fast-process bottlenecks in suitable cases.

  • Takeaways & Limitations

    Reliable KMC practice requires sensitivity analysis, uncertainty quantification, careful steady-state identification, and explicit attention to timescale disparity and lateral interactions.

Abstract

from arXiv · show

This review article is intended as a practical guide for newcomers to the field of kinetic Monte Carlo (KMC) simulations, and specifically to lattice KMC simulations as prevalently used for surface and interface applications. We will provide worked out examples using the kmos code, where we highlight the central approximations made in implementing a KMC model as well as possible pitfalls. This includes the mapping of the problem onto a lattice and the derivation of rate constant expressions for various elementary processes. Example KMC models will be presented within the application areas surface diffusion, crystal growth and heterogeneous catalysis, covering both transient and steady-state kinetics as well as the preparation of various initial states of the system. We highlight the sensitivity of KMC models to the elementary processes included, as well as to possible errors in the rate constants. For catalysis models in particular, a recurrent challenge is the occurrence of processes at very different timescales, e.g. fast diffusion processes and slow chemical reactions. We demonstrate how to overcome this timescale disparity problem using recently developed acceleration algorithms. Finally, we will discuss how to account for lateral interactions between the species adsorbed to the lattice, which can play an important role in all application areas covered here.

I. INTRODUCTION

KMC coarse-grains rare atomistic events into stochastic trajectories, enabling simulations of surface diffusion, crystal growth, and catalysis across relevant scales. The review focuses on practical lattice-KMC modeling, including lattice mapping, rate constants, approximations, and simulation pitfalls.

  • KMC describes transport, material evolution, and equilibrium or nonequilibrium chemistry by statistically propagating discrete elementary events.
  • Surface processes such as adsorption, desorption, diffusion, bond breaking, and bond forming can be represented as elementary events on a catalyst surface.
  • High activation barriers make relevant surface transitions rare, so meaningful conclusions may require simulations extending from picoseconds to seconds or longer.
  • The master equation becomes infeasible for surface systems because configurations grow exponentially with the number of sites.
  • KMC avoids explicit matrix solution by generating stochastic state-to-state trajectories, selecting processes from their rate constants and advancing simulation time.

B. Mapping the problem onto a lattice

Lattice KMC maps potential-energy-surface minima onto periodic lattice positions, reducing the process space through crystalline symmetry and local-environment approximations. These efficiencies require a suitably chosen ordered lattice and can omit reconstruction or other lattice transformations.

  • Lattice KMC represents system states through adsorbate distributions over lattice positions mapped from relevant potential-energy-surface minima.
  • Crystalline translational symmetry allows equivalent elementary processes to reuse the same computed rate constants.
  • Nearest-neighbor lateral interactions require five distinct diffusion rates for the possible occupations of four neighboring sites, while broader interactions require more calculations.
  • Rate catalogs precompute local process rates, whereas on-the-fly evaluation calculates cluster-expansion corrections during each KMC step.
  • An unsuitable ordered lattice may omit important energy minima and cannot capture reaction-induced reconstruction, morphological transitions, or loading-induced lattice transformations.

C. Mean-field approximation

Mean-field approximation simplifies the master equation by replacing spatially resolved configurations with average coverages and related rate equations. This enables standard numerical or sometimes analytical solutions, but sacrifices spatial correlations and remains an approximation.

  • Mean-field approximation replaces detailed spatial resolution over the lattice with mean coverages for each species and site type.
  • For second-order reactions, mean-field treatment uses products of averaged coverages instead of explicit pair probabilities.
  • Geometrical connectivity factors account for the arrangement of sites involved in the initial and final states.
  • The approximation produces coupled differential equations solvable by standard algorithms and can support analytical reaction-rate expressions with rate-determining-step assumptions.
  • Mean-field models simplify the problem substantially but should be treated as approximations rather than exact descriptions.

III. GETTING PRACTICAL: ALGORITHMS AND INPUT DATA

The BKL or n-fold-way KMC algorithm samples among possible elementary processes according to their rates and advances time using the total escape rate. Randomized trajectory generation replaces explicit master-equation solution while preserving the appropriate stochastic evolution.

  • The BKL algorithm selects an executable process from all possible processes using rate-weighted random sampling.
  • KMC time increments arise from exponentially distributed first-escape statistics rather than simply using the inverse rate of one selected process.
  • The time to escape a state depends on the total rate constant, whereas the selected process determines the next state.
  • Independent random seeds generate statistically independent trajectories for repeated simulations.
  • At each step, the algorithm lists possible processes, extracts one, executes it, and updates the clock using t → t −ln(ρ2)/ktot.

IV. RATE CONSTANTS FROM FIRST PRINCIPLES: TRANSITION STATE THEORY

Transition State Theory supplies KMC rate constants from activation barriers and state/transition-state partition functions, but common prefactor approximations have important limits, especially for non-activated adsorption.

  • KMC surface models predominantly obtain elementary-process rate constants through Transition State Theory under assumptions including no recrossing and classical barrier crossing.
  • The activation barrier ΔEij comes from the potential-energy surface, making it accessible to first-principles calculations, while the prefactor can in principle be calculated.
  • Harmonic TST estimates partition functions from vibrational modes at the initial state and transition state.
  • 1012–1013 s−1: studies often approximate the prefactor as k0 ≃ 1–10 to avoid the substantial cost of vibrational calculations.
  • Non-activated adsorption requires a different rate expression because the prefactor must account for strong entropy reduction from the gas phase to the surface-bound state.

A. Master equation and detailed balance

At steady state, transition flows into and out of each state balance; thermodynamic equilibrium further requires microscopic forward and reverse rates to balance. KMC models therefore need reversible processes and thermodynamically consistent rate expressions and state free energies.

  • At steady state, the Master equation requires the time derivative of state probabilities to vanish.
  • This conservation law equates the summed rates of transitions into and out of every state.
  • At thermodynamic equilibrium, microscopic reversibility and detailed balance require each microscopic process to balance its reverse process on average.
  • KMC models should define a reverse process for every microscopic process and use forward and reverse rate expressions satisfying Eq. 15.
  • Forward and reverse free energies must use the same numerical approximations; otherwise, the kinetic model can become thermodynamically inconsistent with drastic consequences.

B. Obtaining rate constants: Transition state search

Transition-state search strategy depends on whether the mechanism and endpoint structures are known. Lattice KMC normally predefines elementary processes and endpoints, enabling interpolation methods, whereas adaptive KMC may require local exploration from the current state.

  • Lattice KMC assumes or derives a mechanism and compiles its elementary processes before rate constants are determined.
  • Because lattice-KMC initial and final states are predetermined, interpolation methods can be used to locate transition states.
  • Adaptive KMC often lacks a prior mechanism assumption and may need to explore possible, probable escape pathways from the current state.
  • Local methods are mandatory for adaptive searches when only information about the initial state is available.

1. Interpolation methods

Interpolation methods connect known initial and final states to locate transition states, ranging from coordinate-driving searches to chain-of-state methods such as NEB. Their accuracy and efficiency depend strongly on force treatment and technical parameters, while local methods avoid endpoint information but impose other requirements.

  • Interpolation methods: Coordinate-driving methods constrain selected internal coordinates between initial and final structures while optimizing the remaining degrees of freedom.
  • Interpolation methods: Ridge and other chain-of-state methods refine multiple structures distributed along an initial path between known endpoints.
  • Interpolation methods: NEB minimizes image energies plus a spring penalty that distributes images along the reaction path.
  • Interpolation methods: Nudging uses the spring-force component parallel to the path and the energy-force component perpendicular to it to reduce corner-cutting and down-sliding.
  • Interpolation methods: CI-NEB drives the highest-energy image toward the saddle point by removing its spring force and reversing the parallel potential-force component.
  • Interpolation methods: NEB can be computationally expensive and sensitive to image count and spring constant, motivating adaptive-resolution approaches such as AutoNEB.
  • Local methods: Local methods use PES values and derivatives without endpoint geometries, but Newton–Raphson requires a sufficiently close starting point and a full Hessian.
  • Local methods: The dimer method avoids explicit Hessian calculation by using two symmetrically displaced replicas to determine the saddle-search direction.

3. BEP and scaling relations

Approximate energy and scaling relations reduce the cost of generating first-principles rate constants, but their transferability is limited. KMC results remain sensitive to the completeness and accuracy of the supplied processes and rates.

  • BEP and scaling relations: The Brønsted-Evans-Polanyi relation estimates activation barriers through linear relationships based on initial- and final-state energies.Its fitted parameters must be obtained from suitable first-principles calculations.
  • BEP and scaling relations: BEP parameters are generally transferable only to site types similar to those used for fitting.The passage gives fcc(111) sites as an example of a similar site type.
  • BEP and scaling relations: BEP and binding-energy scaling relations can substantially reduce the computational cost of obtaining first-principles rate constants for catalyst screening.Scaling relations correlate molecular binding energies with those of a small set of base elements.
  • Input-data limitations: KMC model predictions are limited by the quality of their input processes and rate constants under the garbage in–garbage out principle.First-principles KMC therefore depends critically on the reliability of electronic-structure inputs.
  • Input-data limitations: Advanced KMC approaches seek to identify relevant processes automatically using accelerated molecular dynamics and adaptive on-the-fly process searches.The passage names hyperdynamics, temperature-accelerated dynamics, replica-exchange dynamics, dimer methods, and high-temperature searches.
  • Implementation: The worked examples use the kmos software package and are available in the supplemental data.

A. Adatom diffusion on Au(100)

The Au(100) diffusion example compares hopping with exchange diffusion in a lattice KMC model. Including the lower-barrier exchange mechanism changes the predicted diffusion coefficient by more than three orders of magnitude, while rate errors matter most for dominant processes.

  • Hopping diffusion: The basic model represents Au adatoms hopping among four neighboring hollow sites on a (20 × 20) square lattice.Neglecting lateral interactions, all hops use the same DFT-LDA diffusion barrier of 0.83 eV.
  • Hopping diffusion: The diffusion coefficient is obtained by tracking adatom mean squared displacements over simulation time.The coefficient uses lattice dimension d = 2 and averages over adatoms and 25 independently initialized simulations.
  • Exchange diffusion: 4.7 nm2/s is the diffusion coefficient after adding exchange diffusion, more than three orders of magnitude above the hopping-only result.Exchange diffusion replaces a surface atom and has a 0.65 eV barrier, lower than hopping’s 0.83 eV barrier.
  • Sensitivity and model completeness: Lowering the hopping barrier changes the result substantially only when it approaches the exchange barrier, whereas lowering the dominant exchange barrier by 0.05 eV produces a seven-times-higher diffusion constant.The contrast illustrates the exponential sensitivity of rates to barriers and the importance of including relevant processes.
  • Sensitivity and model completeness: Comprehensive KMC simulations can identify rate-determining processes in complex catalytic cycles rather than assuming them.

A. CO oxidation on RuO2(110)

The RuO2(110) CO-oxidation model uses alternating bridge and coordinately unsaturated sites, with 26 adsorption, desorption, diffusion, and reaction processes. Its evolution demonstrates that apparent steady states can precede the true steady state by orders of magnitude in time and activity.

  • Model construction: The model represents RuO2(110) as alternating bridge and coordinately unsaturated site columns, each occupied by O, CO, or a vacancy.The 26 processes include CO and O2 adsorption/desorption, O and CO diffusion, and CO2 formation.
  • Transient evolution: Starting from an empty lattice at 450 K and 1 bar CO and O2, surface coverage builds rapidly on nanosecond timescales.O coverage initially grows roughly twice as fast as CO coverage because each O2 adsorption event produces two O atoms.
  • Transient evolution: After 20–30 ns, nearly constant coverages and TOF can create a premature apparent steady state, with approximately two-thirds O and one-third CO coverage.These early coverages reflect the relative impingement situation rather than the true long-time steady state.
  • Steady-state behavior: On second timescales, bridge-site O decreases to about one-third, cus sites become essentially CO-covered, and TOF falls by more than three orders of magnitude.Strongly bound, relatively unreactive bridge oxygen slows the longer-term transformation, while CO replaces reactive cus oxygen.
  • Steady-state verification: The same steady state is reached from a distinct initialization with oxygen-covered cus sites and vacant bridge sites, despite a different nanosecond-scale quasi-steady state.Multiple initial conditions are therefore useful for testing convergence, although automatic fool-proof steady-state detection is unavailable.
  • Steady-state verification: Once steady state and ergodicity are established, coverages, elementary-step frequencies, and TOFs can be computed as time averages weighted by KMC escape times.The simulation must run long enough to reduce statistical error, and trajectory fractions must exceed the decorrelation time to be independent.

VII. SENSITIVITY ANALYSIS AND UNCERTAINTY QUANTIFICATION

KMC models can be highly sensitive to uncertain rate constants, while disparate process timescales can make important slow dynamics computationally inaccessible. The review describes sensitivity-analysis practices and acceleration schemes, including their assumptions and failure modes.

  • Sensitivity analysis and uncertainty quantification: Sensitivity analysis identifies which rate constants require accurate calculation and which conclusions remain reliable despite input errors.It also supports design and optimization questions about performance limits and how material changes should alter rate constants.
  • Sensitivity analysis and uncertainty quantification: Campbell’s degree of rate control measures how the average product-formation rate responds locally to changing one rate constant while holding other rate constants and its equilibrium constant fixed.Positive, negative, and zero values indicate increasing, decreasing, or insensitive reaction rates, respectively.
  • Sensitivity analysis and uncertainty quantification: DRC results are local because they use linear response theory, whereas nonlinear kinetic models can change substantially under reaction-barrier errors of ±0.2 eV.Global sensitivity analysis and UQ address broader parameter uncertainty, including correlated errors from DFT functionals and scaling relations.
  • Sensitivity analysis and uncertainty quantification: A practical baseline is to vary key rate constants manually, identify critical dependencies, and escalate to more sophisticated analysis when sensitivity is high.The review notes that sophisticated sensitivity analysis is not yet used by most practitioners.
  • Timescale disparity problem: Fast processes can consume nearly all CPU time in timescale-disparate KMC models, leaving slow but important dynamics poorly sampled, especially for surface reactions on metals.The disparity commonly involves fast diffusion and slow chemical reactions.
  • Timescale disparity problem: Acceleration schemes scale down the fastest processes, with manual verification for simple networks and algorithms that automate scaling for more complex systems.The τ-leap method fires multiple processes simultaneously but requires approximately constant surface populations and therefore does not suit microscopic lattices with abrupt site-population changes.
  • Timescale disparity problem: Superbasin-based acceleration scales fast, quasi-equilibrated processes while retaining slow escape processes, and can reduce the time needed to leave the current superbasin.Related approaches have reproduced KMC simulations at significantly reduced cost and improved sensitivity-measure sampling.
  • Timescale disparity problem: Acceleration can fail when two independently produced, low-coverage species must meet, because overly aggressive scaling may undersample neighboring configurations.A nearest-neighbor correction works for directly neighboring production but not when species form at distant sites and require diffusion before reacting.

IX. LATERAL INTERACTIONS

Lateral interactions between adsorbed species can be attractive or repulsive and often matter for correctly describing surface dynamics. Lattice KMC represents them through configuration-dependent rates or cluster-expansion energetics, but increasing interaction complexity can become computationally expensive.

  • IX. LATERAL INTERACTIONS: Lateral interactions are attractive or repulsive interactions between species adsorbed to a KMC lattice.Recent studies identify them as essential for correctly describing system dynamics.
  • IX. LATERAL INTERACTIONS: In lattice KMC, lateral interactions can be represented by assigning hopping-rate constants to each neighboring configuration.Repulsive interactions make particles more likely to leave configurations with occupied neighboring sites than empty ones.
  • IX. LATERAL INTERACTIONS: Cluster expansion expresses lattice energy as discrete pairwise, three-body, and higher-order cluster interactions, while transition-state effects can be approximated through a Brønsted-Evans-Polanyi relation.The approach is available in ZACROS for adsorbates.
  • IX. LATERAL INTERACTIONS: Increasing the number of clusters from 3 to 1224 increased KMC computational cost by about 5 orders of magnitude in one benchmark.The review notes that the exact cost depends on implementation and algorithm.
  • IX. LATERAL INTERACTIONS: Despite added cost, lateral interactions often matter for surface diffusion, crystal growth, and heterogeneous catalysis.The review therefore presents crystal-growth and catalysis examples that include lateral interactions.

A. kmos models with lateral interactions

The kmos examples show how lateral interactions can improve model realism while creating substantial computational costs that depend on backend design. The review also frames broader KMC practice around rate-constant accuracy, acceleration of timescale disparities, and careful lattice-model construction.

  • Backend choices: The local smart backend uses a pre-calculated rate catalog to prioritize runtime performance, whereas the on-the-fly backend calculates rate constants during simulation.The on-the-fly approach is designed for complex lateral interaction models and keeps the total number of processes constant with respect to interaction count.
  • Crystal growth: The crystal-growth model includes four same-height neighboring sites whose occupations modify desorption rates, producing 16 inequivalent desorption process types.The site below the desorbing species is always occupied and is included in the low-coverage desorption barrier.
  • Crystal growth: At 450 K, the grown crystal structure is smoother than at 350 K because unfavorable adsorption sites are more likely to desorb.The model uses attractive pairwise lateral interactions and begins from a fixed substrate layer.
  • CO oxidation: For CO oxidation, the local smart backend is effectively limited to four pairwise interactions because its rate-catalog cost grows exponentially with interaction count.The number of processes in the rate catalog grows exponentially as more interactions are included.
  • Summary and outlook: Lateral interactions can increase KMC accuracy through cluster expansions but also increase the number of required rate constants and computational demands.This trade-off is presented as a limitation relative to simpler mean-field treatments.
  • Summary and outlook: The review emphasizes lattice-approximation pitfalls, transition-state bottlenecks, rate-constant sensitivity, and acceleration algorithms for disparate timescales.Acceleration methods can identify fast, quasi-equilibrated processes, although some cases still cause them to break down.

XII. AUTHOR CONTRIBUTIONS

The author-contributions section records that all listed authors made substantial intellectual contributions and approved the publication, alongside funding acknowledgements and supplemental-data availability.

  • All listed authors made substantial, direct, and intellectual contributions to the work and approved it for publication.
  • The project received Horizon 2020 funding, with additional DFG support and publication support from DFG and TUM.
  • Files required to set up and run the described kmos models, together with additional explanations, are provided as supplemental data.
Loading 1904.02561v1…