## Full-dimensional (12D) variational vibrational states of CH 4 · F : -interplay of anharmonicity and tunneling ∗

Gustavo Avila 1, † and Edit Matyus 1, ‡

1 Institute of Chemistry, ELTE, E¨ otv¨s Lor´nd University, o a P´ azm´ny a P´ eter s´t´ny 1/A, 1117 Budapest, e a Hungary

(Dated: January 14, 2020)

## Abstract

The complex of a methane molecule and a fluoride anion represents a 12-dimensional (12D), fourwell vibrational problem with multiple large-amplitude motions, which has challenged the quantum dynamics community for years. The present work reports vibrational band origins and tunneling splittings obtained in a full-dimensional variational vibrational computation using the GENIUSH program and the Smolyak quadrature scheme. The converged 12D vibrational band origins and tunneling splittings confirm complementary aspects of the earlier full- and reduced-dimensionality studies: (1) the tunneling splittings are smaller than 0 02 . cm -1 ; (2) a single-well treatment is not sufficient (except perhaps the zero-point vibration) due to a significant anharmonicity over the wells; and thus, (3) a full-dimensional treatment appears to be necessary. The present computations extend to a higher energy range than earlier work, show that the tunneling splittings increase upon vibrational excitation of the complex, and indicate non-negligible 'heavy-atom' tunneling.

∗ Dedicated to Professor Attila Cs´sz´r on the occasion of his 60th birthday. a a

† Electronic address: Gustavo˙Avila@telefonica.net

‡ Electronic address: matyuse@caesar.elte.hu

## I. INTRODUCTION

The CH 4 · F -complex has been subject of experimental interest over the past decades. Its infrared spectrum has been recorded and studied in the methane's stretching region [15], and it has been used as a precursor in anion photoelectron spectroscopy to probe the transition state region of the F + CH 4 reaction [6, 7]. This experimental activity motivated the computational (ro)vibrational quantum dynamics study of the complex [8-10]. This complex also serves as a good prototype for molecular interactions with relatively large monomer distortions and strong binding.

CH 4 · F -has turned out to be challenging for the current (ro)vibrational methodologies, due to its high vibrational dimensionality and multi-well, highly anisotropic potential energy landscape. The vibrational states from Refs. [8], [9], and [10], using the MULTIMODE [11], the MCTDH [12, 13], and the GENIUSH [14, 15] quantum dynamics program packages, respectively, show several (tens of) wavenumbers (dis)agreement. In the present work, we aim to resolve this controversy.

There is currently a single, full-dimensional potential energy surface (PES) available for the complex developed by Czak´, Braams, and Bowman in 2008 [8], which we will refer to o as 'CBB08 PES'. The CBB08 PES was obtained by fitting permutationally invariant (up to 6th-order) polynomials to 6547 (plus 3000) electronic energy points of the interaction (plus fragment) region computed at the frozen-core CCSD(T)/aug-cc-pVTZ level of the ab initio theory. The root-mean-square deviation (rmsd) of the fitting, within the energy range below 22 000 cm -1 , was reported to be 42 cm -1 , and in practice, the PES describes the intermolecular region well up to moderate ion-molecule separations.

The equilibrium structure of the complex has C 3v point-group (PG) symmetry with the fluoride binding to one of the apexes of the methane tetrahedron. Since the F -anion can bind to any of the four hydrogens of methane, there are four equivalent minima on the PES, which are separated by 'surmountable' barriers, and thus, the molecular symmetry (MS) group is T d (M). The complex is bound by D e = 2434 cm -1 on the CBB08 PES [8], which corresponds to D 0 = 2316 cm -1 including the zero-point vibrational energy correction, and we have found the lowest barrier connecting the equivalent wells to be V = 1104 cm -1 (vide infra) . It is interesting to compare these values with the similar parameters of prototypical systems of hydrogen bonding and tunneling. The prototype of strong hydrogen bond, (HF) 2

features a dissociation energy of D 0 ≈ 1050 cm -1 [16], which is less than half of the binding energy of CH 4 · F . The prototype for (double-well) tunneling, malonaldehyde has a barrier -of 1410 cm -1 [17], which is ca. 30% higher than the barrier in CH 4 · F , but of course, for -estimating the tunneling splitting one has to consider the effective mass and also the shape (width) of the barrier.

The strong interaction of the methane and the fluoride in CH 4 · F -is accompanied by a relatively large distortion of the methane fragment. For the interaction (int) region, Table 3 of Ref. [8] reports the equilibrium structure with an elongated C-H bond, r int eq (C-H ) = b 1 112 ˚ , for the H which binds (b) to the F . A -, while for the other three hydrogens, r int eq (C-H) = 1 095 A. The corresponding distorted tetrahedral structure . ˚ is characterized by the α (HC-H ) = 110 46 b . o angle. In the CH +F 4 -channel, the practically isolated (isol) methane molecule is a regular tetrahedron with a C-H equilibrium distance, r isol eq (C-H) = 1 090 ˚ . . A

In the forthcoming sections, we briefly summarize the quantum dynamics methodology used in this work (Sec. II), explain the symmetry analysis and assignment of the vibrational states (Sec. III), and report the vibrational energies obtained in the fulland reduceddimensionality treatments (Sec. IV). After assessment of the convergence of the vibrational energies, a detailed comparison is provided with the vibrational energies reported in earlier studies [8-10]. The article is concluded (Sec. V) with the computation of tunneling splittings for excited vibrations of the complex up to and slightly beyond the energy range of the barrier separating the equivalent wells.

## II. THEORETICAL AND COMPUTATIONAL DETAILS

The present work is among the first applications of the GENIUSH-Smolyak algorithm and computer program [18]. The GENIUSH-Smolyak approach combines the non-product grid (Smolyak) method of Ref. [19], which has been used for several high-dimensional, semi-rigid molecules [20, 21] as well as molecules with a single large-amplitude motion [22, 23], and the numerical kinetic energy operator (numerical KEO) approach of Ref. [14] implemented in the GENIUSH program [14, 15], which includes by now dozens of vibrational-coordinate definitions for floppy systems [10, 14, 15, 24-30].

Concerning the coordinate definition for the fluoride-methane complex, we used the ( R, cos θ, φ ) spherical polar coordinates to describe the relative orientation of the methane fragment and the fluoride ion, and the nine normal coordinates, q , q 1 2 , . . . , q 9 , of methane (the coordinate definition is provided in the Supplementary Material) to describe its internal vibrations. In the full-dimensional computations, we used the KEO given in Eq. (50) of Ref. [18], which reads for the ( ξ , ξ 1 2 , ξ 3 , . . . , ξ D ) general coordinates with the special ξ 2 = c choice

/negationslash

/negationslash

$$\hat { T } ^ { \nu } & = - \frac { 1 } { 2 } \sum _ { j = 1 } ^ { D } \frac { \partial } { \partial c } G _ { c, j } \frac { \partial } { \partial \xi _ { j } } - \frac { 1 } { 2 } \sum _ { i = 1, i \neq 2 } ^ { D } \sum _ { j = 1 } ^ { D } G _ { i, j } \frac { \partial } { \partial \xi _ { i } } \frac { \partial } { \partial \xi _ { j } } - \frac { 1 } { 2 } \sum _ { i = 1 } ^ { D } B _ { i } \frac { \partial } { \partial \xi _ { i } } + U, \\ & B _ { i } & = \sum _ { k = 1, k \neq 2 } ^ { D } \frac { \partial } { \partial \xi _ { k } } G _ { k, i } \,, \\, & \quad - \sup ( N + \lambda \nu ( N + \lambda ) \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \, \overline { \omega } \quad - 1 \ \tilde { \omega } \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \end{array}$$

where g ∈ R ( D +3) × ( D +3) is the mass-weighted metric tensor, G = g -1 , ˜ g = det g , the extrapotential term,

$$U = \frac { 1 } { 3 2 } \sum _ { k = 1 } ^ { D } \sum _ { l = 1 } ^ { D } \left [ \frac { G _ { k l } } { \tilde { g } ^ { 2 } } \frac { \partial \tilde { g } } { \partial \xi _ { k } } \frac { } { \partial \xi _ { l } } + 4 \frac { \partial } { \partial \xi _ { k } } \left ( \frac { G _ { k l } } { \tilde { g } } \frac { \partial \tilde { g } } { \partial \xi _ { l } } \right ) \right ],$$

and the volume element is d V = d ξ 1 d c . . . d ξ D . In the full-dimensional treatment of CH 4 · F , -D = 12 and the coordinates are ξ 1 = R ξ , 2 = c = cos θ , ξ 3 = φ , ξ 3+ i = q i ( i = 1 2 , , . . . 9). We treat c = cos θ differently from the other coordinates in Eq. (1) in order to avoid a non-symmetric finite basis representation of the Hamiltonian due to inaccurate integration caused by singular terms in the KEO [18]. The Hamiltonian matrix was constructed using a finite basis representation (FBR) for all coordinates except c , for which the sin-cot discrete variable representation (DVR) [31] was used as it is explained in Sec. IV.E of Ref. [18]. The reduced-dimensionality computations have been carried out with the original GENIUSH

program [14], using the Podolsky form of the KEO (constructed in an automated way for the imposed geometrical constraints) and the Hamiltonian matrix was constructed using DVR [32]. The lowest eigenvalues and eigenfunctions of the Hamiltonian matrix were computed with an iterative Lanczos eigensolver.

Concerning the full-dimensional computations, it is necessary to reiterate some methodological details from Ref. [18] and to specify them for the case of the fluoride-methane complex. First of all, full-dimensional (12D) computations were possible for this complex because we used normal coordinates for the methane fragment together with harmonic oscillator basis functions, which provide a good zeroth-order description. Hence, the 9D product basis set of the methane fragment can be pruned [33], i.e., we can discard high-energy basis functions. We used the simple

$$\sum _ { i = 1 } ^ { 9 } n _ { q _ { i } } \leq b$$

pruning condition for the harmonic oscillator indexes, n q i . Since several basis functions are discarded from the methane basis set complying with this condition, it is possible to substantially reduce also the number of quadrature points which are used to calculate the overlap and low-order polynomial integrals with the retained basis functions. Pruning the grid following this observation was first realized by Avila and Carrington [19, 20] in vibrational computations using the Smolyak algorithm.

Concerning the 3-dimensional ion-molecule 'intermolecular' part, described by the R , cos θ , and φ coordinates, we retained the full direct-product basis and grid. For R , we used a Morse tridiagonal basis set constructed similarly to Ref. [18], but using the D 0 = 1975 27 cm . -1 , α = 0 9, and . γ = 18 parameter values which correspond to the 1D cut of the current PES (all other coordinates fixed at their equilibrium value). The cos θ degree of freedom was described with sin-cot-DVR basis functions and quadrature points [31], while we used Fourier functions for the φ angle.

We used large basis sets and grids for the intermolecular degrees of freedom, ( R, cos θ, φ ), both in the 3D and in the 12D computations, which is necessary to ensure that the degeneracies (some of them obtained numerically, only) and tunneling splittings are well converged (vide infra). In order to converge the energies with respect to the intramolecular (methane) part of the basis and grid, we have carried out computations with increasing values of the b parameter, b = 2 3 , , and 4, in the basis pruning condition, Eq. (3), and determined a

Smolyak grid which integrates exactly the overlap and fifth-order polynomials with all basis functions retained in the pruned basis.

Table I reports the computed vibrational energies in comparison with literature values. We estimate the vibrational excitation energies from the largest 12D computation, in column B C , to be converged within 1 cm -1 and the (apparently very small) tunneling splittings are converged with an uncertainty lower than 0.05 cm -1 . Note that for b = 2 and 3 we used only 23 sin-cot-DVR basis functions, B A and B B , respectively, which results in a 0.2-0.3 cm -1 (artificial) splitting for some higher excited states, but this splitting is reduced to less than 0 05 cm . -1 upon the increase of the basis set, which is reported in column B C of the table (25 sin-cot-DVR basis functions and b = 4 basis-pruning parameter).

Finally, we mention that we were able to put together a 'fitted' 3D model which reproduced the 12D GENIUSH-Smolyak vibrational band origins with an rmsd of 1.9 cm -1 (' G fit ' in the Table I). ' G fit ' was obtained by fine-tuning the regular tetrahedral methane structure used in the KEO and in the PES. In this 'fitting' procedure, we obtained r PES (C-H) = 2 143 624 bohr to define the PES cut (which is slightly different . from the value used in Ref. [10] and reproduces slightly better the 12D stretching excitations of the complex), but we had to use a (much!) larger value in the KEO, r KEO (C-H) = 2 518 620 bohr . (which corresponds to a much smaller effective rotational constant for the methane) to have a good agreement with the vibrational band origins especially for the bending excitations.

Note that '3D( G fit )' is merely a fitted model, which was designed to reproduce the 12D GENIUSH-Smolyak vibrational band origins, and for which we deliberately used two adjustable parameters, the C-H distance in the KEO and in the PES. It differs in this twoparameter adjustment scheme from the 'rigorous' 3D reduced-dimensionality treatment (for which only a single structure is selected) used in the first 3D computation of CH 4 · F -in Ref. [10].

## III. SYMMETRY ANALYSIS AND ASSIGNMENT OF THE VIBRATIONAL STATES

First of all, we assigned the computed states to irreducible representations (irreps) of the MS group of the complex. Then, we identified single molecular vibrations classified according to the PG symmetry of the equilibrium structure as groups of states close in energy (slightly)

split due the interaction (spread) of the wave function over the equivalent wells.

At the equilibrium structure of C 3v PG symmetry, one of the hydrogens of the methane binds to the fluoride anion. Any of the four hydrogen atoms of the methane can bind to the fluoride, which gives rise to the four equivalent wells and these wells are connected with 'surmountable' barriers. Thereby, the MS group of the complex is T d (M), for which the symmetry analysis of (the global minimum of) CH 4 · Ar [34] is almost verbatim adapted.

In order to assign irrep labels to the CH 4 · F -vibrational states computed in the present work, we analyzed the wave function of the 3D fitted model computations ('3D G fit ' column in Table I) and the labels were transferred to the 12D results based on the energy ordering (direct analysis of the 12D hybrid DVR-FBR computation would have been prohibitively expensive). We assigned T d (M) molecular symmetry labels to the 3D wave functions by computing their overlap with 2D coupled-rotor (CR) functions, labelled with [ j, j ] 00 ( j = 0 1 , , . . . ) [28, 34], where j is the angular momentum quantum number of the methane and the diatom (corresponding to the relative motion of the center of mass of the methane and the fluoride), coupled to a zero total angular momentum state. The characters and the irrep decomposition of the Γ CR ( j ) representation spanned by the [ j, j ] 00 coupled-rotor functions in T d (M) are [34]:

$$\mathbb { A } \text{;} \\ \Gamma ^ { \text{CR} } ( 0 ) & = A _ { 1 } \,, \\ \Gamma ^ { \text{CR} } ( 1 ) & = F _ { 2 } \,, \\ \Gamma ^ { \text{CR} } ( 2 ) & = E \oplus F _ { 2 } \,, \\ \Gamma ^ { \text{CR} } ( 3 ) & = A _ { 1 } \oplus F _ { 1 } \oplus F _ { 2 } \,, \\ \Gamma ^ { \text{CR} } ( 4 ) & = A _ { 1 } \oplus E \oplus F _ { 1 } \oplus F _ { 2 } \,, \\ \Gamma ^ { \text{CR} } ( 5 ) & = E \oplus F _ { 1 } \oplus F _ { 2 } \,, \\ \Gamma ^ { \text{CR} } ( 6 ) & = A _ { 1 } \oplus A _ { 2 } \oplus E \oplus F _ { 1 } \oplus 2 F _ { 2 } \,, \dots \quad \quad ( 4 ) \\ \text{ed} \text{ relative rotation of the molecule and the ion over the four wells cycles rise}$$

The (hindered) relative rotation of the molecule and the ion over the four wells gives rise to tunneling splittings of the vibrations, which can be classified by C 3v point-group labels of the symmetry of the local minima. (If there was no interaction between the wells, each vibrational state would be 4-fold degenerate due to this feature.) The MS group species within the tunneling manifold of the vibrational modes classified by the PG symmetries

(irreps) are [34]

$$\Gamma ( A _ { 1 } ^ { C 3 v } ) & = A _ { 1 } \oplus F _ { 2 } \,, \\ \Gamma ( A _ { 2 } ^ { C 3 v } ) & = A _ { 2 } \oplus F _ { 1 } \,, \\ \Gamma ( E ^ { C 3 v } ) & = E \oplus F _ { 1 } \oplus F _ { 2 } \,.$$

Note that we use the C 3v superscript for the PG irreps, i.e., A C 3v 1 , A C 3v 2 , and E C 3v , in order to distinguish them from the MS group irreps, which are labelled with A , A , E, F , and 1 2 1 F . 2 The result of this analysis for the computed vibrational wave functions is summarized in the 'Γ(MS)' and 'Γ(PG)' columns of Table I, respectively.

The fourth column, ' n R ', of the table gives the index of the wave function of the 1D model with active R (all other coordinates fixed at their equilibrium value) for which the 3D wave function has the largest overlap. Hence, ' n R ' is an index for the excitation along the ion-molecule separation, which we were able to unambiguously assign for all states listed in the table (due to the weak coupling of the the radial and angular degrees of freedom in this complex).

We also note that due to the small tunneling splittings, identification of the PG vibrations as a set of states of similar character and close in energy comprising the appropriate MS group species, Eq. (5), was possible without ambiguities for all states listed in the table. (This may be contrasted with the floppy CH 4 · Ar complex, for which unambiguous assignment of the PG vibrations beyond the zero-point state of the global minimum was hardly possible [34].) Once, the complete tunneling manifold was assigned and the PG symmetry was found, we attached a qualitative description to the states (listed in column 'Label' in Table I) based on the nodal structure along the ion-molecule separation coordinate R , i.e., n th 'stretching' excitation, labelled with v s , 2 v s , etc. Excitations different from pure stretching excitation were termed 'bending' in this qualitative description, and were labelled with v b , 2 v b , etc., or combinations of stretching and bending, v s + v b , etc. Similar qualitative labels had been provided in the earlier studies [9, 10], which we used to compare with our full-dimensional vibrational energies.

## IV. COMPUTED VIBRATIONAL STATES

The 12D vibrational energies computed in the present work are reported in Table I (column '12D, GENIUSH-Smolyak'), and the best results are listed in column ' B C : (25,4)' of the table. The level of convergence and comparison with earlier work is addressed in the following paragraphs.

We observe a monotonic decrease of the zero-point vibrational energy (ZPVE) upon the increase of the b basis pruning parameter, which indicates that the computations are almost variational (the sin-cot-DVR basis and grid is large). Our best ZPVE value is 9791.6 cm -1 . It is 5.1 cm -1 larger than the 12D MCTDH result, which the authors of Ref. [9] claim to be a variational upper bound to the exact ZPVE. The 12D MULTIMODE ZPVE [8] is 3.1 cm -1 (8.2 cm -1 ) larger than the lowest value obtained by us (by Ref. [9] using MCTDH). The totally symmetric ZPV state (A C 3v 1 ) is split by the relative rotation of the methane and the fluoride to an A 1 and an F 2 symmetry species, Eq. (5). The B C 12D computation (numerically) reproduces the degeneracy of the F 2 state within 0.001 cm -1 and predicts a tunneling splitting (much) smaller than 0.05 cm -1 .

At the same time, the MCTDH result for the ZPV manifold gives four states with energies differing by 0.1-0.9 cm -1 (up to 1.2-3.4 cm -1 , depending on the basis set), which suggests that the (intermolecular) basis or integration grid used in Ref. [9] was too small. For higher excited states the (artificial) splittings increase to 2 cm -1 [9], whereas our computations, using a large intermolecular basis set and grid, indicate that the tunneling splittings are smaller than 0.05 cm -1 , and we obtain the triple degeneracies converged (numerically) better than 0.001 cm -1 .

The small tunneling splittings obtained in the present work are in agreement with the earlier, 3D computations including the intermolecular, ( R, cos θ, φ ) coordinates as active vibrational degrees of freedom using the GENIUSH program [10]. Although our results confirm the small splittings obtained in the 3D computation of Ref. [10], we observe larger deviations for the vibrational excitation energies (band origins) from the 3D results of Ref. [10]. The first 7 vibrational excitations have a 23.5 cm -1 root-mean-square deviation (rmsd) from our best 12D energies, which suggests that monomer (methane) flexibility effects are important.

A better agreement, with an 11.7 cm -1 rmsd, is observed for the first 6 vibrational excitations in comparison with the 12D MCTDH result [9], which accounts for the flexibility

of the methane fragment and the motion of the fluoride over the four wells, but it is affected by incomplete convergence (manifested in the artificial splittings).

The three lowest-energy vibrational fundamentals obtained within a 12D but single-well treatment with MULTIMODE [8] has a (surprisingly) large, 20.3 cm -1 rmsd, which suggests that in spite of the small tunneling splittings, a multi-well treatment is necessary for which the normal-coordinate representation is inadequate. It should be noted that the largest deviation, 32.3 cm -1 , is observed for the 299.9 cm -1 vibration, and this is the state which had the worse convergence properties in the MULTIMODE computation (Table 5 of Ref. [8]).

These observations suggest that although tunneling (and the corresponding splittings of the vibrational bands) is almost negligible up to 725 cm -1 above ZPVE under an energy resolution of 0.05 cm -1 , there is a non-negligible anharmonicity due to the quantum mechanical motion over the multiple wells. For these reasons, a well converged, full-dimensional variational treatment, carried out in the present work, appears to be necessary to capture all quantum dynamical features of this strongly bound ion-molecule complex.

TABLE I: Vibrational states, ˜ in cm ν -1 , of CH 4 · F -up to 730 cm -1 above the zero-point vibrational energy (ZPVE), corresponding to the full-dimensional PES of Ref. [8]. Vibrational energies computed in the present work and values taken from literature are shown together for comparison. The largest computed splitting for each vibrational manifold is given in parenthesis after the vibrational energy value. The most accurate vibrational band origins (from this work) are in column ' B C : (25,4)'.

|       |               |          |       |               |             | GENIUSH-Smolyak   | [this            | GENIUSH-Smolyak   | GENIUSH-Smolyak        | 12D MCTDH B3/B4 [9]   | 12D MULTIMODE [8]   | 3D ∗ GENIUSH [10]   | 3D GENIUSH [this work]   |
|-------|---------------|----------|-------|---------------|-------------|-------------------|------------------|-------------------|------------------------|-----------------------|---------------------|---------------------|--------------------------|
|       |               |          |       |               | B A : (23 , | B A : (23 ,       | 2) e,f B B : (23 | 2) e,f B B : (23  | e,f B C : (25 , 4) e,f |                       |                     |                     |                          |
| #     | Γ(MS) a       | Γ(PG) b  | n R c | Label d       | [3 . 9 ·    | 10 5 ] e          | [1 . 6           | · 10 6 ] e        | [5 . 6 · 10 6 ] e      |                       |                     |                     | G fit f,g                |
| 0-3   | A 1 ⊕ F 2     | A C 3v 1 | 0     | ZPVE          | 9845.6      | (0.0)             | 9799.1           | (0.0)             | 9791.6 (0.0)           | 9786.5(0.5)           | 9794.7              | 461.0 (0.0)         | 378.8 (0.0)              |
| 4-7   | A 1 ⊕ F 2     | A C 3v 1 | 1     | [ v s ]       | 193.2       | (0.0)             | 193.4            | (0.0)             | 193.6 (0.0)            | 194.4 (0.1)           | 201.1               | 182.5 (0.0)         | 194.5 (0.0)              |
| 8-15  | E ⊕ F 1 ⊕ F 2 | E C 3v   | 0     | [ v b ]       | 266.3       | (0.0)             | 268.3            | (0.0)             | 267.6 (0.0)            | 271.7 (1.1)           | 299.9               | 284.5 (0.0)         | 267.7 (0.0)              |
| 16-19 | A 1 ⊕ F 2     | A C 3v 1 | 2     | [2 v s ]      | 378.0       | (0.0)             | 378.7            | (0.0)             | 379.2 (0.0)            | 380.6 (0.1)           | 391                 | 355.8 (0.0)         | 380.3 (0.0)              |
| 20-27 | E ⊕ F 1 ⊕ F 2 | E C 3v   | 1     | [ v s + v b ] | 452.7       | (0.0)             | 454.8            | (0.0)             | 454.3 (0.0)            | 460.2 (1.3)           |                     | 458.8 (0.0)         | 455.2 (0.0)              |
| 28-31 | A 1 ⊕ F 2     | A C 3v 1 | 0     | [2 v b ]      | 506.7       | (0.0)             | 509.8            | (0.0)             | 509.1 (0.0)            | 528.7 (0.4)           |                     | 533.0 (0.1)         | 509.8 (0.0)              |
| 32-39 | E ⊕ F 1 ⊕ F 2 | E C 3v   | 0     | [2 v b ]      | 523.5       | (0.2) h           | 527.0            | (0.2) h           | 526.0 (0.0)            | 545.7 (2.1)           |                     | 555.3 (0.1)         | 525.6 (0.0)              |
| 40-43 | A 1 ⊕ F 2     | A C 3v 1 | 3     | [3 v s ]      | 555.0       | (0.0)             | 556.4            | (0.0)             | 557.5 (0.0)            |                       |                     | 519.2 (0.0)         | 556.8 (0.0)              |
| 44-51 | E ⊕ F 1 ⊕ F 2 | E C 3v   | 2     | [2 v s + v b  | ] 630.3     | (0.0)             | 632.8            | (0.0)             | 632.7 (0.0)            |                       |                     |                     | 633.6 (0.0)              |
| 52-55 | A 1 ⊕ F 2     | A C 3v 1 | 1     | [ v s +2 v b  | ] 685.7     | (0.0)             | 689.0            | (0.0)             | 688.5 (0.0)            |                       |                     |                     | 690.2 (0.0)              |
| 56-63 | E ⊕ F 1 ⊕ F 2 | E C 3v   | 1     | [ v s +2 v b  | ] 702.6     | (0.2) h           | 706.2            | (0.2) h           | 705.5 (0.0)            |                       |                     |                     | 705.9 (0.0)              |
| 64-67 | A 1 ⊕ F 2     | A C 3v 1 | 4     | [4 v s ]      | 724.2       | (0.0)             | 726.6            | (0.0)             | 728.5 (0.0)            |                       |                     |                     | 722.9 (0.0)              |
| rmsd  | i             |          |       |               | 2.4         |                   | 0.8              |                   | 0                      | 11.7                  | 20.3                | 23.5                | 1.9                      |

## Footnotes to Table I:

- a Symmetry assignment (irrep decomposition) in the T d (M) molecular symmetry group of the complex.
- b Symmetry assignment (irrep) within the C 3v point group of the equilibrium structure.
- c Dominant overlap of the wave function with the n R th state of a 1-dimensional vibrational model along the R degree of freedom (all other coordinates are fixed at their equilibrium value). n R = 0 1 2 , , , . . . labels the states of this 1D model in an increasing energy order.
- d Qualitative description based on the nodal structure, overlap with lower-dimensional models, and symmetry assignment. These labels are used to compare the vibrational
- energies computed in the present work with earlier results [8-10].
- e ( N , b c ): short label for indicating the basis set size. The N c value gives the number of the sin-cot-DVR functions used for cos θ and b is the basis pruning parameter, ∑ 9 i =1 n q i ≤ b . The R and φ degrees of freedom are described by N R = 8 Morse tridiagonal and N φ = 39 Fourier basis functions, respectively. The overall basis set size, [ N N N R c φ ( b +9)! /b / ! 9!], is also shown.
- f In the vibrational computations we used atomic masses, m (H) = 1 007 825 032 23 u, . m (C) = 12 u, and m (F) = 18 998 403 162 73 u [35]. .
- g Reduced-dimensionality model with active ( R, cos θ, φ ) degrees of freedom fitted to reproduce the 12D GENIUSH-Smolyak B C result. The methane was treated as a regular tetrahedron and we used r PES (C-H) = 2 143 624 bohr in the PES (to reproduce well the . stretches) and r KEO (C-H) = 2 518 620 bohr in the KEO, which was necessary to obtain . good bending energies.

h These splittings disappear upon increase of the intermolecular basis set.

i Root-mean-square deviation of the vibrational excitation energies (ZPVE not included) listed in the table (without considering the splittings) from the 12D GENIUSH-Smolyak B C result.

∗ Higher energy vibrations, including the 12 fundamental vibrations, are reported in Ref. [8], but they are beyond the (energy) range of the present GENIUSH-Smolyak computations.

## V. EXCITED-STATE TUNNELING MANIFOLDS

For the ZPV and the lower-energy vibrations, we obtained tiny tunneling splittings ( &lt; 0 05 cm . -1 , Sec. IV). The question arises whether the energy splitting due to tunneling of the heavy fragments (heavier than the hydrogen atom) becomes more significant upon the vibrational excitation of the complex. The symmetry species in the tunneling manifold are specified in Eq. (5), but of course, the symmetry analysis by itself does not provide any information about the level energies and the size of the splittings.

The height of the barrier, separating the energy wells, can be indicative for the energy range above which we may expect larger splittings. Ref. [9] estimated the barrier height connecting the equivalent wells to be ca. 1270 cm -1 , and we have pinpointed a lower barrier height, V barrier = 1104 cm -1 at ( θ, φ ) = (90 o , 45 o ) (see Figure 2 of Ref. [10]), R = 5 745 bohr, . and ( q , q 1 2 , . . . , q 9 ) = (0 , 0 0 0 , , , -0 523 0 0 0 . , , , , -0 524). It is an interesting question what the . 'most appropriate' value for the barrier height is (which is not an experimental observable) for comparison with our variational (fully anharmonic, multi-well) vibrational energies. In the present work, we have not attempted to include quantum nuclear corrections in the height and will continue using the purely electronic value. So, we expect the appearance of 'sizeable' splittings at around ca. 1100 cm -1 measured from the PES minima, which is ca. 700 cm -1 above the intermolecular ZPVE (note that the different intermolecular models have somewhat different ZPVEs). Furthermore, it is necessary to remember that the methane's bending vibration becomes important in the energy range near the isolated methane's value at 1311 cm -1 (above the ZPVE), which is an obvious limitation for the applicability of the 3D models and approaching this energy range will require a more efficient 12D treatment.

In order to explore the tunneling splittings and their dependence on the vibrational excitation, we have carried out extensive 3D and more limited (up to only lower vibrational excitations) 12D computations with large, (near) saturated angular basis sets (vide infra). The tunneling splittings obtained from the different computations are visualized in Figures 1 and 2, and the energy lists are provided in the Supplementary Material.

We used the two types of 3D vibrational models to explore the excited-state tunneling manifold. First, we have carried out computations with the 3D reduced-dimensionality model used already in Ref. [10], cited in Table I, which corresponds to imposing rigorous

geometrical constraints in the model through the Lagrangian in GENIUSH [14]. The fixed methane regular tetrahedral structure used in this 3D(rig) model was r (C-H) = 1 104 ˚ . A [10], a value close to the isolated methane's 〈 r (C-H) 〉 0 value. Second, we have carried out computations using the 3D( G fit ) model designed ('fitted') to minimize the rmsd of the vibrational band origins with respect to the 12D GENIUSH-Smolyak results up to 730 cm -1 (Table I). In the largest 3D computations, we used 31 PO-DVR functions [36-38] for R obtained from L ( α ) n generalized Laguerre basis functions (with α = 2 and n = 0 1 , , . . . , 399) and grid points scaled to the R ∈ [2 , 5] bohr interval; 61 sin-cot-DVR functions for the cos θ degree of freedom, which was constructed by extending the cos( mθ ) ( m = 0 , ..., 58) basis set with the sin( θ ) and the sin(2 θ ) functions [18, 31]; and 121 Fourier (sine and cosine) functions for the φ degree of freedom. This basis is highly saturated and represents a space comparable with the set of the Y m j ( θ, φ ) spherical harmonic functions with j = 0 , .., 60 and m = -j, ...., j . These computations (for both 3D models) reproduce the triple degeneracies, which are spatial symmetry features not covered by the DVR grid but obtained numerically, within 0 000 1 cm . -1 for all computed states.

In the 12D computation of the excited state splittings, we used a large angular representation (relative methane-ion rotation) to tightly converge the splittings but we had to keep the methane basis small (due to the finite computational resources), which was sufficient to capture methane's flexibility effects on the splittings. So, we chose b = 2 in the pruning condition of the methane basis (9D), we used 8 tridiagonal-Morse basis functions and 10 GaussLaguerre-Morse points for R ; and a large basis and grid for the angular motion including 31 sin-cot-DVR basis functions and points for cos θ ; and 45 Fourier basis functions and 48 trapezoidal points for φ . So, this 12D basis set designed for the computation of the tunneling splittings (and henceforth labelled with B T ) included (2+9)! / 2! / 9! × × 8 31 × 45 = 613 800 functions and the 12D Smolyak grid ( H = 16, D = 12 dimensions) contained 3 481(for the CH 4 monomer) × 10 × 31 × 48 = 51 797 280 points to calculate the integrals. The angular part of B T is more than 40 % larger than that of B C in Table I, which provides the current benchmark for the vibrational band origins up to 730 cm -1 . Although the uncertainty of the vibrational band origins for B T is 3-4 cm -1 (due to the small b value), the 12D tunneling splittings are converged with an uncertainty better than 0.02 cm -1 up to 856 cm -1 , including 96 vibrational states.

The 12D computations show that the tunneling splittings increase upon the vibrational

excitation of the complex and we find several splittings larger the 0.02 cm -1 uncertainty threshold above ca. 680 cm -1 measured from the ZPVE (Figure 1). It is interesting to note that this value is comparable with the value of the electronic barrier separating the equivalent wells, which is 1104 cm -1 from the PES minimum and ca. 700 cm -1 measured from the intermolecular ZPVE (different models have different ZPVEs).

Figure 2 presents the 12D and the 3D splittings. The 3D( G fit ) model, which reproduces the known 12D vibrational band origins very well (Table I), predicts tunneling splittings orders of magnitude smaller than the 12D computation (notice the logarithmic scale in Figure 2). At the same time, the 3D(rig) model (with rigorous geometrical constraints) appears to get the splittings right (currently, we have only a few points to compare) while it has a large error in predicting the vibrational centers (Table I) which may hinder the comparison of the 3D(rig) and 12D splittings at higher energies.

The differences between the different 3D treatments can be rationalized in terms of the difference of the 'effective' rotational constant of methane in the computation (the small differences in the PES cut does not affect these results). In the 3D treatments, a 'smaller' effective methane rotational constant appears to be necessary to have good vibrational centers, 3D( G fit ), and a 'larger' methane rotational constant is necessary to have good tunneling splittings, 3D(rig). Altogether, these observations suggest that in order to have a correct overall description, it is necessary to explicitly account for the methane's vibrations together with the relative rotation of the methane and the fluoride ion.

FIG. 1: 12D vibrational energies measured from the zero-point vibration. Vibrational levels with tunneling splittings, ∆˜, larger than 0.02 cm ν -1 are highlighted in red. The position of the vibrational band origins may have an overall uncertainty as large as 3-4 cm -1 , but the tunneling splittings within a band are converged better than 0.02 cm -1 (see Figure 2).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00809v2-figure-6c453656-fef8-4624-8cb0-228ad11e811f.jpg)

FIG. 2: Logarithm of the energy splitting, ∆˜ in cm ν -1 , among the different symmetry species of the tunneling manifold, Eq. (5), shown with respect to the ˜ energy of the vibrational state ν measured from the zero-point vibrational energy (ZPVE). The splittings from the rigorous 3D(rig) model, in green triangles, and the fitted 3D( G fit ) model, in blue squares, are converged with an uncertainty of 0.000 1 cm -1 . The 12D tunneling splittings, in red circles, are converged better than 0.02 cm -1 (the splittings smaller than this value may be affected by numerical artifacts and are labelled with ∗ in the figure). The lowest electronic barrier on the CBB08 PES, 1104 cm -1 above the minimum, measured from the intermolecular ZPVE is indicated with the grey shaded area (note that the intermolecular ZPVEs are different for the different models). The full list of the computed energies is available in the Supplementary Material.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00809v2-figure-fa1279c2-1355-47ee-b1ad-ad627181af81.jpg)

## VI. SUMMARY AND CONCLUSIONS

Full-dimensional (12D), near-variational, vibrational states are reported for the strongly bound complex of the methane molecule and the fluoride anion. This is the first application of the recently developed GENIUSH-Smolyak algorithm and computer program [18] with a fully coupled, high-dimensional potential energy surface.

Benchmark-quality vibrational band origins are computed with energies and tunneling splittings converged better than 1 cm -1 and 0.05 cm -1 , respectively. These computations confirm complementary aspects of earlier studies [8-10], which relied on different assumptions or approximations about the dynamics of this complex. Regarding controversial aspects of earlier work, we can confirm that tunneling splittings (up to 730 cm -1 above the zero-point energy) are small [10], &lt; 0 05 cm . -1 , but due to the strong binding and significant monomer distortions, a full-dimensional (12D) treatment is necessary [8, 9]. Although the tunneling splittings are small, a single-well (normal coordinate) description [8] is not sufficient due to the significant anharmonicity of the multi-well potential energy landscape [9]. Even so, it is necessary to use a large basis set (and grid) to properly describe the relative, hindered rotation of the methane fragment and the fluoride anion [10].

We also show from 12D computations that the tunneling splittings, increase with the vibrational excitation of the complex. Both the 12D and approximate 3D rigid-monomer computations indicate that sizeable tunneling splittings, &gt; 0 1 cm . -1 , appear near the 'top of the barrier' which separates the equivalent wells. Interestingly, a 'fitted' 3D model, which gives excellent vibrational band origins, predicts orders of magnitude smaller splittings than the 12D computation, while the rigorous, reduced-dimensionality 3D model gives reasonable splittings, but fails to reproduce the 12D vibrational band energies.

With further progress of the quantum dynamics methodology reported in the present work and a potential energy surface with a broader coverage of the nuclear coordinates (especially, for the ion-molecule separation), it will become possible to study the effect of the methane vibrational excitation on the tunneling dynamics and the intra- to inter-molecular energy transfer under predissociation.

## Supplementary Material

Tunneling splittings and vibrational energies obtained from 12D and 3D computations and

the definition of normal coordinates of methane are provided in the Supplementary Material.

## Acknowledgment

Financial support of the Swiss National Science Foundation through a PROMYS Grant (no. IZ11Z0 166525) is gratefully acknowledged. We also thank NIIFI for providing us computer time at the Miskolc node of the Hungarian Computing Infrastructure. EM is thankful to ETH Z¨ urich for supporting a stay as visiting professor during 2019 and the Laboratory of Physical Chemistry for their hospitality, where part of this work has been completed.

- [1] D. A. Wild, Z. M. Loh, and E. Bieske, Int. J. Mass. Spectrom. 220 , 273 (2002).
- [2] Z. M. Loh, R. L. Wilson, D. A. Wild, E. J. Bieske, and M. S. Gordon, Aust. J. Chem. 57 , 1157 (2004).
- [3] Z. M. Loh, L. Wilson, D. A. Wild, J. Bieske, J. M. Lisy, B. Njegic, and M. S. Gordon, J. Phys. Chem. A 110 , 13736 (2006).
- [4] D. M. Neumark, J. Phys. Chem. A 112 , 13287 (2008).
- [5] T. I. Yacovitch, E. Garand, J. B. Kim, C. Hock, T. Theis, and D. M. Neumark, Faraday Discuss. 157 , 399 (2012).
- [6] M. Cheng, Y. Feng, Y. Du, Q. Zhu, W. Zheng, G. Czak´, and J. M. Bowman, J. Chem. Phys. o 134 , 191102 (2011).
- [7] T. Westermann, J. B. Kim, M. L. Weichman, C. Hock, T. I. Yacovitch, J. Palma, D. M. Neumark, and U. Manthe, Angew. Chem. Int. Ed.. 53 , 1122 (2014).
- [8] G. Czak´ o, B. J. Braams, and J. M. Bowman, J. Phys. Chem. A 112 , 7466 (2008).
- [9] R. Wodraszka, J. Palma, and U. Manthe, J. Phys. Chem. A 116 , 11249 (2012).
- [10] C. F´bri, A. G. Cs´sz´r, and G. Czak´, J. Phys. Chem. A a a a o 117 , 6975 (2013).
- [11] J. M. Bowman, S. Carter, and X. Huang, International Reviews in Physical Chemistry 22 , 533 (2003).
- [12] H.-D. Meyer, F. Gatti, and G. A. Worth, MCTDH for Density Operator (Wiley-Blackwell, 2009), chap. 7, pp. 57-62, ISBN 9783527627400.
- [13] M. Beck, A. Jackle, G. Worth, and H.-D. Meyer, Physics Reports 324 , 1 (2000), ISSN

0370-1573.

- [14] E. M´ atyus, G. Czak´, and A. G. Cs´sz´r, J. Chem. Phys. o a a 130 , 134112 (2009).
- [15] C. F´bri, E. M´tyus, and A. G. Cs´sz´r, J. Chem. Phys. a a a a 134 , 074105 (2011).
- [16] M. Quack and M. A. Suhm, J. Chem. Phys. 95 , 28 (1991).
- [17] W. Mizukami, S. Habershon, and D. P. Tew, J. Chem. Phys. 141 , 144310 (2014).
- [18] G. Avila and E. M´tyus, J. Chem. Phys. a 150 , 174107 (2019).
- [19] G. Avila and J. T. Carrington, J. Chem. Phys. 131 , 174103 (2009).
- [20] G. Avila and J. T. Carrington, J. Chem. Phys. 134 , 054126 (2011).
- [21] G. Avila and J. T. Carrington, J. Chem. Phys. 134 , 064101 (2011).
- [22] D. Lauvergnat and A. Nauts, 119 , 18 (2014), ISSN 1386-1425.
- [23] A. Nauts and D. Lauvergnat, Mol. Phys. 116 , 3701 (2018).
- [24] C. F´bri, E. M´tyus, and A. G. Cs´sz´r, Spectrochim. Acta a a a a 119 , 84 (2014).
- [25] C. F´bri, J. Sarka, and A. G. Cs´sz´r, J. Chem. Phys. a a a 140 , 051101 (2014).
- [26] J. Sarka and A. G. Cs´sz´r, J. Chem. Phys. a a 144 , 154309 (2016).
- [27] J. Sarka, A. G. Cs´sz´r, S. C. Althorpe, D. J. Wales, and E. M´tyus, Phys. Chem. Chem. a a a Phys. 18 , 22816 (2016).
- [28] J. Sarka, A. G. Cs´sz´r, and E. M´tyus, Phys. Chem. Chem. Phys. a a a 19 , 15335 (2017).
- [29] C. F´bri, M. Quack, and A. G. Cs´sz´r, J. Chem. Phys. a a a 147 , 134101 (2017).
- [30] I. Simk´ o, T. Szidarovszky, and A. G. Cs´sz´r, J. Chem. Theory Comput. a a 15 , 4156 (2019).
- [31] G. Schiffel and U. Manthe, Chem. Phys. 374 , 118 (2010), ISSN 0301-0104.
- [32] J. C. Light and T. Carrington Jr., Discrete-Variable Representations and their Utilization (John Wiley &amp; Sons, Ltd, 2007), pp. 263-310.
- [33] R. J. Whitehead and N. C. Handy, J. Mol. Spectrosc. 55 , 356 (1975).
- [34] D. Ferenc and E. M´tyus, Mol. Phys. a 117 , 1694 (2019).
- [35] J. S. Coursey, D. J. Schwab, J. J. Tsai, and R. A. Dragoset, Atomic Weights and Isotopic Compositions (version 4.1): http://physics.nist.gov/Comp [last accessed on 12 May 2018]. National Institute of Standards and Technology, Gaithersburg, MD. (2015).
- [36] H. Wei and T. Carrington, Jr., J. Chem. Phys. 97 , 3029 (1992).
- [37] J. Echave and D. C. Clary, Chem. Phys. Lett. 190 , 225 (1992).
- [38] V. Szalay, G. Czak´ o, A. Nagy, T. Furtenbacher, and A. G. Cs´ asz´r, a J. Chem. Phys. 119 , 10512 (2003).
Gustavo Avila Edit Matyus
Full-dimensional (12D) variational vibrational states of CH$_4\cdot$F$^-$: interplay of anharmonicity and tunneling
The complex of a methane molecule and a fluoride anion represents a
12-dimensional (12D), four-well vibrational problem with multiple
large-amplitude motions, which has challenged the quantum dynamics community
for years. The present work reports vibrational band origins and tunneling
splittings obtained in a full-dimensional variational vibrational computation
using the GENIUSH program and the Smolyak quadrature scheme. The converged 12D
vibrational band origins and tunneling splittings confirm complementary aspects
of the earlier full- and reduced-dimensionality studies: (1) the tunneling
splittings are smaller than 0.02 cm$^{-1}$; (2) a single-well treatment is not
sufficient (except perhaps the zero-point vibration) due to a significant
anharmonicity over the wells; and thus, (3) a full-dimensional treatment
appears to be necessary. The present computations extend to a higher energy
range than earlier work, show that the tunneling splittings increase upon
vibrational excitation of the complex, and indicate non-negligible `heavy-atom'
tunneling.