Density functional theory (DFT) is the method of choice for predicting structures and reaction energies of molecular systems. However, it remains a daunting task to systematically improve the accuracy of an approximate density functional. The recently proposed many-pair expansion (MPE) [ Phys. Rev. B 2016, 93, 201108] is a density functional hierarchy that systematically corrects any deficiencies of an approximate functional to converge to the exact energy, and was shown to give accurate results for lattice models. In this work, we extend MPE to molecular systems and implement it using Gaussian basis sets. The self-attractive Hartree (SAH) decomposition [ J. Chem. Theory Comput. 2018, 14, 92-103] is employed to generate localized v-representable pair densities for performing MPE calculations. We demonstrate that MPE at the second order (MPE2) already predicts accurate molecular and reaction energies for a series of small molecules and hydrogen chains, with the EXX functional as its starting point. We also show that MPE correctly describes the symmetric bond breaking in hydrogen rings, indicating its ability to remove strong correlation errors. MPE thus provides a promising framework to systematically improve density functional calculations of molecules.