Using hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster calculations, we investigate the stabilization of silicon and oxygen dopants in GaN. Formation energies of Si on a Ga site and O on an N site are calculated at two levels of theory using conventional thermochemical and kinetic exchange and correlation density functionals (B97-2 and BB1k). We confirm the shallow donor nature of these substitutional defects. We find that the 0/1+ transition levels for both Si and O species lie well above the bottom of the conduction band, in agreement with previous supercell-based simulations. The origin of this artifact is discussed in the context of relevant experimental results and we show how correct in-gap shallow levels can be ascertained in good agreement with experiment.