A first principles-based methodology for efficiently and accurately finding thermodynamically stable and metastable atomic structures is introduced and benchmarked. The approach is demonstrated for gas-phase metal-oxide clusters in thermodynamic equilibrium with a reactive (oxygen) atmosphere at finite pressure and temperature. It consists of two steps. First, the potential-energy surface is scanned by means of a global-optimization technique, i.e., a massive-parallel first-principles cascade genetic algorithm for which the choice of all parameters is validated against higher-level methods. In particular, we validate (a) the criteria for selection and combination of structures used for the assemblage of new candidate structures, and (b) the choice of the exchange-correlation functional. The selection criteria are validated against a fully unbiased method: replica-exchange molecular dynamics. Our choice of exchange-correlation functional, the van der Waals-corrected PBE0 hybrid functional, is justified by comparisons up to the highest level currently achievable within density-functional theory, i.e., the renormalized second-order perturbation theory. In the second step, the low-energy structures are analyzed by means of ab initio atomistic thermodynamics in order to determine compositions and structures that minimize the Gibbs free energy at given temperature and pressure of the reactive atmosphere.