Central to this goal is the adequate description of the potential energy surface (PES). Wavefunction and density functional theory (DFT)-based methods are well-established approaches to explore reaction mechanisms by characterizing critical points along a reaction path, such as the reactant state (RS), transition state (TS), and product state (PS). However, these static methods have limitations as, in addition to being computationally costly, they do not account for dynamics, which is essential for exploring bifurcating PESs, flexible molecules and solvent effects. Dynamics simulations driven by ab initio methods, where energies and forces are computed "on the fly" by solving the Schrödinger equation at each time step, such as ab initio molecular dynamics (AIMD) and quasiclassical trajectories, enable realistic modelling of reaction mechanisms. These approaches have been applied, for example, in the identification of entropic intermediates, the prediction of product ratios for reactions exhibiting bifurcating surfaces, and the modelling of reactions in solution. 12-14 Despite these successes, the high computational cost of the underlying ab initio method limits the size, complexity, and timescale of the systems that can be studied. Moreover, these methods often require trade-offs between accuracy (level of theory used) and simulation time, leading to insufficient sampling and failure to achieve a converged free energy surface (FES).