# Vibrational transitions of H 2 Cl +: Potential energy surface and anharmonic computations

Kokou Afansounoudji, Rabiou Issa, Komi Sodoga, David Lauvergnat

► To cite this version:

Kokou Afansounoudji, Rabiou Issa, Komi Sodoga, David Lauvergnat. Vibrational transitions of H 2 Cl +: Potential energy surface and anharmonic computations. Chemical Physics Letters, 2023, 833, pp.140914. 10.1016/j.cplett.2023.140914. hal-04267773

HAL Id: hal-04267773

https://cnrs.hal.science/hal-04267773v1

Submitted on 2 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Vibrational transitions of $H_2Cl^+$: potential energy surface and anharmonic computations

Kokou M. Robert Afansounoudji<fnref n="1" />*, Rabiou Issa<fnref n="1" />†, Komi Sodoga<fnref n="1" />‡and David Lauvergnat<fnref n="2" />§

## Abstract

A new three-dimensional potential energy surface of the electronic ground state of the chloronium ion, $H_2Cl^+$, based on the explicitly correlated coupled cluster method with a triple zeta basis set adapted to this method has been expanded in an analytical representation. This potential energy surface is later incorporated into our home-made Fortran code to compute variationally the vibrational levels, zero-point ground average structural parameters and the rotational constants of the chloronium ion and several isotopologues. Our results show a good agreement with experimental data and that our results will help to detect $H_2Cl^+$ isotopologues in the interstellar medium.

## 1 Introduction

Recently, the chloronium ion $H_2Cl^+$ was detected for the first time in the interstellar medium (ISM) by the Herschel Space Observatory using the Heterodyne Instrument for the far-Infrared (HIFI) [1, 2]. The dominant ionization state of each element in diffuse interstellar gas clouds is determined by its ionization potential. Chlorine has a slightly lower ionization potential than that of a hydrogen molecule, and is therefore mostly ionized individually in the form of $Cl^+$. The latter ion can react with the dominant constituent of the interstellar medium $H_2$ to form the ion $HCl^+$ [1] which, in turn, reacts with hydrogen to form the chloronium ion $H_2Cl^+$. Chemical modeling suggests that the chloronium ion could be very abundant around the diffuse ISM [3]. However, prior to the launch of Herschel, $H^{35}Cl$ and $H^{37}Cl$ were the only $Cl^-$bearing molecules detected in and around the ISM [4, 5, 6, 7]. This can be explained by the fact that the $H_2Cl^+$ ion transforms into hydrogen chloride, $HCl$ by dissociative recombination or by proton transfer to carbon monoxide CO. In addition to the detection by Herschel, another study reported by Neufeld et al [8], also detected absorption and emission of $H_2Cl^+$ from two interstellar cloud sources. The chemical reactions that give rise to the deuterated forms $HDCl^+$ and $D_2Cl^+$ are similar to those giving the chloronium ion, $H_2Cl^+$, with the same orders of magnitude of the reaction rate coefficient. It follows that, like the chloronium ion, its deuterated isotopologues should be abundant in the ISM. Surprisingly, to date, the existence of these deuterated isotopologues ($HDCl^+$ and $D_2Cl^+$) has not been observed in the ISM.

The availability of accurate spectroscopic data on chloronium and its deuterated forms could greatly facilitate their identification in interstellar space. Despite their importance, very few spectroscopic data on this ion and even less on its deuterated forms $HDCl^+$ and $D_2Cl^+$ are available in the literature. Experimentally, the first structural information on the chloronium ion was reported in 1986 by Kawaguchi and Hirota [9], who studied its $v_2$ band with infrared diode laser spectroscopy. Saito et al. [10] detected microwave spectral lines of $H_2^{35}Cl^+$ in 1988. In the same year, Lee et al [11] worked on the difference-frequency laser spectroscopy of the chloronium. These latter authors have observed the gas phase of the $v_1$ and $v_2$ fundamental bands of $H_2Cl^+$ and determined its equilibrium structure. The rotational spectra of $H_2Cl^+$ and $HDCl^+$ were measured by Araki et al [12] using a submillimeter-wave spectrometer. Furthermore, theoretically, Botschwina [13] computed a potential energy surface (PES) at coupled electron-pair approximation (CEPA) level and the vibrational frequencies for $H_2Cl^+$ and its isotopologues $HDCl^+$ and $D_2Cl^+$. More recently, after the detection of the chloronium ion, Majumdar et al [14] explored the possibility of finding the isotopologues of $H_2Cl^+$ in and around the ISM. These authors used Moller-Plesset perturbation theory (MP2) along with a quadruple-zeta

<footnote n="1">¹Laboratoire de Physique des Matériaux et des Composants à Semi-conducteurs, Département de Physique, Université de Lomé, Bd Eyadema, 01BP1515, Lomé, Togo</footnote>

<footnote n="2">²Université Paris-Saclay, CNRS, Institut de Chimie Physique, UMR8000, 91405 Orsay Cedex, France</footnote>

<footnote>*email: afankokourobert@gmail.com</footnote>

<footnote>†email: issagoudouya@gmail.com</footnote>

<footnote>‡email: ksodoga@univ-lome.tg</footnote>

<footnote>§Corresponding author: David.Lauvergnat@universite-paris-saclay.fr</footnote>
correlation consistent (aug-cc-pVQZ) basis set to optimize the geometry of chloronium and computed the harmonic vibrational frequencies of $H_2Cl^+$ and its two deuterated isotopologues $HDCl^+$ and $D_2Cl^+$. Within their chemical modeling, they found that $HDCl^+$ could be efficiently formed in the gas phase and that the emission line should be strong enough to be observed.

The lack of availability of laboratory data, in particular in the case of reactive species (transition molecules, radicals, ions) which remain more difficult to study in the laboratory than their stable counterparts, is the main factor currently limiting astronomical molecular detection. In order to facilitate the detection of the chloronium isotopologues $HDCl^+$ and $D_2Cl^+$, in this work, we provide information on the spectroscopic properties of these isotopologues. We use the explicitly correlated F12, coupled-cluster method, CCSD(T)-F12b, approach to optimize the geometric parameters of the ion and to generate the 3D-potential energy surface. This allows the computation of vibrational energy levels of chloronium $H_2Cl^+$, its two isotopologues $HDCl^+$ and $D_2Cl^+$ (with 35Cl and 37Cl chlorine isotopes), using the variational approach implemented in our home-made Fortran code. We compare our theoretical spectra with those available in the literature, experimental ones with Lee et al [11] and theoretical ones with Botschwina [13]. Our results show good agreements with those of the literature.

The present paper is organized as follows. In Section 2, the methods for the determination of the kinetic energy operator, the electronic structure, the potential energy surface and the computational details are presented. The results are discussed and concluding remarks are given in Section 3. In the Appendix, the fitting procedure of PES is given.

# 2 Methods

## 2.1 Coordinates and kinetic energy operator

The vibrational motion of N atom molecular system, moving on a potential energy surface, $V(\mathbf{q})$ can be described by $n = 3N - 6$ internal coordinates, $\mathbf{q} = [q^1, q^2, \dots, q^n]$. Considering the body-fixed (BF) Cartesian coordinates expressed as functions of the internal coordinates, $\mathbf{x} = \mathbf{x}(\mathbf{q})$, i.e $x^a = x^a(q^1, \dots, q^n)$, $\alpha = 1, \dots, 3N$, the standard expression of the deformation part of the kinetic energy operator (KEO) in terms of the contravariant components of the metric tensor $G^{ij}(\mathbf{q})$ and with the Euclidean volume element, $dv^E$, reads as [15]:

$$ \hat{T}_{def}^{E}(\mathbf{q}, \partial_{\mathbf{q}}) = -\frac{\hbar^2}{2} \sum_{i,j=1}^{n,n} J(\mathbf{q})^{-1} \partial_i J(\mathbf{q}) G^{ij}(\mathbf{q}) \partial_j, \quad \partial_i \equiv \partial/\partial q^i, \qquad (1) $$

where $J(\mathbf{q}) = \sqrt{\det(\mathbf{g})}$, $[G^{ij}(\mathbf{q})] = [g_{ij}(\mathbf{q})]^{-1}$. The covariant components of the metric tensor, $g_{ij}(\mathbf{q})$, are split in several contributions, $\mathbf{g} = \begin{bmatrix} \mathbf{S}_{n \times n} & \mathbf{C}_{n \times 3}^t \\ \mathbf{C}_{3 \times n} & \mathbf{I}_{3 \times 3} \end{bmatrix}$ where $\mathbf{S}_{n \times n}$, $\mathbf{I}_{3 \times 3}$ and $\mathbf{C}_{3 \times n}$ are the deformation part, the rotation part or the inertia tensor and the Coriolis part of the metric tensor, respectively. Furthermore, $dv_{def}^E = J(\mathbf{q}) dq^1 \cdot dq^2 \cdots dq^n$ is the deformation part of the Euclidean volume element.

By expanding the kinetic energy operator (1) and changing the volume element, it can be rewritten in the form:

$$ \hat{T}_{def}^{\rho}(\mathbf{q}, \partial_{\mathbf{q}}) = \sum_{i \le j} f_{2}^{ij}(\mathbf{q}) \partial_{ij}^{2} + \sum_{i} f_{1}^{i}(\mathbf{q}) \partial_{i} + V_{e}(\mathbf{q}), \quad \partial_{ij}^{2} \equiv \partial^{2}/\partial q^{i} \partial q^{j}, \qquad (2) $$

where $V_e(\mathbf{q})$ is the extrapolation term [16, 17, 18] resulting from the use of a non-Euclidean volume element $dv^\rho = \rho(\mathbf{q}) dq^1 dq^2 \cdots dq^n$. $\rho(\mathbf{q})$ is a real positive function which has to be adapted so that the basis set be orthonormal. The quantities $f_2^{ij}(\mathbf{q})$, $f_1^i(\mathbf{q})$ and $V_e(\mathbf{q})$ in eq. (2) are expressed as follows:

$$ f_2^{ii}(\mathbf{q}) = -\frac{\hbar^2}{2} G^{ii}(\mathbf{q}), \qquad (3) $$

$$ f_2^{ij}(\mathbf{q}) = -\hbar^2 G^{ij}(\mathbf{q}) \quad (i \neq j), \qquad (4) $$

$$ f_1^i(\mathbf{q}) = -\frac{\hbar^2}{2} \sum_j [G^{ij} \partial_j \ln \rho(\mathbf{q}) + \partial_j \ln G^{ij}(\mathbf{q})], \qquad (5) $$

$$ V_e(\mathbf{q}) = \left[ \sqrt{\frac{J(\mathbf{q})}{\rho(\mathbf{q})}} \hat{T}_{def}^{E}(\mathbf{q}, \partial_{\mathbf{q}}) \sqrt{\frac{\rho(\mathbf{q})}{J(\mathbf{q})}} \right]. \qquad (6) $$
The formalism of the derivation of the kinetic energy operator (KEO) Eq. (1) is well known [15, 19, 20, 21]. The first numerical implementation of the KEO operator dates back to the work of Laane et al. [22, 23]. Later, other groups proposed numerical approaches to determine the kinetic energy operator in the context of quantum dynamics [24, 25, 26]. M. Ndong and coworkers [27, 28] developed the Fortran program TANA which provides an automatic computation of the kinetic energy operator in analytical form for polyspherical coordinates for any number of atoms. In this approach, the expression of the metric tensor is obtained analytically and it can be used in Eq. (1) or in Eq. (2). Furthermore, these expressions can be exported in several formats (MCTDH,[29, 30] MidasCpp,[31] Fortran).

The 3D-vibrational motion of the chloronium ion can be described in terms of the three valence coordinates $\mathbf{q} = \{q^1, q^2, q^3\} \equiv \{R_1, R_2, \theta\}$, where $R_1$ and $R_2$ are the distances between the chlorine atom and the two hydrogen atoms and $\theta$ is the angle between the two $Cl-H$ bonds. The application of the general KEO formula Eqs. (2, 7, 8) with the Wilson normalization convention, $\rho(\mathbf{q}) = 1$, ($dv^1 = dR_1 \cdot dR_2 \cdot d\theta$) yields the following KEOs of the $H_2Cl^+$ and its isotopologues:

$$ \begin{aligned} \hat{T}(R_1, R_2, \theta) = & F_2^{11} \frac{\partial^2}{\partial R_1^2} + F_2^{22} \frac{\partial^2}{\partial R_2^2} + F_2^{33} \frac{\partial^2}{\partial \theta^2} + F_2^{12} \frac{\partial^2}{\partial R_1 \partial R_2} + F_2^{13} \frac{\partial^2}{\partial R_1 \partial \theta} + \\ & F_2^{23} \frac{\partial^2}{\partial R_2 \partial \theta} + F_1^1 \frac{\partial}{\partial R_1} + F_1^2 \frac{\partial}{\partial R_2} + F_1^3 \frac{\partial}{\partial \theta} + V_{ep}, \end{aligned} \quad (7) $$

where the functions, $F_2^{ij}$, $F_1^i$ and $V_{ep}$, are given below:

$$ \begin{aligned} F_2^{11} &= -\frac{\hbar^2}{2\mu_{R_1}}, & F_2^{22} &= -\frac{\hbar^2}{2\mu_{R_2}}, & F_2^{12} &= -\hbar^2 \frac{\cos(\theta)}{m_{Cl}}, \\ F_2^{33} &= -\frac{\hbar^2}{2} \left( \frac{1}{\mu_{R_1} \cdot R_1^2} + \frac{1}{\mu_{R_2} \cdot R_2^2} - \frac{2 \cos(\theta)}{m_{Cl} \cdot R_1 \cdot R_2} \right), \\ F_2^{13} &= \frac{\hbar^2}{2} \frac{1}{m_{Cl} \cdot R_2} \left( \sin(\theta) + \frac{1}{\sin(\theta)} - \frac{\cos^2(\theta)}{\sin(\theta)}, \right), \\ F_2^{23} &= \frac{\hbar^2}{2} \frac{1}{m_{Cl} \cdot R_1} \left( \sin(\theta) + \frac{1}{\sin(\theta)} - \frac{\cos^2(\theta)}{\sin(\theta)}, \right), \end{aligned} \quad (8) $$

$$ \begin{aligned} F_1^1 &= \frac{\hbar^2}{4} \frac{1}{m_{Cl} \cdot R_2} \left( 3 \cos(\theta) - \frac{\cos(\theta)}{\sin^2(\theta)} - \frac{\cos^3(\theta)}{\sin^2(\theta)}, \right), \\ F_1^2 &= \frac{\hbar^2}{4} \frac{1}{m_{Cl} \cdot R_1} \left( 3 \cos(\theta) - \frac{\cos(\theta)}{\sin^2(\theta)} - \frac{\cos^3(\theta)}{\sin^2(\theta)}, \right), \\ F_1^3 &= \frac{\hbar^2}{m_{Cl} \cdot R_1 \cdot R_2} \left( \frac{\cos^2(\theta)}{\sin(\theta)} - \frac{1}{\sin(\theta)}, \right), \\ V_e(\mathbf{q}) &= -\frac{\hbar^2}{4} \frac{1}{m_{Cl} \cdot R_1 \cdot R_2} \left( 2 \cos(\theta) + \frac{\cos^3(\theta)}{\sin^2(\theta)} - 2 \frac{\cos(\theta)}{\sin^2(\theta)} \right) - \\ & \frac{\hbar^2}{8} \frac{1}{\mu_{R_1} \cdot R_1^2} \left( 2 + \frac{\cos^2(\theta)}{\sin^2(\theta)} \right) - \frac{\hbar^2}{8} \frac{1}{\mu_{R_2} \cdot R_2^2} \left( 2 + \frac{\cos^2(\theta)}{\sin^2(\theta)} \right), \end{aligned} \quad (9) $$

where $\mu_{R_1}$ and $\mu_{R_2}$ are reduced masses: $\frac{1}{\mu_{R_1}} = \frac{1}{m_{Cl}} + \frac{1}{m_{H^1}}$ and $\frac{1}{\mu_{R_2}} = \frac{1}{m_{Cl}} + \frac{1}{m_{H^2}}$. Furthermore, $m_{Cl}$, $m_{H^1}$ and $m_{H^2}$ are, respectively, the isotopic masses of the chlorine, the first hydrogen and the second hydrogen atoms. It is worth noting, that the TANA code has been used to obtain the analytical expression and the Fortran subroutines of the KEO. Furthermore, some simplifications are possible in the previous expressions (Eq. (8-9)). For example, $(\sin(\theta) + \frac{1}{\sin(\theta)} - \frac{\cos^2(\theta)}{\sin(\theta)})$ can be simplified in $(2 \sin(\theta))$. However, we prefer to present the unsimplified expressions, which are used in our quantum dynamics code and come directly from the TANA code [27]. The present triatomic KEO has been published in the past by several group [32, 33, 34].
<table><caption>**Table 1:** Isotopic masses used in this study and obtained from the 2018 NIST table.</caption><thead><tr><th>Isotope</th><th>masses (g.mol<sup>-1</sup>)</th><th>masses (au)</th></tr></thead><tbody><tr><td><sup>35</sup>Cl</td><td>34.968852682</td><td>63744.31895</td></tr><tr><td><sup>37</sup>Cl</td><td>36.965902602</td><td>67384.71826</td></tr><tr><td>H</td><td>1.00782503223</td><td>1837.152648</td></tr><tr><td>D</td><td>2.01410177812</td><td>3671.482943</td></tr></tbody></table>

## 2.2 Electronic structure and potential energy surface

The electronic structure has been carried out with the explicitly correlated, F12, coupled cluster method, CCSD(T)-F12b [35], with the F12 adapted triple zeta basis set [36], cc-pVTZ-F12 (CCSD(T)-F12b/VTZ-F12) using molpro2010 code [37]. As noted in several studies, this level of computation provides a good compromise between quality and computational efficiency [38]. The geometry was optimized at CCSD(T)-F12b/VTZ-F12 and the optimized H-Cl distance and the HClH angle are 1.306 Å and 94.3°, respectively. These values are close to the Botschwina CEPA-1 values (1.302 Å and 94.2°) [13] and the experimental ones (1.304 Å and 94.24° for both $H_2^{35}Cl^+$ or $H_2^{37}Cl^+$ isotopologues) [11].

The $H_2Cl^+$ potential energy surface (PES), $V(Q)$, was obtained by fitting 145 ab initio points at CCSD(T)-F12b/VTZ-F12 level using symmetrized coordinate displacements ($\theta - \theta_{eq}$, $1/2(R_1 + R_2) - R_{eq}$ and $1/2(R_1 - R_2)$). Along these displacements, the grid ranges in bohr are [-0.90 : 1.03], [-0.41 : 0.65] and [0.00 : 0.48], respectively. The full grid covers energies up to about 10000 $cm^{-1}$ and it has been fitted as linear combinations of coordinate displacement products. The fitted potential reproduces almost perfectly the ab initio points. Indeed, the root-mean-square error and the maximum error are about 0.1 $cm^{-1}$ and 0.5 $cm^{-1}$, respectively. Furthermore, Fig. 1 shows the ab initio energy points (red crosses) and the fitted potential (solid blue line) along the three symmetrized coordinates ($\theta$, $R_+$ and $R_-$), which clearly shows the almost perfect agreement.

The grid generation and the fitting procedure are explained in detail in the Appendix. This PES has been implemented in the Quantum Model Library, freely available on github [39].

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Ab initio energy point (red cross) and fitted potential (full blue line) along the three symmetrized coordinates. The energy is in Hartree, $\theta$ in radian and $R_+$ and $R_-$ are in bohr.</figcaption></figure>

Fig. 2 shows the two iso-contours of the PES cuts. The left panel displays the contour as function of the two HCl distances with the angle fixed around its equilibrium value (1.64 radian), while the right panel displays the contour as function of the angle and one HCl distance with the other distance fixed around its equilibrium value (2.46 bohr).
<figure><img src="image_3.png" /><figcaption>**Figure 2:** 2D-contour plots of $H_2Cl^+$ potential energy surface at CCSD(T)F12b/cc-pVTZ-F12. Left panel: $R_1 - R_2$ contour. Right panel: $R_1 - \theta$ contour.</figcaption></figure>

## 2.3 Vibrational computational approach

The computation of the vibrational levels is relatively standard and it is related on early works on variational methods [32, 40, 41, 42]. For the present study, it has been performed with the help of our home-made pseudo-spectral Fortran code, *TI_Schrod* [43]. It allows us to solve the time-independent Schrödinger equation with a Hamiltonian expressed in terms of the three internal coordinates, **q** = [$R_1, R_2, \theta$]. The 3D-quantum vibrational bound states are expanded on a direct product of one-dimensional harmonic oscillator basis sets, and thus, $\rho(\mathbf{q})$ associated with the non-Euclidean volume element is equal to one ($\rho(\mathbf{q}) = 1$). Furthermore, the integration procedure is performed using the direct product of Gauss-Hermite quadrature grids.

More precisely, each eigenstate is developed as follows:

$$ \Psi(R_1, R_2, \theta) = \sum_{i_1, i_2, i_3} C_{i_1, i_2, i_3} HO_{i_1}(R_1) HO_{i_2}(R_2) HO_{i_3}(\theta), \quad (10) $$

where, for the coordinate, $q^j$, the $HO_i(q^j)$ are the normalized Harmonic oscillator basis functions, $HO_i(q^j) = N_i^{sc} \cdot H_{i-1}(x^j(q^j)) \cdot \exp[-\frac{1}{2}x^j(q^j)^2] \cdot x^j(q^j)$ is the scaled and shifted coordinate, $x^j(q^j) = sc^j \cdot (q^j - q_{ref}^j)$, $N_i^{scj}$ is a normalization factor and $H_{i-1}(x)$ are the Hermite polynomials.

<table><caption>**Table 2:** Basis set parameters (in atomic unit) for the isotopologues.</caption><thead><tr><th>q<sup>j</sup></th><th>1 (R<sub>1</sub>)</th><th>2 (R<sub>2</sub>)</th><th>3 (θ)</th></tr></thead><tbody><tr><td>q<sup>j</sup><sub>ref</sub></td><td>2.467</td><td>2.467</td><td>1.646</td></tr><tr><td>k<sup>j</sup></td><td>0.282</td><td>0.282</td><td>0.168</td></tr><tr><td rowspan="7">1/G<sup>ij</sup></td><td>H<sub>2</sub><sup>35</sup>Cl<sup>+</sup></td><td>1785.69</td><td>1785.69</td><td>5424.03</td></tr><tr><td>HD<sup>35</sup>Cl<sup>+</sup></td><td>1785.69</td><td>3471.53</td><td>7158.55</td></tr><tr><td>D<sub>2</sub><sup>35</sup>Cl<sup>+</sup></td><td>3471.53</td><td>3471.53</td><td>10523.94</td></tr><tr><td>H<sub>2</sub><sup>37</sup>Cl<sup>+</sup></td><td>1788.39</td><td>1788.39</td><td>5432.85</td></tr><tr><td>HD<sup>37</sup>Cl<sup>+</sup></td><td>1788.39</td><td>3481.78</td><td>7173.92</td></tr><tr><td>D<sub>2</sub><sup>37</sup>Cl<sup>+</sup></td><td>3481.78</td><td>3481.78</td><td>10557.19</td></tr><tr><td>H<sub>2</sub><sup>35</sup>Cl<sup>+</sup></td><td>4.74</td><td>4.74</td><td>5.49</td></tr><tr><td>HD<sup>35</sup>Cl<sup>+</sup></td><td>4.74</td><td>5.59</td><td>5.89</td></tr><tr><td rowspan="5">sc<sup>j</sup></td><td>D<sub>2</sub><sup>35</sup>Cl<sup>+</sup></td><td>5.59</td><td>5.59</td><td>6.48</td></tr><tr><td>H<sub>2</sub><sup>37</sup>Cl<sup>+</sup></td><td>4.74</td><td>4.74</td><td>5.49</td></tr><tr><td>HD<sup>37</sup>Cl<sup>+</sup></td><td>4.74</td><td>5.60</td><td>5.89</td></tr><tr><td>D<sub>2</sub><sup>37</sup>Cl<sup>+</sup></td><td>5.60</td><td>5.60</td><td>6.49</td></tr></tbody></table>
$q_{ref}^j$ is the value of $q^j$ at the equilibrium geometry and the expression of $sc^j$ is given by $\sqrt[4]{k^j/G^{jj}}$ where $k^j$ and $G^{jj}$ are, respectively, the curvature and the value of contravariant component of the metric tensor at the $q_{ref}$ geometry [44]. The values of the $q_{ref}^j$, $sc^j$ and $G^{jj}$ are given in Table [2] for all studied isotopologues.

Before discussing in more detail the vibrational spectra of $H_2Cl^+$ and its isotopologues, we have checked the convergence of the energy levels with respect to the basis set and the grid sizes. In our tests, the three 1D-basis sets have the same number of basis functions, $nb$ which increases from 9 to 12 and the grid point number, $nq$, is fixed to 12. These tests have shown, that, the largest error between the energy levels up to 4000-5000 cm$^{-1}$ calculated with $nb = 11$ and $nb = 12$ is smaller than $10^{-2}$ cm$^{-1}$. Therefore, we have chosen 12 basis functions for each Harmonic oscillator basis and for all isotopologues.

## 3 Discussion and conclusion

The table 3 lists the computed and experimental anharmonic vibrational frequencies, expressed in cm$^{-1}$, for the $H_2^{35}Cl^+$ and $H_2^{37}Cl^+$ isotopologues. The differences between the experimental and the computed values are 0.5 cm$^{-1}$ for the bending mode ($v_2$) and about 5 cm$^{-1}$ for the two stretching modes ($v_1$ and $v_3$). This shows a good agreement between our calculated values and the experimental ones. Furthermore, these values show a clear improvement with respect to those computed by Botschwina with his unscaled potential [13], in particular, for the bending mode. Indeed, Botschwina developed two potentials: (i) the first one (unscaled) is obtained after the fitting of the *ab-initio* calculations; (ii) the second one is derived from the first one, but some parameters are scaled to reproduce the experimental fundamental transitions. Therefore, we compare our full *ab-initio* potential with the first Botschwina potential, and our results clearly show the quality of our potential. Table 4 gives the vibrational frequencies for the other isotopologues of this ion, namely $HD^{35}Cl^+$, $HD^{37}Cl^+$, $D_2^{35}Cl^+$, $D_2^{37}Cl^+$.

<table><caption>**Table 3:** computed and experimental anharmonic vibrational transitions (in cm$^{-1}$) for $H_2^{35}Cl^+$ and $H_2^{37}Cl^+$ isotopologues.</caption><thead><tr><th rowspan="2">Band</th><th colspan="3">H<sub>2</sub><sup>35</sup>Cl<sup>+</sup></th><th colspan="2">H<sub>2</sub><sup>37</sup>Cl<sup>+</sup></th></tr><tr><th>This work</th><th>Botschwina[13]</th><th>Exp.</th><th>This work</th><th>Exp.</th></tr></thead><tbody><tr><td><i>v</i><sub>2</sub></td><td>1184.6</td><td>1206.8</td><td>1184.1<sup>a</sup></td><td>1183.7</td><td>1183.2<sup>a</sup></td></tr><tr><td>2<i>v</i><sub>2</sub></td><td>2356.6</td><td>2398.1</td><td></td><td>2354.8</td><td></td></tr><tr><td><i>v</i><sub>3</sub></td><td>2634.9</td><td>2637.4</td><td>2630.1<sup>b</sup></td><td>2632.8</td><td>2628.1<sup>b</sup></td></tr><tr><td><i>v</i><sub>1</sub></td><td>2648.2</td><td>2650.0</td><td>2643.2<sup>b</sup></td><td>2646.5</td><td>2641.5<sup>b</sup></td></tr><tr><td>3<i>v</i><sub>2</sub></td><td>3515.3</td><td>3565.7</td><td></td><td>3512.6</td><td></td></tr><tr><td><i>v</i><sub>2</sub> + <i>v</i><sub>3</sub></td><td>3794.0</td><td>3816.5</td><td></td><td>3791.1</td><td></td></tr><tr><td><i>v</i><sub>1</sub> + <i>v</i><sub>2</sub></td><td>3813.1</td><td>3838.3</td><td></td><td>3810.5</td><td></td></tr><tr><td>4<i>v</i><sub>2</sub></td><td>4660.2</td><td>4690.5</td><td></td><td>4656.6</td><td></td></tr><tr><td>2<i>v</i><sub>2</sub> + <i>v</i><sub>3</sub></td><td>4940.7</td><td>4979.8</td><td></td><td>4937.0</td><td></td></tr><tr><td>2<i>v</i><sub>2</sub> + <i>v</i><sub>1</sub></td><td>4965.5</td><td>5011.1</td><td></td><td>4962.1</td><td></td></tr><tr><td>ZPE</td><td>3320.9</td><td>3335.1</td><td></td><td>3318.4</td><td></td></tr></tbody></table>

a: Ref [9]

b: Ref [11]
<table><caption>**Table 4:** Computed anharmonic vibrational frequencies (in cm⁻¹) for HD³⁵Cl⁺, HD³⁷Cl⁺, D₂³⁵Cl⁺ and D₂³⁷Cl⁺ isotopologues.</caption><thead><tr><th>Band</th><th>HD<sup>35</sup>Cl<sup>+</sup></th><th>HD<sup>37</sup>Cl<sup>+</sup></th><th>D<sub>2</sub><sup>35</sup>Cl<sup>+</sup></th><th>D<sub>2</sub><sup>37</sup>Cl<sup>+</sup></th></tr></thead><tbody><tr><td>ν<sub>2</sub></td><td>1034.6</td><td>1033.5</td><td>856.9</td><td>855.7</td></tr><tr><td>2ν<sub>2</sub></td><td>2059.4</td><td>2057.4</td><td>1707.4</td><td>1705.0</td></tr><tr><td>ν<sub>3</sub></td><td>2641.3</td><td>2639.4</td><td>1915.3</td><td>1912.4</td></tr><tr><td>ν<sub>1</sub></td><td>1917.5</td><td>1914.8</td><td>1920.0</td><td>1917.5</td></tr><tr><td>3ν<sub>2</sub></td><td>3074.2</td><td>3071.2</td><td>2551.4</td><td>2547.8</td></tr><tr><td>ν<sub>2</sub> + ν<sub>3</sub></td><td>3653.5</td><td>3650.6</td><td>2759.0</td><td>2754.9</td></tr><tr><td>ν<sub>1</sub> + ν<sub>2</sub></td><td>2939.9</td><td>2936.2</td><td>2766.7</td><td>2763.0</td></tr><tr><td>4ν<sub>2</sub></td><td>4078.5</td><td>4074.4</td><td>3388.5</td><td>3383.8</td></tr><tr><td>2ν<sub>2</sub> + ν<sub>3</sub></td><td>4656.2</td><td>4652.4</td><td>3596.5</td><td>3591.3</td></tr><tr><td>2ν<sub>2</sub> + ν<sub>1</sub></td><td>3952.7</td><td>3948.0</td><td>3607.0</td><td>3602.2</td></tr><tr><td>ZPE</td><td>2862.9</td><td>2860.0</td><td>2391.0</td><td>2387.6</td></tr></tbody></table>

The equilibrium molecular structure was therefore determined by taking the expectation values of each coordinate (⟨R₁⟩, ⟨R₂⟩ and ⟨Θ⟩) for the vibrational ground-state. The values resulting from these calculations are reported in the table 5. These values obtained for the isotopologues, H₂³⁵Cl⁺, HD³⁵Cl⁺ and H₂³⁷Cl⁺ agree very well with the experimental values (the distances and the angle are extracted from the experimental rotational constants) [12, 9]. Indeed, for the three isotopologues, the angle values are very close to the experimental ones and the errors are less than 0.1° [12, 9]. With respect to the Cl-H distance of H₂³⁵Cl⁺ obtained from ref [9], the difference with our calculation is about 0.1 Å, which is surprisingly large. However, we are confident with our values. Indeed, our equilibrium values agree very well with Botschwina ones[13] and also from the experimental ones [12, 9]. Furthermore, with respect to a more recent study of Araki et al. [12], the experimental Cl-H and Cl-D distances (associated to the zero point average structure[12]) are in good agreement with our calculated values since the largest difference with respect to the experimental value is only 0.005 Å with respect to ref [12]. Furthermore, our calculations nicely reproduce the shortening of the Cl-D bond length with respect to the Cl-H one for the HD³⁵Cl⁺ ion [12]. In the Table 5, we also give the expectation values for other isotopologues, D₂³⁵Cl⁺, HD³⁷Cl⁺ and D₂³⁷Cl⁺.

<table><caption>**Table 5:** Experimental and computed expectation values of coordinates and rotationnal constants for several isotopologues with the distances in Å, the angles in ° and the rotationnal constants in GHz.</caption><thead><tr><th rowspan="2"></th><th colspan="2">H₂³⁵Cl⁺</th><th colspan="2">HD³⁵Cl⁺</th><th colspan="2">H₂³⁷Cl⁺</th></tr><tr><th>This work</th><th>Ref [12]</th><th>This work</th><th>Ref[12]</th><th>This work</th><th>Ref[12]</th></tr></thead><tbody><tr><td>⟨R₁⟩</td><td>1.325</td><td>1.3206</td><td>1.325</td><td>1.3207</td><td>1.325</td><td>1.3206</td></tr><tr><td>⟨R₂⟩</td><td>1.325</td><td>1.3206</td><td>1.319</td><td>1.3165</td><td>1.325</td><td>1.3206</td></tr><tr><td>⟨Θ⟩</td><td>94.33</td><td>94.24</td><td>94.32</td><td></td><td>94.33</td><td>94.24</td></tr><tr><td>A</td><td>336.9</td><td>337.352</td><td>308.8</td><td>309.476</td><td>336.0</td><td>336.360</td></tr><tr><td>B</td><td>272.5</td><td>273.587</td><td>152.7</td><td>153.168</td><td>272.5</td><td>273.587</td></tr><tr><td>C</td><td>147.9</td><td>148.100</td><td>100.7</td><td>100.837</td><td>147.7</td><td>147.907</td></tr><tr><td></td><th colspan="2">HD³⁷Cl⁺</th><th colspan="2">D₂³⁵Cl⁺</th><th colspan="2">D₂³⁷Cl⁺</th></tr><tr><td></td><th>This work</th><th></th><th>This work</th><th></th><th>This work</th><th></th></tr><tr><td>⟨R₁⟩</td><td>1.325</td><td></td><td>1.320</td><td></td><td>1.319</td><td></td></tr><tr><td>⟨R₂⟩</td><td>1.319</td><td></td><td>1.320</td><td></td><td>1.319</td><td></td></tr><tr><td>⟨Θ⟩</td><td>94.32</td><td></td><td>94.31</td><td></td><td>94.31</td><td></td></tr><tr><td>A</td><td>308.2</td><td></td><td>177.6</td><td></td><td>176.7</td><td></td></tr><tr><td>B</td><td>152.4</td><td></td><td>136.5</td><td></td><td>136.5</td><td></td></tr><tr><td>C</td><td>100.5</td><td></td><td>76.2</td><td></td><td>76.0</td><td></td></tr></tbody></table>
In addition, we were able to compute the three rotational constants associated with the vibrational ground state. However, this has to be done with care. Indeed, since we did not perform a full ro-vibrational calculation, the rotational Eckart conditions[45] must be applied to minimize the Coriolis coupling. The Cartesian equilibrium geometry (obtained from $R_1 = R_2 = 2.46729$ bohr and $\theta = 1.64588$ rad) is used as the reference geometry to apply the Eckart conditions. They have been implemented in Tnum[28] as described in ref [46]. The rotational constants have been computed with the Fortran codes ElVibRot[47] coupled with Tnum using the following procedure: (i) Calculate the average values ($\langle \mu_{\alpha,\beta} \rangle$, with $\alpha$ and $\beta$ being x, y or z) of the contravariant rotational components of the metric tensor. This is done with the numerical conditions (basis set and grid parameters) used to compute the energy levels. (ii) Diagonalize the previous 3 × 3 matrix to get the rotational constants as its eigenvalues. The comparison between the computed and the experimental values given in Table 5 shows a nice agreement. Indeed, the largest error is about 1. GHz. Furthermore, the rotational constant shifts due to the isotopic substitution (from $H_2^{35}Cl^+$ to $H_2^{37}Cl^+$) are well reproduced.

In this computational work, we have generated the 3D-potential energy surface of the ground state of the chloronium ion using a state-of-the-art electronic structure level, namely the F12-explicitly correlated coupled cluster method with a F12-triple zeta atomic basis set. This fitted polynomial expansion was incorporated into a home-made variational Fortran code to compute the energy levels, the zero-point average coordinate values and the rotational constants of several isotopologues, $H_2^{35}Cl^+$, $H_2^{37}Cl^+$, $HD^{35}Cl^+$, $HD^{37}Cl^+$, $D_2^{35}Cl^+$ and $D_2^{37}Cl^+$. Our results are in good agreement with the available experimental data.

## 4 Acknowledgement

DL gratefully acknowledges the hospitality of the Laboratoire de Physique des Matériaux et des Composants à Semi-conducteurs, Département de Physique, Université de Lomé, where part of this work was carried out. RA gratefully thanks the ASESMANET network for giving me the opportunity to spend a month in the Institut de Chimie Physique of the Université Paris-Saclay.

## A Potential energy surface fitting procedure

The potential, $V(Q)$, has been expanded in terms of displacements, $\Delta Q_i = Q_i - Q_i^{ref}$, ($i = 1,2,3$)) of symmetrized coordinates, $Q$ ($[\theta, R_+, R_-]$) with $Q^{ref} = [1.645916$ rad, $2.467341$ bohr, $0. bohr]$ and expressed as a sum of four n-mode contributions, $V_0(Q)$, $\Delta V_1(Q)$, $\Delta V_2(Q)$ and $\Delta V_3(Q)$, which depend on the coupling between modes. The first contribution, $V_0(Q)$, is just the reference energy, i.e. the energy of the minimum. The other contributions are fitted from the energies generated similarly to the n-mode expansion:

$$ \Delta V_K(Q) = \sum_{\mathcal{R}_K^V(i_1,i_2,i_3)} C_{i_1,i_2,i_3} \cdot \Delta Q_1^{i_1} \Delta Q_2^{i_2} \Delta Q_3^{i_3}. \quad (11) $$

The $C_{i_1,i_2,i_3}$ are the coefficients obtained by a least square procedure of the multimode expansion ab initio values, $\Delta V_K^{ab}(Q)$. Each $\mathcal{R}_K^V(i_1,i_2,i_3)$ defines a constraint on the $[i_1,i_2,i_3]$ and its general expression is $(i_1 + i_2 + i_3) \le L_K^V$ and with $i_1 \ge 0, i_2 \ge 0$ and $i_3 \ge 0$. Furthermore, Since the potential is symmetric with respect to $\Delta Q_3$, only the even values of $i_3$ are possible. $L_K^V$ is a parameter which controls the number of terms or functions in the expansion. More precisely, for each value of K, $\mathcal{R}_K$ is defined as follows:

(i) One mode expansion (K=1), $\mathcal{R}_1^V$:

The possible values of the indices are $[i_1,0,0]$ or $[0,i_2,0]$ or $[0,0,i_3]$ (with $i_1 \ge 1, i_2 \ge 1$ and $i_3 \ge 2$) and $L_1^V = 9$, the number of functions is 22 (9 functions for $[i_1,0,0]$ and $[0,i_2,0]$ and only 4 functions for $[0,0,i_3]$ since $i_3$ is even.

(ii) Two mode expansion (K=2), $\mathcal{R}_2^V$:

The possible values of the indices are $[i_1,i_2,0]$ or $[i_1,0,i_3]$ or $[0,i_2,i_3]$ (with $i_1 \ge 1, i_2 \ge 1$ and $i_3 \ge 2$). With $L_2^V = 8$, the number of functions is 52 (28 functions for $[i_1,i_2,0]$ and 12 functions for $[i_1,0,i_3]$ or $[0,i_2,i_3]$.

(iii) Three mode expansion (K=3), $\mathcal{R}_3^V$:

Here, all indices are larger than 0 and with $L_3^V = 6$, the number of functions is 7.

Those three contributions are fitted independently of each other using well-defined grids generated in terms of displacement with respect to a reference geometry, $Q^{ref}$. A least-square procedure is used with the
$\Delta V_K^{ab}(\mathbf{Q})$ values on the grid points. Those values are related to the *ab initio* energy values, $V^{ab}(\mathbf{Q})$. The full 3D-grid is split in four parts, $G_0$, $G_1$, $G_2$ and $G_3$ depending on the coupling between modes. The first grid, $G_0$, contains a single point, the reference geometry and the corresponding ab initio energy, $E_0$. The other points are generated with constraints, $\mathcal{R}_K^G$, on the indices $[u_1, u_2, u_3]$; $u_1 + u_2 + u_3 \le L_K^G$, where $L_K^G$ control the size of the grid. The constraints for $K = 1, 2, 3$ are the following:

(i) One mode expansion (K=1), $\mathcal{R}_1^G$:

The possible values of the indices are $[u_1, 0, 0]$ or $[0, u_2, 0]$ or $[0, 0, u_3]$ (with $u_1 \ge 1$, $u_2 \ge 1$ and $u_3 \ge 1$). With $L_1^G = 5$, the number of grid points is 25 (10 for $[u_1, 0, 0]$ or $[0, u_2, 0]$ and only 5 for $[0, 0, u_3]$). Then, the $\Delta V_1^{ab}(\mathbf{Q})$ values are obtained as follows: $\Delta V_1^{ab}(\mathbf{Q}) = V^{ab}(\mathbf{Q}) - E_0$.

(ii) Two mode expansion (K=2), $\mathcal{R}_2^G$:

The possible values of the indices are $[u_1, u_2, 0]$ or $[u_1, 0, u_3]$ or $[0, u_2, u_3]$ (with $u_1 \ge 1$, $u_2 \ge 1$ and $u_3 \ge 2$). With $L_2^G = 5$, the number of grid points is 80 (40 for $[u_1, y_2, 0]$ and 20 $[u_1, 0, u_3]$ or $[0, u_2, u_3]$). The $\Delta V_2^{ab}(\mathbf{Q})$ is given by:

$$ \begin{aligned} \Delta V_2^{ab}(\mathbf{Q}) = V^{ab}(\mathbf{Q}) - E_0 & - \Delta V_1^{ab}(Q_1, Q_2^{\text{ref}}, Q_3^{\text{ref}}) \\ & - \Delta V_1^{ab}(Q_1^{\text{ref}}, Q_2, Q_3^{\text{ref}}) - \Delta V_1^{ab}(Q_1^{\text{ref}}, Q_2^{\text{ref}}, Q_3). \end{aligned} \quad (12) $$

(iii) Three mode expansion (K=3), $\mathcal{R}_3^G$:

Here, all index absolute values are larger than 0 and with $L_3^G = 5$, the number of grid points is 40. The $\Delta V_3^{ab}(\mathbf{Q})$ is given by:

$$ \begin{aligned} \Delta V_3^{ab}(\mathbf{Q}) = V^{ab}(\mathbf{Q}) - E_0 & - \Delta V_1^{ab}(Q_1, Q_2^{\text{ref}}, Q_3^{\text{ref}}) \\ & - \Delta V_1^{ab}(Q_1^{\text{ref}}, Q_2, Q_3^{\text{ref}}) - \Delta V_1^{ab}(Q_1^{\text{ref}}, Q_2^{\text{ref}}, Q_3) \\ & - \Delta V_2^{ab}(Q_1, Q_2, Q_3^{\text{ref}}) - \Delta V_2^{ab}(Q_1, Q_2^{\text{ref}}, Q_3) \\ & - \Delta V_2^{ab}(Q_1^{\text{ref}}, Q_2, Q_3) \end{aligned} \quad (13) $$

In those three grids ($G_K$), for each point the $k^{th}$-coordinate reads:

$$ Q_k = Q_k^{\text{ref}} + \Delta Q_k(u_k) \quad (14) $$

where the $\Delta Q_k(u_k)$ ($k = 1, 2, 3$) are given in table 6.

<table><caption>**Table 6:** Values of the $\Delta Q_k(u_k)$ as function of $u_k$ ($k = 1, 2, 3$) to define the grid points.</caption><thead><tr><th>u<sub>k</sub></th><th>0</th><th>1</th><th>-1</th><th>2</th><th>-2</th><th>3</th><th>-3</th><th>4</th><th>-4</th><th>5</th><th>-5</th></tr></thead><tbody><tr><th>ΔQ<sub>1</sub></th><td>0.00</td><td>0.21</td><td>-0.18</td><td>0.41</td><td>-0.36</td><td>0.62</td><td>-0.54</td><td>0.83</td><td>-0.72</td><td>1.03</td><td>-0.90</td></tr><tr><th>ΔQ<sub>2</sub></th><td>0.00</td><td>0.13</td><td>-0.08</td><td>0.26</td><td>-0.17</td><td>0.39</td><td>-0.25</td><td>0.52</td><td>-0.33</td><td>0.65</td><td>-0.41</td></tr><tr><th>ΔQ<sub>3</sub></th><td>0.00</td><td>0.10</td><td></td><td>0.19</td><td></td><td>0.29</td><td></td><td>0.38</td><td></td><td>0.48</td><td></td></tr></tbody></table>

## References

[1] D. C. Lis, J. C. Pearson, D. A. Neufeld, P. Schilke, H. S. P. Müller, H. Gupta, T. A. Bell, C. Comito, T. G. Phillips, E. A. Bergin, et al., Herschel/HIFI discovery of interstellar chloronium $H_2Cl^+$, Astronomy & Astrophysics 521 (2010) L9.

[2] G. L. Pilbratt, J. R. Riedinger, T. Passvogel, et al., Herschel Space Observatory-An ESA facility for far-infrared and submillimetre astronomy, Astronomy & Astrophysics 518 (2010) L1.

[3] D. A. Neufeld, M. G. Wolfire, The chemistry of interstellar molecules containing the halogen elements, The Astrophysical Journal 706 (2) (2009) 1594.

[4] G. A. Blake, J. Keene, T. G. Phillips, Chlorine in dense interstellar clouds-the abundance of HCl in OMC-1, Astrophysical Journal 295 (2) (1985) 501-506.
[5] J. Zmuidzinas, G. A. Blake, J. Carlstrom, J. Keene, D. Miller, HCl absorption toward sagittarius B2, The Astrophysical Journal 447 (2) (1995) L125.

[6] P. Schilke, T. Phillips, N. Wang, Hydrogen chloride in OMC-1, The Astrophysical Journal 441 (1995) 334–342.

[7] M. Salez, M. A. Frerking, W. D. Langer, Hydrogen Chloride Detection in Orion A and Monoceros R2 and Derivation of the H³⁵Cl/H³⁷Cl Isotopic Ratio, The Astrophysical Journal 467 (1996) 708.

[8] D. A. Neufeld, E. Roueff, R. L. Snell, et al., Herschel observations of interstellar chloronium, The Astrophysical Journal 748 (1) (2012) 37.

[9] K. Kawaguchi, E. Hirota, Magnetic field modulated infrared laser spectroscopy of the chloronium ClH₂⁺ ion ν₂ band, The Journal of chemical physics 85 (12) (1986) 6910–6913.

[10] S. Saito, S. Yamamoto, K. Kawaguchi, The microwave spectrum of the H₂Cl⁺ ion, The Journal of Chemical Physics 88 (4) (1988) 2281–2283.

[11] S. K. Lee, T. Amano, K. Kawaguchi, M. Oldani, Difference-frequency laser spectroscopy of the ν₁ and ν₃ fundamental bands of H₂Cl⁺: Determination of the equilibrium molecular structure, Journal of Molecular Spectroscopy 130 (1) (1988) 1–19.

[12] M. Araki, T. Furuya, S. Saito, Submillimeter-wave spectra of H₂Cl⁺ and its isotopic species: Molecular structure, Journal of Molecular Spectroscopy 210 (1) (2001) 132–136.

[13] P. Botschwina, Anharmonic potential-energy surfaces, vibrational frequencies and infrared intensities calculated from highly correlated wavefunctions, Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 84 (9) (1988) 1263–1276.

[14] L. Majumdar, A. Das, S. K. Chakrabarti, Formation of different isotopomers of chloronium in the interstellar medium, The Astrophysical Journal 782 (2) (2014) 73.

[15] A. Nauts, X. Chapuisat, Momentum, quasi-momentum and hamiltonian operators in terms of arbitrary curvilinear coordinates, with special emphasis on molecular hamiltonians, Mol. Phys. 55 (6) (1985) 1287–1318.
doi:10.1080/00268978500102031.

[16] D. Lauvergnat, A. Nauts, Exact numerical computation of a kinetic energy operator in curvilinear coordinates, J. Chem. Phys 116 (2002) 8560.
doi:10.1080/00268978500102031.

[17] E. Mátyus, G. Czakó, A. G. Császár, Toward black-box-type full-and reduced-dimensional variational (ro) vibrational computations, The J. of Chem. Phys 130 (13) (2009) 134112.

[18] X. Chapuisat, A. Belafhal, A. Nauts, N-body quantum-mechanical hamiltonians: Extrapotential terms, J. of Mol. Spec. 149 (1) (1991) 274–304.
doi:https://doi.org/10.1016/0022-2852(91)90159-8.

[19] R. Meyer, H. H. Günthard, General internal motion of molecules, classical and quantum-mechanical hamiltonian, The J. Chem. Phys. 49 (4) (1968) 1510–1520.

[20] E. B. Wilson, J. C. Decius, P. C. Cross, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra, McGraw-Hil Edition, no. 9, New York, 1955.
doi:10.1149/1.2430134.

[21] G. Brocks, A. Van Der Avoird, B. Sutcliffe, J. Tennyson, Quantum dynamics of non-rigid systems comprising two polyatomic fragments, Molecular Physics 50 (5) (1983) 1025–1043.

[22] J. Laane, M. A. Harthcock, P. Killough, L. Bauman, J. Cooke, Vector representation of large-amplitude vibrations for the determination of kinetic energy functions, J. of Mol. Spec. 91 (2) (1982) 286–299.

[23] M. A. Harthcock, J. Laane, Calculation of kinetic energy terms for the vibrational hamiltonian: Application to large-amplitude vibrations using one-, two-, and three-dimensional models, J. of Mol. Spec. 91 (2) (1982) 300–324.

[24] Y. G. Smeyers, F. J. Melendez, M. L. Senent, Ab initio theoretical study of the methyl and phosphine torsion modes in ethylphosphine, J. Chem. Phy 106 (5) (1997) 1709–1717.
[25] M. L. Senent, Determination of the kinetic energy parameters of non-rigid molecules, Chem. phys. lett. 296 (3-4) (1998) 299-306.

[26] D. J. Rush, K. B. Wiberg, Ab initio CBS-QCI calculations of the inversion mode of ammonia, J. Phys. Chem. A 101 (17) (1997) 3143-3151.

[27] M. Ndong, A. Nauts, L. Joubert-Doriol, H.-D. Meyer, F. Gatti, D. Lauvergnat, Automatic computer procedure for generating exact and analytical kinetic energy operators based on the polyspherical approach: General formulation and removal of singularities, J. Chem. Phys. 139 (20) (2013) 034107.

[28] D. Lauvergnat, Tnum-Tana: a Fortran code to deal with curvilinear coordinates and kinetic energy operator (2002).
URL https://github.com/lauvergn/Tnum-Tana

[29] F. Richter, F. Thaunay, D. Lauvergnat, P. Carbonnière, Anharmonic Vibrational Treatment Exclusively in Curvilinear Valence Coordinates: The Case of Formamide, The Journal of Physical Chemistry A 119 (48) (2015) 11719-11728.
doi:10.1021/acs.jpca.5b08482
URL http://pubs.acs.org/doi/10.1021/acs.jpca.5b08482

[30] M. Schröder, F. Gatti, D. Lauvergnat, H.-D. Meyer, O. Vendrell, The coupling of the hydrated proton to its first solvation shell, Nature Communications 13 (1) (2022) 6170. arXiv:2204.03744, doi:10.1038/s41467-022-33650-w
URL https://www.nature.com/articles/s41467-022-33650-w

[31] E. L. Klinting, D. Lauvergnat, O. Christiansen, Vibrational Coupled Cluster Computations in Polyspherical Coordinates with the Exact Analytical Kinetic Energy Operator, Journal of Chemical Theory and Computation (2020) acs.jctc.0c00261doi:10.1021/acs.jctc.0c00261.
URL https://pubs.acs.org/doi/10.1021/acs.jctc.0c00261

[32] G. D. Carney, L. L. Sprandel, C. W. Kern, Variational approaches to vibration-rotation spectroscopy for polyatomic molecules, Advances in Chemical Physics 37 (1978) 305-379.

[33] S. Carter, N. Handy, A variational method for the calculation of vibrational levels of any triatomic molecule, Molecular Physics 47 (6) (1982) 1445-1455.

[34] B. Sutcliffe, A comment on a recent proposal for the calculation of vibrational energies in the general triatomic molecule, Molecular Physics 48 (3) (1983) 561-566.

[35] T. B. Adler, G. Knizia, H.-J. Werner, A simple and efficient CCSD(T)-F12 approximation, The Journal of Chemical Physics 127 (22) (2007) 221106. doi:10.1063/1.2817618.
URL https://pubs.aip.org/aip/jcp/article/918519

[36] K. A. Peterson, T. B. Adler, H.-J. Werner, Systematically convergent basis sets for explicitly correlated wavefunctions: The atoms H, He, B-Ne, and Al-Ar, The Journal of Chemical Physics 128 (8) (feb 2008).
doi:10.1063/1.2831537
URL https://pubs.aip.org/jcp/article/128/8/084102/71213/protect/discretionary{\char\hyphenchar\font}{}{}Systematically-convergent-basis-sets-for

[37] H.-J. Werner, P. J. Knowles, MOLPRO, A package of ab initio programs (see http://www.molpro.net for more details) (2010).

[38] J. M. L. Martin, M. K. Kesharwani, Assessment of CCSD(T)-F12 Approximations and Basis Sets for Harmonic Vibrational Frequencies, Journal of Chemical Theory and Computation 10 (5) (2014) 2085-2090.
doi:10.1021/ct500174q
URL https://pubs.acs.org/doi/10.1021/ct500174q

[39] D. Lauvergnat, Quantum Model Library: a library of potentials.
URL https://github.com/lauvergn/QuantumModelLib

[40] R. Whitehead, N. Handy, Variational calculation of vibration-rotation energy levels for triatomic molecules, Journal of Molecular Spectroscopy 55 (1-3) (1975) 356-373.
[41] P. Botschwina, Vibrational frequencies from anharmonic ab initio/empirical potential energy functions. I. Method and application to H₂O, HNO, HOF and HOCl, Chemical Physics 40 (1-2) (1979) 33–44.

[42] J. Tennyson, B. T. Sutcliffe, The ab initio calculation of the vibrational-rotational spectrum of triatomic systems in the close-coupling approach, with KCN and H₂Ne as examples, The Journal of Chemical Physics 77 (8) (1982) 4061–4072.

[43] R. Afansounoudji, TI_Schrod: a Fortran code to solve the Schrödinger equation.
URL https://github.com/asodoga/TI_{}Schrod_{}Robert

[44] E. Marsili, M. Olivucci, D. Lauvergnat, F. Agostini, Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model, Journal of Chemical Theory and Computation 16 (10) (2020) 6032–6048. doi:10.1021/acs.jctc.0c00679.

[45] C. Eckart, Some Studies Concerning Rotating Axes and Polyatomic Molecules, Physical Review 47 (1935) 552–558.

[46] D. Lauvergnat, J. M. Luis, B. Kirtman, H. Reis, A. Nauts, Numerical and exact kinetic energy operator using Eckart conditions with one or several reference geometries: Application to HONO, The Journal of Chemical Physics 144 (8) (2016) 084116. doi:10.1063/1.4942172.

[47] D. Lauvergnat, ElVibRot: a Fortran Quantum Dynamics code. (2006).
URL https://github.com/lauvergn/ElVibRot