The semiclassical (SC) initial value representation (IVR) provides a way for including quantum effects into classical molecular dynamics simulations. Implementation of the SC-IVR to the thermal rate constant calculation, based on the reactive flux correlation function formalism, has two major obstacles: (1) the SC integrand may be highly oscillatory with respect to the initial phase space variables; and (2) matrix elements of the Boltzmannized flux operator, which are crucial in generating the initial (or final) distribution for the SC trajectories, are generally not available in analytic forms. In this paper, we present practical ways of overcoming these two barriers for the SC calculation of thermal rate constants. For the first problem, we show that use of a symmetric flux-flux correlation function, together with the generalized Filinov transformation technique, can significantly smooth the corresponding SC integrand and make the calculation practical for quite large systems. For the second problem, we propose a general method for evaluating matrix elements of the Boltzmannized flux operator "on-the-fly," based on the combination of the imaginary-time path integral technique with the Metropolis random walk algorithm. Using these approaches, it is shown that thermal rate constants can be obtained for systems with more than 100 degrees of freedom, as well as for reactions in the deep tunneling regimes where quantum effects are significant.