The noncovalent interactions of nucleobases and hydrogen-bonded (Watson–Crick) base-pairs on graphene are investigated with the DFT-D method, i.e., all-electron density functional theory (DFT) in generalized gradient approximation (GGA) combined with an empirical correction for dispersion (van der Waals) interactions. Full geometry optimization is performed for complexes with graphene sheet models of increasing size (up to C₁₅₀H₃₀). Large Gaussian basis sets of at least polarized triple-ζ quality are employed. The interaction energies are extrapolated to infinite lateral size of the sheets. Comparisons are made with B2PLYP-D and SCS-MP2 single point energies for coronene and C₅₄H₁₈ substrates. The contributions to the binding (Pauli exchange repulsion, electrostatic and induction, dispersion) are analyzed. At a frozen inter-fragment distance, the interaction energy surface of the rigid C₉₆H₂₄ and base monomers is explored in three dimensions (two lateral and one rotational). Methodologically and also regarding the results of an energy decomposition analysis, the complexes behave like other π-stacked systems examined previously. The sequence obtained for the interaction energy of bases with graphene (G > A > T > C > U) is the same for all methods and supports recent experimental findings. The absolute values are rather large (about −20 to −25 kcal mol⁻¹) but in the expected range for π-systems of that size. The relatively short equilibrium inter-plane distance (about 3 Å) is consistent with atomic force microscopy results of monolayer guanine and adenine on graphite. In graphene⋯Watson–Crick pair complexes, the bases lie differently from their isolated energy minima leading to geometrical anti-cooperativity. Together with an electronic contribution of 2 and 6%, this adds up to total binding anti-cooperativities of 7 and 12% for AT and CG, respectively, on C₉₆H₂₄. Hydrogen bonds themselves are merely affected by binding on graphene.