Full-dimensional (12D) variational vibrational states of CH$_4\cdot$F$^-$: interplay of anharmonicity and tunneling

TL;DR

This study presents full-dimensional variational vibrational calculations of CH$_4$·F$^-$, revealing small tunneling splittings, significant anharmonicity, and the necessity of 12D treatment for accurate vibrational states and tunneling phenomena.

Contribution

It provides the first full-dimensional (12D) vibrational analysis of CH$_4$·F$^-$, demonstrating the importance of anharmonicity and tunneling effects using advanced computational methods.

Findings

Tunneling splittings are smaller than 0.02 cm$^{-1}$.

Single-well models are insufficient due to anharmonicity.

Tunneling splittings increase with vibrational excitation.

Abstract

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.

Tables

Table 1. Table 1: Vibrational states, ν ~ ~ 𝜈 \tilde{\nu} in cm − 1 superscript cm 1 \text{cm}^{-1} , of CH ⋅ 4 {}_{4}\cdot F - up to 730  cm − 1 superscript cm 1 \text{cm}^{-1}  above the zero-point vibrational energy (ZPVE), corresponding to the full-dimensional PES of Ref.  Czakó et al. ( 2008 ) . 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 ‘ ℬ C subscript ℬ 𝐶 \mathcal{B}_{C} : (25,4)’.

					12D						12D		12D	3D		3D
					GENIUSH–Smolyak						MCTDH		MULTIMODE∗	GENIUSH		GENIUSH
					[this work]						B3/B4 Wodraszka et al. (2012)		Czakó et al. (2008)	Fábri et al. (2013)		[this work]
					${ℬ}_{\text{A}}$: $(23,2)$${}^{\text{e,f}}$		${ℬ}_{\text{B}}$: $(23,3)$${}^{\text{e,f}}$		${ℬ}_{\text{C}}$: $(\mathrm{𝟐𝟓},\mathrm{𝟒})$${}^{\text{e,f}}$
#	$\Gamma \left(\text{MS}\right)$${}^{\text{a}}$	$\Gamma \left(\text{PG}\right)$${}^{\text{b}}$	$n_R$${}^{\text{c}}$	Label${}^{\text{d}}$	$\left[3.9\cdot {10}^5\right]$${}^{\text{e}}$		$\left[1.6\cdot {10}^6\right]$${}^{\text{e}}$		$\left[5.6\cdot {10}^6\right]$${}^{\text{e}}$							$G_{\text{fit}}$${}^{\text{f,g}}$
0–3	A1$\oplus$F2	A${}_{}{}^{C_{\text{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	A1$\oplus$F2	A${}_{}{}^{C_{\text{3v}}}{}_1{}^{}$	1	$\left[v_{\text{s}}\right]$	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$\oplus$F1$\oplus$F2	E${}^{C_{\text{3v}}}$	0	$\left[v_{\text{b}}\right]$	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	A1$\oplus$F2	A${}_{}{}^{C_{\text{3v}}}{}_1{}^{}$	2	$\left[2v_{\text{s}}\right]$	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$\oplus$F1$\oplus$F2	E${}^{C_{\text{3v}}}$	1	$\left[v_{\text{s}}+v_{\text{b}}\right]$	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	A1$\oplus$F2	A${}_{}{}^{C_{\text{3v}}}{}_1{}^{}$	0	$\left[2v_{\text{b}}\right]$	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$\oplus$F1$\oplus$F2	E${}^{C_{\text{3v}}}$	0	$\left[2v_{\text{b}}\right]$	523.5	(0.2)${}^{\text{h}}$	527.0	(0.2)${}^{\text{h}}$	526.0	(0.0)	545.7	(2.1)		555.3	(0.1)	525.6	(0.0)
40–43	A1$\oplus$F2	A${}_{}{}^{C_{\text{3v}}}{}_1{}^{}$	3	$\left[3v_{\text{s}}\right]$	555.0	(0.0)	556.4	(0.0)	557.5	(0.0)				519.2	(0.0)	556.8	(0.0)
44–51	E$\oplus$F1$\oplus$F2	E${}^{C_{\text{3v}}}$	2	$\left[2v_{\text{s}}+v_{\text{b}}\right]$	630.3	(0.0)	632.8	(0.0)	632.7	(0.0)						633.6	(0.0)
52–55	A1$\oplus$F2	A${}_{}{}^{C_{\text{3v}}}{}_1{}^{}$	1	$\left[v_{\text{s}}+2v_{\text{b}}\right]$	685.7	(0.0)	689.0	(0.0)	688.5	(0.0)						690.2	(0.0)
56–63	E$\oplus$F1$\oplus$F2	E${}^{C_{\text{3v}}}$	1	$\left[v_{\text{s}}+2v_{\text{b}}\right]$	702.6	(0.2)${}^{\text{h}}$	706.2	(0.2)${}^{\text{h}}$	705.5	(0.0)						705.9	(0.0)
64–67	A1$\oplus$F2	A${}_{}{}^{C_{\text{3v}}}{}_1{}^{}$	4	$\left[4v_{\text{s}}\right]$	724.2	(0.0)	726.6	(0.0)	728.5	(0.0)						722.9	(0.0)
rmsd${}^{\text{i}}$					2.4		0.8		0		11.7		20.3	23.5		1.9

Equations

$${\hat{T}}^{\mathrm{v}}$$

$$B_i$$

$$U$$

$$\sum _{i=1}^9{n}_{q_i}\le b$$

$${\Gamma }^{\text{CR}}(0)$$

$${\Gamma }^{\text{CR}}(1)$$

$${\Gamma }^{\text{CR}}(2)$$

$${\Gamma }^{\text{CR}}(3)$$

$${\Gamma }^{\text{CR}}(4)$$

$${\Gamma }^{\text{CR}}(5)$$

$${\Gamma }^{\text{CR}}(6)$$

$$\Gamma ({\text{A}}_1^{C_{3\text{v}}})$$

$$\Gamma ({\text{A}}_2^{C_{3\text{v}}})$$

$$\Gamma ({\text{E}}^{C_{3\text{v}}})$$

Full text

I Introduction

The CH${}_{4}\cdot$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 Wild et al. (2002); Loh et al. (2004, 2006); Neumark (2008); Yacovitch et al. (2012), and it has been used as a precursor in anion photoelectron spectroscopy to probe the transition state region of the $\text{F}+\text{CH}_{4}$ reaction Cheng et al. (2011); Westermann et al. (2014). This experimental activity motivated the computational (ro)vibrational quantum dynamics study of the complex Czakó et al. (2008); Wodraszka et al. (2012); Fábri et al. (2013). This complex also serves as a good prototype for molecular interactions with relatively large monomer distortions and strong binding.

CH${}_{4}\cdot\text{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. Czakó et al. (2008), Wodraszka et al. (2012), and Fábri et al. (2013), using the MULTIMODE Bowman et al. (2003), the MCTDH Meyer et al. (2009); Beck et al. (2000), and the GENIUSH Mátyus et al. (2009); Fábri et al. (2011) 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 Czakó et al. (2008), which we will refer to 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 $\text{cm}^{-1}$, was reported to be 42 $\text{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_{3\text{v}}$ 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_{\text{d}}$(M). The complex is bound by $D_{\text{e}}=2434$ $\text{cm}^{-1}$ on the CBB08 PES Czakó et al. (2008), which corresponds to $D_{0}=2316$ $\text{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}\approx 1050$ $\text{cm}^{-1}$ Quack and Suhm (1991), which is less than half of the binding energy of CH${}_{4}\cdot\text{F}^{-}$. The prototype for (double-well) tunneling, malonaldehyde has a barrier of 1410 $\text{cm}^{-1}$ Mizukami et al. (2014), which is ca. 30% higher than the barrier in CH${}_{4}\cdot$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}\cdot\text{F}^{-}$ is accompanied by a relatively large distortion of the methane fragment. For the interaction (int) region, Table 3 of Ref. Czakó et al. (2008) reports the equilibrium structure with an elongated C–H bond, $r^{\text{int}}_{\text{eq}}(\text{C--H}_{\text{b}})=1.112$ Å, for the H which binds (b) to the F-, while for the other three hydrogens, $r^{\text{int}}_{\text{eq}}(\text{C--H})=1.095$ Å. The corresponding distorted tetrahedral structure is characterized by the $\alpha($H–C–H${}_{\text{b}})=110.46^{\text{o}}$ angle. In the CH${}_{4}+$F- channel, the practically isolated (isol) methane molecule is a regular tetrahedron with a C–H equilibrium distance, $r^{\text{isol}}_{\text{eq}}(\text{C--H})=1.090$ Å.

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 full- and reduced-dimensionality 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 Czakó et al. (2008); Wodraszka et al. (2012); Fábri et al. (2013). 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 Avila and Mátyus (2019). The GENIUSH–Smolyak approach combines the non-product grid (Smolyak) method of Ref. Avila and T. Carrington (2009), which has been used for several high-dimensional, semi-rigid molecules Avila and T. Carrington (2011a, b) as well as molecules with a single large-amplitude motion Lauvergnat and Nauts (2014); Nauts and Lauvergnat (2018), and the numerical kinetic energy operator (numerical KEO) approach of Ref. Mátyus et al. (2009) implemented in the GENIUSH program Mátyus et al. (2009); Fábri et al. (2011), which includes by now dozens of vibrational-coordinate definitions for floppy systems Mátyus et al. (2009); Fábri et al. (2011, 2013); Fábri et al. (2014a, b); Sarka and Császár (2016); Sarka et al. (2016, 2017); Fábri et al. (2017); Simkó et al. (2019).

Concerning the coordinate definition for the fluoride-methane complex, we used the $(R,\cos\theta,\phi)$ spherical polar coordinates to describe the relative orientation of the methane fragment and the fluoride ion, and the nine normal coordinates, $q_{1},q_{2},\ldots,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. Avila and Mátyus (2019), which reads for the $(\xi_{1},\xi_{2},\xi_{3},\ldots,\xi_{D})$ general coordinates with the special $\xi_{2}=c$ choice

[TABLE]

where $\boldsymbol{g}\in{R}^{(D+3)\times(D+3)}$ is the mass-weighted metric tensor, $\boldsymbol{G}=\boldsymbol{g}^{-1}$, $\tilde{g}=\text{det}\boldsymbol{g}$, the extrapotential term,

[TABLE]

and the volume element is $\text{d}V=\text{d}\xi_{1}\text{d}c\ldots\text{d}\xi_{D}$. In the full-dimensional treatment of CH${}_{4}\cdot\text{F}^{-}$, $D=12$ and the coordinates are $\xi_{1}=R$, $\xi_{2}=c=\cos\theta$, $\xi_{3}=\phi$, $\xi_{3+i}=q_{i}$ ($i=1,2,\ldots 9$). We treat $c=\cos\theta$ 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 Avila and Mátyus (2019). 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) Schiffel and Manthe (2010) was used as it is explained in Sec. IV.E of Ref. Avila and Mátyus (2019). The reduced-dimensionality computations have been carried out with the original GENIUSH program Mátyus et al. (2009), 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 Light and Carrington Jr. (2007). 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. Avila and Mátyus (2019) 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 Whitehead and Handy (1975), i.e., we can discard high-energy basis functions. We used the simple

[TABLE]

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 Avila and T. Carrington (2009, 2011a) in vibrational computations using the Smolyak algorithm.

Concerning the 3-dimensional ion-molecule ‘intermolecular’ part, described by the $R$, $\cos\theta$, and $\phi$ coordinates, we retained the full direct-product basis and grid. For $R$, we used a Morse tridiagonal basis set constructed similarly to Ref. Avila and Mátyus (2019), but using the $D_{0}=1975.27$ cm*-1*, $\alpha=0.9$, and $\gamma=18$ parameter values which correspond to the 1D cut of the current PES (all other coordinates fixed at their equilibrium value). The $\cos\theta$ degree of freedom was described with sin-cot-DVR basis functions and quadrature points Schiffel and Manthe (2010), while we used Fourier functions for the $\phi$ angle.

We used large basis sets and grids for the intermolecular degrees of freedom, $(R,\cos\theta,\phi)$, 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 1 reports the computed vibrational energies in comparison with literature values. We estimate the vibrational excitation energies from the largest 12D computation, in column $\mathcal{B}_{C}$, to be converged within 1 $\text{cm}^{-1}$ and the (apparently very small) tunneling splittings are converged with an uncertainty lower than 0.05 $\text{cm}^{-1}$. Note that for $b=2$ and $3$ we used only 23 sin-cot-DVR basis functions, $\mathcal{B}_{A}$ and $\mathcal{B}_{B}$, respectively, which results in a 0.2–0.3 $\text{cm}^{-1}$ (artificial) splitting for some higher excited states, but this splitting is reduced to less than $0.05$ $\text{cm}^{-1}$ upon the increase of the basis set, which is reported in column $\mathcal{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 $\text{cm}^{-1}$ (‘$G_{\text{fit}}$’ in the Table 1). ‘$G_{\text{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_{\text{PES}}(\text{C--H})=2.143\,624$ bohr to define the PES cut (which is slightly different from the value used in Ref. Fábri et al. (2013) and reproduces slightly better the 12D stretching excitations of the complex), but we had to use a (much!) larger value in the KEO, $r_{\text{KEO}}(\text{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_{\text{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 two-parameter 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}\cdot$F- in Ref. Fábri et al. (2013).

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_{3\text{v}}$ 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_{\text{d}}(\text{M})$, for which the symmetry analysis of (the global minimum of) CH${}_{4}\cdot$Ar Ferenc and Mátyus (2019) is almost verbatim adapted.

In order to assign irrep labels to the CH${}_{4}\cdot$F- vibrational states computed in the present work, we analyzed the wave function of the 3D fitted model computations (‘3D $G_{\text{fit}}$’ column in Table 1) 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_{\text{d}}(\text{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,\ldots$) Sarka et al. (2017); Ferenc and Mátyus (2019), 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 $\Gamma^{\text{CR}}(j)$ representation spanned by the $[j,j]_{00}$ coupled-rotor functions in $T_{\text{d}}(\text{M})$ are Ferenc and Mátyus (2019):

[TABLE]

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_{3\text{v}}$ 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 Ferenc and Mátyus (2019)

[TABLE]

Note that we use the $C_{\text{3v}}$ superscript for the PG irreps, i.e., $\text{A}_{1}^{C_{\text{3v}}}$, $\text{A}_{2}^{C_{\text{3v}}}$, and $\text{E}^{C_{\text{3v}}}$, in order to distinguish them from the MS group irreps, which are labelled with $\text{A}_{1}$, $\text{A}_{2}$, E, $\text{F}_{1}$, and $\text{F}_{2}$. The result of this analysis for the computed vibrational wave functions is summarized in the ‘$\Gamma(\text{MS})$’ and ‘$\Gamma(\text{PG})$’ columns of Table 1, 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}\cdot$Ar complex, for which unambiguous assignment of the PG vibrations beyond the zero-point state of the global minimum was hardly possible Ferenc and Mátyus (2019).) 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 1) based on the nodal structure along the ion-molecule separation coordinate $R$, i.e., $n$th ‘stretching’ excitation, labelled with $v_{\text{s}}$, $2v_{\text{s}}$, etc. Excitations different from pure stretching excitation were termed ‘bending’ in this qualitative description, and were labelled with $v_{\text{b}}$, $2v_{\text{b}}$, etc., or combinations of stretching and bending, $v_{\text{s}}+v_{\text{b}}$, etc. Similar qualitative labels had been provided in the earlier studies Wodraszka et al. (2012); Fábri et al. (2013), 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 1 (column ‘12D, GENIUSH–Smolyak’), and the best results are listed in column ‘$\mathcal{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 $\text{cm}^{-1}$. It is 5.1 $\text{cm}^{-1}$ larger than the 12D MCTDH result, which the authors of Ref. Wodraszka et al. (2012) claim to be a variational upper bound to the exact ZPVE. The 12D MULTIMODE ZPVE Czakó et al. (2008) is 3.1 $\text{cm}^{-1}$ (8.2 $\text{cm}^{-1}$) larger than the lowest value obtained by us (by Ref. Wodraszka et al. (2012) using MCTDH). The totally symmetric ZPV state (A${}_{1}^{C_{\text{3v}}}$) is split by the relative rotation of the methane and the fluoride to an A1 and an F2 symmetry species, Eq. (5). The $\mathcal{B}_{C}$ 12D computation (numerically) reproduces the degeneracy of the F2 state within 0.001 $\text{cm}^{-1}$ and predicts a tunneling splitting (much) smaller than 0.05 $\text{cm}^{-1}$.

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

The small tunneling splittings obtained in the present work are in agreement with the earlier, 3D computations including the intermolecular, $(R,\cos\theta,\phi)$ coordinates as active vibrational degrees of freedom using the GENIUSH program Fábri et al. (2013). Although our results confirm the small splittings obtained in the 3D computation of Ref. Fábri et al. (2013), we observe larger deviations for the vibrational excitation energies (band origins) from the 3D results of Ref. Fábri et al. (2013). The first 7 vibrational excitations have a 23.5 $\text{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 $\text{cm}^{-1}$ rmsd, is observed for the first 6 vibrational excitations in comparison with the 12D MCTDH result Wodraszka et al. (2012), 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 Czakó et al. (2008) has a (surprisingly) large, 20.3 $\text{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 $\text{cm}^{-1}$, is observed for the 299.9 $\text{cm}^{-1}$ vibration, and this is the state which had the worse convergence properties in the MULTIMODE computation (Table 5 of Ref. Czakó et al. (2008)).

These observations suggest that although tunneling (and the corresponding splittings of the vibrational bands) is almost negligible up to 725 $\text{cm}^{-1}$ above ZPVE under an energy resolution of 0.05 $\text{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.

Footnotes to Table 1:

${}^{\text{a}}$ Symmetry assignment (irrep decomposition) in the $T_{\text{d}}$(M) molecular symmetry group of the complex.

${}^{\text{b}}$ Symmetry assignment (irrep) within the $C_{\text{3v}}$ point group of the equilibrium structure.

${}^{\text{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,\ldots$ labels the states of this 1D model in an increasing energy order.

${}^{\text{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 Czakó et al. (2008); Wodraszka et al. (2012); Fábri et al. (2013).

${}^{\text{e}}$ $(N_{c},b)$: short label for indicating the basis set size. The $N_{c}$ value gives the number of the sin-cot-DVR functions used for $\cos\theta$ and $b$ is the basis pruning parameter, $\sum_{i=1}^{9}n_{q_{i}}\leq b$. The $R$ and $\phi$ degrees of freedom are described by $N_{R}=8$ Morse tridiagonal and $N_{\phi}=39$ Fourier basis functions, respectively. The overall basis set size, $[N_{R}N_{c}N_{\phi}(b+9)!/b!/9!]$, is also shown.

${}^{\text{f}}$ In the vibrational computations we used atomic masses, $m(\text{H})=1.007\,825\,032\,23$ u, $m(\text{C})=12$ u, and $m(\text{F})=18.998\,403\,162\,73$ u Coursey et al. (2015).

${}^{\text{g}}$ Reduced-dimensionality model with active $(R,\cos\theta,\phi)$ degrees of freedom fitted to reproduce the 12D GENIUSH-Smolyak $\mathcal{B}_{C}$ result. The methane was treated as a regular tetrahedron and we used $r_{\text{PES}}(\text{C--H})=2.143\,624$ bohr in the PES (to reproduce well the stretches) and $r_{\text{KEO}}(\text{C--H})=2.518\,620$ bohr in the KEO, which was necessary to obtain good bending energies.

${}^{\text{h}}$ These splittings disappear upon increase of the intermolecular basis set.

${}^{\text{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 $\mathcal{B}_{C}$ result.

∗ Higher energy vibrations, including the 12 fundamental vibrations, are reported in Ref. Czakó et al. (2008), 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 ($<0.05$ $\text{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. Wodraszka et al. (2012) estimated the barrier height connecting the equivalent wells to be ca. 1270 $\text{cm}^{-1}$, and we have pinpointed a lower barrier height, $V_{\text{barrier}}=1104$ $\text{cm}^{-1}$ at $(\theta,\phi)=(90\ ^{\text{o}},45\ ^{\text{o}})$ (see Figure 2 of Ref. Fábri et al. (2013)), $R=5.745$ bohr, and $(q_{1},q_{2},\ldots,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 $\text{cm}^{-1}$ measured from the PES minima, which is ca. 700 $\text{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 $\text{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. Fábri et al. (2013), cited in Table 1, which corresponds to imposing rigorous geometrical constraints in the model through the Lagrangian in GENIUSH Mátyus et al. (2009). The fixed methane regular tetrahedral structure used in this 3D(rig) model was $r(\text{C--H})=1.104$ Å Fábri et al. (2013), a value close to the isolated methane’s $\langle r(\text{C--H})\rangle_{0}$ value. Second, we have carried out computations using the 3D$(G_{\text{fit}})$ model designed (‘fitted’) to minimize the rmsd of the vibrational band origins with respect to the 12D GENIUSH–Smolyak results up to 730 $\text{cm}^{-1}$ (Table 1). In the largest 3D computations, we used 31 PO-DVR functions Wei and Carrington, Jr. (1992); Echave and Clary (1992); Szalay et al. (2003) for $R$ obtained from $\mathcal{L}_{n}^{(\alpha)}$ generalized Laguerre basis functions (with $\alpha=2$ and $n=0,1,\ldots,399$) and grid points scaled to the $R\in[2,5]$ bohr interval; 61 sin-cot-DVR functions for the $\cos\theta$ degree of freedom, which was constructed by extending the $\cos(m\theta)$ ($m=0,...,58$) basis set with the $\sin(\theta)$ and the $\sin(2\theta)$ functions Schiffel and Manthe (2010); Avila and Mátyus (2019); and 121 Fourier (sine and cosine) functions for the $\phi$ degree of freedom. This basis is highly saturated and represents a space comparable with the set of the $Y_{j}^{m}(\theta,\phi)$ 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$ $\text{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 Gauss–Laguerre–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\theta$; and 45 Fourier basis functions and 48 trapezoidal points for $\phi$. So, this 12D basis set designed for the computation of the tunneling splittings (and henceforth labelled with $\mathcal{B}_{\text{T}}$) included $(2+9)!/2!/9!\times 8\times 31\times 45=613\,800$ functions and the 12D Smolyak grid ($H=16$, $D=12$ dimensions) contained $3\,481$(for the CH4 monomer)$\times 10\times 31\times 48=51\,797\,280$ points to calculate the integrals. The angular part of $\mathcal{B}_{T}$ is more than 40 % larger than that of $\mathcal{B}_{C}$ in Table 1, which provides the current benchmark for the vibrational band origins up to 730 $\text{cm}^{-1}$. Although the uncertainty of the vibrational band origins for $\mathcal{B}_{T}$ is 3–4 $\text{cm}^{-1}$ (due to the small $b$ value), the 12D tunneling splittings are converged with an uncertainty better than 0.02 $\text{cm}^{-1}$ up to 856 $\text{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 $\text{cm}^{-1}$ uncertainty threshold above ca. 680 $\text{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 $\text{cm}^{-1}$ from the PES minimum and ca. 700 $\text{cm}^{-1}$ measured from the intermolecular ZPVE (different models have different ZPVEs).

Figure 2 presents the 12D and the 3D splittings. The 3D$(G_{\text{fit}})$ model, which reproduces the known 12D vibrational band origins very well (Table 1), 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 1) 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_{\text{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.

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 Avila and Mátyus (2019) 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 $\text{cm}^{-1}$ and 0.05 $\text{cm}^{-1}$, respectively. These computations confirm complementary aspects of earlier studies Czakó et al. (2008); Wodraszka et al. (2012); Fábri et al. (2013), 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 $\text{cm}^{-1}$ above the zero-point energy) are small Fábri et al. (2013), $<0.05$ $\text{cm}^{-1}$, but due to the strong binding and significant monomer distortions, a full-dimensional (12D) treatment is necessary Czakó et al. (2008); Wodraszka et al. (2012). Although the tunneling splittings are small, a single-well (normal coordinate) description Czakó et al. (2008) is not sufficient due to the significant anharmonicity of the multi-well potential energy landscape Wodraszka et al. (2012). 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 Fábri et al. (2013).

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, $>0.1$ $\text{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ürich 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.