Multiple Coherent States for First-Principles Semiclassical Initial Value Representation Molecular Dynamics

A multiple coherent states implementation of the semiclassical approximation is introduced and employed to obtain the power spectra with a few classical trajectories. The method is integrated with the time-averaging semiclassical initial value representation to successfully reproduce anharmonicity and Fermi resonance splittings at a level of accuracy comparable to semiclassical simulations of thousands of trajectories. The method is tested on two different model systems with analytical potentials and implemented in conjunction with the ﬁrst-principles molecular dynamics scheme to obtain the power spectrum for the carbon dioxide molecule. © 2009 American Institute of Physics . (cid:3) DOI: 10.1063/1.3155062 (cid:4)


I. INTRODUCTION
One of the goals of quantum nuclear molecular dynamics simulations is that of simulating the interplay between electronic properties and quantum nuclear effects. Electronic properties are local in nuclear configuration state, whereas nuclear wave functions are nonlocal for the same coordinates. In the first-principles ͑FP͒ approach to quantum molecular dynamics, nonlocal properties of the nuclear Schrödinger equation have been to some extent localized and are calculated directly by quantum chemistry methods. The main advantage of this method is that the interplay between the electronic structures and nuclear motion is reproduced within the limits of the Born-Oppenheimer ͑BO͒ approximation. From the computational standpoint, the use of FP methods allows avoiding the use of precomputed potential energy surfaces ͑PESs͒. Fitting PESs to empirical or to calculated electronic structure data is usually a tradeoff between the desired accuracy and the human and computational effort required to construct them and it is often biased by the functional forms chosen. In FP simulations, the interatomic potentials and derivatives are calculated directly from electronic wave functions, i.e., on-the-fly, thus avoiding any artifact resulting from fitting as well as the cumbersome fitting process, which gets formidable as the degrees of freedom of the system increase. On-the-fly simulations are sometimes advantageous when describing rapidly changing regions of the PES, such as bond-breaking processes, as fitting errors can be larger in these cases. Simulations of the dynamics of floppy molecules can benefit from this approach as well. FP, on-the-fly simulations have been employed very successfully for studying conical intersections. 1 To date, there have been several implementations of FP approaches for the treatment of nuclear molecular dynamics with or without quantum effects. In a typical FP molecular dynamics simulation following Newton's equations and therefore not taking quantum effects into consideration, a velocity autocorrelation function is calculated and it is Fourier transformed to obtain the vibrational frequencies. In this classical dynamics approach, anharmonic effects are properly accounted for, but the zero point energy and tunneling effects are not taken into account in the classical formalism. Furthermore, long simulations are necessary in order to obtain a good frequency resolution. [2][3][4] Several techniques have been developed to speed up classical simulations. For example, in the BO molecular dynamics approach, the electronic structure calculations for a given simulation step are converged based on previous step information. However, this approach can lead to systematic energy drifts. 5 Alternatively, extended Lagrangian molecular dynamics approaches ͑ELMD͒ [6][7][8][9] involve the classical propagation of nuclear and electronic degrees of freedom on an equal footing by assigning classical variables to the electronic ones. The electronic degrees of freedom, in turn, are expanded in terms of plane waves, 6 Gaussian functions, 8 or real-space grids. 9 ELMD propagation is computationally efficient, however, the resultant energy surface is not always close to the actual BO PES, as typically these dynamics oscillate around the BO PESs.
Of concern are the dependencies that have been found on the fictitious electronic masses. 8,10 In recent years, several avenues have been successfully pursued in order to include quantum effects in a FP approach to quantum molecular dynamics. 7,[11][12][13] One can cast the path-integral expression for the quantum propagator into an extended Lagrangian formalism, 11 or adopt the differential approach of the variational multiconfiguration Gaussian wave packet method, 12 where the potential is approximated to be locally harmonic.
Other approaches include a mean field approximation for the a͒ time-dependent self-consistent quantum propagation 14,15 and multiple spawning methods developed by Ben-Nun and Martinez. 13 Semiclassical methods based on classical trajectories [16][17][18][19][20][21][22][23][24][25][26][27] are naturally suitable for carrying out on-thefly molecular dynamics simulations. It is well known from the SC developments in the early 1970s 28,29 that these kinds of approximations include all quantum effects at least qualitatively and often proved to be quantitative in accurately describing molecular systems. 17,18,25,[30][31][32][33][34][35][36][37][38][39] More recently we implemented the semiclassical initial value representation ͑SC-IVR͒ of the quantum propagator into a FP molecular dynamics scheme 40 and showed how SC-IVR can be coupled with FP electronic structure approaches for classical molecular dynamics. A very recent study 41 has shown similar results for vibrational absorption spectra. By only employing a single FP classical trajectory, we obtained the vibrational power spectrum for the CO 2 molecule with 1 cm −1 spectrum resolution. However, the accuracy of this simulation worsens for the reproduction of anharmonic shifts for high-energy vibrational levels. The main goal of this paper is to develop a method that requires a small number of trajectories such that it can be still employed for FP simulations and, at the same time, can better reproduce anharmonic effects. For this purpose, we suggest the use of a fixed coherent states basis set determined by information gathered at the minimum of the PES.
In Sec. II, we briefly review the SC-IVR formulation for obtaining power spectra. In Sec. III, the fixed coherent state basis is introduced; we write the related power spectrum expression in terms of coherent states basis. In Sec. IV, we describe our tests for the accuracy of the method and its validity by carrying out simulations on well-known, realistic, analytical potentials, where SC-IVR and exact grid calculations have been performed. In Sec. V the method is applied to the FP simulation of a molecular system. Conclusions are drawn in Sec. VI.

II. FIRST-PRINCIPLES SC-IVR
In the SC-IVR method for quantum propagators calculation, the original sum over all classical trajectories at given end positions is replaced by a phase space integral in F dimensions, ϫe iS t ͑p͑0͒,q͑0͒͒/ប ͉p͑t͒,q͑t͒͗͘p͑0͒,q͑0͉͒, ͑1͒ where ͑p͑t͒ , q͑t͒͒ are the set of classically evolved phase space coordinates, S t is the classical action, and C t is a preexponential factor. One of the advantages of the IVR formulation is that it is amenable to Monte Carlo integration. In the Heller-Herman-Kluk- Kay 25,42 version of the SC-IVR, the pre-exponential factor involves mixed phase space derivatives, and by introducing a 2F ϫ 2F symplectic ͑monodromy or stability͒ matrix M͑t͒ϵ‫͑͑ץ‬p t , q t ͒ / ‫͑ץ‬p 0 , q 0 ͒͒, one calculates the prefactor of Eq. ͑2͒ as its determinant from blocks of size F ϫ F. Simplectic properties of the monodromy matrix are a good indicator of the accuracy of the classical approximate propagation and these are monitored by checking the deviation of this determinant from unity. A more restricted check can be done by monitoring the deviation of the determinant of the positive-definite matrix M T M ͑Ref. 43͒ from unity. In this work, we monitored this deviation and it never became greater than a strict threshold of 10 −6 . We found that a time step of 10 a.u. satisfied the strict monodromy matrix restrictions even for the light hydrogen atoms. The semiclassical propagator of Eq. ͑1͒ is usually calculated employing the direct product of one dimensional coherent reference states, ͗q͉p͑t͒,q͑t͒͘

͑3͒
In the previous expression, the coherent states have a fixed width, ␥ i . For bound systems, a reasonable choice for the width parameter is that of employing the width of the harmonic oscillator approximation to the wave function at the minimum of the PES where the dynamics will be carried out. In our previous experience and that of others, we found no significant dependency upon variation in the width parameter. 19,40,42 In this paper, the main application of FP-SC-IVR molecular dynamics will be the calculation of spectral density, where ͉͘ is some reference state and ͕͉ n ͖͘ and ͕E n ͖ are the collections of exact eigenstates and eigenvalues of the Hamiltonian Ĥ . The spectrum of Eq. ͑4͒ embodies all the information coming from the eigenvalues and eigenvectors; the peaks are located at the positions of the eigenvalues and their intensity is proportional to the overlap between the trial state ͉͘ and the actual eigenvectors associated with the PES. Replacing the Dirac delta function of Eq. ͑4͒ by its Fourier representation, one can more conveniently achieve the same spectrum with the following time dependent representation: 26 Here, the spectrum is written in terms of the Fourier trans-form of the survival probability. The SC-IVR approximation for the survival probability, 17,42 ͉͗e −iĤ t/ប ͉͘ = 1 ϫe iS t ͑p͑0͒,q͑0͒͒/ប ͉͗p͑t͒,q͑t͒͗͘p͑0͒,q͑0͉͒͘,

͑6͒
suggests a suitable choice for the phase-space reference state ͉͘ = ͉p eq , q eq ͘. Usually, a coherent wave packet, which mimics the ground state wave packet and its motion following Eq. ͑3͒, is chosen. The purpose of this paper, however, is to suggest an alternative expression to the single harmonic reference state that improves the accuracy of the results for excited states, while maintaining a reasonable computational effort.
In order to reduce the number of trajectories required for Monte Carlo phase space sampling in Eq. ͑6͒, one can introduce a time average integral at the cost of a longer simulation time. Insertion of the time average integral is not an approximation per se and it does not imply any ergodicity property. This is in principle exact, by virtue of Liouville's theorem. Using this approach, Kaledin and Miller 44 and Ceotto 45 obtained the time-averaging 46,47 ͑TA͒ SC-IVR approximation for the spectral density, Here, the p͑t 1 ͒ , q͑t 1 ͒ and p͑t 2 ͒ , q͑t 2 ͒ position and momentum variables evolve from the same initial conditions, but to different times, and T is the total simulation time. The first time integral is the Fourier transform of Eq. ͑5͒, while the second time integration corresponds to the time-averaging filter. This double-time integration implies that for any given trajectory from time 0 to time T, all possible intervals from t 1 to t 2 are considered. The prefactor of Eq. ͑2͒ for each interval now depends on two times, i.e., C t 2 ͑p͑t 1 ͒ , q͑t 1 ͒͒. This approach would be really advantageous for FP applications if the number of trajectories can be reduced to a few. In the following, we will introduce and describe an implementation of the reference state representation aimed to achieve this goal. In order to make Eq. ͑7͒ even less computationally demanding, its prefactor ͓Eq. ͑7͔͒ is approximated as phase difference, C t 2 ͑p͑t 1 ͒ , q͑t 1 ͒͒ = exp͓i͑͑t 2 ͒ − ͑t 1 ͒͒ / ប͔, where ͑t͒ = phase͓C t ͑p͑0͒ , q͑0͔͒͒. 44 With this "separable approximation" Eq. ͑7͒ becomes leading to a simplification of the double-time integration to a single-and positive-definite time integral. This approximation offers also the advantage of reducing the instabilities of the complex phase space integrand by means of the timeaveraging filter. The PES, gradients, and Hessians are computed at each nuclear configuration directly from a Gaussian-based density functional theory. 48 The Kohn-Sham orbitals were expanded on a nonorthogonal Gaussian basis. The computationally intense part of our calculation is the evaluation of the potential together with its higher derivatives at each time step. The computer time required to complete this calculation on a standard desktop ͑3.2 GHz Intel CPU͒ was usually a few hours for a 1 cm −1 spectra resolution for the small molecules studied here. The classical nuclear dynamics forces, have been calculated on the BO PES and the classical propagation is performed according to the velocity-Verlet algorithm, as implemented in the Q-CHEM electronic structure package. 49

III. MULTIPLE COHERENT STATES REPRESENTATION
By inspection of the definition of the power spectrum given in Eq. ͑4͒, one can develop a new strategy for SC-IVR power spectra calculation that would enhance its precision. According to this equation, the reference state should significantly overlap with each eigenstate in order to have a well defined intensity peak. In our previous work, we showed that when a single ground harmonic state is chosen as the reference, the spectral accuracy is poor for the excited states. 40 Instead, by constructing a reference state that is a combination of states that mimic the spectral eigenfunctions with the right eigenvalue spacings, one should obtain more intense spectral peaks. This is of relevance because we found by numerical simulation that usually the more intense the peaks are, the higher the accuracy obtained for a fixed-time simulation.
For illustrative purposes, a representation of the multiple reference state coherent basis is shown on Fig. 1. In Fig.  1͑a͒, the representation of a potential well and its vibrational energy levels are depicted. The power spectrum for this potential is shown on Fig. 1͑b͒. With the goal of obtaining high peak intensities and therefore improved accuracy, a set of carefully, yet systematically chosen, classical trajectories are sought. Figure 1 shows representative classical trajectories that can be associated with each vibrational level. These "eigentrajectories" are chosen to have turning points located at the values of potential energy equal to the peaks of the power spectrum. On panel Fig. 1͑c͒, the phase portrait for the eigentrajectories is presented. In light of these considerations and inspired by previous coherent state grids calculations by Davis and Heller, 50 we have chosen the reference states ͉͘ to be a combination of coherent states, located at either the black or gray circles on Fig. 1͑c͒. In the first case ͑black circles͒, the coherent states are located at the equilibrium molecular position q eq and the initial momentum is prescribed in such a way that the kinetic energy is just enough to reach the classical turning point of the targeted vibrational energy level. The harmonic approximation is employed to set these initial momenta, namely, p eq 2 / 2m = ប͑n +1/ 2͒. In the second case ͑gray circles͒, the positions of the initial coherent states are placed at the turning points q eq and the momenta are set to be a small value close to zero.
There have been several examples of coherent state grids for semiclassical applications, where the basis set is updated at each infinitesimal time step and it follows coherently the quantum wave packet. 51 Instead, in the approach proposed above, the coherent states are placed in a nonlocal fashion and their centers are kept fixed all along the simulation time.
To gain numerical insight about the multiple coherent ͑MC͒ states approach presented above, one can look at the survival probability expression in Eq. ͑6͒ in terms of position and momenta for a single coherent state, The previous equation shows that the signal for the Fourier transform is significant only if the classical trajectories happen to be on the neighborhood of ͑p eq, i , q eq i ͒. Choosing trajectories around this neighborhood is important since one can facilitate phase space Monte Carlo integration if the coherent basis states and the classical trajectories are chosen to be close to the modes associated with the desired vibrational energy peaks. As stated above, in the case of calculations of bound states, a good guess can be obtained from the evaluation of the nuclear Hessian at the equilibrium geometry of the potential well of interest. By inserting Eq. ͑10͒ into Eq. ͑7͒ one obtains the final expression for the MC states power spectra calculations where the phase space integral has been reduced to a sum over the coherent states centers. An analogous expression holds when the separable approximation is applied. In the Appendix we show how one can calculate the full harmonic spectrum starting from Eq. ͑7͒ and show how the center of the momentum distribution is associated with the spectral peak intensities.

IV. NUMERICAL TESTS ON MODEL POTENTIALS
In this section, we describe two applications of the proposed method employing model potentials. The systems studied are the carbon monoxide molecule adsorbed in the copper ͑100͒ surface and the vibrational spectra of the water molecule. These two applications allow for testing our method in well-known potentials, for which results are readily available.

A. Carbon monoxide adsorbed on the copper "100… surface
To test our ideas on a simple, yet physically relevant model system, we performed SC-IVR simulations to model the interactions of a carbon monoxide ͑CO͒ molecule chemisorbed on the copper Cu͑100͒ surface. We chose this model due to the availability of a well-tested analytical potential in the form of 52 where r C and r O are, respectively, the carbon and the oxygen Cartesian coordinates and r i are that of the ith coordinates of the copper atom. The intramolecular term V CO is a standard Morse potential, while the interaction of CO with each copper atom is described as The first term of Eq. ͑14͒ describes the oxygen-copper repulsion. The remainder of the terms in Eq. ͑14͒ represents the carbon-copper interactions. The molecule orientation is taken into account by the angle between the O-C and C-Cu bonds. The metal is represented by three layers of 36 ͑6 ϫ 6͒ copper atoms arranged according to fcc lattice structure. The molecular axis is fixed perpendicular to the surface and the resulting molecular motion is described by two degrees of freedom, the internal CO mode and the stretching mode perpendicular to the surface, as depicted in Fig. 2.
The potential parameters were adjusted to obtain the first normal mode at 2084 cm −1 and the second one at 353 cm −1 . 53 Since there is an order of magnitude difference between mode frequencies, the time-averaging filter will more likely perform better on the fast mode, while the accuracy of this filter on the slow and floppy mode remains to be verified.
We performed single and MC states calculations of the CO-Cu system to investigate the accuracy of the MC states approach to SC-IVR dynamics proposed in this work. A single-trajectory simulation using the time-averaging filter yields the data shown in the panel ͑0,0͒ of Fig. 3. The corresponding vibrational eigenvalues are reported on Table I in the column denoted as MC state ͓1͔, ͑MC͓1͔͒. For comparison, we show the result obtained employing the separable approximation of Eq. ͑8͒ as MC state/separable, MC͓1͔/S. Although these frequencies have values very close to the exact numerical ones, a closer inspection of the vibrational eigenvalue spacings shows that single-trajectory simulations take into account of anharmonicity only for the first quantum excitation; higher vibrational energy levels show a harmoniclike spacing as reported in Table II.
In order to find the source of this deficiency, we calculated a single trajectory power spectrum using different reference states. The chosen state is denoted by the quantum numbers in parenthesis in each one of the panel of Fig. 3. This figure highlights how the intensities of the spectral peaks are highly dependent on the initial state. This is of importance because a higher spectral peak usually results in higher accuracy. The value and the spacing of the eigenvalues corresponding to the vibrational reference state ͑highlighted in each panel of Fig. 3͒ are reported in Tables I and II under the columns peak-centered ͑PC͒, and peak-centered/ separable approximation PC/S. By choosing various reference states, the vibrational anharmonicity is well reproduced and the spacings obtained are almost exact. The percentage error of the calculated vibrational level spacings for these calculations was always less than 1%. Obviously it is desirable to reproduce these data with the same level of accuracy in a single spectrum. Hence, we employed the MC-TA-SC-IVR approach using the following procedure. For every desired state, a trajectory which starts at the equilibrium coordinates ͑p eq , q eq ͒ and has the appropriate initial momentum associated with the harmonic approximation for the desired mode is propagated. The results of the simulations using this approach are shown in Tables I and II under Fig. 2. Note that the peak intensity is maximized for the frequency ͑shown by arrows͒ corresponding to that of the harmonic mode associated with each initial coherent state as described in Fig. 1.

234113-5
First-principles semiclassical dynamics J. Chem. Phys. 130, 234113 ͑2009͒ 6000 cm −1 was obtained with errors of the order of 2-3 cm −1 . It is also worth noting that the MC-SC-IVR approach is able to reproduce vibrational states that are not included into the harmonic combinations of Eq. ͑10͒ with the same level of accuracy.

B. Water molecule
A more challenging test for the time-averaging MC-SC-IVR is the calculation of the vibrational spectra of the water molecule. This system has strongly coupled modes and strong anharmonic effects in the higher vibrational levels. For this study, we employed the realistic potential surface of Ref. 55.
We summarize the results of our calculations in Table  III. The spectroscopic vibrational terms and their associated symmetry have been reported on the first column. For comparison purposes, the harmonic levels are reported on the second column. The next three columns are TA-SC-IVR calculations performed with the separable approximation of Eq. ͑8͒ but at a different computational level. Under the MC͓1͔/S column, we show the calculated frequencies for a single 8000 time-step trajectory. Under the MC͓8͔/S label, the results of eight classical trajectories sampled according to the harmonic approximation of the first eight vibrational states and the MC states basis of Eq. ͑10͒ is employed. On the column labeled "separable," the results from a simulation employing 32 000 trajectories using the separable approximation and a single ground harmonic reference state are reported. The data shown in this column are converged at the semiclassical level and we employ this column as a reference for comparison. The purpose of this table is that of comparing the results of our few-trajectory MC states approach to the considerably computationally more expensive 32 000 trajectories simulation at the same level of semiclassical approximation. The results of Table III show that the MC͓8͔/S simulation reproduces with a certain accuracy the semiclassically converged vibrational levels. The MC͓8͔/S simulations correctly describe anharmonic effects, a deficiency observed for single-trajectory calculations. As in the case of the CO/Cu͑100͒ simulation shown previously, even if the ZPE is reproduced accurately at all levels of semiclassical dynamics, the vibrational energy levels spacings of higher-energy states TABLE I. Vibrational energy levels for the different internal and external stretching excitations: On the second column the exact grid calculation results. The MC states MC͓1͔/S column represents a single trajectory timeaveraging semiclassical calculation carried out with the separable approximation. The same calculation without the separable approximation is listed under the MC͓1͔ column. The peak centered with separable approximation ͑PC/S͒ and peak-centered ͑PC͒ columns correspond to the initially chosen coherent states shown in the panels of Fig. 3. The MC states simulations of six trajectories are shown in columns MC͓6͔/S and MC͓6͔ with and without the separable approximation, respectively. The six-trajectory calculations are able to reproduce the anharmonic spacings associated with each mode, as shown in Table II are different for each type of semiclassical simulation. The single-trajectory MC͓1͔/S calculation gives a harmonic bending spacing of 1607 cm −1 . In contrast, the eighttrajectory MC͓8͔/S simulation provides an anharmonic spacing of 1599, 1560, and 1537 cm −1 . This is to be compared to the spacing of 1588, 1552, and 1528 cm −1 resulting from the 32 000 trajectories calculation. A similar trend is observed for the symmetric and asymmetric stretch modes. For example, the ͑1,1,1͒ vibrational state has an energy of 13 789 cm −1 under the MC͓1͔/S approximation. The same state has an energy of 13 705 cm −1 for the more accurate eight-trajectory case MC͓8͔/S and 13 690 cm −1 for the 32000 trajectories simulation. The corresponding energy spacings from the ground state are, respectively, 9154, 9068, and 9046 cm −1 . As in the case of the CO/Cu͑100͒ simulation, the single-trajectory semiclassical calculation does not show any anharmonic contribution. The anharmonic spacing is recovered using the proposed MC͓8͔/S approach, as compared to the 32 000 trajectories case under the column labeled separable. On the last two columns, the same calculations carried out using Eq. ͑7͒ are reported. Similar conclusions can be deduced as in the case of employing the separable approximation. Once obtained the power spectrum by using a harmonic basis, one can use these results to center a new coherent state basis exactly at the peak location found. However, by doing so we did not observe any significant improvement. This suggests that an accurate choice of the coherent basis set initial momenta is not crucial for the systems studied and a reasonable guess is enough to achieve an improvement over a single reference state. These numerical findings suggest that for bound harmonic-like systems like the ones explored in this work, the choice of the MC states basis within the semiclassical framework can be employed to obtain results similar to simulations of thousands of trajectories at a fraction of the computational cost.

V. FIRST-PRINCIPLES SC-IVR APPLICATIONS
A challenging test for the MC states FP-SC-IVR method is the calculation of the full dimensional vibrational power spectrum of the CO 2 molecule. A successful method should reproduce spectral features such as degenerate bending modes, strong intermodal couplings, as well as the Fermi resonances. 56 To evaluate the accuracy of the MC-FP-SC-IVR method, the vibrational spectrum of the CO 2 molecule will be compared with numerically exact discrete variable representation ͑DVR͒ eigenvalue calculations.
These calculations were performed on a potential fitted to a set of FP points obtained at the same level of theory and are reported on Table IV under the DVR heading. 40 On the first column of the same table, we report the experimental values, so one can also evaluate how realistic is the description offered by the DFT functional. In these calculations, Cartesian coordinates have been used to describe the classical molecular motion of Eq. ͑9͒ on a PES generated by a B3LYP density functional 48 with the cc-pVDZ basis set. 58 On the second column of Table IV one can find the spectroscopic label of each vibrational state and on the third the harmonic estimate for the vibrational energy levels. The main results on Table IV are from the application of the MC states FP-SC-IVR method. The columns labeled MC͓8͔ and MC͓8͔/S list the data obtained from calculations where the first eight states are chosen as the MC states reference basis. The MC͓8͔ column shows results employing Eq. ͑7͒ with only eight on-the-fly classical trajectories. The subsequent TABLE III. Water vibrational eigenvalues at different computational levels. On the first column, the spectroscopic vibrational terms are shown. The second column lists the second vibrational energy levels using normal mode analysis. The column labeled MC͓1͔/S shows semiclassical results obtained with a single trajectory of 8000 time steps. The column labeled MC͓8͔/S lists the results obtained from an eight-trajectory calculation of 4000 time steps each using the MC states approach of Eq. ͑10͒. The column labeled "separable" shows a semiclassically converged result of 32 000 trajectories using the TA-SC-IVR approach resulting from Monte Carlo integration of Eq. ͑8͒ with a standard sampling width and 2000 time steps for each trajectory. The remainder of the columns shows results employing the full TA-SC-IVR approach of Eq. ͑7͒. Uncertain peaks are marked with ‫ء‬ . column shows results from the same type of calculation but employing the separable approximation of Eq. ͑8͒. Both columns provide evidence for the ability to account for anharmonicity and change in the energy levels due to Fermi resonances. These findings confirmed previous results on reduced dimensionality model potentials. 59 A graphical comparison of our results is shown in Fig. 4 to illustrate the accuracy of the different numerical approaches. In this figure, the energy scale is divided into three parts, each one corresponding to groups of vibrational levels associated with their Fermi resonances.
Under label ͑a͒ the harmonic levels are shown, while ͑b͒ and ͑c͒ are the MC states results with eight trajectories, respectively, for the FP-SC-IVR propagator ͑MC͓8͔͒ and for the separable approximation ͑MC͓8͔/S͒. On the last column the DVR levels are drawn. Clearly both FP-SCI-IVR cannot account for the full amount of quantum anharmonicity and Fermi repulsion, however, any harmonic degeneracy is removed in a way to mimic the correct quantum behavior. The physical structure of these splittings is also well reproduced. A detailed analysis of the vibrational energy level spacings reveals that the MC states approach account very well for anharmonic spacing in comparison with the single trajectory results of our previous work reported on Ref. 40.

VI. CONCLUSIONS
In this paper, a novel semiclassical implementation for power spectrum calculation is presented. The results of this method are comparable with those of Monte Carlo phase space integration employing thousands of trajectories. Our results are valid for bound harmonic-like potentials such as those of molecules with a few degrees of freedom. The proposed MC states, MC͓n͔ method, is a strategy for employing a basis set of a small number of coherent states, which resemble a linear combination of approximate eigenfunctions of the system. Each one of the trajectories for our basis set is chosen to resemble a target vibrational state. In this way, a multiple vibrational components signal is successfully reproduced and results in a well-defined power spectrum. We found that a harmonic MC͓n͔ ansatz provides a great improvement over a single ground harmonic reference state simulation. The method has been developed by analytical considerations that show the importance of the reference state choice in deriving the harmonic spectrum and tested on both precomputed and on-the-fly FP PESs. The proposed methodology can be extended to the use of additional trajectories. A possible extension is the use of trajectories sampled from a Gaussian distribution around the Harmonic approximation to the chosen target vibrational modes. Further investigations along this line are currently underway in our laboratories.

234113-9
First-principles semiclassical dynamics J. Chem. Phys. 130, 234113 ͑2009͒ leads to obtaining a spectral density composed of a single peak corresponding to the ground-state energy, I͑E͒ = ␦͑E − / 2͒ / 2. In summary, we have shown that if the ͉0,0͘ reference state is used, the single spectral peak obtained will be that associated with the ground-state energy. In the case of employing an arbitrary reference state for the momentum and by employing a series expansion for the relevant part of the exponent of Eq. ͑A9͒, one can derive the final expression for a harmonic power spectrum, spans the entire power spectrum and it has been obtained with a reference state of nonzero momenta. The same reasoning can be applied to reference states of the form ͉0,q eq ͘ as described in Fig. 1͑c͒. In other words, the arbitrary reference states of the type ͉p eq , q eq ͘, ͉0,q eq ͘, or ͉p eq ,0͘ exhibit a full power spectrum since they can be expanded in terms of harmonic oscillator eigenstates. In conclusion, the analytical considerations presented above highlight the relevance of the ansatz suggested in Eq. ͑10͒, where the reference state is given in terms of a set of carefully chosen coherent states.