We present an anharmonic approach for molecule-surface vibrations that employs rigid body coordinates, based on the Rodrigues rotation formula, to describe curvilinear displacements (rotations and translations) of the molecule along normal modes. These displacements are used to calculate energy data points from which one-dimensional polynomial potentials are fitted using cubic splines. In these potentials, for each of the six rigid body modes separately, one-dimensional Schrödinger equations are solved with harmonic oscillator or Fourier functions (for pure rotations) as basis sets. The anharmonic vibrational energies obtained are used to calculate partition functions and from them enthalpies, entropies, and Gibbs free energies of adsorption. Our numerical implementation has been successfully tested for Morse and cosine potentials with known analytical solutions. The methods have been applied to adsorption of CH4 on the hydroxyl group of the proton form of the chabazite zeolite (H-CHA), as well as to adsorption of CH4 and CO on the Mg2+ ions of the metal-organic framework (MOF) Mg2(dobdc). To obtain the best estimates for thermodynamic functions, we include the coupling between molecule-surface vibrations and intrasystem vibrations at the harmonic level. The calculated Gibbs free energies show that the coupling is small for CH4/H-CHA and CO/MOF (between -0.7 and +0.1 kJ/mol) but substantial for CH4/MOF (-3.4 kJ/mol). The predicted anharmonic effect on the Gibbs free energy of adsorption for CH4/H-CHA, CH4/MOF, and CO/MOF is -4.7, 0.3 ± 0.7, and -2.4 ± 0.6 kJ/mol, respectively, which results in +4.2, +0.9 ± 0.7, and -0.4 ± 0.6 kJ/mol, respectively, for the deviation from experiment. This is well within chemical accuracy limits (±4.2 kJ/mol) for the adsorption of CH4 and CO in the MOF. The larger deviation for CH4/H-CHA, at the edge of the chemical accuracy range, is most likely due to contributions from soft zeolite modes which are neglected in our approach.