Population analyses are part of the theoretical chemist’s toolbox. They provide means to extract information about the repartition of the electronic density among molecules or solids. The values of atomic multipoles in a molecule can shed light on its electrostatic properties and may help to predict how different molecules could interact or to rationalize chemical reactivity for instance. Not being physical observables to which a quantum mechanical operator can be associated, atomic charges and higher order atomic multipoles cannot be defined unambiguously in a molecule, and therefore, several population schemes (PS) have been devised in the last decades. In the context of density functional theory (DFT), PS based on the electron density seem to be best grounded. In particular, some groups have proposed various iterative schemes the outcomes of which are very encouraging. Modern implementations of DFT that are for example based on density fitting techniques permit the investigation of molecular systems comprising of hundreds of atoms. However, population analyses following iterative schemes may become very CPU time consuming for such large systems. In this article, we investigate if the computationally less expensive analyses of the variationally fitted electronic densities can be safely carried out instead of the Kohn-Sham density. It is shown that as long as flexible auxiliary function sets including f and g functions are used, the multipoles extracted from the fitted densities are extremely close to those obtained from the KS density. We further assess if the multipoles obtained through the Hirshfeld’s approach, in its standard or iterative form, can be a useful approach to calculate interaction energies in non-covalent complexes. Relative energies computed with the AMOEBA polarizable forced field combined to iterative Hirshfeld multipoles are encouraging.