Iterative training is essential to construct robust MLPs and can be summarised as follows: i) developing a training set, ii) performing the regression, iii) repeating the process until the desired accuracy is obtained. Automated approaches to construct large training sets have been previously reported; however, their application remain limited to small systems. Moreover, the need to perform several thousands of energy and gradient evaluations on a reference potential energy surface (PES) makes it difficult to go beyond density functional theory (DFT) and employ accurate but computationally expensive wavefunction (WF)-based methods. Exceptions are rare and limited to systems with ≲ 10 atoms. Finally, achieving the desired accuracy when training reactive MLPs is challenging because the energy scale required for accurate dynamics is much larger than for non-reactive MLPs (where sampling around the minima is often sufficient). As such, more training data is required to reach the same accuracy. This is particularly relevant when sampling areas around the transition states (TSs), which often require WF-based methods, such as coupled-cluster (CC). Here, we show that recently developed MLP methods can be used to generate accurate potentials for modestly sized reactions (≈ 50 atoms) in an automated fashion and outline some of the enabled insights. 2 Results and Discussion