The large error of the CNEO-DFT harmonic analysis for the bound O-H stretch seems to indicate a suboptimal performancein this particular case, but one should note a special feature of H 4 O 2 : its bound O-H stretch strongly couples with its low-frequency modes, making the system highly anharmonic. Therefore, although CNEO-DFT harmonic analysis incorporates quantum nuclear delocalization effects and thus some single mode anharmonicity, the absence of mode coupling effects beyond the harmonic picture could be the reason for its lower accuracy in describing the bound O-H stretch. Since MD simulations provide an efficient way of incorporating mode coupling effects through finite temperature thermal motions, we performed CNEO-MD with the direct NVE -eqp method to investigate the importance of these coupling effects. Note that the water dimer system is more complex than the methanol molecule and different trajectories tend to give significantly different spectra (Figure ). Therefore, instead of using 5 trajectories as in the methanol case, we used 50 trajectories to obtain the converged IR spectrum, although we found that in fact 30 Figure shows the calculated IR spectra from CNEO-MD and AIMD with the total energy corresponding to 300 K. For comparison, we also show their harmonic analysis results with artificial peak broadening as well as the experimental peak position references. We can see that in spectra obtained from both CNEO-MD and AIMD, the bound O-H stretch blueshifts about 100 cm -1 compared with their corresponding harmonic results. This shift Here we emphasize that the simulation temperatures/energies do not mean that experiments are carried out at this temperature, instead, they should be purely perceived as a parameter that can modulate the mode coupling strength. We found that the 300 K parameter works well for both the water dimer and the protonated water clusters studied here, and it also seems to work well for other hydrogen bonded systems. with experimental references. In experimental spectra, the dominant feature is the broad progression of bands from 1800 to 2400 cm -1 with the maximum peak located at 1858 cm -1 for the bare cluster and 1878 cm -1 for the D 2 -tagged cluster. VSCF/VCI analysis assigned the most intense peak to the hydronium H 3 O + asymmetric stretch with strong coupling to some low-frequency modes, and the broad tail extended to 2400 cm -1 is attributed to many overtones and combination bands. CNEO harmonic analysis overestimates the asymmetric stretch of bound O-H of H 3 O + by 200 cm -1 (2086 cm -1 ), although it is already much better than the DFT harmonic frequency (2349 cm -1 ). With the incorporation of mode coupling effects, CNEO-MD gives further improved results over CNEO harmonic analysis with an intense peak positioned around 1900 cm -1 and shows a progression of bands extended to 2400 cm -1 . Although the relative intensity contrast between the main peak and the band progression is not as sharp as that in the experimental spectra, this significant improvement again demonstrates that CNEO-MD overall successfully reproduces most features presented in experiment and captures the strong mode coupling effects between bound O-H stretches and other low-frequency modes. In contrast, conventional AIMD predicts an overly broadened progression from 1800 to 2800 cm -1 with its most intense peak appearing around 2100 cm -1 . For the free O-H stretches around 3600 cm -1 that have limited coupling with other modes, MD barely changes relative to the harmonic analysis results, and both CNEO harmonic analysis and CNEO-MD can accurately predict their vibrational frequencies. This observation is again consistent with the water dimer case as well as a series of molecular systems our group has previously studied. 44 spectra for H 9 O + 4 is the band around 2650 cm -1 , which is assigned to the hydronium H 3 O + asymmetric stretch coupled to low-frequency modes. 12 CNEO harmonic analysis predicts the frequency of this mode to be 2585 cm -1 , which is 65 cm -1 lower than experiment. With CNEO-MD, this peak blueshifts to 2600 cm -1 , slightly closer to the experimental results.