Vibrational Spectra of Highly Anharmonic Water Clusters: Molecular Dynamics and Harmonic Analysis Revisited with Constrained Nuclear-Electronic Orbital Methods

: Vibrational spectroscopy is widely used to gain insights into structural and dynamic properties of chemical, biological, and materials systems. Thus, an efficient and accurate method to simulate vibrational spectra is desired. In this paper, we justify and employ a microcanonical molecular simulation scheme to calculate the vibrational spectra of three 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+ ). We find that with the accurate description of quantum nuclear delocalization effects through the constrained nuclear-electronic orbital framework, including vibrational mode coupling effects through molecular dynamics simulations can additionally improve the vibrational spectrum calculations. In contrast, without the quantum nuclear delocalization picture, conventional ab initio molecular dynamics may even lead to less accurate results than harmonic analysis.


INTRODUCTION
−4 Theoretical predictions and simulations of vibrational spectra can provide valuable insights into both static and dynamic properties of a system, including structural information, interactions with the surrounding environment, reaction pathways, and dynamic evolutions.
Many theoretical methods have been developed to compute vibrational spectra.Based on precalculated high-quality potential energy surfaces (PESs), there are highly accurate methods such as vibrational self-consistent field/configuration interaction (VSCF/VCI), 5,6 multiconfiguration time-dependent Hartree (MCTDH), 7 and diffusion Monte Carlo (DMC). 8,9−13 In addition, there are path-integral-based methods that can calculate vibrational spectra with reasonable accuracy and less computational cost.−20 Combined with accurate force fields or on-the-fly ab initio calculations, these path-integral-based methods have been applied to gas phase and condensed phase systems and produced vibrational spectra that agree well with experimental results, 17,25,26 although in general, the direct combination with ab initio calculations remains expensive.A less expensive method is the vibrational second-order perturbation theory (VPT2). 27VPT2 takes into account anharmonicities and mode coupling effects through a perturbative expansion around the local stationary geometry and has been very successful in predicting the vibrational frequencies of many molecular systems.However, VPT2 may face challenges in highly anharmonic systems such as shared proton systems. 28n addition to the aforementioned methods, currently in the field, the most widely used methods are harmonic analysis and molecular dynamics (MD) simulations based on either force fields or ab initio calculations.In harmonic analysis, the PES is assumed to be harmonic, and vibrational frequencies are directly obtained from diagonalizing the mass-weighted Hessian matrix at a stationary geometry.In practical calculations, empirical scaling parameters are often applied to the harmonic frequencies to account for the lack of anharmonic and nuclear quantum effects, thereby reducing deviation from experimental results.Despite the popularity of empirical scaling, the choice of scaling parameters is ad hoc and can therefore vary significantly depending on systems, vibrational modes, and underlying electronic structure methods.Compared with harmonic analysis, MD simulations, especially ab initio MD (AIMD) simulations, can often give better vibrational spectra with more detailed information.They can sample the anharmonic region of PESs as well as incorporate coupling effects between different vibrational modes through thermal motions. 29 commonly used procedure to obtain vibrational spectra with AIMD is first to thermally equilibrate the system with a canonical ensemble (NVT) simulation, subsequently perform many microcanonical (NVE) simulations starting from uncorrelated frames picked from the equilibrated NVT trajectory, and perform a Fourier transform on obtained NVE autocorrelation functions before taking an average to obtain the spectrum.30−32 Unfortunately, this procedure, which we refer to as the NVT−NVE scheme hereafter, requires numerous NVE simulation runs to reach sufficient Boltzmann sampling and obtain converged spectra.Later, Kirchner and co-workers proposed the direct use of equilibrated NVT trajectories to calculate vibrational spectra, which we will refer to as the direct NVT scheme hereafter.33,34 This direct NVT scheme greatly reduced the simulation cost and worked well for bulk systems, although for gas phase small molecules with small couplings between a limited number of degrees of freedom, it often needs an impractically long time for different vibrational modes to exchange energy and to reach equilibrium.34 Furthermore, as has been shown in model systems, the direct NVT scheme may introduce artifacts caused by thermostats into vibrational spectra.34 In addition to the NVT−NVE and direct NVT schemes, there also exist simulations that directly use a microcanonical ensemble to calculate vibrational spectra.35,36 In these NVE simulations, the initial mode energies were either chosen to be zero-point energies (NVE-zpe) or obtained through normal mode sampling (NVE-nms) 37 with a total thermal energy corresponding to the simulation temperature.These treatments, which we refer to as the direct NVE scheme, were initially used for simulating reaction dynamics 38 but later used for vibrational spectrum calculations in an ad hoc manner.35,36 Furthermore, the quasiclassical NVE-zpe simulations that use zero-point energies as initial mode energies often lead to the zero-point leakage problem, causing the system to break apart due to mode energy transfer in the simulations.39,40 Recently, our group developed an MD framework based on constrained nuclear-electronic orbital density functional theory (CNEO-DFT) 41,42 to efficiently incorporate nuclear quantum effects, especially zero-point effects due to quantum delocalization, into molecular simulations.43,44 In CNEO-DFT, both electrons and key nuclei are treated quantum mechanically, but the classical molecular geometry picture is retained by using the expectation values of position operators for quantum nuclei together with the positions of classical nuclei.This treatment is justified by the physical intuition that quantum nuclei are still relatively localized in space in most chemical systems of interest.Therefore, in CNEO-DFT, constraints on the expectation values of quantum nuclear positions are imposed, which lead to effective PESs that inherently incorporate nuclear quantum effects.On these effective PESs, MD, termed CNEO-MD, 44 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 the Supporting Information of ref 28), CNEO-MD significantly outperforms conventional AIMD in describing vibrational modes with a large hydrogen-motion character.28,44 However, like conventional AIMD, CNEO-MD may still suffer from the insufficient Boltzmann sampling problem from either 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 simulations of vibrational spectra. Wefirst provide a brief theoretical justification on the use of the direct NVE scheme 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.We find that in these systems, CNEO-MD can further improve over CNEO-DFT harmonic analysis in vibrational calculations due to the incorporation of strong mode coupling effects on top of quantum nuclear delocalization effects.In contrast, conventional DFT-based AIMD may not necessarily perform better than DFT harmonic analysis because of the lack of quantum nuclear delocalization effects.

Review of Calculating Infrared Spectra with
Classical Simulations.To introduce and justify the direct NVE scheme, we first briefly review infrared (IR) spectrum calculations with classical simulations.IR spectroscopy measures the frequency-dependent absorbance (A(ω)) of the IR waves by a sample.The measured absorbance is directly proportional to the sample's absorption cross section (σ(ω)) and thus the rate of energy loss from the radiation E ( ( )) where l is the length and V is the volume of the sample under radiation and S(ω) is the energy flux or the time-averaged amplitude of the Poynting vector.With a continuous and perturbative electromagnetic field E(t) = E 0 cos ωt, the radiation energy loss rate can be computed using Fermi's golden rule (1 e ) d e (0) ( )

Journal of Chemical Theory and Computation
where R(ω) and R(−ω) are the excitation and de-excitation transition rates, β = 1/k B T, with k B the Boltzmann constant and T the temperature of the system, and t (0) ( ) qm −47 The energy flux can be calculated by classical electromagnetic theory, and in the end, the IR spectra can be obtained with the Fourier transform of the quantum dipole autocorrelation function where n r r = is the frequency-dependent refractive index, c 1/ 0 0 = is the speed of light in vacuum, and the relative permeability μ r is often approximated to be unity.
In practice, because quantum correlation functions are usually difficult to calculate, classical autocorrelation functions obtained from classical simulations are more often used for spectrum computation.These two types of autocorrelation functions are not rigorously equal, but for the harmonic oscillator model, if the quantum and classical oscillators both reach equilibrium at the temperature T within a canonical ensemble, the Fourier transforms of quantum and classical position autocorrelation functions only differ by a frequencydependent prefactor Although this prefactor is derived from the harmonic oscillator model, it is often used to approximately account for the differences between quantum and classical position and dipole autocorrelation functions in real systems.Thus, with this harmonic approximation, the IR absorbance A(ω) can be expressed with classical dipole autocorrelation functions where B( ) . For further simplification, integration by parts can be performed, and the absorbance can be expressed with the autocorrelation function of the time derivative of the dipole vector This equation implies that in principle, in order to obtain the IR absorption spectra, one needs to generate enough sampling from the canonical ensemble, which is usually achieved by running an NVT trajectory at a given temperature and picking uncorrelated configurations to initialize subsequent NVE simulations.The final absorption spectra can be obtained by taking the average of the Fourier transformed dipole derivative autocorrelation functions with a prefactor that is essentially frequency-independent if the refractive index n(ω) does not significantly vary with frequency.This equation is the theoretical foundation for the calculations of vibrational spectra with the NVT−NVE scheme.The NVT−NVE scheme may sound easy to perform, but to achieve enough sampling from the canonical ensemble, the number of configurations needed as well as subsequent NVE trajectories is very large even in model systems. 34n order to reduce the computational cost, Kirchner and coworkers introduced the direct NVT scheme for spectrum calculations. 33,34Instead of an NVT simulation followed by many NVE simulations, the direct NVT scheme uses only one equilibrated NVT trajectory to calculate the ensemble averaged dipole derivative autocorrelation function.In principle, this scheme is not theoretically rigorous because the thermostat will affect the evolution of the system and thus bring artifacts into the spectrum.However, in practical simulations and especially those in the condensed phase, if a weak thermostat is used, the artifacts are negligible. 33.2.Direct NVE Scheme for Calculating Vibrational Spectra.In addition to the NVT−NVE and direct NVT schemes, the direct NVE scheme has also been employed occasionally for calculating vibrational spectra, 35,36,48 although it was more used in an ad hoc manner with no theoretical justification yet.Here, we provide a brief theoretical foundation for the direct NVE scheme by revisiting the autocorrelation function of the harmonic model within a microcanonical ensemble.This scheme is different from the previous two schemes that are built on calculating the dipole derivative autocorrelation functions for a canonical ensemble equilibrated at a given temperature T.
For a microcanonical ensemble of a harmonic oscillator with energy E, the classical position autocorrelation function is This autocorrelation function is very similar to the classical position autocorrelation function within the canonical ensemble at temperature T These two equations indicate that within the harmonic oscillator model, the classical position autocorrelation functions of the microcanonical ensemble and the canonical ensemble differ only by a prefactor k B T/E. Since classical autocorrelation functions from a canonical ensemble have been widely used to approximate quantum autocorrelation functions, when we take into account the subtle difference in the prefactor, the IR spectra could also be calculated with This equation seems trivially similar to eq 6, but there is a key difference: the autocorrelation function in eq 6 needs to be obtained from a canonical ensemble, whereas here, it only needs to be obtained from a microcanonical ensemble.This difference is beneficial because, for a given system, the phase space of a microcanonical ensemble with a fixed energy E is notably smaller than that of a canonical ensemble with a fixed temperature T, especially for gas phase systems with limited degrees of freedoms.As such, the sampling need for a microcanonical ensemble simulation is much smaller.However, this brief justification is for a single 1D harmonic oscillator only.For multiple harmonic oscillators, the assignment of mode energies E will lead to different variants of the direct NVE scheme.Specifically, the NVE-zpe method assigns the corresponding harmonic zero-point energy to each mode, and therefore, if the intermode energy transfer is neglected and there is no zero-point leakage issue, the intensity of each mode can be obtained by performing a Fourier transform on the dipole derivative autocorrelation function and then dividing by a prefactor based on the respective mode energy.The NVEnms method assigns energies randomly to each mode using the normal mode sampling technique. 37With sufficient trajectory sampling, this method satisfies the energy equipartition condition among each mode, and therefore, the intensities can be obtained directly by Fourier transforming the dipole derivative autocorrelation functions (see the Supporting Information for a detailed justification).However, the challenge is that the NVE-nms method may still require a large number of trajectories to reach energy equipartition, thus leading to a suboptimal sampling efficiency.Here, we propose a slightly different NVE variant in which we start the NVE simulations from configurations with energy equipartition preset (NVE-eqp), i.e., all modes begin with the same mode energy.Admittedly, this initial condition is less physical than the initial condition obtained from normal mode sampling and will artificially favor those configurations close to energy equipartition during the dynamics simulations, but we will show later that it presents higher computational efficiency, and the simulated spectra are reasonably accurate compared to experimental results in a series of highly challenging water cluster systems.
We note that in conventional NVT simulations, the simulation temperature can affect the final spectra, including both peak positions and peak intensities.This is in contrast to quantum simulations, in which only peak intensities change with temperature.The temperature dependence of peak positions in classical simulations is a known artifact, 33,34 although temperature elevation is often used to account for anharmonicities and mode coupling effects, which are absent in harmonic analysis.Similarly, since the direct NVE scheme is also based on classical simulations, the final spectra will depend on the simulation energy E. A finite energy E will naturally lead to the incorporation of anharmonic and mode coupling effects, which we will show later to be important for accurate spectrum prediction in challenging water cluster systems.Nevertheless, this finite energy does not need to be high in CNEO-MD as compared to conventional quasiclassial MD (NVE-zpe) because zero-point energies are already incorporated into the underlying effective PESs in CNEO-MD.Therefore, CNEO-MD should not be used in combination with the NVE-zpe method and thus will not suffer from the zero-point leakage problem.

Brief Review of CNEO-MD Working Equations.
In this paper, we will combine the direct NVE-eqp method with both conventional DFT-based AIMD and CNEO-MD.The CNEO-MD framework was recently developed by our group to incorporate nuclear quantum effects, especially zero-point effects due to quantum delocalization, in practical molecular simulations. 44Its working equations are 43,44,49 x p m t and where ⟨x⟩ and ⟨p⟩ are expectation values of nuclear position and momentum operators, respectively, and V CNEO (⟨x⟩) is the CNEO-DFT energy surface.Note that CNEO-MD is essentially an AIMD method using the energies and gradients of CNEO-DFT computed on-the-fly at each MD step.
Compared with conventional DFT-based AIMD, CNEO-MD only leads to a small increase in computational cost (<15% with recent implementation) but much more accurate vibrational spectra with quantum nuclear delocalization effects incorporated. 28,44,50 COMPUTATIONAL DETAILS The gas phase methanol molecule was used as an example to compare the performance of different MD sampling schemes and to demonstrate the strength of the direct NVE-eqp method.For the scheme comparison, we only performed CNEO-MD since the results from conventional AIMD will be highly similar.To further show the power of the direct NVEeqp method along with CNEO-MD, as well as the importance of vibrational mode coupling effects, we chose three highly challenging water clusters, including 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 + ) for study.It is known that the choice of electronic density functional could significantly influence the calculation of vibrational spectra, 44,51 and here we adopted the PBE0 52,53 functional for both CNEO-DFT and DFT calculations as it was previously found that PBE0 is the best functional for free X−H stretches as benchmarked against the CCSD(T) method 54,55 for a series of small organic molecules, 44 where X is a heavy atom such as carbon, nitrogen, and oxygen.However, because hydrogen bonds are not present in the small organic molecules that were benchmarked before, here, for water cluster systems, we additionally benchmarked several commonly used functionals against CCSD(T) on the O−H stretch modes in H 4 O 2 .We found that the range-separated hybrid meta-GGA functional ωB97MV 56 with VV10 dispersion correction 57 has the best performance and PBE0 is a close second, with less than 25 cm −1 deviation from the CCSD(T) reference (Table S1).Because of the relatively low computational cost and relatively high accuracy, we continued using PBE0 in this paper.We used the def2-TZVP 58 electronic basis set for methanol calculations and def2-TZVPPD 58,59 for water cluster calculations.As for nuclei, we only treated protons quantum mechanically because they feature the most significant nuclear quantum effects, but we note that a full quantum treatment for all nuclei is possible. 42The PB4D basis set 60 was used for protons.Because protons are relatively localized with little overlap between each other, they were treated as distinguishable particles, and no proton−proton correlation was included.We did not include electron−proton correlation (epc), either, because we have previously found that without an epc functional, 61−63 CNEO-DFT and CNEO-MD are already sufficiently accurate for vibrational calculations. 28,44,64We leave the investigation of the epc effects for future studies.Both DFT and CNEO-DFT harmonic analyses were performed using a locally modified version of PySCF, 65,66 which is available through our group GitHub. 67The MD simulations were performed with atomic simulation environment (ASE). 68or the direct NVE scheme, unless otherwise specified, 50 2 ps simulations were performed for simulating spectra.To initialize each trajectory, the starting geometry was the optimized geometry and a kinetic energy E = k B T was added to each

Journal of Chemical Theory and Computation
mode with random positive or negative velocity directions. 69or the direct NVT scheme, following ref 33, we employed the Nose−Hoover chain thermostat, 70−72 which has been implemented in an in-house version of ASE.First, a 10 ps simulation was performed for the system to reach equilibrium, and then, another 10 ps simulation was performed for the production run.It is known that the Nose−Hoover chain thermostat may suffer from ergodicity issues in small or strongly harmonic systems; 73 therefore, for the NVT−NVE scheme, we used the Andersen thermostat 74 for NVT sampling.A 20 ps NVT simulation was first performed, and then uncorrelated configurations from the later 15 ps were picked as initial configurations for subsequent 2 ps NVE simulations.All the MD simulations were performed under 300 K or equivalent of total energy.

Comparison of Simulation Schemes
. Figure 1 presents the power and IR spectra of a single methanol molecule at 300 K or its equivalent energy by CNEO-MD with direct NVE (including both NVE-eqp and NVE-nms variants), direct NVT, and direct NVT−NVE schemes.In the direct NVE scheme, we used five and 15 2 ps trajectories to obtain converged spectra for eqp and nms variants, respectively.As power spectra reflect the energy distribution among different vibrational modes, we can see from Figure 1a that energy is about equally distributed among different modes from the direct NVE scheme, which leads to good agreement between their IR spectra and the experimental spectrum.For the direct NVT scheme, from the power spectrum, we can see that with the same amount of total simulation time as the NVE-eqp method, energy equipartition in the direct NVT scheme is far from being achieved.Specifically for this trajectory, energies of the C−H and O−H stretch modes above 2500 cm −1 are much lower than those of the remaining modes, which leads to weak IR intensities for the highly IR active C−H and O−H stretch modes.This observation agrees well with previous literature that similar unequal energy distributions across different vibrational modes were observed for a single methanol molecule in conventional AIMD simulations. 33As to the NVT−NVE scheme, here, the spectra were obtained from an average of 30 2 ps NVE trajectories, whose initial configurations were sampled every 500 fs from an NVT trajectory.This choice of 500 fs sampling time gap is justified by the fast decay of autocorrelation functions of the NVT trajectory (Figure S1).From the power and IR spectra, it can be seen that the results from the NVT−NVE scheme are close to those of the direct NVE scheme and in good agreement with the experimental spectrum as well.
Note that here we only ran 5 trajectories for NVE-eqp and 15 for NVE-nms methods, but 30 trajectories for the NVT− NVE scheme.This is because the spectra calculated from the direct NVE scheme need fewer trajectories to converge.Figure 2 shows the convergence test results for power spectra as a function of the number of trajectories used in different methods.We can see that for the methanol molecule, at least 30 trajectories are needed for convergence within the NVT− NVE scheme.In contrast, 15 trajectories are sufficient to achieve convergence with the NVE-nms variant and only 1−5 trajectories are needed with the direct NVE-eqp variant.Although later in more complicated water clusters, more trajectories will be needed for spectrum convergence, the general trend is always that the direct NVE scheme needs significantly fewer trajectories than the NVT−NVE scheme to reach convergence, and the NVE-eqp variant is the more efficient variant.
Another difference between the direct NVE and the NVT− NVE schemes is that the peaks are generally broader in the NVT−NVE scheme.It is understandable that with the NVT− NVE scheme, there is a total energy distribution with more fluctuation, whereas the total energy is fixed in the direct NVE scheme.However, this does not mean that energy will not flow with the direct NVE scheme.Instead, through mode coupling effects, energy can exchange among different modes and also broaden their peaks.Therefore, both schemes are able to capture the mode coupling effects through a finite temperature or energy, and it is simply the fact that the fluctuation is larger in the NVT−NVE scheme that makes the peaks broader.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 artifacts.
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 the direct NVE scheme is significantly more efficient in spectrum convergence with a minimal simulation time.The direct NVT scheme may be comparable to the NVE-eqp method in computational efficiency, but it may suffer from difficulty in effectively equilibrating different vibrational modes in a single trajectory due to relatively 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, equipartition will be easily achieved and thus the direct NVT scheme will be more successful in that paradigm. 33,34

H 4 O 2 .
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 3 (coordinates are given in the Supporting Information).−83 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), PBE0, and CNEO-PBE0 harmonic analyses, along with experimental references, are presented in Table 1.Within the harmonic approximation, CCSD(T) overestimates all O− H stretch frequencies by over 150 cm −1 and PBE0 gives similar results.In contrast, although the harmonic approximation is still used, because of the incorporation of quantum nuclear delocalization effects, CNEO-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-PBE0 harmonic analysis gives a larger error and underestimates the frequency by 107 cm −1 , but it still slightly outperforms PBE0, which overestimates the frequency by 130 cm −1 .
The large error of the CNEO-DFT harmonic analysis for the bound O−H stretch seems to indicate a suboptimal performance in 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. 84,85Therefore, 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, 33 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 S3).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 trajectories are already sufficient to reach good spectrum convergence (see Figures S4−S6 for convergence test results for H 4 O 2 , H 7 O 3 + , and H 9 O 4 + , respectively). Figure 4 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 makes CNEO-MD accurately capture the frequency of the bound O−H stretch but makes the AIMD results worse in that the underlying DFT already    S2 for the temperature dependence test).Here, we emphasize that the simulation temperatures/energies do not mean that experiments are carried out at this temperature; instead, they should be 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. 86 , whose vibrational spectra have been studied extensively with high-level and computationally intensive methods such as VSCF/VCI, 10−12 MCTDH, 13 and DMC. 87Relatively inexpensive AIMD and TRMPD have also been applied to these two systems, but the spectra of the former are blue-shifted significantly and the latter suffers from artificial broadening. 88igure 5 shows the calculated IR spectra of H 7 O 3 + from AIMD and CNEO-MD along 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. 10VSCF/VCI analysis assigned the most intense peak to the hydronium H 3 O + asymmetric stretch with strong coupling to some lowfrequency modes, and the broad tail extended to 2400 cm −1 was attributed to many overtones and combination bands. 10NEO 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 .Compared to experiments 10 and the highly accurate VSCF/VCI spectra, 88 the relative intensity contrast between the main peak and the band progression is not sharp enough in CNEO-MD, which might be due to the rapid energy transfer out of the bound O−H stretch modes.However, CNEO-MD is successful overall in reproducing most features presented in  the experiment with the capture of strong mode coupling effects between bound O−H stretches and other lowfrequency 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. 44igure 6 shows the calculated IR spectra of H 9 O 4 + from AIMD and CNEO-MD along with experimental references.For this cluster, there were controversies in the past about the existence of the linear-chain Zundel isomer in the observed gas phase IR spectrum, 32 but later joint experimental and theoretical efforts concluded that only the Eigen isomer exists in the experimental spectrum. 11Therefore, we performed only CNEO-MD and conventional AIMD simulations on the Eigen isomer.The prominent feature of the experimental IR 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 lowfrequency modes. 12CNEO harmonic analysis predicts the frequency of this mode to be 2585 cm −1 , which is 65 cm −1 lower than that of experiment.With CNEO-MD, this peak blue-shifts to 2600 cm −1 , slightly closer to the experimental results.In contrast, DFT harmonic analysis overestimates the frequency by about 200 cm −1 and AIMD makes the overestimation slightly worse.Despite the great success of CNEO-MD in reproducing the main bands of the experimental IR spectra, we point out a small feature that seems absent in CNEO-MD: a diffuse feature that lies about 380 cm −1 above the dominant peak in both experimental spectra. 12,89With several high-level methods, this feature has been attributed to various overtones and combination bands, especially the combination band involving the O−O stretch mode. 12,90IMD seems to have fortuitously captured some diffuse feature at around 3200 cm −1 , but far from the experimental reference, and for CNEO-MD, this feature is not observed in the spectrum.The reason is that although CNEO-MD accounts for quantum nuclear delocalization effects, it is still based on classical MD, and it is known in the past that the accurate description of overtones, combination bands, and Fermi resonances has been particularly challenging for classical MD. 33 Therefore, these features should not be anticipated from either AIMD or CNEO-MD, although some Fermi resonance has been fortunately captured by CNEO-MD in the past. 44The frequencies of free O−H stretches can be accurately described with both CNEO harmonic analysis and CNEO-MD, which is consistent with the results from all previous investigations.
Finally, we comment that despite their geometric similarity, the two protonated water clusters have very different vibrational features, especially in the dominant bands featuring bound O−H stretches.The H 7 O 3 + spectrum is relatively more complex, and CNEO harmonic analysis greatly overestimates its bound O−H stretch frequency, whereas the H 9 O 4 + spectrum features an intense band associated with the asymmetric bound O−H stretches, and CNEO harmonic analysis gives a reasonable prediction of the frequency.We found that CNEO-MD improves much over the poor CNEO harmonic results for H 7 O 3 + but makes only a small correction on top of the already good CNEO harmonic results for H 9 O 4 + .Thus, we can see that CNEO-MD can correctly incorporate the coupling effects between different vibrational modes, whether they are strong or weak.

SUMMARY
In summary, we provided a brief theoretical justification for the direct use of microcanonical ensemble simulations for efficient computations of vibrational spectra.To further improve the efficiency, we proposed an NVE-eqp method in which all modes start with the same amount of energy.With the gas phase methanol molecule as an example, we showed that the NVE-eqp method can reach spectrum convergence faster with less total simulation time than other methods and thus is much more efficient for calculating vibrational spectra of gas phase clusters.In combination with the NVE-eqp method, we demonstrated again the power of CNEO-MD for accurate , as examples, we showed that CNEO-DFT harmonic analysis improves over DFT harmonic analysis with the incorporation of quantum nuclear delocalization effects, and CNEO-MD further improves the results with the incorporation of mode coupling effects.In contrast, DFT-based AIMD may not necessarily improve over DFT harmonic analysis due to the lack of quantum nuclear delocalization effects.Therefore, both quantum nuclear delocalization effects and mode coupling effects are major reasons for the excellent performance of CNEO-MD.As a cost-effective method based on classical simulations, CNEO-MD in combination with the NVE-eqp method holds significant promise as a robust tool for conducting vibrational spectrum simulations in more challenging systems.

Figure 1 .
Figure 1. and IR spectra of a single methanol molecule by different MD schemes as indicated.Dashed vertical lines are harmonic analysis results from CNEO-DFT.The experimental IR spectrum of gas phase methanol is obtained from the National Institute of Standards and Technology WebBook.

Figure 2 .
Figure 2. Methanol molecule power spectra as a function of the number of trajectories averaged for direct NVE and NVT−NVE schemes.The NVE-eqp method within the direct NVE scheme needs significantly fewer trajectories than does the NVT−NVE scheme to reach spectrum convergence.

Figure 4 .
Figure 4. Calculated IR spectra of H 4 O 2 by (a) DFT harmonic analysis, (b) AIMD, (c) CNEO harmonic analysis, and (d) CNEO-MD.Harmonic spectra are obtained by adding 10 cm −1 Gaussian width to each transition.The experimental references (dashed lines) are from ref 75.The orange peaks represent the bound O−H stretch, and the green peaks represent the free O−H stretch modes.

Figure 5 .
Figure 5. IR spectra of H 7 O 3 + obtained by the indicated methods.Experimental IR spectra are from ref 10.

Figure 6 .
Figure 6.IR spectra of H 9 O 4 + obtained by the indicated methods.The experimental IR spectrum of bare cluster is from ref 12 and the experimental IR spectrum of D 2 -tagged cluster is from ref 89.

Table 1 .
Vibrational Frequencies of the O−H Stretch Modes in H 4 O 2 (in cm −1 ) from Harmonic Analysis and Experiments

■ ASSOCIATED CONTENT * sı Supporting Information The
Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.3c01037.Detailed derivation and justification of the direct NVE scheme; structures of water clusters; benchmark of different electronic functionals on H 4 O 2 ; temperature dependence of H 4 O 2 spectra; spectrum convergence test; and full-range IR spectra of water clusters (PDF)