QUBEKit 13 interfaces with the Gaussian09 47 and ONETEP QM software packages to perform bond, angle, torsion, charge and Lennard-Jones parameter derivation. BOSS/MCPRO style z-matrices, and the corresponding PDB files, of the 22 NNRTIs were generated using the LigParGen web server. Gaussian09 input files were prepared using QUBEKit. Structural optimizations and Hessian matrix calculations were performed with the ωB97X-D functional and a 6-311++G(d,p) basis set. Harmonic bond stretching and angle bending parameters were derived using the modified Seminario method. Non-bonded parameter derivation was performed with the linear-scaling density functional theory code, ONETEP, 20 using previously reported protocols. Density derived electrostatic and chemical (DDEC) electron density partitioning, as implemented in ONETEP, was used to assign atom-centered point charges and atomic volumes. Lennard-Jones parameters were derived by QUBEKit from the atomic volumes using the Tkatchenko-Scheffler method. Torsion parameter fitting on core fragments of the ligand sets followed the general methods used in previous studies, though here we employ an interface between QUBEKit and the TorsionDrive package, which improves the quality of both QM and MM torsion scans through its recursive wavefront propagation algorithm. Final force fields were output in xml format. Further descriptions of parameters used in force field derivation may be found in the Supporting Methods.