Vibrational Tunneling Spectra of Molecules with Asymmetric Wells: A Combined Vibrational Configuration Interaction and Instanton Approach.

A combined approach that uses the vibrational configuration interaction (VCI) and semiclassical instanton theory was developed to study vibrational tunneling spectra of molecules with multiple wells in full dimensionality. The method can be applied to calculate low-lying vibrational states in the systems with an arbitrary number of minima, which are not necessarily equal in energy or shape. It was tested on a two-dimensional double-well model system and on malonaldehyde, and the calculations reproduced the exact quantum mechanical (QM) results with high accuracy. The method was subsequently applied to calculate the vibrational spectrum of the asymmetrically deuterated malonaldehyde with nondegenerate vibrational frequencies in the two wells. The spectrum is obtained at a cost of single-well VCI calculations used to calculate the local energies. The interactions between states of different wells are computed semiclassically using the instanton theory at a comparatively negligible computational cost. The method is particularly suited to systems in which the wells are separated by large potential barriers and tunneling splittings are small, for example, in some water clusters, when the exact QM methods come at a prohibitive computational cost.

Introduction

Physical systems with multiple energetically stable minima are ubiquitous in chemistry and physics.[1] Bound states that are localized in such wells, separated by potential barriers, interact via quantum tunneling, which results in observable shifts of their energies.[2,3] For equivalent, symmetry-related wells, the states that would be degenerate in the absence of tunneling produce a splitting pattern of energy levels.

Molecules and molecular complexes with two or more equivalent stable configurations are multidimensional systems that display these effects in their vibrational spectrum. The inversion of ammonia,[4] proton tunneling in malonaldehyde,[5] double proton transfer in porphycene,[6] or bond rotation in the vinyl radical[7] are examples of symmetric double-well systems that produce measurable tunneling splittings (TSs) of their vibrational state energies. Water clusters are prototype multiwell systems that exhibit nontrivial splitting patterns caused by tunneling rearrangements between many stable configurations of the cluster.[8]

The asymmetric systems, which have nonequivalent wells, have been less studied. When the state energies of different wells are in resonance, the tunneling dynamics will again cause the delocalization of the wavefunction across the wells and the energy shifts in the spectrum.[3] Away from the resonance, the states remain localized in one well. The asymmetry can be induced in symmetric molecular systems by asymmetric isotopic substitutions.[9] The normal modes and vibrational frequencies in equivalent symmetry-related potential wells then differ, and the correspondence of the vibrational wavefunctions of different wells is not preserved in general. As an example, the malonaldehyde molecule deuterated at the D7/D9 position (see Figure ) thus has an asymmetric level structure with the localized states and those that are delocalized across the two minima.[9] A mixing angle between the left–right ground vibrational states of partially deuterated (PD) malonaldehyde has been determined experimentally.[10] Further examples of the mixing have been studied in the HF–DF dimer[11] and PD vinyl radical,[12] CHD–CH, using full-dimensional exact calculations. The bifurcation splitting patterns in PD water trimers HDO(H2O)2 and D2O(H2O)2 have been determined in experiments[13] and by us using the instanton theory.[14]

Figure 4

Annotated equilibrium geometry of malonaldehyde and schematic representation of four lowest-frequency normal modes.

The asymmetry in molecules can also be found in some tautomers. In this case, the potential energy surface (PES) does not possess a symmetry relating the wells and their shapes, and the minimum energies are different. A possible candidate belonging to this class is 2-hydroxy-1-naphthaldehyde, shown in Figure . The hydroxyl proton forms a hydrogen bond with the oxygen atom of the carbonyl group and can tunnel to it to form a tautomer, which is a local minimum.

Figure 1

Tunneling tautomers of 2-hydroxy-1-naphthaldehyde.

Thiomalonaldeyde has two nearly degenerate minima in the form of enol and enethiol tautomers, shown in Figure . Enethiol is about 70 cm–1 more stable,[15] with the barrier height to interconversion slightly lower than in the malonaldehyde. This implies that the TS is similar in magnitude to the energy asymmetry of the wells and the states in different wells that lie below the barrier are expected to interact. Interestingly, it has been suggested[15] that the replacement of hydrogen, shared by the hydrogen bonds OH–S and SH–O, by deuterium reverses the stability order of tautomers due to the zero-point energy effect.

Figure 2

Tunneling tautomers of thiomalonaldehyde.

The asymmetry can also be caused by the environment. Molecules in rare gas matrices can have energy asymmetry between the wells comparable to their TS in isolation. Delocalization of the tunneling hydrogen was observed[16] in 9-hydroxyphenalone embedded in a neon matrix. Molecules in crystals in the vicinity of a suitable guest molecule can also have comparable energies of the splitting and energy asymmetry of the wells.[17]

Quantum tunneling has also been observed in macroscopic systems. Tunneling of Bose–Einstein condensates,[18] electron spin tunneling in the nanomagnetic molecules,[19] or tunneling of the magnetic flux in superconducting circuits based on Josephson junctions[20] are some recent examples. In a collective macroscopic variable, these processes can be described by a double well with externally controllable parameters that can induce asymmetry between the wells.

Calculation of the TS in moderately large molecules is prohibitively costly. Exact variational methods for determining the bound states of molecules scale exponentially with the basis set size, while large basis sets are often required.[21] Basis functions need to span over two or more wells sufficiently densely to obtain enough resolution to extract the splittings from the difference of the energies in their spectrum. The asymmetry of the wells also suggests that the symmetry cannot be used to reduce the size of the problem. Full-dimensional studies of malonaldehyde using multiconfigurational time-dependent Hartree[22,23] (MCTDH) and variational calculations on the HF dimer[21] or H2O dimer[24] represent the state-of-the-art calculations of the vibrational levels using formally exact methods.

A direct calculation of the TS in larger systems can be performed using a recently developed path integral molecular dynamics method[25] based on the potential sampling around the minimum action paths (MAPs) connecting different wells. The multiwell splitting patterns of the water trimer and hexamer[26] were obtained in this way using a matrix model of Hamiltonian in the basis of local vibrational states. The tunneling matrix (TM) elements are extracted from the zero-temperature limit of the partition function, which means that the method only works for the vibrational ground state in symmetric well systems.

Alternatively, the TM elements can be estimated using semiclassical methods. From that class, the instanton method, which comes in several forms,[27−30] has some particularly appealing features. It can be applied in Cartesian coordinates[29,31] to any molecule without modification. Numerically, it relies on the optimization of the MAP that connects the symmetry-related minima[32] and requires the potential and Hessians of the potential along the MAP to evaluate the splittings. It thus relies on a modest number of potential and gradient evaluations[33] in comparison with the exact quantum mechanical (QM) methods. This allows one to perform calculations in full dimensionality or in combination with on-the-fly evaluation of the electronic potential. Additionally, its accuracy is higher for large barriers and small splittings. Precisely in this regime, the exact variational methods become inefficient and resource-intensive.

The first derivation of the multidimensional instanton theory was accomplished by means of Jacobi fields integration (JFI).[30] The JFI method has been used to determine TSs for a range of symmetric double-well systems, such as malonaldehyde,[30,32,33] the vinyl radical,[34] and the formic acid dimer.[35] The instanton method was later rederived in the ring polymer form (RPI),[29] which could treat asymmetric potentials along MAPs and multiple wells. The RPI was used to calculate and interpret experimental ground-state splitting patterns of water clusters in terms of their rearrangement dynamics[8] for the dimer,[29,36,37] trimer,[29,38] hexamer,[39] and octamer.[40] We recently generalized the JFI method[31] to treat the multiwell systems and used it to explain the ground-state splitting pattern of 320 states in the water pentamer in terms of five dominant rearrangement pathways.[41] The extension of the method to low-lying vibrational states[42] is based on the work of Mil’nikov and Nakamura[43] and forms the groundwork of calculating the TM elements between local vibrational states of different wells in the present study below.

Weakly biased double-well systems have been considered in previous studies by several authors. Analytical results in one dimension have been obtained using the semiclassical Wentzel–Kramers–Brillouin (WKB) and instanton methods. Garg has demonstrated[44] that the instanton method and the WKB method with the Herring formula[45] give equivalent results for the TS in symmetric systems. Cesi et al.[46] considered a one-dimensional (1D) double-well with the shape asymmetry and no energy asymmetry using instantons and obtained an expression for the ground-state TS. An approximate solution for a 1D double well with a weak bias was also obtained by Mugnai and Ranfagni,[47] using instantons based on the MAP that does not fully connect the minima of the two wells. Leggett et al. obtained a solution[48] by adding a parabolic correction potential to remove the asymmetry between the wells, the contribution of which was then subsequently subtracted from the action integral. Dekker[49] derived the ground-state TS from the quantization condition via asymptotic matching of the semiclassical wavefunction in the barrier to the parabolic cylinder wavefunctions of harmonic oscillators in the two wells. Song[50,51] extended Dekker’s method[49] (as Halataie and Leggett[52] have performed independently) to obtain the TS in vibrationally excited states of asymmetric 1D potentials with an arbitrarily large shape and energy asymmetry. Song also showed[51] that the instanton wavefunctions with the Herring formula in a 2 × 2 matrix model give equivalent results to those obtained using Dekker’s method.[49]

In multidimensional systems, tunneling can be assisted or supressed by the excitation of transversal vibrational modes.[43,53] In the presence of asymmetry, the excited states of one well can be in a resonance with the states of another well with a different set of local quantum numbers, which results in a delocalization of the wavefunction across these wells.[51] Benderskii et al. devised a multidimensional perturbative instanton method[54] in which they treat the asymmetry of the potential in an analytic two-dimensional (2D) model as a correction of first order in ℏ, the same as that for energy. In this way, the MAP remains symmetric, and the asymmetry is moved to the transport equation along with energy. They also show that the equivalent expressions for the TS are obtained using the instanton quantization condition of Dekker[49] and using the instanton or WKB wavefunctions with the Herring formula[45] in one dimension. The method was applied to calculate TSs in excited vibrational states of malonaldehyde[55] with the asymmetric isotopic substitutions using a fit of model potential parameters to quantum-chemical data. The extensions of the RPI and JFI methods to the ground-states of the asymmetric systems with a weak bias have recently been derived and applied to PD malonaldehyde[9] and the water trimer,[14] respectively.

The object of this paper is to propose a method for calculating the vibrational tunneling spectrum of multiwell systems of mid-sized molecules that are outside the reach of the exact quantum methods. For this purpose, we extend the usual 2 × 2 matrix model to the ∑N × ∑N model, which represents the molecular Hamiltonian in the basis of all N local vibrational states of each well m. We rederive a generalized Herring formula[45,54] in order to calculate the off-diagonal TM elements that represent the interaction of local vibrational states of different wells. The semiclassical wavefunctions at the dividing plane, in the barrier that separates the wells, are obtained using the recently generalized JFI method.[31,43] The JFI wavefunctions are thus used to calculate the couplings between states that have different energies and normal-mode excitations for the first time. The diagonal energies of the local vibrational states can be calculated using any accurate quantum method with a basis set that spans only one well. Vibrational configuration interaction[56−58] (VCI) is used in this work. The effect of rotations is neglected.

The method, presented in Section , allows one to study the vibrational structure in asymmetric systems with multiple wells, separated by large potential barriers, in an approximate manner. The accuracy of the method is tested on a 2D double-well model in Section . In Section , it is applied to the (symmetric) malonaldehyde molecule, which tests the accuracy of the matrix model using a combination of VCI and JFI matrix elements on a realistic PES, in vibrationally excited states, against the exact MCTDH calculations. The vibrational tunneling spectrum of PD malonaldehyde is calculated in Section , which features the mixing of inequivalent well states due to tunneling. The paper concludes in Section .

Tunneling Matrix

Without the loss of generality, we start by considering a system with two minima separated by a large potential barrier. The minima, denoted as “left” (L) and “right” (R), are not necessarily symmetric either in shape or energy. For low-energy spectra, the vibrational Hamiltonian can be represented in the basis of states that are localized in the wells, {ϕ(L),ϕ(R)}, asSquare blocks H(L/R) are formed using basis functions of the same minimum and are not necessarily of equal size. Their off-diagonal elements describe the interaction between different basis functions localized in the same minimum and can be made small by a suitable choice of the basis. In the instanton theory of TSs, the usual presumption is that the local vibrational wavefunctions are harmonic oscillator states. In that case, the off-diagonal terms describe anharmonic contributions that originate from the difference between the actual potential and the harmonic potential.

In our approach here, we replace the harmonic surface of each well by an n-mode representation[59,60] of the well potential and calculate local eigenfunctions and eigenvalues using the vibrational self-consistent field (VSCF) and VCI methods.[56−58] The technical details of the calculations are described in Section 2 of the Supporting Information. Using the more accurate local wavefunctions as a basis reduces the magnitude of the off-diagonal matrix elements in H(L/R), which we then neglect. The matrices H(L/R) become diagonal, and the diagonal matrix elements are referred to as the local vibrational energies of the left/right (L/R) well. For symmetric wells, the local energies are doubly degenerate.

The block h in matrix (1) contains the TM elements that describe the interaction of local wavefunctions of the left and right minima. The exact quantal calculation of these elements requires a large basis set that can accurately represent the form of the wavefunction inside the barrier. Instead, we obtain them by means of the Herring formula[45] in combination with the semiclassical wavefunctions from the instanton theory.[42,43] Since the only effect of matrix (1) is to mix local wavefunctions of different minima via tunneling, we refer to it as the TM.[29]

We now derive the Herring formula without the usual assumptions of the two-state model and the L/R symmetry. Rather, we consider the Schrödinger equation with Hamiltonian (or tunneling) matrix (1), from which it follows thatwhere E(L/R) are the local vibrational energies and the summation over repeated index k is assumed. We proceed in mass-scaled Cartesian coordinates and define a dividing plane inside the barrier via the implicit equation fD(x) = 0, which separates the left minimum from the right minimum. Equation are multiplied by ϕ(R) and ϕ(L), respectively, subtracted, and integrated over the left part of the domain (i.e., over the space on the “left” side of the dividing plane). The local wavefunctions ϕ(L/R), either harmonic or VCI, have been obtained as eigenfunctions of a Hermitian matrix and are therefore taken to be orthonormal. For a sufficiently high barrier, the wavefunctions ϕ(L/R) can be considered small in the R/L domain, respectively. We thus neglect the integrals involving the like products ϕ(R)ϕ(R) in the L volume and extend the integrals involving ϕ(L)ϕ(L) over the entire domain to produce δ. The integrals involving the mixed products ϕ(L)ϕ(R) have also been neglected. The error introduced by the neglect of these terms outside the resonance, that is, for E(L) ≠ E(R), is analyzed in the Appendix. The TM element is then expressed aswhere we exploited the use of mass-scaled Cartesian coordinates and, in the last step, used the divergence theorem to turn the spatial integration into the integral over the dividing plane. S in eq denotes the coordinate that describes an orthogonal shift from the dividing plane.

Local wavefunctions, which we designed to calculate the local vibrational energies on the diagonal of matrix (1), are constructed using the VSCF/VCI on an approximate PES (see Section 2 in the Supporting Information), and their accuracy drops inside the barrier that separates the wells. In order to evaluate the surface integral in the Herring formula, eq , inside the barrier, we employ the JFI wavefunctions instead, which we recently derived in ref (42). These are based on the WKB method in which the energy is treated as a term of order ℏ1 and is moved to the transport equation, leaving the Hamilton–Jacobi equation energy independent. It was shown that this approach gives equivalent results to the standard WKB method in one dimension.[44] Moreover, the ground-state TS obtained from the Herring formula using the ground-state JFI wavefunctions[42] is identical to the standard instanton result derived from the steepest descent approximation of the partition function in the path integral formulation.[31]

The characteristic of the Hamilton–Jacobi equation that connects the minimum of a well to a point in the configuration space obeys the equation[42]and represents a classical trajectory x(τ) on the inverted PES, parameterized by the “imaginary” time τ. In order to represent the quantities in the neighborhood of the characteristic, N local coordinates (S, Δx) are defined,[42,43] where S is the mass-scaled arc length distance from the minimum along the characteristic and Δx is the orthogonal shift from the nearest point on the characteristic. The classical momentum is defined asand S can be used, instead of τ, to reparameterize the characteristic. Local wavefunctions in the harmonic vicinity of the characteristic are obtained by integrating the Hamilton–Jacobi and transport equations on the characteristic[42] asFor vibrationally excited states, label ν in eq represents the number of quanta in the excited vibrational mode of frequency ωe. Matrices A(L/R) are Gaussian widths of the wavefunction in the directions orthogonal to the characteristic and are obtained fromH(S) in eq is the Hessian of the potential at S, which is used to approximate the potential up to quadratic terms in the neighborhood of the characteristic. The initial condition for eq at the minimum is A0(L/R) = (H0(L/R))1/2, where H0(L/R) is the Hessian at the L/R minimum. For vibrationally excited states, the prefactor in the parentheses in eq contains terms F(S) and U(S), which are defined via equationsF(L/R) terms account for the change in the amplitude of the excited-state wavefunction along the characteristic, while the U(L/R) term describes the nodal plane. The initial condition for F(L/R) is found by matching the instanton wavefunction to that of the harmonic oscillator at a small distance S = ε from the minimum as F(L/R) (ε) = U0(L/R),T(x(ε) – xmin(L/R)). The U0(L/R) is the excited-state normal mode, that is, the eigenvector of H0 having frequency ωe, and serves as the initial condition for U in eq .

The local instanton wavefunctions, eq , for the left and right minima are next inserted into the Herring formula given in eq without the previous assumption[42] that ϕ(L) and ϕ(R), each, refer to the excitation of the same normal mode and the same number of quanta ν. For this purpose, a connection point x(Scp) is chosen on the dividing surface fD(x) inside the barrier, and characteristics are determined, which connect it to the minima on both sides of the dividing surface. The shape of the characteristic between two points in the configuration space is determined by minimizing the Jacobi action.[32] The surface integral in eq can then be computed analytically.[42]

This approach yields the best results if the connection point is chosen so that both wavefunctions are near their maxima in the dividing plane at the connection point. This can be obtained by the minimization of the sum of action integrals . For minima of the same energy, this procedure yields the MAP that connects the minima and any point on that path is a suitable candidate for the connection point. The dividing surface can then be chosen as the plane orthogonal to the MAP at the connection point. If the minima do not have the same energies but differ by the amount d, this procedure is equivalent to determining the MAP on the modified PES Ṽ(x) = V(x) – Θ(S – Scp)d, where Scp is the position of the connection point on the characteristic and Θ is the Heaviside step function. In this case, the position of the connection point has to be given a priori, and the resulting path will depend on its position. The safest choice is to pick the connection point in the middle of the MAP, which is expected to be near the maximum of the potential energy barrier. For minima at different energies, the resulting MAP is going to have a tangent dicontinuity at the connection point as p0(L) ≠ p0(R) at Scp. The tangent direction at the connection point is then defined as the average tangent of its L and R limits at Scp. Again, the dividing plane is taken to be orthogonal to the MAP, and the surface integral in eq is solved analytically. The TM element then becomeswith all quantities in the brackets evaluated at S = Scp. In eq , Stot is the total length of the MAP, , and the symbol ⊥ means that the tangent direction to the MAP was explicitly projected out. The det′ in eq denotes the product of all non-zero eigenvalues. Matrices A0 have zero eigenvalues associated with the overal translations and rotations, while A̅ has an additional zero eigenvalue associated with the tangent to the MAP. For energy-equivalent minima, the tangent vector is an eigenvector of A̅ with zero eigenvalue, and the explicit projection to the orthogonal space is not needed. The TM element in eq is valid for ν, ν′ = 0–1. For ν > 1 and multiple excitations in different modes, the TM element can still be evaluated using the Herring formula and wavefunctions of form (6) using analytical integrals, but we have only implemented it numerically, without trying to write down the explicit form. We also remark here that the wavefunction in eq for the multiply excited normal modes, ν > 1, does not correspond to the harmonic oscillator wavefunction near the minimum as the prefactor in eq is not a Hermite polynomial. We further note that the TM element, (9), is not invariant with respect to the position of the connection point Scp unless the two local states are in resonance, as shown in the Appendix.

Numerical Tests

Numerical tests were carried out on a model 2D PES and on the malonaldehyde molecule with some atoms substituted with heavier isotopes. MAPs that connect the minima were determined using the string method.[32,33] The criterion for convergence was chosen to be the largest component of the gradient of Jacobi action perpendicular to the path and was set to be 10–6 au. The number of beads used to discretize the string was 301 for model potential, which is much larger than necessary for convergence but was used to ensure that all results obtained using different parameters of the potential are sufficiently converged. For the potential with minima at different energies, the dividing plane was set to pass through the central bead and perpendicular to the MAP. In the case of malonaldehyde, the number of beads was 201, and the minima were oriented toward the first neighboring bead in each step of the optimization to minimize the root-mean-square distance between their geometries.[32] After optimization, Hessians of the potential were determined at each bead on the MAP. Translations and rotations were explicitly projected out from Hessians.[31] Geometries along the path in mass-scaled Cartesian coordinates, potential, and Hessian matrix elements were parameterized by the arc length S along the MAP and interpolated using natural cubic splines. Matrices A(L/R) in eq were propagated using the previously described approach,[42] with the initial “jump” at ε = 0.1 au for model potential and ε = 0.25 au for malonaldehyde. The fourth order Runge–Kutta method was used for integration of eq . Matrices A(L/R) (S) were saved at each bead, and their matrix elements were interpolated using natural cubic splines, as implemented for the Hessians above. The interpolant was then used to propagate F(L/R) and U(L/R) in eq from minima up to the dividing plane.

The particular implementation of the VSCF/VCI method that is employed in our calculation here is described in Section 2 of the Supporting Information. We determined the 1-mode and 2-mode terms of the PES and neglected the terms beyond. In each normal mode, the potential was evaluated at Gauss–Hermite discrete variable representation (DVR) points, which correspond to the zeroes of Hermite polynomials. We used eight DVR points for the 2D model potential and 11 DVR points for malonaldehyde. This approach utilizes the natural lengthscales of the harmonic oscillators in each normal mode, which gives a balanced description of the potential at different minima. The 1-mode terms were then fitted to the eighth-order polynomials using linear regression. For 2-mode terms, the potential was computed on a rectangular grid of DVR points determined above, and a fit was performed analogously. For each 1-mode potential, a quick QM calculation was performed using sine DVR basis with 100 basis functions. The difference in the lowest two energies from that calculation was used as a frequency for the harmonic oscillator basis set, which was used to solve the VSCF equations. This approach provides a better basis for determining the 1-mode potentials that quickly deviate from the harmonic curve, reducing the number of basis functions needed to describe the 1-mode functions in the VSCF. The basis sets of 7 and 16 harmonic oscillator functions in each normal mode in the 2D model potential and malonaldehyde, respectively, were needed to converge the energies (to less than 1 cm−1 for the VSCF ground state). The existence of a plateau over which the energies of interest do not change appreciably with the basis set size, in comparison with the size of the TM elements, indicates that the approximate separability is possible, which is a necessary condition for the applicability of the proposed combined approach. A larger basis should not be used as functions corresponding to larger energies start to penetrate into unphysical parts of the fitted potential, which can cause the appearance of intruder states and worse energies. After VSCF calculation, the computed 1-mode functions were used for VCISD calculation, where the highest excitation in each mode was limited to six in both the 2D model system and malonaldehyde.

2D Model Potential

The 2D model potential with two minima, which we use in our test calculations below, is defined by the following equationswhere x are not mass-scaled. Minima are located at x(L/R). Coefficients γ1 and γ2 are chosen so that in the vicinity of the left minimum, the potential is approximately harmonic and equals V ≈ V(L), while in the vicinity of the right minimum, the potential is approximately harmonic and shifted in energy by d, that is, V ≈ V(R) + d. α1,L/R and α2,L/R are eigenvalues of the Hessian, while U0(L/R) are normal modes. Parameter θ denotes the angle of inclination of the normal mode to the x axis. The mass of the system was taken to be m = 3.5 in both dimensions so that the harmonic frequencies are given by .

The above form of the potential can be used to independently vary harmonic frequencies ω1/2(R), by changing parameters α1/2,R, or the shift d without affecting the other parameters of either the left minimum or the right minimum. In this paper, the parameters of the left minimum were α1,L = 1.6 and α2,L = 4.0. The parameters of the right minimum were the same as the parameters of the left, for the symmetric case with d = 0. To obtain the asymmetric potentials below, one of the three parameters was varied, with α1,R going from 1.6 to 36, parameter α2,R going from 4 to 49 and d going from 0 to 1.1. Positions of the minima were set with β = 2.0 and angle θ = π/12. This angle corresponds to approximately equal contributions of F(L/R) and U(L/R) in the TM elements.[42]Figure shows the model potential for a selection of parameters α1,R, α2,R, and d.

Figure 3

PES of the 2D model given by eq . Top left panel corresponds to the symmetric potential, top right to ω1(R) > ω1(L), bottom left to ω2(R) > ω2(L), and bottom right to d > 0, with other parameters set equal to the symmetric case.

Frequency ω1 is the lower frequency, and the MAP enters the minima along the corresponding normal mode. Consequently, ω1 does not contribute toward the zero-point energy in the plane orthogonal to the MAP. The effective barrier for the tunneling motion from the ground state in the left minimum, corrected by the zero-point motion contribution, can be defined asVmax in eq is the maximum of the potential V(Smax) along the MAP. λ(L) is the non-zero eigenvalue of the A⊥ = PAP matrix, where P projects out the tangent direction to the MAP at S = Smax. The effective barrier can be defined for other states similarly. Figure (in the second column panels) shows that for the symmetric case, ω1(R) = ω1(L), the JFI theory provides accurate TSs in the ground state, and in the second excited state, which corresponds to the excitation of the transversal frequency ω2. In the first excited state, the JFI theory slightly overestimates the TS. In that state, the effective barrier is much smaller and equals Veff = 0.545, in contrast with the barriers of 1.221 and 0.976 for the ground and second excited states, respectively. This overestimation is a known property of the instanton method.[42]In the symmetric case, the only contribution to the splitting comes from the off-diagonal matrix elements so that the harmonic and VCI energies yield the same results. However, it can be observed (from the first column panels in Figure ) that the harmonic vibrational energies overestimate the exact QM energies by 3–5%.

Figure 5

Dependence of vibrational energies of the lowest six states in the double-well potential given by eq on ω1(R). Circles represent QM values, blue lines are obtained using the instanton method with harmonic energies, and red lines are obtained using a combined VCI/instanton approach. Frames I–III in the top panel are shown magnified in the left column panels below, and the dependence of the associated TSs on ω1(R) is shown in the right column panels below.

As the frequency ω1(R) is increased and the difference in the local L/R energies begins to contribute to the overall splitting, the TSs computed using harmonic energies quickly begin to deviate from the QM values due to the neglect of anharmonicities, which no longer cancel out. When the difference in the lower frequency, ω1(R) – ω1(L), is only 0.02, which corresponds to the asymmetry (Δω1/ω1(L)) of 3%, the error in the TS of the ground state is 24%, whereas it is 90% for the transversal mode (ω2) excitation. A larger error in the excited state reflects the fact that the local excited-state wavefunction penetrates deeper into the barrier, where anharmonicity is larger. However, the VCI energies correctly account for the anharmonicity and provide an excellent agreement, as can be observed in Figure , both in the absolute energies (first column panels in Figure ) and in the TSs (second column panels).

With a further increase in the frequency ω1(R), different local vibrational states of the left and right minima enter into resonance, and vibrational energies exhibit avoided crossings, as shown in frames IV–VI of the top panel in Figure . Harmonic energies do not provide accurate positions of these avoided crossings, as seen in Figure , due to errors in the local energies. In the case of the avoided crossing between the higher-frequency ω2(L)-excited state of the left minimum and the ω1(R)-excited state of the right minimum, as shown in frame IV of the top panel in Figure , the error in the position of the avoided crossing (in ω1(R) – ω1(L)) using harmonic energies is 16% (and falls outside frame IV in Figure ). VCI energies, as shown magnified in Figure together with the exact QM energies, provide a significantly more accurate position with the error of only 0.4%. The small discrepancy can be attributed to the fact that as frequency ω1(R) is increased, the local wavefunction in the right minimum penetrates deeper into the barrier. In this region, the approximate n-mode representation of the potential used in the VCI calculations begins to deviate from the actual potential, which introduces an error in the local energies.

Figure 6

Dependence of vibrational energies and TSs in the 2D model potential given by eq on the frequency ω1(R) in the region of the avoided crossing between the first (ω1(R)-) excited state in the right minimum and the second (ω2(L)-) excited state in the left minimum, shown in frame IV in the top panel of Figure . Circles represent QM values, while red lines represent values obtained using a combined VCI/instanton approach.

The TS in the avoided crossing is reproduced with great accuracy, shown as the minima in the lower panel in Figure , with the error of 5%. Errors in the positions of other avoided crossings (IV–VI in Figure ), namely, between the ground state of the right minimum and the ω1(L)- and ω2(L)-excited states in the left minimum (frames V and VI in Figure , respectively) become larger even using VCI energies. The local wavefunction of the right minimum has a larger energy and penetrates deeper into the region where the n-mode representation of the potential becomes unreliable. Nevertheless, the TSs in the avoided crossings are again reproduced accurately, which indicates that the JFI method can indeed give reliable TM elements between different vibrational states of L/R minima and, in combination with the VCI energies, is a useful tool for computing vibrational tunneling spectra. Similar results were observed with the frequency ω2(R) varied (shown in Figures S3 and S4 in the Supporting Information).

Figure shows the dependence of energy levels on the variation in the depth d of the right minimum. Overall, the introduction of the energy asymmetry between the wells results in a similar energy level pattern to that observed above. A notable difference is that in this case, the TSs obtained using harmonic energies are much closer to the exact QM values. This is an artefact of the construction of the PES, in which the frequencies in the left and right minima are the same. As a result, the shapes of the local potentials in both minima are similar, and a large part of the error introduced by the anharmonic terms cancels out. However, in realistic applications, it is unlikely that systems with minima of different energies have the same L/R frequencies. The error in the position of the avoided crossing IV is also much smaller for the harmonic energies (≈2%), while it is further reduced using VCI energies (0.7%), as shown in Figure . The error in the TS in the avoided crossing is 12%, which is comparable to the error in the case of the frequency variation, as shown in Figure .

Figure 7

Dependence of vibrational energies of the lowest six states in the double-well potential given by eq on the energy shift d of the right well. Circles represent QM values, blue lines are obtained using the instanton method with harmonic energies, and red lines are obtained using a combined VCI/instanton approach. Frames I–III in the top panel are shown magnified in the left column panels below, and the dependence of the associated TSs on d is shown in the right column panels below.

Figure 8

Dependence of vibrational energies and TSs in the 2D model potential given by eq on the energy shift d of the right well in the region of the avoided crossing between the first (ω1(R)-) excited state in the right minimum and the second (ω2(L)-) excited state in the left minimum, shown in frame IV in the top panel of Figure . Circles represent QM values, blue lines are obtained using the instanton method with harmonic energies, and red lines represent values obtained using a combined VCI/instanton approach.

Malonaldehyde

We next employ our combined approach to study the symmetric, homoisotopic malonaldehyde on the PES developed by Wang et al.[61] The molecule is shown labeled in the top panel in Figure . It has two equivalent wells with hydrogen 6 attached to either oxygen 1 or 5. We study below the effect of adding additional states in the TM. For this purpose, vibrational energies are computed either from a 2 × 2 matrix involving corresponding states in the two wells, an 8 × 8 matrix involving four local states at both sides of the barrier, or a 16 × 16 matrix model. Thereby, we again calculate the local single-well states using the VSCF/VCI, while the TM matrix elements are computed using the recently developed JFI method.[42]

Malonaldehyde has been extensively studied in the past[62] and presents a benchmark system for the development of quantum dynamical methods. Most recent calculations on the same PES using exact quantum methods were obtained using MCTDH by Hammer and Manthe[23] and Schröder and Meyer[22] and show a good level of agreement with experimental results.[63] We use the results of ref (23) for comparison as they report TSs for a number of vibrationally excited states having a large transition dipole moment and are believed to be more accurate.[22]

Local harmonic and VSCF/VCI energies, calculated in a 2-mode representation of the single-well potential, as described in Section 2 of the Supporting Information, for the lowest 8 vibrational states, that we consider below, are shown in Table . The ground state is labeled GS, while the excited states are labeled by the frequencies ν of the excited normal modes, numbered in the order of increasing frequency in the subscript and separated by a “+” sign for multiple excitations. A noticeable shift can be observed between all harmonic and VCI energies in Table due to anharmonicity, but the order in energies remains unchanged. The lowest four normal modes that can get excited in the lowest eight local vibrational states and that play a role in our calculations below are depicted in Figure . Higher vibrational states become more densely spaced in energy and start to mix vibrational modes at minima. Our approach relies on being able to uniquely define the excited normal modes at minima for each local vibrational state considered because the instanton wavefunctions that are used to calculate the TM elements that connect these states tend to harmonic oscillator eigenstates at minima. Moreover, a higher density of states at higher energies would require the inclusion of many additional states in the TM, which are not known as precisely as for the low-lying states and would thus degrade the accuracy. We limit ourselves, therefore, to the lowest eight local states in the studies of tunneling spectra of malonaldehyde below.

Table 1

Harmonic and VCI Energies in Inverse Centimeters (cm–1) of the First Eight Local Vibrational States of Malonaldehyde Labeled by the Excited Normal Mode Frequenciesa

state	harmonic	VCI
GS	14,950.11 (0.00)	14,682.46 (0.00)
ν1	15,218.68 (268.57)	15,012.65 (330.19)
ν2	15,245.53 (295.42)	15,042.51 (360.05)
ν3	15,333.29 (383.17)	15,133.95 (451.49)
ν1 + ν1	15,487.25 (537.14)	15,262.05 (579.59)
ν4	15,472.20 (522.08)	15,281.89 (599.43)
ν2 + ν2	15,540.95 (590.83)	15,318.85 (636.40)
ν1 + ν2	15,514.10 (563.99)	15,336.16 (653.70)

Energies relative to the local ground state are given in parentheses.

The TM elements in the h matrix that connect the two sets of local states in the L and R wells are calculated using the JFI method and listed in Table . Both minima of malonaldehyde belong to the C symmetry group, and its local vibrational states can be classified according to the irreducible representation of the excited normal mode ν at the minimum. The C symmetry is preserved along the MAP so that the TM elements that connect normal modes of different symmetries vanish exactly, as seen in Table .

Table 2

TM Elements Connecting the First Eight Local Vibrational States of Different Minima in Malonaldehyde in Inverse Centimeters (cm–1)

	GS(R)	ν1(R)	ν2(R)	ν3(R)	(ν1 + ν1) (R)	ν4(R)	(ν2 + ν2) (R)	(ν1 + ν2) (R)
GS(L)	–12.30	0.00	–21.94	0.00	–4.98	–4.62	–25.53	0.00
ν1(L)	0.00	–6.70	0.00	6.85	0.00	0.00	0.00	–11.95
ν2(L)	–21.94	0.00	–44.20	0.00	–8.87	–7.54	–55.97	0.00
ν3(L)	0.00	6.86	0.00	8.53	0.00	0.00	0.00	12.22
(ν1 + ν1) (L)	–4.98	0.00	–8.88	0.00	–4.84	–1.87	–9.14	0.00
ν4(L)	–4.61	0.00	–7.53	0.00	–1.87	7.82	–8.14	0.00
(ν2 + ν2) (L)	–25.53	0.00	–55.97	0.00	–9.14	–8.15	–75.86	0.00
(ν1 + ν2) (L)	0.00	–11.95	0.00	12.22	0.00	0.00	0.00	–21.31

In a 2 × 2 matrix model, only the diagonal elements of the h matrix are used, and the degenerate vibrational states of L/R wells are split into doublets. Equivalent results are obtained using the first-order perturbation theory for degenerate states, yielding the TS of Δ = 2h. Energies of the GS and the first three excited states obtained in this manner already show a good agreement with the MCTDH results of ref (23), as can be seen in Table (from the second column and the last column). The vibrational states are numbered in the order of increasing energy in Table . The wavefunction content, obtained from the eigenvectors of the TM, is listed in Table and can be used to identify states in Table in terms of the excited normal modes.

Table 3

Vibrational Energy Levels of Malonaldehyde in Inverse Centimeters (cm–1) Obtained Using a Combined VCI/Instanton Approacha

no.	E(pairs)	E(4)	E(8)	E(MCTDH)
1	14,670.15	14,668.69	14,667.08	14,671.3
2	14,694.76	14,693.54	14,692.76	14,694.8
3	14,998.31	14,999.77	14,987.74	14,941.5
4	15,005.95	15,005.60	15,005.09	15,008.2
5	15,019.35	15,018.91	15,018.54	15,014.9
6	15,086.70	15,087.92	15,077.14	15,005.4
7	15,125.42	15,125.86	15,125.14	15,108.3
8	15,142.47	15,142.82	15,142.04	15,124.6
9	15,243.00		15,249.04
10	15,257.21		15,263.41
11	15,266.89		15,266.12
12	15,274.07		15,273.84	15,249.6
13	15,289.71		15,291.11	15,268.4
14	15,314.85		15,316.14
15	15,357.47		15,358.55
16	15,394.71		15,407.27

E(pairs), E(4), and E(8) are energies obtained from the 2 × 2, 8 × 8, and 16 × 16 matrix models, respectively, as explained in the text. E(MCTDH) are MCTDH energies from ref (23).

Table 4

Dominant Configurations of Vibrational States of Malonaldehyde, Obtained as the Eigenvectors of the TM in the 2 × 2 (Pairs) and 16 × 16 (Eight-State) Models, as Described in the Text

The TSs for the GS and the singly excited modes ν1–4 in the 2 × 2 TM model are obtained as Δ(GS) = 24.60 cm–1, Δ(ν1) = 13.40 cm–1, Δ(ν2) = 88.40 cm–1, Δ(ν3) = 17.06 cm–1, and Δ(ν4) = 15.64 cm–1. The MCTDH results[23] for the TSs in the same states are Δ(GS) = 23.5 cm–1, Δ(ν1) = 6.7 cm–1, Δ(ν2) = 69.9 cm–1, Δ(ν3) = 16.3 cm–1, and Δ(ν4) = 18.8 cm–1. Differences in TSs, apart from the ν1- and ν2-excited modes, are well within the estimated error of the MCTDH calculations, which validates the accuracy of our approach. The ν2 mode corresponds to the longitudinal mode as it lies parallel to the MAP at minima. The excitation of this mode effectively lowers the barrier of the tunneling motion, and the instanton theory is known to overestimate TSs in the shallow tunneling regime.[29,42] The wavefunction also penetrates deeper into the barrier where the anharmonic effects are larger, and the VCI energies degrade as a result. Thus, the accuracy in absolute energies in Table is also expected to be affected for these states. The large increase in the TS for the excitation of the longitudinal mode is, however, expected,[42] as confirmed by our results. The TS for the ν1 mode is overestimated by a factor of two. This is most likely due to the anharmonicity along this normal mode, indicated by the large difference between the harmonic (268.57 cm–1) and VCI (330.19 cm–1) energies. Since the TS for the pair of states is significantly suppressed compared to the GS, the frequency and energy in its direction change substantially along the MAP. Therefore, if the anharmonicity also changes significantly, it could cause the observed discrepancy. As an aside, we also note here that the other TSs computed using MCTDH in ref (23), which do not result in the mixture of normal modes at minima, are Δ(ν5) = 21.1 cm–1, Δ(ν7) = 33.3 cm–1, Δ(ν8) = 14.6 cm–1, and Δ(ν11) = 19.5 cm–1 and are in good agreement with the values we obtained using the JFI theory as Δ(ν5) = 24.4 cm–1, Δ(ν7) = 39.5 cm–1, Δ(ν8) = 15.6 cm–1, and Δ(ν11) = 22.1 cm–1.

We next consider the TM using four local states in each well. This takes into account interactions between the doublets considered above, whereby only the states of the same symmetry interact. If the states of the same symmetry are well-separated with respect to the size of their TM element, the shift in energy can also be computed using the second-order perturbation theory. When four local states are taken into account in the 8 × 8 TM model, slight shifts are observed in the GS and ν2-doublets in Figure (left-side spectrum). The absolute energies change by 1.22–1.46 cm–1, while the perturbation theory gives the shift of 1.34 cm–1. However, the change in the TS is negligible.

Figure 9

Vibrational tunneling spectrum of the lowest 8 (left panel) and 16 (right panel) states of malonaldehyde. Green and blue lines represent VCI energies of local wavefunctions in the D7 and D9 minima, respectively. Dashed red lines are obtained using a 2 × 2 TM model. Black lines in the left panel represent energies from an 8 × 8 TM model, and in the right panel, they represent energies from a 16 × 16 model. See the text for details.

In the 16 × 16 TM model, consisting of eight local states in each well, a strong interaction with the doubly excited (ν2 + ν2) mode causes a significant shift in the energies of the GS and the ν2-excited doublets as well as their splittings. The TSs change from 24.85 to 25.68 cm–1 and from 88.15 to 89.4 cm–1, which can clearly be observed in Figure (right-side spectrum). A particularly strong mixing also occurs between the doubly excited (ν2 + ν2) mode and the doubly excited (ν1 + ν1) mode, for which the lower levels in the doublets are very close in energy (14.21 cm–1), and they interact strongly (h = 9.14 cm–1 in Table ). The mixing results in visible changes in the dominant coefficients of TM eigenvectors in Table and leads to observable energy shifts. Furthermore, the singly excited ν4 mode interacts and mixes with the doubly excited (ν1 + ν1) mode, which results in the change of its TS from 15.64 to 17.27 cm–1, which is in closer agreement with the MCTDH value of 18.8 cm–1. Finally, we remark that the TS of the doubly excited (ν1 + ν2) state amounts to 42.62 cm–1, which is in good agreement with 49.5 cm–1 obtained by Schröder and Meyer.[22]

The above results clearly show that the interactions of different vibrational states can have a non-negligible effect, both on the absolute values of the vibrational energies and on the values of the TSs. This effect is especially pronounced if two or more states of the same symmetry are close in energy and if the TM elements that connect them are large. This scenario is expected to play a significant role in the higher vibrationally excited states, where the density of states becomes larger and the interactions increase due to the presence of multiple excitations.

PD Malonaldehyde

In the previous subsection, we have learned what accuracy one might expect in the calculation of the tunneling spectra of malonaldehyde through comparison with the exact QM results. We now consider the PD malonaldehyde, where hydrogen in position 7/9 is substituted by deuterium (see Figure ) and the system is no longer symmetric. Since deuterium is not placed in equivalent positions in the two minima, their local vibrational frequencies and energies are no longer equal, even though the PES remains unchanged. The particular choice of deuteration was chosen for our study because the mixing angle in its GS was determined experimentally by Baughcum et al.[10] and the TS by Jahr et al.[9] using the RPI method. Furthermore, the size of the relative energy shifts between the left and right minima is comparable to the size of the TM elements, which makes the system interesting in that both the VCI energies and the instanton TM elements are expected to make a significant contribution to the TSs in this system.

In the PD malonaldehyde, the isotopic substitution causes a significant lowering of the zero-point energy, given in Table , from 14,682.45 to 13,978.19 cm–1 for the D7 minimum and to 14,013.04 cm–1 for the D9 minimum. Additionally, the excitation energies for the first seven excited states decrease as well, by up to 40 cm–1. As a result, the vibrational states are more closely spaced, see Figure , and larger interstate L/R mixings are expected.

Table 5

Harmonic and VCI Energies in Inverse Centimeters (cm–1) of the First Eight Local Vibrational States of PD Malonaldehyde Labeled by the Excited Normal Mode Frequenciesa

	harmonic		VCI
state	D7	D9	D7	D9
GS	14,228.18 (0.00)	14,253.67 (25.49)	13,978.19 (0.00)	14,013.04 (34.85)
ν1	14,492.55 (264.37)	14,492.10 (263.92)	14,298.70 (320.51)	14,311.75 (333.56)
ν2	14,522.65 (294.47)	14,547.57 (319.39)	14,327.02 (348.83)	14,361.38 (383.19)
ν3	14,568.85 (340.67)	14,626.56 (398.38)	14,384.49 (406.30)	14,444.96 (466.76)
ν1 + ν1	14,756.92 (528.74)	14,546.67 (502.34)	14,950.11 (568.48)	14,543.15 (564.96)
ν4	14,744.31 (516.13)	14,769.16 (540.98)	14,562.35 (584.16)	14,595.52 (617.32)
ν2 + ν2	14,817.13 (588.95)	14,841.47 (613.29)	14,601.27 (623.08)	14,637.01 (658.82)
ν1 + ν2	14,787.02 (558.84)	14,786.00 (557.82)	14,614.31 (636.12)	14,624.21 (646.02)

Energies relative to the local ground state of the D7 minimum are given in parentheses.

Figure 10

Vibrational tunneling spectrum of the lowest 8 (left panel) and 16 (right panel) states of PD malonaldehyde. Green and blue lines represent VCI energies of local wavefunctions in the D7 and D9 minima, respectively. Dashed red lines are obtained using a 2 × 2 TM model. Black lines in the left panel represent energies from an 8 × 8 TM model, and in the right panel, they represent energies from a 16 × 16 model. See the text for details.

The normal modes in PD malonaldehyde are qualitatively similar to the homoisotopic malonaldehyde, depicted in Figure . The ordering of local single-well states, labeled by the excited normal mode at the minimum, is also preserved upon deuteration, with the exception of the |(ν1 + ν2)(D9)⟩ and |(ν2 + ν2)(D9)|⟩ states, which exchange order. The TM elements, shown in Table , are remarkably similar to the homoisotopic malonaldehyde, which indicates that the wavefunctions in the barrier region are not significantly affected by the asymmetry. The error estimates due to the variation of the position of the dividing plane are shown in parentheses in Table and are discussed in more detail in the Appendix.

Table 6

TM Elements Connecting the First Eight Local Vibrational States of Different Minima in PD Malonaldehyde in Inverse Centimeters (cm–1)a

	GS(D9)	ν1(D9)	ν2(D9)	ν3(D9)	(ν1 + ν1)(D9)	ν4(D9)	(ν2 + ν2)(D9)	(ν1 + ν2)(D9)
GS(D7)	–12.32 (0.005)	0.00 (0.00)	–21.95 (0.133)	0.00 (0.00)	–5.04 (0.042)	–4.63 (0.034)	–25.54 (0.343)	0.00 (0.00)
ν1(D7)	0.00 (0.00)	–6.86 (0.001)	0.00 (0.00)	6.52 (0.014)	0.00 (0.00)	0.00 (0.00)	0.00 (0.00)	–12.22 (0.077)
ν2(D7)	–21.89 (0.109)	0.00 (0.00)	–44.12 (0.031)	0.00 (0.00)	–8.95 (0.037)	–7.52 (0.029)	–55.89 (0.452)	0.00 (0.00)
ν3(D7)	0.00 (0.00)	7.37 (0.007)	0.00 (0.00)	8.64 (0.007)	0.00 (0.00)	0.00 (0.00)	0.00 (0.00)	13.14 (0.054)
(ν1 + ν1)(D7)	–5.16 (0.042)	0.00 (0.00)	–9.19 (0.034)	0.00 (0.00)	–4.58 (0.000)	–1.94 (0.001)	–9.45 (0.024)	0.00 (0.00)
ν4(D7)	–4.44 (0.029)	0.00 (0.00)	–7.25 (0.021)	0.00 (0.00)	–1.82 (0.001)	7.97 (0.004)	–7.84 (0.010)	0.00 (0.00)
(ν2 + ν2)(D7)	–25.42 (0.304)	0.00 (0.00)	–55.79 (0.349)	0.00 (0.00)	–9.18 (0.015)	–8.10 (0.001)	–75.67 (0.090)	0.00 (0.00)
(ν1 + ν2)(D7)	0.00 (0.00)	–12.19 (0.072)	0.00 (0.133)	11.59 (0.037)	0.00 (0.00)	0.00 (0.00)	0.00 (0.00)	–21.72 (0.006)

Values in parentheses refer to the estimated error introduced by the neglect of the overlap between L/R local states, as explained in the Appendix.

We again consider pairwise interactions of the corresponding states in a 2 × 2 TM model. This is possible since the normal modes at both minima can approximately be mapped to one another using a symmetry operation. The pairs of states are no longer degenerate in this case, and the first-order perturbation theory cannot be used to estimate the TSs. Instead, the TS, obtained from the eigenvalues of the TM, is seen to be equal to the local energy difference corrected by the second-order perturbative termswhere, in the last line of eq , we assumed that the TM element |h| ≪ |E(D7) – E(D9)|. This assumption is certainly violated if there are other local states that are energetically close and coupled by the TM elements that are comparable in size.

The TM element for the GS is determined to be 12.32 cm–1 using the JFI method, which is in excellent agreement with the 12.4 cm–1 obtained by Jahr et al.[9] using the RPI. The mixing angle for the GS was estimated experimentally by Baughcum et al.[10] to be ϕ = 41° from the dipole moment measurements. The dipole moment of the D7/D9 isomer was modeled as a superposition of the dipole moments of D7 and D9 minima. These were, in turn, approximated by the dipole moments of D6D7D8 and D6D8D9 isomers, taken as two separate species, with the tunneling assumed to be suppressed. Reference (9) estimates the angle to be ϕ = 44° using local harmonic energies. Using VCI energies, we estimate the mixing angle to be ϕ = 35.3°, which indicates that the anharmonicity is indeed responsible for a decrease in its value, as speculated by Jahr et al.[9] We were also able to estimate the effect of the inclusion of other local vibrational states on the mixing angle from the components of the TM eigenvectors in Table aswhich gives ϕ = 36.8°. It thus appears that the inclusion of additional interactions corrects the mixing angle toward the experimental value, which was obtained indirectly, as stated above, and may carry a considerable error.

Table 8

Dominant Configurations of Vibrational States of PD Malonaldehyde, Obtained as the Eigenvectors of the TM in the 2 × 2 (Pairs) and 16 × 16 (Eight-State) Models, as Described in the Text

Changes in the vibrational levels of the excited states in the 8 × 8 and the 16 × 16 matrix models, listed in Table , are qualitatively similar to the homoisotopic malonaldehyde due to the similarity in their TM elements. The vibrational tunneling spectrum is shown graphically in Figure . One significant difference here is that some doublet states change the order of their components after the inclusion of additional vibrational states in the model due to their proximity in energy after deuteration, as seen in Figure . Another difference is the absence of symmetry in the wavefunctions with respect to the symmetry operation that connects the minima in the homoisotopic case. As a result, the extensions of the 2 × 2 model to higher-dimensionality matrix models will mix both the lower and higher components of doublets, with all other doublet states. Finally, due to the proximity of vibrational states, the lower components of the (ν1 + ν1), ν4, and (ν2 + ν2) doublets are significantly mixed, as can be seen in Table . This mixing between the states changes their energies, but it is also expected to affect the intensity of the transition to the 11th state as its ν4 component (see Table ) has a higher transition dipole moment, being the singly excited state.

Table 7

Vibrational Energy Levels of PD Malonaldehyde in Inverse Centimeters (cm–1) Obtained Using a Combined VCI/Instanton Approacha

no.	E(pairs)	E(4)	E(8)
1	13,974.27	13,972.91	13,971.48
2	14,016.96	14,015.52	14,014.49
3	14,296.86	14,298.42	14,286.79
4	14,295.75	14,295.44	14,294.91
5	14,314.69	14,313.95	14,313.54
6	14,391.55	14,392.78	14,381.17
7	14,383.28	14,384.05	14,383.32
8	14,446.16	14,446.45	14,445.64
9	14,540.00		14,537.31
10	14,549.82		14,549.32
11	14,541.38		14,556.57
12	14,560.53		14,561.08
13	14,597.33		14,598.21
14	14,596.99		14,598.33
15	14,641.54		14,642.66
16	14,696.90		14,709.18

E(pairs), E(4), and E(8) are energies obtained from the 2 × 2, 8 × 8, and 16 × 16 matrix models, respectively, as explained in the text.

Conclusions

We applied a combination of VCI and the instanton theory to calculate vibrational tunneling spectra of some exemplary double-well systems in full dimensionality at a much reduced computational cost in comparison with the exact QM methods. The VCI method was used to compute the single-well vibrational spectra, while the recently developed instanton method was used to determine the wavefunctions inside the barrier that separates the wells at a comparatively negligible computational cost. The interaction between the states of different wells was obtained from the Herring formula evaluated at a dividing surface inside the barrier. The Herring formula was rederived in an extended N × N matrix model (N > 2) , and the size of the associated leading error term was analyzed.

The accuracy of our approach was first tested on a model 2D system. It was shown that the JFI method can be used to compute TM elements that connect states in inequivalent wells and that have excitations in different normal modes. The energy levels of an asymmetric system exhibit avoided crossings with the variation of frequency or depth of one well relative to that of the other. The VCI calculation of local energies proved to be necessary in order to reproduce the exact QM results with high accuracy. The method was then tested on malonadehyde in full dimensionality, where good agreement was achieved with the exact MCTDH results in the absolute energies and the splittings. It was shown that the extension of the standard 2 × 2 model to include more states can influence the vibrational energies. The results are not affected dramatically in the case of malonaldehyde, but it was shown that the model is able to accommodate the additional vibrational states in the systems where they lie close in energy.

Finally, the method was used to calculate the vibrational spectrum of the low-lying states in PD malonaldehyde, which is near the computational limit of the presently available exact QM methods. The ground-state mixing angle was compared to the experimental value, and the influence of including additional vibrational states was shown to affect the angle and the order of some states in the spectrum.

The method is expected to perform well for mid-sized molecules, where rotational motion, which is neglected in this work, can be separated from the tunneling dynamics, and for moderately anharmonic systems with high barriers and, consequently, small TSs. It is exactly in these circumstances that the exact QM methods come at a prohibitive computational cost. The developed combined approach can be used to calculate the low-lying vibrational spectra in systems with an arbitrary number of wells, which are not necessarily related by symmetry. This makes the method particularly suitable to the studies of clusters, for example, for the assignment of spectra in water clusters, which feature multiple minima and high barriers in their bifurcation dynamics (where hydrogen bonds are broken and reformed). The computational cost of our approach is concentrated in solving the single-well spectra separately. The instanton theory can also be combined with other high-level methods, instead of VCI, and the combined, dual-level, approach can be used to calculate the tunneling spectra in general multidimensional asymmetric well systems, beyond molecular applications and chemistry.