Markov State Models for Tracking Reaction Dynamics on Catalytic Nanoparticles

TL;DR

This work develops and applies Markov State Models to dynamically coarse-grain MD trajectories of hydrogen interacting with rhodium surfaces and nanoparticles, enabling automatic identification of slow kinetic modes without predefined reactive events. By combining a local, atom-centered MSM representation with a machine-learned ACE interatomic potential trained on DFT data, the study covers multiple geometries and hydrogen coverages to reveal how surface features and crowding shape reactivity. Key findings include that nanoparticle corners and edges can trap H2, slowing association/dissociation, while diffusion and facet exchange constitute distinct slow modes; hydrogen association rates show nonmonotonic dependence on coverage, highlighting dynamical effects not captured by standard transition-state theory. Overall, the MSM framework provides mechanistic insight into operando catalysis and offers guidance for tailoring nanoparticle geometry and surface occupancy to optimize activity.

Abstract

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.

Figures (12)

• Figure 1: Markov State Models for hydrogen dissociation/association on rhodium. A) Simulation data is generated for hydrogen atoms interacting with rhodium surfaces with four different geometries (2 nm nanoparticle, 5 nm nanoparticle, (100) slab, (111) slab) and various H concentrations. B) The local environment around each hydrogen atom $i$ is represented by a feature vector $\bm \phi^i({\bf x}_t)$. The feature values, $\phi_j^i({\bf x}_t),$ change as a function of time reflecting changes in local coordination and composition. C) Density plot of TICA projection of the high dimensional ($d_\phi = 131$) feature space into the five dynamically most relevant dimensions (tIC 1-5); Kmeans++ clustering is performed within this lower dimensional space to discretize the states, white points indicate cluster centers.
• Figure 2: The slowest dynamical processes of hydrogen atoms interacting on a large rhodium nanoparticle. The surface coverage of hydrogen in this example is about 50%. A) Simulation snapshot where hydrogen atoms (small spheres) are colored by their geometric state label and rhodium atoms are represented white, large spheres. H2 on corners/edges/facets are red/orange/yellow and H on (100)/(111)/edges/top sites are blue/green/light purple/dark purple. The local concentration of H on surfaces is also indicated by dark/light being concentrated/dilute. B) Labeled data projected into TICA space (sampled every ps, top graph) and a bar plot of the equilibrium distribution obtained from the values of the first eigenvector and integrated over the corresponding microstates. C) Projected flux of the two slowest dynamical modes and the states that contribute to these modes. The slowest process is the diffusion from a (111) facet to a (100) facet. The second slowest process is the association/dissociation reaction of H$\leftrightarrow$$\frac{1}{2}$H2. D) Spectral plot of the eigenvalues of the transition matrix and the implied timescale convergence of the two slowest processes.
• Figure 3: Configurations from the simulation data described in Fig. \ref{['fig:flux']} projected into TICA space and colored by A) the nearest distance to a neighboring hydrogen atom and B) the committor value for the dissociation reaction.
• Figure 4: A) Timescales of the association/dissociation reactions on rhodium slabs and nanoparticles as a function of hydrogen concentration. B) Rates of the association (H$\rightarrow\frac{1}{2}$H2, left graph) and dissociation ($\frac{1}{2}$H2$\rightarrow$H, right graph) of hydrogen as a function of concentration. Dashed lines show TST rates on (100) and (111) surfaces.
• Figure S1: Parity plot of the ACE predicted and DFT calculated energies for the entire set of MLIP training data. The error is more pronounced for the very high energy structures that were added in active learning steps.