We have theoretically designed five different m-phenylene coupled high-spin bis-heteroverdazyl diradicals and their analogous p-phenylene coupled low-spin positional isomers. The geometry-based aromaticity index, harmonic oscillator model of aromaticity (HOMA) values for both the couplers (local HOMA), and the whole diradicals (global HOMA) have been calculated for all the diradicals. We also qualitatively relate these HOMA values with the intramolecular magnetic exchange coupling constants (J), calculated using a broken symmetry approach within unrestricted density functional theory framework. Structural aromaticity index HOMA of linkage specific benzene rings in our designed diradical systems shows that the aromatic character depends on the planarity of the molecule and it controls the sign and magnitude of J. The predicted J values are explained on the basis of spin polarization maps, average dihedral angles, and magnetic orbitals. The effect of the spin leakage phenomenon on magnetic exchange coupling constant and that on HOMA values of certain phosphaverdazyl systems has been explicitly discussed. In addition, a similar comparison is made between the calculated exchange coupling constants and corresponding HOMA values. The main novelty of this work stands on the consideration of the aromatic behavior by means of the geometrical index HOMA. We also estimate another aromaticity index, nucleus independent chemical shift (NICS) values for the phenylene coupler in each diradical to measure aromaticity and compare its value with that of HOMA. The ground state stabilities of these diradicals have also been compared.