Reference energies and forces were calculated with DFT using the PBE exchange correlation functional with the D3(BJ) dispersion correction (PBE-D3). These were used to train NNPs of the Behler-Parinello high-dimensional NNP form as implemented in the n2p2 code. The input to the NNPs are vectors of radial and angular atomcentered symmetry functions (ACSFs). ACSFs were selected from a larger set by a CUR decomposition following the procedure detailed by Imbalzano et al., 54 which offers an improvement over a general set of symmetry functions (Table ). For oxalic acid, 64 symmetry functions per element were selected, while for TTBI and resorcinol 128 symmetry functions were selected. In all cases a maximum radial cutoff of 8.0 Å was used. Further details of the reference calculations and NNP models are provided in the supporting information.