Vibrational spectroscopy has been a powerful tool to obtain structural information and uncover reaction dynamics in both gas and condensed phase . Compared to conventional linear vibrational spectroscopy, nonlinear vibrational spectroscopy can additionally provide more information by probing two or more photon processes. For example, two-dimensional infrared spectroscopy (2DIR) can go beyond the vibrational ground state and monitor the vibrational excited state dynamics and thus has been used to investigate reaction dynamics in a variety of materials and biological systems There have been many theoretical methods to model vibrational spectroscopy. A most widely used method to estimate vibrational frequencies is to diagonalize the mass-weighted Hessian matrices obtained from either density functional theory (DFT) or wave function theories. This method invokes the harmonic approximation and can generally provide qualitatively correct vibrational modes results; however, when describing systems with strong anharmonicities, it often needs to rely on empirical scaling factors to obtain quantitatively right results. Going beyond the harmonic approximation, vibrational second order perturbation theory (VPT2) is often used to obtain anharmonicity-corrected vibrational frequencies, which utilizes information of local higher-order derivatives. VPT2 is currently considered to have reached a good balance between accuracy and efficiency, although it may face challenges when local higher-order derivatives are not sufficient to capture the behavior of the whole potential energy surface (PES), especially in some shared proton systems with double-well PESs. Vibrational self-consistent field theory (VSCF) is currently one of the most accurate methods for obtaining vibrational spectra. It utilizes the high-dimensional PES to compute vibrational ground and excited states. However, due to the high computational cost of constructing the high-dimensitional PES, VSCF is often limited to small systems.