Hybrid QM/MM is now the method of choice for the studies of enzymatic reactions and chemical reactions in the condensed phase. However, the application of QM/MM methods is often plagued with extremely poor computational scaling of the ab initio QM methods, as well as the long time scales of molecular dynamics propagation that are usually required before any essential dynamic processes for the degrees-of-freedom (DoF) orthogonal to the boosted one can be observed. Furthermore, exploring the nuclear quantum effect (NQE), such as tunneling, nuclear delocalization phenomena, and zero-point energy of molecules, gains in popularity over the years. One of the most appealing approaches to incorporating NQE is path integral molecular dynamics (PIMD). PIMD is based on the isomorphism to approximately transform the representation of quantum effect in single-particle systems into a classical system with discrete finite number of pseudo-particles (or beads) along a cyclic path. Being one of the most commonly used simulation methods, PIMD can be used to obtain thermostatistical equilibrium properties and dynamic properties, with free energy profiles of reactions and reaction rate coefficients as two representative examples. However, applications of PIMD simulations in the studies of chemical reactions of complex systems are extremely demanding. On the one hand, PIMD requires P times more computational time than that of the classical molecular dynamics (MD), where P is the number of beads, onto which each quantum particle is mapped. P determines the convergence of PIMD simulations, and it needs to satisfy P > hω max /k B T , where ω max is the maximal vibrational frequency of the system. P surges as ω max increases or T decreases. It is thus computational expensive to use PIMD to study a reaction in a sufficiently long timescale at an acceptable accuracy with convergence-ensured number of beads. To reduce the computational cost, some progresses have been achieved. Tuckerman et al combined a noncanonical transformation on the quadratic part of the action with multiple time scale integration techniques, and proposed a very efficient algorithm. Markland and Manolopoulos proposed an approach to reduce the computational effort by separating the short-range interaction from the long-range one. Cheng et al proposed the multiple-timestep molecular dynamics (MTS-MD) algorithm to accelerate the propagation by dividing the force into slowly oscillating part and quickly varying component, but the former factor is still the efficiency bottleneck. Large integration time step can also be made possible by applying proper transformations, so the analytic integration over the high-frequency modes is allowed. After the transformation, the number of the beads can be reduced via ring-polymer contraction, in which some highest normal modes can be ignored. Marsalek and Markland proposed a ring polymer contraction approach that can dramatically reduce the computational cost of PIMD at an ab initio level. On the other hand, the accessible simulation time is further limited when it is necessary to use high-level electronic structure methods such as post-Hartree-Fock theories and density functional theory (DFT). Approximated levels of theory, such as semi-empirical (SE) methods, can improve computational efficiency and allow sufficient exploration in the phase space for systems in condensed phase. However, these methods are crude in the treatment of electronic structures, which may lead to inexact results for certain physical properties.