A multilevel approach is presented to assess the ability of several popular dispersion corrected density functionals (M06-2X, CAM-B3LYP-D3, BLYP-D3, and B3LYP-D3) to reliably describe two-body interaction potential energy surfaces (IPESs). To this end, the automated Picky procedure ( Cacelli et al. J. Comput. Chem. 2012 , 33 , 1055 ) was exploited, which consists in parametrizing specific intermolecular force fields through an iterative approach, based on the comparison with quantum mechanical data. For each of the tested functionals, the resulting force field was employed in classical Monte Carlo and Molecular Dynamics simulations, performed on systems of up to 1000 molecules in ambient conditions, to calculate a number of condensed phase properties. The comparison of the resulting structural and dynamic properties with experimental data allows us to assess the quality of each IPES and, consequently, even the quality of the DFT functionals. The methodology is tested on the benzene dimer, commonly used as a benchmark molecule, a prototype of aromatic interactions. The best results were obtained with the CAM-B3LYP-D3 functional. Besides assessing the reliability of DFT functionals in describing aromatic IPESs, this work provides a further step toward a robust protocol for the derivation of sound force field parameters from quantum mechanical data. This method can be relevant in all those cases where standard force fields fail in giving accurate predictions.