During optimization, all atoms were allowed to relax, and the optimization criteria was set to 0.01 eV/Å. Normal mode frequencies were calculated using the harmonic approximation and the finite displacement method were used to get the force constants using the Phonopy software . During frequency calculations, the atoms below the 2 nd topmost layer were frozen. The normal mode frequencies were then used to calculate the vibrational part of the partition function and used to get thermodynamic properties for the surface system. Calculations on gas phase molecules were performed by placing the molecule in a cubic cell with side 10 Å to ensure separation from the neighboring cells. For gaseous molecules, modes corresponding to translational and rotational movement were removed from the vibrations and the translational and rotational contributions to the partition function were calculated separately.