A perturbative extension of the diatomics-in-ionic-systems (DIIS) is formulated as a practical method for describing global many-body potential energy surfaces with accuracy and economy. The method is applied to HF clusters, generalizing the prior accurate DIIS treatment of the dimer [Grigorenko, Nemukhin, and Apkarian, J. Chem. Phys. 108, 4413 (1998)] to arbitrary numbers of HF molecules. The calculated geometries, energetics, and harmonic frequencies of (HF)(n), n = 2 - 6 clusters agree in detail with the available data on this well-studied system. The formulation is based on treating intermolecular interactions within perturbation theory. It is shown that second-order perturbation, which includes bimolecular excitations, is necessary and sufficient in describing the many-body potential energy surfaces with spectroscopic accuracy. The approach allows the analysis of H-bonding and its nonadditive induction and dispersion forces in terms of mixings and exchange between ground-and excited states of dimers including intra- and intermolecular charge-transfer states as well as molecular triplet states. The speed of evaluation of the potential scales is the cube of the number of molecules, providing a practical method for dynamical simulations of extended hydrogen-bonded networks.