Accurate thermal rate constants for the D + H2 reactions are determined for the Liu-Siegbahn-Truhlar-Horowitz potential energy surface over the temperature range 300-1500 K. We evaluate the rate constants via the quatum flux-flux autocorrelation function formulation of Miller [J. Chem. Phys. 61, 1823 (1974)] using the adiabatically adjusted principal axis hyperspherical coordinates of Pack [Chem. Phys. Lett. 108, 333 (1984)] and a symmetry adapted discrete variable representation used earlier for the H + H2 reaction [T. J. Park and J. C. Light, J. Chem. Phys. 91, 974 (1989)]. The initial L2 basis of approximately 15 000 functions is sequentially diagonalized and truncated, with a final reduction to approximately 420 accurate eigenvectors of the symmetry adapted (C2-upsilon) Hamiltonians for J = O. Direct products of these functions with symmetry adapted rotation functions are then used as the basis for the J > O Hamiltonians. Nuclear spin symmetries are also included. For J > O, the individual J, K(J) blocks of the Hamiltonian are diagonalized, the Coriolis coupling is neglected, and the K(J) +/- 2 coupling is included by perturbation theory. The thermal rate constants are evaluated for each total angular momentum from the flux-flux autocorrelation function. Angular momenta up to J = 25 are required to converge the rate constants at 1500 K to approximately 5%. Thermal rate constants as functions of T (and J) are presented for the D + H2 reaction and compared with experiment and other calculations. Agreement with experiment for D + H2 is excellent up to about 1000 K and remains within a factor of 2 of the experimental rate constant up to 1500 K. Thus agreement of the rates over more than four orders of magnitude is quite reasonable.