However, we cannot argue which peak width is more physical since peak broadening can be observed even in simple harmonic and Morse oscillator models, in which quantum references give sharp peaks and thus the broad peaks obtained from classical simulations are purely From the spectrum comparison and spectrum convergence results, we can conclude that both the direct NVE scheme and the NVT-NVE scheme can produce reliable spectra, but the NVE -eqp variant within direct NVE scheme is significantly more efficient in spectrum convergence with minimal simulation time. The direct NVT scheme may be comparable to the NVE -eqp method, but it may suffer from difficulty in effectively equilibrating different vibrational modes in a single trajectory due to relative weak mode couplings between a limited number of degrees of freedom in gas phase clusters. However, we note that in condensed phase systems where modes couple more strongly and interact with the surroundings more frequently, reaching equipartition will be easily achieved and thus the direct NVT scheme will be more successful in that paradigm. To further show the capability of the combination of the direct NVE scheme and CNEO-MD, we calculated the IR spectra of three water cluster systems. The structures of three water clusters to be studied are shown in Figure (Structure files are given in the Supporting Information). To begin, the water dimer H 4 O 2 is the simplest water cluster and its vibrational spectra have been studied experimentally and theoretically for decades. In H 4 O 2 , the two water molecules are bound by the intermolecular hydrogen bond, and there are three free O-H stretch modes and one bound O-H stretch mode. The vibrational frequencies of these four O-H stretch modes obtained from CCSD(T), DFT, and CNEO-DFT harmonic analysis, along with experimental references, are presented in Table . Within the harmonic approximation, CCSD(T) overestimates all O-H stretch frequencies by over 150 cm -1 and DFT with PBE0 gives similar results. In contrast, although the harmonic approximation is still used, because of the incorporation of quantum nuclear delocalization effects, CNEO-DFT with PBE0 accurately describes the vibrational frequencies of the free O-H stretches, whose errors are within 30 cm -1 of the experimental references. Interestingly, for the bound O-H stretch, CNEO-DFT harmonic analysis gives a larger error and underestimates the frequency by 107 cm -1 , but it still slightly outperforms DFT, which overestimates the frequency by 130 cm -1 .