We have used path-integral Monte Carlo ͑PIMC͒ to simulate molecular hydrogen on graphite at submonolayer coverage. First we use a flat substrate and we study the first layer for various values of the coverage up to layer completion. We found that the first layer has a solid-gas coexistence phase at low densities and a triangular solid phase at and above the equilibrium density 0 ϭ0.0705 Å Ϫ2 . We also determine that the first layer promotion coverage is at 0.094 Å Ϫ2 in agreement with experiment. Second we introduce the full H 2 -graphite interaction, i.e., we include the effects of substrate corrugations. In this case we carry our PIMC simulations on a variety of systems at and below the 1/3 coverage. We calculate the energy as a function of coverage, contour plots of the molecule probability distribution, the pair distribution function, the static structure function and the specific heat. When the substrate corrugation part of the interaction is included we find that at 1/3 coverage the system is in a ͱ3ϫͱ3 commensurate solid phase. At coverages below that and at low enough temperature the system exists in solid clusters surrounded by vapor. At coverages below a critical density, defining a tricrical point, as the system is heated up these clusters melt into a uniform fluid phase. We find that below the commensurate density and above the tricritical point, as the clusters are heated up, first they undergo a transition into a phase where the vapor phase disappears and a commensurate phase with vacancies arises. This commensurate solid melts at higher temperature into a uniform fluid phase.