The stable structures and melting behavior of Hgn clusters, 2 ≤ n < 60, have been theoretically investigated using an updated diatomics-in-molecules (DIM) model initially proposed by Kitamura [Chem. Phys. Lett.2006, 425, 2056]. Global optimization and sampling at finite temperature are achieved on the basis of hierarchical and nested Markov chain Monte Carlo methods, respectively. The DIM model predicts highly symmetric icosahedral global minima that are generally similar to the standard van der Waals atomic clusters, without any indication of distorted or low-coordinated geometries, but also at variance with the global minima found with the pairwise Hg2 potential. The combined influences of surface and many-body effects due to s–p mixing are considerable on the melting point: although the model predicts a bulk melting temperature in fair agreement with experimental results, it is found to decrease with increasing cluster size.