  Aiming towards materials design for methane activation, we study temperature
($T$), pressure ($p$) dependence of the composition, structure, and stability
of metal oxide clusters in a reactive atmosphere using a prototypical model
catalyst having wide applications: free transition metal (Ni) clusters in a
combined oxygen and methane atmosphere. A robust methodological approach is
employed, to show that the conventional harmonic approximation miserably fails
for this class of materials and capturing anharmonic effects to the vibration
free energy contribution is indispensable. To incorporate the anharmonicity in
the vibrational free energy, we evaluate the excess free energy of the clusters
numerically by thermodynamic integration method with hybrid density functional
theory and {\em ab initio} molecular dynamics simulation inputs. We find that
the anharmonic effect has a significant impact in detecting the activation of
C$-$H bond, whereas the harmonic infrared spectrum completely fails due to the
wrong prediction of the vibrational modes.
