In this paper, we use molecular dynamics (MD) simulations to investigate the origin of strain hardening in glassy polymers, with a particular focus on the influences of strain rate and temperature. We demonstrate that strain hardening in uniaxial tension arises from bond stretching and is relaxed by bond rotation, characterized by the evolution of dihedral angles. Based on this rotational relaxation, strain rate plays a role analogous to temperature, resulting in a steady state behavior in the hardening region identifiable by strain rate and temperature. Following these observations, we propose a constitutive model to describe the rate- and temperature-dependent strain hardening behavior of glassy polymers using a viscous potential as a function of strain rate. To capture the elastic and yield behavior at small strains, we integrate this viscous potential with a conventional elasto-viscoplastic (EVP) model, giving rise to our V-EVP model. The V-EVP model is shown to be thermodynamically consistent and is validated by comparison to MD simulation results of glassy polymers undergoing uniaxial tension tests across a broad range of temperatures and strain rates.