It is desirable to perform the optimisation of eq 13 directly in QUBEKit, but storage and imports of the full atomic electron density is inefficient. Instead, V i QM is now reconstructed from the QM atomic multipole moments, up to quadrupole order, which are readily extracted from either Chargemol (using the DDEC AIM method) or the PSI4 software (using the MBIS AIM method). For atoms with isotropic electron density, the atomic multipoles are low, and the QM electrostatic potential is well-described by an atom-centred point charge (monopole). For atoms with anisotropic electron density (defined here as Φ > 1 kcal/mol), eq 13 is minimised with respect to the virtual site charge and position, using the Scipy python package. In previous works, new non-bonded force field parameters tend to be investigated in conjunction with the bonded parameters of standard transferable force fields, which may lead to incompatibility given the close inter-dependency between parameter types. Here, the benefit of assembling the force fields using the QUBEKit package, is that torsion parameters for rotatable bonds may be re-fit in a simple, automated manner for compatibility with the rest of the bonded and non-bonded parameters. Using a new interface between QUBEKit, ForceBalance and TorsionDrive, dihedral parameters of any freely rotatable bonds were optimised via a least squares minimisation of the root mean square error (RMSE) between QM and MM torsional potential energy surfaces. We use the same procedure for parameter fitting as the recent Open Force Field Parsley force field, and full details are provided in Supporting Information Section S3.1. A new feature in this work is that we have modified the ForceBalance code to allow parameter fitting in the presence of local coordinate virtual sites. These virtual sites then inherit their exceptions and exclusions from the parent atom meaning they interact with the same rules as the parent. The code may be obtained from conda-forge or github (), and is available from version 1.9.0 onward. Figure (c) demonstrates an example of bespoke dihedral parameter fitting for the molecule acetic anhydride, which has two virtual sites on the central oxygen atom. With the initial parameter set (which is incompatible with the new Lennard-Jones parameters and virtual sites), the RMS error is more than 3 kcal/mol, but following fitting it falls to just 0.12 kcal/mol. In fact, the average RMSE across all molecules in the 51 molecule test set (a total of 117 scans) after fitting was just 0.13 kcal/mol.