In this study we investigate π-stacking interactions of a variety of aromatic heterocycles with benzene using dispersion corrected density functional theory. We calculate extensive potential energy surfaces for parallel-displaced interaction geometries. We find that dispersion contributes significantly to the interaction energy and is complemented by a varying degree of electrostatic interactions. We identify geometric preferences and minimum interaction energies for a set of 13 5- and 6-membered aromatic heterocycles frequently encountered in small drug-like molecules. We demonstrate, that the electrostatic properties of these systems are a key determinant for their orientational preferences. The results of this study can be applied in lead optimization for the improvement of stacking interactions, as it provides detailed energy landscapes for a wide range of coplanar heteroaromatic geometries. These energy landscapes can serve as a guide for ring replacement in structure-based drug design.