In this work, a constitutive theory is developed, within the context of continuum mechanics, to describe the creep deformation of single crystal superalloys. The con- stitutive model that is developed here is based on the fact that as bodies deform the stress free state that corresponds to the current configuration (referred to as the "natural configuration", i.e., the configuration that the body would attain on the removal of the external stimuli) evolves. It is assumed that the material possesses an infinity of natural (or stress-free) configurations, the underlying natural configuration of the body changing during the deformation process, with the response of the body being elastic from these evolving natural configurations. It is also assumed that the evolution of the natural configurations is determined by the tendency of the body to undergo a process that maximizes the rate of dissipation. Central to the theory is the prescription of the forms for the stored energy and rate of dissipation functions. The stored energy reflects the fact that the elastic response exhibits cubic symmetry. Consistent with experiments, the elastic response from the natural configuration is assumed to be linearly elastic and the model also takes into account the fact that the symmetry of single crystals does not change with inelastic deformation. An ap- propriate form for the inelastic stored energy (the energy that is `trapped' within dislocation networks) is also utilized based on simple ideas of dislocation motion. In lieu of the absence of any experimental data to corroborate with, the form for the inelastic stored energy is assumed to be isotropic. The rate of dissipation function is chosen to be anisotropic, in that it reflects invariance to transformations that belong to the cubic symmetry group. The rate of dissipation is assumed to be proportional to the density of mobile dislocations and another term that takes into account the damage accumulation due to creep. The model developed herein is used to simulate uniaxial creep of <001>, <111> and <011> oriented single crystal nickel based su- peralloys for a range of temperatures. The predictions of the theory match well with the available experimental data for CMSX-4. The constitutive model is also imple- mented as a User Material (UMAT) in commercial finite element software ABAQUS to enable the analysis of more general problems. The UMAT is validated for simple problems and the numerical scheme based on an implicit backward difference formula works well in that the results match closely with those obtained using a semi-inverse approach.