The reactive chemistry of molecular hydrogen at surfaces, notably dissociative sticking and hydrogen evolution, plays a crucial role in energy storage and fuel cells. Theoretical studies can help to decipher underlying mechanisms and reaction design, but studying dynamics at surfaces is computationally challenging due to the complex electronic structure at interfaces and the high sensitivity of dynamics to reaction barriers. In addition, ab initio molecular dynamics, based on density functional theory, is too computationally demanding to accurately predict reactive sticking or desorption probabilities, as it requires averaging over tens of thousands of initial conditions. High-dimensional machine learning-based interatomic potentials are starting to be more commonly used in gas-surface dynamics, yet robust approaches to generate reliable training data and assess how model uncertainty affects the prediction of dynamic observables are not well established. Here, we employ ensemble learning to adaptively generate training data while assessing model performance with full uncertainty quantification (UQ) for reaction probabilities of hydrogen scattering on different copper facets. We use this approach to investigate the performance of two message-passing neural networks, SchNet and PaiNN. Ensemble-based UQ and iterative refinement allow us to expose the shortcomings of the invariant pairwise-distance-based feature representation in the SchNet model for gas-surface dynamics.

Simulations of chemical reaction probabilities in gas surface dynamics require the calculation of ensemble averages over many tens of thousands of reaction events to predict dynamical observables that can be compared to experiments. At the same time, the energy landscapes need to be accurately mapped, as small errors in barriers can lead to large deviations in reaction probabilities. This brings a particularly interesting challenge for machine learning interatomic potentials, which are becoming well-established tools to accelerate molecular dynamics simulations. We compare state-of-the-art machine learning interatomic potentials with a particular focus on their inference performance on CPUs and suitability for high throughput simulation of reactive chemistry at surfaces. The considered models include polarizable atom interaction neural networks (PaiNN), recursively embedded atom neural networks (REANN), the MACE equivariant graph neural network, and atomic cluster expansion potentials (ACE). The models are applied to a dataset on reactive molecular hydrogen scattering on low-index surface facets of copper. All models are assessed for their accuracy, time-to-solution, and ability to simulate reactive sticking probabilities as a function of the rovibrational initial state and kinetic incidence energy of the molecule. REANN and MACE models provide the best balance between accuracy and time-to-solution and can be considered the current state-of-the-art in gas-surface dynamics. PaiNN models require many features for the best accuracy, which causes significant losses in computational efficiency. ACE models provide the fastest time-to-solution, however, models trained on the existing dataset were not able to achieve sufficiently accurate predictions in all cases.

Accurate and efficient methods to simulate nonadiabatic and quantum nuclear effects in high-dimensional and dissipative systems are crucial for the prediction of chemical dynamics in the condensed phase. To facilitate effective development, code sharing, and uptake of newly developed dynamics methods, it is important that software implementations can be easily accessed and built upon. Using the Julia programming language, we have developed the NQCDynamics.jl package, which provides a framework for established and emerging methods for performing semiclassical and mixed quantum–classical dynamics in the condensed phase. The code provides several interfaces to existing atomistic simulation frameworks, electronic structure codes, and machine learning representations. In addition to the existing methods, the package provides infrastructure for developing and deploying new dynamics methods, which we hope will benefit reproducibility and code sharing in the field of condensed phase quantum dynamics. Herein, we present our code design choices and the specific Julia programming features from which they benefit. We further demonstrate the capabilities of the package on two examples of chemical dynamics in the condensed phase: the population dynamics of the spin-boson model as described by a wide variety of semiclassical and mixed quantum–classical nonadiabatic methods and the reactive scattering of H₂ on Ag(111) using the molecular dynamics with electronic friction method. Together, they exemplify the broad scope of the package to study effective model Hamiltonians and realistic atomistic systems.

Molecular motion at metallic surfaces is affected by nonadiabatic effects and electron-phonon coupling. The ensuing energy dissipation and dynamical steering effects are not captured by classical molecular dynamics simulations, but can be described with the molecular dynamics with electronic friction method and linear response calculations based on density functional theory. Herein, we present an implementation of electron-phonon response based on an all-electron numeric atomic orbital description in the electronic structure code FHI-aims. After providing details of the underlying approximations and numerical considerations, we present significant scalability and performance improvements of the new code compared to a previous implementation (Maurer et al 2016 Phys. Rev. B 94 115432). We compare convergence behaviour and results of our simulations for exemplary systems such as H₂ adsorption on Cu(111), and CO on Ru(0001) against existing plane wave implementations. We examine different expressions to calculate electronic friction and vibrational lifetimes for their reliability and ease of convergence. Finally, we show the capabilities of the new code by studying the contribution of interband and intraband excitations to the vibrational lifetime of aperiodic adsorbate motion in large, previously unfeasible, periodic surface models.

Dynamics of coarse-grained particle systems derived via the Mori–Zwanzig projection formalism commonly take the form of a (generalized) Langevin equation with configuration-dependent friction tensor and diffusion coefficient matrix. In this article, we introduce a class of equivariant representations of tensor-valued functions based on the Atomic Cluster Expansion framework that allows for efficient learning of such configuration-dependent friction tensors from data. Besides satisfying the correct equivariance properties with respect to the Euclidean group E(3), the resulting heat bath models satisfy a fluctuation-dissipation relation. We demonstrate the capabilities of the model approach by fitting a model of configuration-dependent tensorial electronic friction calculated from first principles that arises during reactive molecular dynamics at metal surfaces.

In the dynamics of atoms and molecules at metal surfaces, electron-hole pair excitations can play a crucial role. In the case of hyperthermal hydrogen atom scattering, they lead to nonadiabatic energy loss and highly inelastic scattering. Molecular dynamics with electronic friction simulation results, based on an isotropic homogeneous electron gas approximation, have previously aligned well with measured kinetic energy loss distributions, indicating that this level of theoretical description is sufficient to describe nonadiabatic effects during scattering. In this study, we demonstrate that friction derived from density functional theory linear response calculations can also describe the experimental energy loss distributions, although agreement is slightly worse than for the simpler isotropic homogeneous electron gas approximation. We show that the apparent agreement of the homogeneous electron gas approximation with experiment arises from a fortuitous cancellation of errors as friction is overestimated close to the surface and the spin transition is neglected. Differences in frictional treatment affect single, double, and multibounce scattering trajectories in distinct ways, altering the shape of low-temperature energy loss distributions. These distinctions are largely absent at room temperature but may be measurable in future low-temperature scattering experiments.