The thermomechanical behavior of solids includes dissipative processes such as plastic deformation and fracture. The relative importance of these processes on the response of energetic materials has been a subject of study for many decades due to their significance on ignition and reaction. However, a constitutive model to simulate the anisotropy of the crack patterns and the effect of plastic deformation due to slip in energetic materials is not yet available. Finite strain thermomechanical constitutive equations that couple crystal plasticity, an equation of state, and an anisotropic phase field damage model are presented. The model is implemented in a multiphysics finite element solver and used to simulate recent experiments on β-HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) by Zaug et al. The simulations reproduce qualitatively the crack pattern and the crystal orientation dependence of the observed damage. Specifically, more damage is observed when the crystal is impacted in the (010) direction, while more plastic deformation is observed when the load is applied in the (110) direction. The present model represents a step forward to understand the interplay between plasticity and fracture in shocked β-HMX single crystals. It can be used to gain insights into temperature increase and hot-spot formation under shock.