For the excited state absorptions, conventional MD cannot distinguish them from fundamental transitions since the underlying theoretical foundation for spectra calculations using MD is a harmonic approximation, in which the fundamental 0 → 1 excitation and all excited state n → n + 1 transitions share the same transition energy. Quasi-or semi-classical methods can possibly describe these excited state absorptions by giving more initial energy to the nuclei. However, their ZPE leakage problems could be more severe because of the even higher initial energy given to the systems. Currently, the commonly used method for describing the vibrational excited state absorptions in the 2DIR spectrum is to first obtain the PES for the mode of interest while fixing all other modes , and then either fit the PES to a Morse potential and approximate the results with the exact quantum solutions for the Morse potential, or utilize the discretized variable representation (DVR) method to solve the Shrödinger equation on this one-mode PES. Both of these two ways can yield qualitatively correct results but the quantitative accuracy can be questionable due to the approximations made during the calculations. Therefore, it remains highly desirable to develop a method that can accurately and efficiently simulate vibrational excited state absorptions.