Thermal rate constants for the H + D2 reaction on the LSTH potential-energy surface are determined quantum mechanically over T = 300-1500 K using the quantum flux-flux autocorrelation function of Miller [J. Chem. Phys. 61, 1823 (1974)]. Following earlier works [T. J. Park and J. C. Light, J. Chem. Phys. 91, 974 (1989); T. J. Park and J. C. Light, ibid. 94, 2946 (1991)], we use the adiabatically adjusted principal axis hyperspherical coordinates of Pack [Chem. Phys. Lett. 108, 333 (1984)] and a direct product C2-upsilon symmetry-adapted discrete variable representation to evaluate the Hamiltonian and flux. The initial representation of the J = 0 Hamiltonian in the L2 basis of approximately 14 000 functions is sequentially diagonalized and truncated to yield approximately 600 accurate eigenvalues and eigenvectors for each symmetry species block. The J > 0 Hamiltonian is evaluated in the direct product basis of truncated J = 0 eigenvectors and parity decoupled Wigner rotation functions. Diagonalization of the J > 0 Hamiltonian is performed separately for each K(J) block by neglecting Coriolis coupling and approximating K coupling by perturbation. Both eigenvalues and eigenvectors are corrected by the perturbation. Thermal rate constants for each J, k(J)(T), are then determined by the flux-flux autocorrelation function considering nuclear spins. Due to the eigenvector corrections, both parity calculations are required to determine k(J)(T). Overall thermal rate constants k(T) are obtained by summing k(J)(T) over J with the weight of 2J + 1 up to J = 30. The results show good agreement with experiments.