Feynman's path integral reformulates the quantum Schrödinger differential equation to be an integral equation. It has been being widely used to compute internuclear quantum-statistical effects on many-body molecular systems. In this Review, the molecular Schrödinger equation will first be introduced, together with the Born-Oppenheimer approximation that decouples electronic and internuclear motions. Some effective semiclassical potentials, e.g., centroid potential, which are all formulated in terms of Feynman's path integral, will be discussed and compared. These semiclassical potentials can be used to directly calculate the quantum canonical partition function without individual Schrödinger's energy eigenvalues. As a result, path integrations are conventionally performed with Monte Carlo and molecular dynamics sampling techniques. To complement these techniques, we will examine how Kleinert's variational perturbation (KP) theory can provide a complete theoretical foundation for developing non-sampling/non-stochastic methods to systematically calculate centroid potential. To enable the powerful KP theory to be practical for many-body molecular systems, we have proposed a new path-integral method: automated integration-free path-integral (AIF-PI) method. Due to the integration-free and computationally inexpensive characteristics of our AIF-PI method, we have used it to perform ab initio path-integral calculations of kinetic isotope effects on proton-transfer and RNA-related phosphoryl-transfer chemical reactions. The computational procedure of using our AIF-PI method, along with the features of our new centroid path-integral theory at the minimum of the absolute-zero energy (AMAZE), are also highlighted in this review.