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.