Using the Burt–Foreman envelope function theory and effective mass approximation, we develop a theoretical model for an arbitrary number of interacting donor atoms embedded in silicon which reproduces the electronic energy spectrum with high computational efficiency, taking into account the effective mass anisotropy and the valley–orbit coupling. We show that the variation of the relative magnitudes of the electronic coupling between the donor atoms with respect to the valley–orbit coupling as a function of the internuclear distance leads to different kinds of spatial interference patterns of the wavefunction. We also report on the impact of the orientation of the diatomic phosphorus donor molecular ion in the crystal lattice on the ionization energy and on the energy separation between the ground state and the lowest excited state.