On these effective PESs, molecular dynamics, termed CNEO-MD, can be performed, and our group has found that with essentially the same computational cost as conventional DFT-based AIMD (less than 15% more expensive as benchmarked, see Supporting Information of ref 28), CNEO-MD significantly outperforms conventional AIMD in describing vibrational modes with large hydrogen-motion character. However, like conventional AIMD, CNEO-MD may still suffer from the insufficient Boltzmann sampling problem either from the NVT-NVE or the direct NVT scheme. Additionally, in the past, CNEO-MD was mostly performed on simple small molecules, and it awaits further tests on more complicated systems with much stronger mode coupling effects. In this paper, we employ a direct NVE scheme for gas-phase vibrational spectra simulations. We first provide a brief theoretical justification on the use of direct NVE schemes in simulating vibrational spectra and propose a slightly different energy equipartition variant for more efficient spectrum convergence, which we denote as the NVE -eqp method. We find that compared to NVT-NVE and direct NVT schemes, the direct NVE -eqp method can obtain converged spectra with much less simulation time. Furthermore, we calculate the vibrational spectra of three highly challenging water clusters: the neutral water dimer (H 4 O 2 ), the protonated water trimer (H 7 O + 3 ), and the protonated water tetramer (H 9 O + 4 ) and demonstrate the excellent performance of CNEO-MD in combination with the direct NVE -eqp method.