We show how the path-integral formulation of quantum statistical mechanics can be used to construct practical ab initio techniques for computing the chemical potential of molecules adsorbed on surfaces, with full inclusion of quantum nuclear effects. The techniques we describe are based on the computation of the potential of mean force on a chosen molecule and generalize the techniques developed recently for classical nuclei. We present practical calculations based on density functional theory with a generalized-gradient exchange-correlation functional for the case of H(2)O on the MgO (001) surface at low coverage. We note that the very high vibrational frequencies of the H(2)O molecule would normally require very large numbers of time slices (beads) in path-integral calculations, but we show that this requirement can be dramatically reduced by employing the idea of thermodynamic integration with respect to the number of beads. The validity and correctness of our path-integral calculations on the H(2)O/MgO(001) system are demonstrated by supporting calculations on a set of simple model systems for which quantum contributions to the free energy are known exactly from analytic arguments.