In addition to the dihedral potential energy surfaces, we also performed vibrational analyses to provide additional data for the bonded parameters. For each dipeptide in both protonation states, an unconstrained RI-MP2/aug-cc-pVTZ energy minimization was carried out for the overall lowest energy structure taken from the grids, followed by a RI-MP2/aug-cc-pVTZ numerical Hessian calculation and vibrational analysis for the vibrational frequencies and normal modes. Scaling factors for the vibrational frequencies were the same as those used for the AMBER-FB15 frequency calculations. Despite our best efforts to scan over the relevant degrees of freedom when generating the QM data, the optimized MM force field often has spurious local energy minima that occur at structures not covered by the QM data set. To address this issue, we generated new QM data energies and gradients from molecular mechanics (MM) energy-minimized structures using the optimized force field and appended them to the training data set, then re-optimized the force field with the expanded data set. After we have optimized the initial force field using the dihedral potential energy surfaces and frequency calculations, each structure on the 2D QM grid is optimized using the new force field, and the resulting structures are clustered using heavy-atom root-mean-square deviation (RMSD) and a cutoff of 0.1 Å. The clustering produces 25-50 geometries per residue / protonation state, from which RI-MP2/CBS singlepoint energies and RI-MP2/aug-cc-pVTZ gradients are computed and added to the dataset.