Later, Kirchner and co-workers proposed the direct use of equilibrated NVT trajectories to calculate vibrational spectra, which we will refer to as the direct NVT scheme hereafter. This direct NVT scheme greatly reduced the simulation cost and worked well for bulk systems, although for gas phase small molecules with small couplings between a limited number of degrees of freedom, it often needs an impractically long time for modes to exchange energy and to reach equilibrium. Furthermore, as has been shown in model systems, the direct NVT scheme may introduce artifacts caused by thermostats into vibrational spectra. In addition to the NVT-NVE and direct NVT schemes, there also exist simulations that directly uses a microcanonical ensemble to calculate vibrational spectra. In these NVE simulations, the initial mode energies were either chosen to be zero-point energies (NVEzpe) or obtained through normal mode sampling (NVE -nms) with a total thermal energy corresponding to the simulation temperature. These treatments, which we refer to as the direct NVE scheme, were initially used for simulating reaction dynamics 38 but later used for vibrational spectrum calculations in an ad hoc manner. Furthermore, the quasiclassical NVE -zpe simulations that use zero-point energies as initial mode energies often lead to the zero-point energy leakage problem, causing the system to break apart due to mode energy transfer in the simulations. Recently, our group developed a molecular dynamics framework based on constrained nuclear-electronic orbital density functional theory (CNEO-DFT) to efficiently incorporate nuclear quantum effects, especially zero-point effects due to quantum delocalization, into molecular simulations. In CNEO-DFT, both electrons and key nuclei are treated quantum mechanically, but the classical molecular geometry picture is retained using the expectation values of position operators for quantum nuclei together with the positions of classical nuclei. This treatment is justified by the physical intuition that quantum nuclei are still relatively localized in space in most chemical systems of interest. Therefore, in CNEO-DFT, constraints on the expectation values of quantum nuclear positions are imposed, which leads to effective PESs that inherently incorporate nuclear quantum effects.