The relative phase stability of VO2 is one of the most fundamental issues concerning the metalinsulator transition in this material but has so far largely unexplored theoretically. We investigate the relative stability of various phases of VO2 using different levels of energy functionals within density functional theory (DFT). It is found that straightforward applications of several popular energy functionals, including the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional, result in a wrong prediction for the ground state of VO2. In particular, although the HSE and DFT+U methods are able to produce a band gap in the M1 phase, they strongly favor the formation of local magnetic moments, a result that clearly disagrees with experiments. We also examine the effect of the occupation and the redistribution of the d derived t2g (i.e., dxz, dyz and d x 2 −y 2 ) orbitals of V atoms on the calculated relative phase stability of VO2. We find that a small change in d-occupation can result in a drastically different theoretical prediction. With the introduction of an orbital-dependent potential, a complete separation between the d x 2 −y 2 derived valence band and dxz and dyz derived conduction bands in the M1 phase is achieved, resulting in a slight redistribution of the d occupation and a more faithful account of the polarization of the t2g orbitals. This slight rearrangement of the d occupation also leads to a relative phase stability of VO2 (including structural and magnetic phases) that agrees well with experiment.