While atom-centered point charges provide a good representation of the QM electrostatic potential (ESP) if the partitioned atomic electron density is spherical, in many cases this simple representation is inadequate. This situation occurs when there is significant anisotropy in the underlying electron distribution, and is common in molecules containing nitrogen, sulfur or halogens. Here, to model electron anisotropy, we employ off-center, "virtual" sites, which have been shown to be competitive with the use of more computationally expensive higher-order multipole electrostatics. Virtual sites are commonly used in water models, such as TIP4P, and various force fields for modelling lone pairs and σ-holes, 66 but the positions and charges of the virtual sites require fitting to experiment. On the other hand, it has recently been shown that virtual site positions may be derived directly from localized QM molecular orbitals, but currently the magnitudes of the charges are derived by fitting to the molecular dipole moment, which may be problematic for extension to larger molecules that contain multiple sites. In keeping with our goal of avoiding fitting FF parameters to experiment and developing methods that scale to biological molecules, we proposed a method that relied on the dipole and quadrupole moments of the partitioned atomic electron density, to optimize the charges and locations of virtual sites. However, the method employed did not consistently converge and resulted in a large number of off-center point charges. Modifications were required to correct these issues and improve the usability of the method in an automated high-throughput scenario.