A flexible and polarizable interatomic potential has been developed to model hydrogen clusters interacting with one hydrogen anion, (H2)nH-, in a broad range of sizes n = 1-54 and parametrized against coupled cluster quantum chemical calculations. Using path-integral molecular dynamics simulations at 1 K initiated from the putative classical global minima, the equilibrium structures are found to generally rely on icosahedral shells with the hydrogen molecules pointing toward the anion, producing geometric magic numbers at sizes n = 12, 32, and 44 that are in agreement with recent mass spectrometry measurements. The energetic stability of the clusters is also connected with the extent of vibrational delocalization, measured here by the fluctuations among inherent structures hidden in the vibrational wave function. As the clusters grow, the outer molecules become increasingly free to rotate, and strong finite size effects are also found between magic numbers, associated with more prominent vibrational delocalization. The effective icosahedral structure of the 44-molecule cluster is found to originate from quantum nuclear effects as well, the classical structure showing no particular symmetry.