A method employing the multi-configurational time-dependent (MCTDH) approach for the direct calculation of the thermal rate constant is presented. It is based on the flux-position correlation function of Miller et al. [J. Chem. Phys. 79, 4889 (1983)]. Eigenvalues of the thermal flux operator (F) over cap(T) = e(-<(H)over cap /2kT>) (F) over cap e(-<(H)over cap /2kT>) are calculated employing an iterative diagonalization scheme suitable for the MCTDH approach. (F) over cap(T), has only a few significant eigenvalues, which can be interpreted as contributions of the ground and excited vibrational states of the activated complex. The rate constant is calculated by propagating the eigenfunctions of (F) over cap(T) in time. As an example, the H+H-2 reaction is studied. Exact results for vanishing total angular momentum (J=0) are given. The total thermal rate constant is calculated by a J shifting scheme which takes the linear geometry of the transition state into account.