π−π interactions in heteroaromatic systems are ubiquitous in biological systems. In the present study, stabilization energies of stacked and hydrogen-bonded dimers of N-heteroaromatic systems (pyridine, pyrazine, sym-triazine, and sym-tetrazine) have been computed using a benchmark quality coupled cluster through the perturbative triples (CCSD(T)) method at the estimated complete basis set (CBS) limit. In the case of stacking, monomer units are found to be stacked in parallel planes with displaced geometries. The stabilization energies for the most stable stacked geometry of pyridine, pyrazine, sym-triazine, and sym-tetrazine dimers are found to be −3.39, −4.14, −4.02, and −3.90 kcal/mol, respectively at the est. CCSD(T)/CBS level of theory, which is clearly larger than the stabilization energy for the most stable geometry of the benzene dimer. In the case of spreading, hydrogen bonded dimers and trimers are stabilized by weak C−H···N interactions. The stabilization energies for the stacked and the spread out complexes are found to be comparable. The stabilization energy for the trimers is computed using the MP2, MP3, and B3LYP-D methods. The present study is aimed at unraveling the basis of preferred conformations of N-heteroaromatic dimers. These model systems explain partly the stability of double helical DNA and RNA structures that are formed by stacking and hydrogen bonding between nucleic acid bases.