Textbook explanations for the associations of multipoint hydrogen-bonded arrays have long hinged on the secondary electrostatic interaction (SEI) model, which suggests that array association strengths depend on the proton donor (D) and acceptor (A) patterns of the interacting units. Here, computational results based on the block-localized wavefunction (BLW) method reveal limitations of the SEI model, demonstrating instead that, in the gas-phase (and in implicit chloroform solvation), the inherent free-energies of associations of multipoint hydrogen-bonded arrays correlate with the degree of "aromaticity gain" (i.e., the amount of increased cyclic π-electron delocalization) in arrays upon complexation. Excellent correlations for 46 triply (r = 0.940) and quadruply (r = 0.959) hydrogen-bonded arrays are presented.