Abstract: Understanding atomic-level processes in surface science and heterogeneous catalysis is complicated by the wide range of time scales and length scales needed for simulations. To accelerate molecular dynamics calculations, we rely on machine learning methods to capture interatomic interactions with quantum accuracy. We then implement and deploy these models on parallel GPUs to reach billions of atoms in size or microseconds in time. One method is a family of equivariant interatomic potential models (NequIP and Allegro) based on symmetry-preserving layer architectures that we use to achieve state-of-the-art accuracy and training efficiency for simulating atomistic. Another method (FLARE) enables autonomous selection of training sets for reactive systems, based on adaptive closed-loop algorithm that constructs accurate and uncertainty-aware Bayesian force fields on-the-fly from a molecular dynamics simulation. We examine the current limitations of machine learning models and highlight the usefulness of ML-accelerated MD simulations to study dislocation dynamics, surface reconstructions, direct heterogeneous reactions, and nanoparticle shape changes.