Elucidation of the mechanism for optical spin initialization of point defects in solids in the context of quantum applications requires an accurate description of the excited electronic states involved. While variational density functional calculations have been successful in describing the ground state of a great variety of systems, doubts have been expressed in the literature regarding the ability of such calculations to describe electronic excitations of point defects. A direct orbital optimization method is used here to perform time-independent, variational density functional calculations of a prototypical defect, the negatively charged nitrogen-vacancy center in diamond. The calculations include up to 512 atoms subject to periodic boundary conditions and the excited state calculations require similar computational effort as ground state calculations. Contrary to some previous reports, the use of local and semi-local density functionals gives the correct ordering of the low-lying triplet and singlet states, namely 3A2 < 1E < 1A1 < 3E. Furthermore, the more advanced meta generalized gradient approximation functionals give results that are in remarkably good agreement with high-level, many-body calculations as well as available experimental estimates, even for the excited singlet state which is often referred to as having multireference character. The lowering of the energy in the triplet excited state as the atom coordinates are optimized in accordance with analytical forces is also close to the experimental estimate and the resulting zero-phonon line triplet excitation energy is underestimated by only 0.15 eV. The approach used here is found to be a promising tool for studying electronic excitations of point defects in, for example, systems relevant for quantum technologies.