We report systematic periodic DFT studies of the ε and ζ phases of solid oxygen between 10 and 140 GPa using several types of exchange-correlation functionals (LDA, GGA, hybrid, and meta-GGA). Enthalpy-driven geometry optimizations with spin-unpolarized periodic DFT are done to study the pressure evolution of the (O2)4 unit cell within the experimental C/2m crystal group. At variance with previous studies using GGA functionals and plane wave basis sets, our pressure-dependent structural and energetic results are in excellent agreement with the most refined experimental data when hybrid functionals are used in conjunction with optimized atomic basis sets. The purely enthalpic description is not enough to explain the relative stability of the phases, therefore phononic calculations for both phases are done at 100 GPa. We show that inclusion of entropic and ZPE vibrational effects are crucial to correctly predict the 95-100 GPa phase coexistence region at 298 K. Our results reproduce the discontinous behavior of the unit cell parameters around the phase transition as well as the evolution of the Raman and Infra-red spectra.