Markov State Models for Tracking Reaction Dynamics on Catalytic Nanoparticles

Markov state models (MSMs) are a powerful tool to analyze and coarse-grain complex dynamical data into interpretable kinetic processes. This capability is particularly important in heterogeneous catalysis, where a medley of reactants and intermediates interact on surfaces that might simultaneously experience structural fluctuations. For these very complex systems, standard transition state theory (TST) approaches are no longer appropriate, motivating alternative approaches that can retain dynamical complexity while providing physical insight. With machine learned interatomic potentials being more and more ubiquitous, directly simulating complex catalytic systems with molecular dynamics (MD) is becoming increasingly feasible. Extending MSMs to dynamically coarse grain MD simulation data of catalytic processes, we analyze hydrogen dynamics on rhodium catalysts with slab and nanoparticle geometries over a range of hydrogen surface concentrations. Somewhat counterintuitively, nanoparticle features, such as corners and edges, effectively slow down the association/dissociation process, and the cooperative behavior of hydrogen-hydrogen interactions leads to a non-monotonic concentration dependence of the rates, which would not be predicted with standard TST.

💡 Research Summary

This paper presents a novel workflow that combines machine‑learned interatomic potentials (MLIPs) with Markov state models (MSMs) to extract mechanistic insight from large‑scale molecular dynamics (MD) simulations of heterogeneous catalysis. Using an atomic‑cluster‑expansion (ACE) potential trained on density‑functional theory data for Rh–H, the authors generate extensive MD trajectories for four distinct surface geometries: two rhodium nanoparticles (≈2 nm and ≈5 nm) and two flat slabs exposing the (100) and (111) facets. Simulations are performed at 450 K over a range of hydrogen surface coverages, providing a rich dataset that captures both chemical reactions (H₂ dissociation/association) and surface restructuring (diffusion, edge and corner dynamics).

The core methodological advance lies in the construction of a local, atom‑centered feature vector (131 components) for each hydrogen atom, describing two‑, three‑, and four‑body interactions within a 7 Å cutoff. These high‑dimensional descriptors are reduced by time‑lagged independent component analysis (TICA) to five slow collective coordinates, preserving the most kinetically relevant information. K‑means++ clustering in this reduced space yields 1 200 discrete microstates, which serve as the basis for estimating transition probability matrices at various lag times. Eigenvalue analysis of the resulting matrices reveals a clear separation between fast vibrational motions and a handful of slow processes, as evidenced by a spectral gap after the fourth eigenvalue.

Four dominant slow modes are identified across all systems: (i) adsorption/desorption equilibria (treated as a fast background process in subsequent analysis), (ii) H₂ association and dissociation, (iii) inter‑facet diffusion between (111) and (100) regions, and (iv) exchange between locally concentrated and dilute hydrogen domains. By projecting the eigenvectors onto the original structural labels (corner, edge, facet, “H‑top”), the authors map each kinetic mode onto physically interpretable pathways. Notably, a previously unreported “H‑top” state emerges on both facets and edges, where a hydrogen atom sits atop a rhodium atom while being trapped by two neighboring hydrogens; this state acts as a kinetic bottleneck for dissociation on nanoparticles.

A striking finding is the non‑monotonic dependence of the association/dissociation rates on hydrogen coverage. At intermediate coverages, cooperative H–H interactions lower the effective barrier, accelerating the reaction, whereas at very low or very high coverages the rates diminish due to insufficient partners or site saturation, respectively. This behavior cannot be captured by conventional transition‑state theory, which assumes isolated events and static surfaces. Moreover, the authors demonstrate that nanoparticle features such as edges and corners actually slow down the overall H₂ turnover compared with ideal flat facets, contrary to the common expectation that low‑coordinated sites are always more reactive.

The paper also discusses practical aspects of MSM construction: sensitivity to hyper‑parameters (number of TICA components, cluster count, lag time), validation through implied timescale convergence, and the use of the MSMbuilder software package for reproducibility. By removing adsorption/desorption from the trajectory analysis, the authors focus on intrinsic surface kinetics under an assumed equilibrium coverage, thereby isolating the intrinsic catalytic steps.

In summary, this work establishes a robust, data‑driven framework for dissecting complex catalytic dynamics on realistic nanostructured surfaces. It shows that MSMs, when combined with accurate MLIPs and careful local feature engineering, can reveal kinetic pathways, hidden metastable states, and cooperative effects that are invisible to traditional static‑energy‑barrier calculations. The approach is broadly applicable to other multi‑component catalytic systems, electrochemical interfaces, and any scenario where surface restructuring occurs on comparable timescales to chemical reactions, opening new avenues for rational catalyst design guided by full dynamical information.