In this paper, we develop a computational model for the anisotropic plastic response of beta-HMX molecular crystals based upon experimental data. While phenomenological models of beta-HMX currently exist, they are inadequate as they cannot account for the orientation dependence of the observed response. The proposed model, on the other hand, successfully predicts the experimentally observed highly anisotropic and orientation dependent yield surfaces of beta-HMX molecular crystals having different crystal orientations. The model also correctly predicts the anisotropic plastic yielding of beta-HMX under uniaxial compression at different temperatures along different crystallographic directions. The model presented in this paper may be extended to predict the yielding and elastic plastic deformation of other energetic molecular crystals.