Many theoretical methods have been developed to compute vibrational spectra. Based on pre-calculated high quality potential energy surfaces (PESs), there are highly accurate methods such as vibrational self-consistent field/configuration interaction (VSCF/VCI), multiconfiguration time-dependent Hartree (MCTDH), and diffusion Monte Carlo (DMC). These methods have been applied to a series of challenging molecular systems and provided accurate predictions of vibrational frequencies and successful identifications of subtle features in vibrational spectra. In addition, there are path-integral based methods that can calculate vibrational spectra with reasonable accuracy and less computational cost. Examples often include centroid molecular dynamics (CMD), ring polymer molecular dynamics (RPMD), thermostatted RPMD (TRPMD) 17 and quasi-centroid molecular dynamics (QCMD). Although CMD may suffer from curvature problems and vibrational spectra calculated by RPMD may be contaminated by internal ring-polymer vibrational modes, recently developed QCMD and TRPMD can mitigate these problems. Combined with accurate force fields or on-the-fly ab initio calculations, these path-integral based methods have been applied to gas phase and condensed phase systems and produced vibrational spectra that agree well with experimental results, although in general the direct combination with ab initio calculations remains expensive. A less expensive method is the vibrational second-order perturbation theory (VPT2). VPT2 takes into account anharmonicities and mode coupling effects through a perturbative expansion around the local stationary geometry, and it has been very successful in predicting the vibrational frequencies of many molecular systems. However, VPT2 may face challenges in highly anharmonic systems such as shared proton systems. In addition to the aforementioned methods, currently in the field, the most widely used methods are harmonic analysis and molecular dynamics (MD) simulations based on either force fields or ab initio calculations. In harmonic analysis, the PES is assumed to be harmonic and vibrational frequencies are directly obtained from diagonalizing the mass-weighted Hessian matrix at a stationary geometry. In practical calculations, empirical scaling parameters are often applied to the harmonic frequencies to account for the lack of anharmonic and nuclear quantum effects, thereby reducing the deviation from experimental results. Despite the popularity of empirical scaling, the choice of scaling parameters is ad hoc and can therefore vary significantly depending on systems, vibrational modes, and underlying electronic structure methods. Compared to harmonic analysis, MD simulations, especially ab initio MD (AIMD) simulations, can often give better vibrational spectra with more detailed information. They can sample the anharmonic region of PESs as well as incorporate coupling effects between different vibrational modes through thermal motions. A commonly used procedure to obtain vibrational spectra with AIMD is first to thermally equilibrate the system with a canonical ensemble (NVT ) simulation, subsequently perform many microcanonical (NVE ) simulations starting from uncorrelated frames picked from the equilibrated NVT trajectory, and perform a Fourier transform on obtained NVE autocorrelation functions before taking an average to obtain the spectrum. Unfortunately, this procedure, which we will refer to as the NVT-NVE scheme hereafter, requires numerous NVE simulation runs to reach sufficient Boltzmann sampling and obtain converged spectra.