Reproducing Deep Tunneling Splittings, Resonances, and Quantum Frequencies in Vibrational Spectra From a Handful of Direct Ab Initio Semiclassical Trajectories

A time-dependent semiclassical approach for vibrational spectra calculations is shown to describe deep tunneling splittings, resonances and quantum frequencies in multidimensional multi-well systems, by propagating a very limited number of classical trajectories. The approach is tested on ammonia by evolving eight trajectories on a full-dimensional PES. Quantum effects are reproduced and results are in good agreement with time-independent quantum calculations. All the features are maintained when ab-initio “on-the-ﬂy” dynamics is adopted, thus demonstrating that pre-computation of the PES can be avoided. The approach overcomes the typical scaling issues of quantum mechanical techniques without introducing any sim-pliﬁcations nor reductions of dimensionality of the problem. The proposed methodology is promising for further applications to systems of major complexity.

Vibrational spectroscopy is a powerful experimental tool for the identification of molecules and characterization of their internal motion. Theoretical simulation of spectra may help interpret experimental results. However, routine spectra calculation still represents a challenging task in quantum mechanics, for at least three reasons: Firstly, accurate potential energy surfaces are difficult and computationally expensive to determine; Secondly, quantum mechanical methods scale exponentially with the number of degrees of freedom of the problem; Finally, standard approximations like a normal mode treatment fail when trying to describe strong anharmonic systems or highenergy vibrational states. Our aim is to demonstrate that it is possible to develop a full-dimensional approach that does not require a pre-computed PES, with a numerical non-exponential scaling, and that provides a reliable approximation even in the case of anharmonic multidimensional multi-well potentials.
Tunneling is one of the main quantum effects related to molecular motion. In fact, quantum vibrational frequencies often present deep tunneling splittings, whose description is still an open problem for time-dependent methods. Using SC-IVR, very accurate splittings and energy levels for a van der Waals complex (the HCl dimer) described by a four-dimensional Hamiltonian were calculated. 1 However, even if this application certainly represents a non-trivial test, it is characterized by a low rotational barrier ( only ≈ 70 cm -1 ) and a relatively large splitting amplitude (≈ 16 cm -1 in the ground state, roughly 1/4 of the barrier height). Just a couple of states are under the barrier and a total of about 3,000 trajectories were needed to converge results, making the approach employed for the HCl dimer not suitable for on-the-fly calculations.
A more challenging problem for any time-dependent approach is the simulation of the vibrational spectrum of ammonia which is characterized by a much higher barrier and that is clearly in the deep-tunneling regime. As far as we are aware, previous semiclassical spectra calculations could not detect such deep tunneling splittings. 2 Approaches based on the instanton approximation are focused on tunneling splitting calculations but cannot evaluate the vibrational levels. [3][4][5] Variational time-dependent quantum mechanical techniques are usually more conveniently employed to short propagation time calculations, e.g. photodissociation or scattering processes. 6 In this manuscript, we report progress towards addressing these issues by developing the multiple coherent states time-averaging semiclassical initial value representation (MC-TA-SC-IVR) method 7-10 for multidimensional multi-well systems. In MC-TA-SC-IVR, the spectrum is computed from the autocorrelation function of a wavepacket evolved "on-the-fly". One of the attractive features of the method is that by careful consideration of initial states, a handful of trajectories is sufficient for convergence. This is promising for application to systems of increased complexity.
Unlike our previous work, 10 in which we resolved the coherent states in a reciprocal space, we resolve them in the direct space. This results in trajectories that originate in one well and are directed towards the other. Furthermore, a multi-reference initial state is introduced to characterize the double well of the ammonia potential.
From a general perspective, semiclassical methods [11][12][13] can be naturally derived by stationary phase approximation (SPA) of the Feynman path integral (PI) propagator representation. 14 They are exact quantum propagators for free particle, linear potential, and harmonic oscillator systems.
The input is a set of quantities stemming from the classical dynamics of the system such as positions, momenta, potential energies, Hessians, and classical actions. The initial value representation of the semiclassical propagator (SC-IVR) 15,16 yields a physical intuition of quantum evolution in terms of coherent states and is performed by a Monte Carlo integration over all possible trajectory initial conditions in phase space. A computationally-cheap version of the SC-IVR propagator 17 has been employed for thermal density matrices calculations of a model monodimensional double well, both isolated 18 and linearly coupled to a harmonic bath 19 , and small Argon clusters. 20 In this manuscript we show that the computational cost of the quantum propagator may be reduced to barely a few classical trajectories one, without losing the essence of the quantum effects necessary to accurately calculate quantum vibrational levels in an anharmonic multidimensional multi-well system. We employ the coherent states representation of the SC-IVR propagator due to Heller, Herman and Kluk (and later re-derived by Kay) 12,21,22 and we look at a pure quantum observable such as the wavepacket survival probability ψ (0) |ψ (t) . This is represented (for a system with F degrees of freedom) in the SC-IVR approximation by the following classical integration, for any given reference state |ψ = p eq , q eq . In eq (1), (p (t) , q (t)) is the set of 2F−dimensional classically-evolved phase space coordinates, S t is the classical action and C t is a pre-exponential factor derived (in part) from the SPA of the PI, i.e. it arises from local second-order fluctuations about the classical paths, and it is given by the square root of the determinant of the combination of the four F × F size blocks of the 2F × 2F monodromy matrix M (t) ≡ (∂ (p (t) , q (t)) /∂ (p (0) , q (0))). The Γ matrix is the coherent state matrix and defines the Gaussian width of the coherent state projection onto the configurational space. We choose Γ diagonal with elements that equal the square root of the Hessian eigenvalues at the equilibrium geometry. The oscillatory behavior of the integrand in eq (1) can be tamed and the number of Monte Carlo trajectories reduced to a few thousands by introducing a time-averaging filter 23 . Thus, the vibrational spectral density, i.e. the Fourier transform of the autocorrelation function defined in eq (1), is written as where (p (t 1 ) , q (t 1 )) and (p (t 2 ) , q (t 2 )) are positions and momenta evolved from the initial phase space point (p 0 , q 0 ) at times t 1 and t 2 , respectively, and T is the total simulation time. One can recognize in eq (3) a Fourier transform integral (t 2 ) (truncated at the simulation time T) and a time averaging one (t 1 ). The computationally-intense nested do-cycles implied by the two timeintegrals can be bypassed if one assumes that the pre-exponential factors is approximated as . This is a reasonable approximation for C t , since it has been demonstrated that it does not introduce significant errors. 23 Then, eq (3) becomes and the single time integral is now positive-definite. The approximation of eq (4) (known as the separable approximation 23 ) has been proved to be much less computationally demanding than eq (3). 2,7,8,[23][24][25] The number of trajectories required for Monte Carlo convergence in eq (3) and eq (4) is of the order of thousands per each degree of freedom. 2 Unfortunately, this amount of computational demand is out of reach for carrying out direct ab initio molecular dynamics.
The main idea behind MC-TA-SC-IVR is to place the coherent reference states composing the wavepacket |ψ in a way to maximize the overlap with the exact quantum eigenfunctions. In the case of ammonia, one can design a successful strategy by examination of a one-dimensional double well model. In Figure Figure 1 We posit that these trajectories are very representative of the actual power spectrum. If the total number of trajectories has to be reduced to a few ab initio ones, then a method able to generate the spectrum exclusively from this handful of trajectories is sought.
In the case of a multidimensional double well, we write the wavepacket in a multi-reference fashion in terms of a combination of two sets of N s /2 coherent states, placed at minima locations q i eq,1 and q i eq,2 . We distribute the momenta in a way to mimic the harmonic approximation of the vibrational spectrum of each separated well, namely (p i eq,1, j ) 2 /2m = hω 1, j n i j + 1/2 for each normal mode frequency ω 1, j of the multidimensional well located at q eq,1 . A similar corresponding procedure is adopted for momenta of the coherent states pertaining the second well. To enhance wavepacket delocalization and quantum interferences, we choose to launch trajectories on each well with opposite momenta, as indicated by the arrows on panel (c) of Figure 1. This SC-IVR choice will allow the time-evolved wavepacket of each well to overlap with the coherent states of the other well at the same evolution time. In this way, a trajectory is generated for every coherent state and, by inserting eq (5) into eq (3), the expression for the MC-TA-SC-IVR spectra calculation (before the separable approximation is introduced) becomes, where p i eq , q i eq is equal to p i eq,1 , q i eq,1 if i is odd and to p i eq,2 , q i eq,2 otherwise. In few words, the integration of eq (3) is reduced to a sum of trajectories starting from the convenient set of coherent state centers pictorially reported in Figure Figure 1. Single-well simulations 8,24,25 showed that the coherent state momenta do not need to be placed at an energy very close to the eigenvalues, because the Gaussian spreading of each coherent state is wide enough to include the peak energy shell 8,24 . The necessary quantum mechanical delocalization is provided by the presence of several coherent states on each well with energy both below and above the barrier.
To start off, we test our MC-TA-SC-IVR approach on the ammonia coupled-cluster potential energy surface of Martin, Lee and Taylor (MLT). 26 The exact quantum values were obtained by direct Hamiltonian Lanczos diagonalization. 27 These quantum results provide the benchmark for our semiclassical method. Below, we show that quantum results may be reproduced with good accuracy and much lighter computational effort. In addition, previous semiclassical results are outperformed. The semiclassical power spectrum represented in Figure Figure 2 is obtained in separable approximation employing eight classical trajectories (as given in Figure Figure 1) evolved for approximately 700 fs. To better identify each spectral peak, we enforce the A 1 symmetry into our coherent states combination of eq (5) by doubling the number of coherent states. 8,23 The A 1 vibrational levels are highlighted in Figure Figure 2 allowing us to prove the presence of tunneling splittings of the order at least of a few wavenumbers. Considering that the barrier height for the  Table Table 1 In our effort to demonstrate that MC-TA-SC-IVR can regain all quantum features and that it is fully suitable for application even when the underlying electronic problem is tackled "on-the-fly", we now turn into B3LYP/cc-pVDZ dynamics using the Q-Chem electronic structure package 28 Figure Figure 3, peaks are slightly more spaced than in Figure Figure 2 , however all vibrational features observed in Figure Figure 2 are preserved in Figure Figure 3. To prove the quantum mechanical nature of the splittings and vibrational couplings (including the Darling-Dennison resonance 29 between the nearly equal stretching   modes turned into three separate states), further calculations have been performed. A simulation of the deuterated ammonia is reported on the bottom panel of Figure Figure 3. This comparison shows that the deuteration turns splittings off, that the amount of ZPE is greatly reduced and that several peaks are shifted at lower frequencies. The deuterated ammonia spectrum resulted also to be in very good agreement with one on an accurate PES. 29 Another simulation consists in evaluat-  ing the effect of the change in the value ofh on tunneling splittings calculated using the MLT PES.
In the classical limit, we expect discretization of energy levels and quantum splittings to disappear.
In Table 2  All these considerations imply that quantum effects observed in our simulations are not an artifact of the procedure adopted and that both on PES and on-the-fly MC-TA-SC-IVR implementation are valuable tools.
Our results are less accurate for energies next to the barrier threshold. This is an expected drawback of semiclassical methods, due to the instability of trajectories at those energies and the consequent loss of accuracy in the separable approximation. The effect is anyway limited to a minor part of the spectrum and the accuracy of splittings for the few involved states is sensitive to the potential in use.
The major achievement of this work is the on-the-fly reproduction of quantum mechanical effects in vibrational spectra, by means of a time-dependent semiclassical technique, from just a handful of selected classical trajectories. Even for the difficult ammonia system, where deep tunneling, potential inversion and vibrational resonances occur, eight trajectories have been demonstrated to be enough. Besides, a comparison to previous semiclassical results for the same problem demonstrates that outcomes are greatly improved, while the computational effort drastically reduced. There are not special symmetry requirements and the approach is appropriately working on-the-fly: in this way, our time-dependent approach may offer an alternative for applications to high dimensional systems to powerful time-independent methods. 6 So far, the advantage of a time-dependent approach to the problem, able to avoid diagonalization of the Hamiltonian matrix, had not been exploited to describe the spectroscopy of fulldimensional multi-well systems, because of the time-delay introduced by the umbrella inversion.
We think that the present approach could be successfully extended to time-dependent quantum wavepacket simulations and could allow time-dependent quantum dynamics to reproduce multiwell tunneling features. In fact, time-dependent quantum wavepacket simulations with low vibrational energy are mainly confined into one well and too a long simulation is needed to consistently tunnel across the umbrella inversion barrier. Time-dependent quantum and semiclassical simulations that are mainly confined into a single well do not provide information about the spectroscopy of umbrella inversion. On the other side, vibrational energy wavepackets with umbrella inverting energy do not properly describe the under-the-barrier vibrational states. The method presented here, instead, yields simply to choose appropriate initial conditions in order to fully describe the quantum properties of the system. A final observation is in order: the name of this multi-configurational method parallels the philosophy of multi-configurational methods in firstprinciples quantum chemistry. 31 In MCSCF and CASSCF methods, it is the careful selection of appropriate molecular orbitals by human intervention that helps reduce the exponential scaling of the configuration interaction method. The MC-TA-SC-IVR method also requires human insight to select the appropriate wavepackets for propagation. As in CASSCF, wrong human choices can lead to bad results. Exploration of compressed sensing approaches to further reduce the computational cost is underway. 32,33