Here, we propose a method for the derivation of virtual site positions and charges directly from the QM electron density in which the virtual sites are positioned so as to reproduce as closely as possible the QM ESP of the partitioned atomic electron density. By determining the virtual site parameters only using atomic properties, the method scales trivially to macromolecules such as proteins. In order to reduce the search space we limit the virtual site positions to those dictated by the symmetry of the atom's bonding environment. Together these improvements allow us to define virtual sites that improve the electrostatic properties of the simulated molecule in an automated manner.