In heterogeneous catalysis, the determination of active phases has been a long-standing challenge, as materials' properties change under operational conditions (i.e. temperature (T) and pressure (p) in an atmosphere of reactive molecules). As a first step towards materials design for methane activation, we study the T and p dependence of the composition, structure, and stability of metal oxide clusters in a reactive atmosphere at thermodynamic equilibrium using a prototypical model catalyst having wide practical applications: free transition metal (Ni) clusters in a combined oxygen and methane atmosphere. A robust methodological approach is employed, where the starting point is systematic scanning of the potential energy surface (PES) to obtain the global minimum structures using a massively parallel cascade genetic algorithm (cGA) at the hybrid density functional level. The low energy clusters are further analyzed to estimate their thermodynamic stability at realistic T, p_O2 and p_CH4 using ab initio atomistic thermodynamics (aiAT). To incorporate the anharmonicity in the vibrational free energy contribution to the configurational entropy, we evaluate the excess free energy of the clusters numerically by a thermodynamic integration method with ab initio molecular dynamics (aiMD) simulation inputs. By analyzing a large dataset, we show that the conventional harmonic approximation miserably fails for this class of materials, and capturing the anharmonic effects on the vibration free energy contribution is indispensable. The latter has a significant impact on detecting the activation of the C–H bond, while the harmonic infrared spectrum fails to capture this, due to the wrong prediction of the vibrational modes.