A method for computing anharmonic thermophysical properties for adsorbates on metal surfaces has been extended to include libration, or frustrated rotation. Classical phase space integration is used with Monte Carlo sampling of the configuration space to obtain the partition function of CO on Pt(111) and CH3OH on Cu(111). A minima-preserving neural network potential energy surrogate is used within the integration routines. Direct state counting using discrete variable representation is used to benchmark the results. We find that the phase space integration approach is in excellent agreement with the direct state counting results. Comparison with standard models such as the harmonic oscillator indicates that anharmonicity contributes significantly to the thermodynamic properties of CH3OH on Cu(111). We find that there is also a considerable difference between the harmonic oscillator and phase space integration for CO on Pt(111), although the discrepancy can largely be attributed to the presence of multiple binding sites within the unit cell. We demonstrate that a multisite harmonic oscillator model might be sufficient for CO–Pt(111). A more thorough description of the potential energy surface, which can be achieved with phase space integration, is necessary for weakly bound adsorbates such as CH3OH. The thermophysical properties were used to calculate free energies of adsorption on the respective metals, and subsequently the equilibrium constants and Langmuir isotherms in relevant temperature ranges. The results show that the choice of model to obtain partition functions greatly affects the resulting surface coverages in kinetic models.