An integro-differential equation, modeling dynamic fractional order viscoelasticity, with a Mittag-Leffler type convolution kernel is considered. A discontinuous Galerkin method, based on piecewise constant polynomials is formulated for temporal semidiscretization of the problem. Stability estimates of the discrete problem are proved, that are used to prove optimal order a priori error estimates. The theory is illustrated by a numerical example.