Despite we have demonstrated encouraging progress to obtain the correct VSC resonance condition ω c = ω 0 with a variety of MQC simulation approaches, the numerical results of k/k 0 overestimate the exact results by 3 to 8 times, depending on the detailed MQC methods. Note that the current quantum and classical subsystem partition is not unique, and one can pick and choose which DOF can be further included inside the quantum subsystem. One can in principle, use the Eq. 2 as the Hamiltonian, and include {q c , pc } inside the quantum subsystem by using Fock state basis. This will in principle bring additional accuracy because it also treats the cavity mode quantum mechanically in MQC simulations. In a photoinduced proton-coupled electron transfer (PI-PCET) model, we have found that this treatment of the proton coordinate indeed brings additional accuracy compared to the classical treatment of the same coordinate. Alternatively, one can quantize qc as a ring polymer (imaginary time path integral) and still keep the key vibrational states described at the state level and use non-adiabatic RPMD to simulate the dynamics. In a polariton mediated electron transfer (PMET) model, we found that this treatment (quantizing qc as ring polymer and describe electron donor and acceptor at the state level) gives very accurate rate constant for PMET rate constant. Other more accurate trajectory-based non-adiabatic dynamics methods, such as spin-mapping variable-based partial-linearized density matrix method (spin-PLDM) or mapping-based surface hopping approach should also be able to generate more accurate results. Further, instead of computing population dynamics and then fitting it to extract the rate constant, one can directly obtain the rate constant by computing the flux-side correlation function using trajectory-based methods. The appealing computational cost of all MQC-based methods will enable large-scale simulations that involve many molecules collectively coupled to many cavity modes. The actual experimental condition involves N ≈ 10 6 ∼ 10 12 molecules collective coupled to many cavity modes that obey the dispersion relation ω k ∝ k 2 ∥ + k 2 ⊥ for a Fabry-Pérot cavity, where k ∥ is the continuous in-plane photonic momentum and k ⊥ is the quantized photonic momentum that is perpendicular to the mirror surface. Experimentally, only when ω k = ω 0 at k ∥ = 0 is satisfied (commonly referred to as the normal incidence) can one observe VSC modification on the rate constant. For a given finite k ∥ , it is possible for ω k = ω 0 , but there will be no apparent VSC effect. Currently, only a few theories can tentatively explain such a resonance condition at the normal incidence. The MQC approaches are thus ideal for simulating this realistic scenario of many molecules coupled to many photon modes inside the Fabry-Pérot cavity and can be used to directly investigate the fundamental mechanism of VSC-modified rate constants.