Torsional constrained optimizations were performed in Gaussian09 61 with the same functional and basis set so as to be consistent with the other bonded terms. The torsional scan optimizations were performed in 15 • increments from 0 • to 360 • . The majority of the dihedral parameter fitting was done using no Boltzmann weighting (corresponding to T=∞) and regularization against OPLS reference values was applied with λ = 0.1. This was only changed in rare cases where it was particularly difficult to recreate the QM energy landscape, in which case λ = 0 and T = 2000K were used as previously suggested. Ground-state electron density calculations for non-bonded parameter derivation were performed using the linear-scaling DFT code ONETEP. Four nonorthogonal generalized Wannier functions (NGWFs), with radii of 10 Bohr, were used for all atoms with the exception of hydrogen, which used one. NGWFs were expanded in a periodic sine (psinc) basis, with a grid size (0.45a o ), corresponding to a plane wave cut-off energy of 1020 eV. The PBE exchange-correlation functional was used with PBE OPIUM norm-conserving pseudopoten-tials. The calculation was carried out in an implicit solvent using a dielectric of 4 to model induction effects. The DDEC module implemented in ONETEP was used to partition the electron density and assign atom-centered point charges and atomic volumes. The charges were assigned with a IH to ISA ratio of 0.02. The ESP error threshold, F thresh , was set to 0.9025 kcal/mol. The additional charges are only added if the decrease in ESP error is larger than F change = 0.0625 kcal/mol. The locations of the virtual sites were restricted using maximum distance cut-offs chosen by element, as virtual sites near the van der Waals radius can be detrimental. The cut-offs were defined as follows: 0.8 Å for N, 1 Å for O, S and F, and 1.5 Å for Cl and Br.