Accurate charge densities are essential for reliable electronic structure calculations because they significantly impact predictions of various chemical properties and, in particular, according to the Hellmann–Feynman theorem, atomic forces. This study examines the accuracy of charge densities obtained from different density functional theory (DFT) exchange–correlation functionals in comparison with coupled cluster calculations with single and double excitations. We find that modern DFT functionals can provide highly accurate charge densities, particularly in case of meta-generalized gradient approximations and hybrid functionals. In connection with Gaussian basis sets, it is necessary to use the largest basis sets available to obtain densities that are nearly free of basis set errors. These findings highlight the importance of selecting appropriate computational methods for generating high-precision charge densities, which are, for instance, needed to generate reference data for training modern machine learned potentials.