Even though a root mean square error of 2.38 meV/atom over the final GAP training set reflects a quite high accuracy, we re-optimized the LMs and TSs for the processes shown in Figures and to ensure a consistent accuracy over all the presented processes. The geometry optimization of the LMs was performed with the conjugate gradient algorithm implemented in VASP, while the optimization to the TS was performed with dimer calculations, starting from the GAP-identified TS structure, along the direction between the initial state and TS structures. Both were performed until residual forces were below 0.01 eV/ Å. Differences in GAP-and DFT-optimized structures and energies are shown in the SI (Figures and). The DFT activation barriers ∆E proc stated below are extracted as the energy difference of the corresponding LM and TS structures, while the corresponding rate constants k proc are determined using the computed local Hessians and harmonic transition state theory. These and all DFT calculations performed for the GAP training employ the Vienna Ab-initio Simulation Package (VASP) version 6.3.2 and the PBE exchangecorrelation functional, with the D3 dispersion correction scheme with full details on the computational settings provided in the SI.