Aromatic stacking interactions between isoalloxazine (ISA) of flavin and three prototypical aromatics (benzene, pyridine, chlorobenzene) were investigated using electronic structure calculations with Monte Carlo simulated annealing. The Effective Fragment Potential (EFP) method was used to locate the low-energy equilibrium configurations for the three dimer systems. These structures were further characterized through DFT (M06-2X) and MP2 calculations. One equilibrium configuration exists for ISA-benzene; characterizing the stacked dimer surface revealed a steep, single-welled potential that funnels benzene directly between rings II and III, positioning a substituent hydrogen adjacent to the redox-active N5. ISA-pyridine and ISA-chlorobenzene minimum-energy structures contain the aromatic ring in very similar position to that in ISA-benzene. However, the added rotational degree of freedom leads to two distinct binding motifs, having approximately antiparallel or parallel dipole moment alignment with ISA. The existence of the latter binding configuration was unexpected but is explained by the shape of the ISA electrostatic potential. Dispersion is the primary noncovalent interaction driving the positioning of aromatic rings above ISA, while electrostatics determine the orientation in dipole-containing substituted benzenes. The interplay of these interactions can be used to tune molecular recognition properties of synthetic redox cofactors, including positioning desired functional groups adjacent to the redox-active N5.