A novel algorithm is introduced for coding all Slater determinants in the covalent space with conserved SZ, the z component of total spin S for a classical valence bond (VB) model. It effectively minimizes the search time and the storing space in the central memory of the computer. In cooperation with symmetry reductions based on molecular point group and spin inversion, the VB calculations have been extended to benzenoid hydrocarbons of up to 28 π-electrons that have 4×10 7 configurations. The low-lying states of benzenoids with 24, 26, and 28 π-electrons have been obtained for 62 species. To rationalize the aromaticity of benzenoids in a VB scheme, the resonance energy per hexagon (REPH) is defined. A linear correlation between the REPH and the energy gap of the ground (singlet) state and the first excited (triplet) state for 89 benzenoids is established.