In the past few years, our group has developed a new method for calculating vibrational spectra based on constrained minimized energy surfaces (CMESs) . The CMES is an effec-tive PES for nuclei, but compared to the conventional PES, the CMES incorporates nuclear quantum effects, especially the zero-point effects, in the effective potential energy surface. It has been shown with a few model systems that MD based on the CMES (CMES-MD) is able to give significantly more accurate fundamental vibrational frequencies than conventional MD, and its performance is comparable to or even better than CMD and RPMD. In real molecular systems, our group has developed constrained nuclear-electronic orbital density functional theory (CNEO-DFT) to approximate the CMES. It was found that CNEO-DFT harmonic frequencies are already comparable to or better than VPT2 , and the vibrational spectra obtained from MD simulations on CNEO-DFT energy surfaces accurately reproduce the experimental spectra . Despite the great success, these past developments were mostly focused on the fundamental 0 → 1 excitations. In order to simulate non-linear vibrational spectra, especially the 1 → 2 excitations, a method that can describe vibrational excited state absorption is essential. One example is that in a 2D-IR spectrum, the 0 → 1 fundamental transitions show up on the diagonal line, whereas the 1 → 2 excited state absorptions can be observed next to them but slightly off the diagonal line, together forming the butterfly-shaped feature of the spectrum.