The analysis and modeling of high-resolution spectra of nonrigid molecules require a specific Hamiltonian and group-theoretical formulation that differs significantly from that of more familiar rigid systems. Within the framework of Hougen–Bunker–Johns (HBJ) theory, this paper is devoted to the construction of a nonrigid Hamiltonian based on a suitable combination of numerical calculations for the nonrigid part in conjunction with the irreducible tensor operator method for the rigid part. For the first time, a variational calculation from *ab initio* potential energy surfaces is performed using the HBJ kinetic energy operator built from vibrational, large-amplitude motion, and rotational tensor operators expressed in terms of curvilinear and normal coordinates. Group theory for nonrigid molecules plays a central role in the characterization of the overall tunneling splittings and is discussed in the present approach. The construction of the dipole moment operator is also examined. Validation tests consisting of a careful convergence study of the energy levels as well as a comparison of results obtained from independent computer codes are given for the nonrigid molecules CH_{2}, CH_{3}, NH_{3}, and H_{2}O_{2}. This work paves the way for the modeling of high-resolution spectra of larger nonrigid systems.

## I. INTRODUCTION

Vibration–rotation spectroscopy remains a formidable tool to understand the spectral features of known molecules or to identify “unknown” molecules in astrophysics, either from laboratory experimental spectra or from accurate theoretical calculations.^{1,2} The standard and historical way of interpreting high-resolution molecular spectra consists of using empirically effective Hamiltonians suitable for data fitting.^{3–6} In turn, some of the resulting fitted parameters may have limited extrapolation power. With the increase in computing capability, particular effort has been made these past two decades in nuclear-motion theory to obtain variational eigenpairs relying on the development of efficient and optimized methods for solving the time-independent Schrödinger equation.^{7–35} Within this context, the derivation of exact quantum kinetic energy operators (KEO) from different coordinate systems as well as the construction of highly accurate *ab initio* potential energy surfaces (PES) remain active fields of research.^{36–49} For semirigid molecules whose PES has a single, well-defined equilibrium configuration, the use of normal coordinates within the framework of the Eckart–Watson Hamiltonian formalism^{50} is well established. Recent papers^{23,51–55} devoted to the construction of accurate line lists for many polyatomic semirigid molecules proved how the Watson Hamiltonian was still of primary importance for the modeling of planetary atmospheres. For nonrigid “floppy” molecules exhibiting one or several large amplitude motions (LAM) connecting multiple potential wells, the use of rectilinear normal coordinates designed for small amplitude vibrations is clearly prohibited. General internal-coordinate Hamiltonians^{19,39,40,47,56–61} were thus constructed in order to study the complex intra-molecular dynamics arising from inversion or internal rotational motions. Indeed, the study of molecules undergoing LAMs is particularly challenging because of their dense and rich spectra due to numerous tunneling splittings.^{62–65} With the spectacular advances in modern experimental techniques, even small splittings can now be resolved, making the development of sophisticated theoretical models necessary to analyze spectra.

Though it does not give any information about the magnitude of the splittings, group theory plays a central role^{66–72} and is inseparable from the theory of nonrigid molecules. It helps to characterize all the tunneling splitting patterns and to provide meaningful labels to the vibration-LAM-rotation states. However, in the absence of a unique equilibrium geometry, it is, of course, obvious that the irreducible representations (irreps) of standard point groups cannot be used to properly label the wave functions of molecules with observable tunneling.^{62,67,68,70,73} Instead, the wave functions of a nonrigid molecule can be classified under the complete nuclear permutation inversion (CNPI) groups.^{66,74} However, when the number of permutation and permutation-inversion elements grows, the use of CNPI becomes rapidly awkward. To avoid manipulating symmetry species associated with tunneling splitting that will never be resolved, it is thus convenient to select only “feasible” operations that bring the molecule into attainable geometries. All these elements form a molecular symmetry (MS) group,^{74} which can be seen as the point group counterpart for floppy molecules and whose irreps help to classify unambiguously the molecular states.

Undoubtedly, curvilinear coordinates are the best choice to describe LAMs, but we can wonder what kind of coordinates can be used to describe the small amplitude vibrations. In 1970, Hougen, Bunker, and Johns^{75,76} (HBJ) derived an exact kinetic energy operator for nonrigid molecules based on a clever combination of normal coordinates for the “rigid” vibrations and curvilinear coordinates for the LAM(s). The other major difference with the standard treatment of rigid molecules lies in the fact that the equilibrium configuration is now replaced by a so-called nonrigid reference configuration depending on the LAM coordinate(s). In the past four decades, the HBJ Hamiltonian has led to the construction of different types of effective Hamiltonians^{77–86} using different types of approximations after applying contact transformation perturbation theory to account for the effects of the small vibrations (more references as well as a non-exhaustive list of computer programs can be found in Sec. III and Table 1 of Ref. 87). To our knowledge, the HBJ Hamiltonian using a KEO expressed in terms of normal coordinates was never used in a full variational calculation to compute vibration–rotation spectra from *ab initio* PESs. Instead, other models, also based on the HBJ philosophy and using a KEO expanded in terms of valence or Morse-like coordinates, were developed and implemented in the MORBID^{56} and TROVE^{57} computer codes.

The aim of this paper is to provide a complete nuclear motion normal-mode HBJ-based Hamiltonian where the LAM coordinate is defined on a grid, combined with advanced group-theoretical methods and optimization tools for solving the Schrödinger equation. The rigid part of the Hamiltonian is treated algebraically with the aid of irreducible tensor operators (ITOs), while the LAM coordinate is treated numerically. A so-called *hybrid* model is thus proposed for the first time. Note that another hybrid model was developed in Ref. 86, but not in that same context. The computation of the vibration-LAM-rotation energy levels is performed variationally using symmetry-adapted primitive basis functions. The construction of a hybrid dipole moment operator is also discussed for line intensity calculation. These past few years, a set of optimized reduction-compression tools was developed to make variational calculations as tractable as possible for semirigid molecules (see, e.g., Refs. 13, and 88–90), so that at this stage, a simple question is raised: *Can the hybrid model benefit from already existing reduction and symmetry tools initially designed for semirigid molecules?* We will try to answer this question in this paper, helped by a series of validation tests and illustrative examples.

In this paper, we will focus on molecules with only one LAM, like inversion or internal rotation, although the extension to several LAMs would be, of course, feasible. Section II recalls the main recipe as well as the basic ingredients that are useful in the HBJ theory, in particular around the choice of the molecule-fixed frame as well as coordinate transformations to conveniently rewrite the potential part. Section III first explains why the HBJ theory is relevant and related to what was previously developed for semirigid molecules. Then, a normally ordered form of the HBJ KEO will be proposed before focusing on symmetry considerations. The construction of the “numerical-ITO” hybrid Hamiltonian and dipole moment operator will be presented, and the choice of the grid points will also be discussed during the computation of the matrix elements. The validation of the present model will be given in Sec. IV, which will also provide a careful examination of the convergence of the energy levels for four nonrigid molecules: CH_{2}, CH_{3}, NH_{3}, and H_{2}O_{2}.

## II. TREATMENT OF A MOLECULE EXHIBITING A LARGE AMPLITUDE MOTION: PRELIMINARIES

### A. Nonrigid reference configuration and choice of the molecule-fixed axis system

In the standard treatment of semirigid molecules, the vibrational displacements of the atomic nuclei are measured with respect to a rigid, equilibrium geometry configuration $aie$ (*i* = 1, …, *N*), with *N* being the number of atoms, and in that case both the reciprocal *μ*_{αβ} tensor contained in the KEO as well as the potential *V* are well approximated by the leading terms in a Taylor series expansion. For nonrigid molecules for which the amplitude of at least one coordinate is not small compared to the equilibrium bond lengths and bond angles, the representation of *μ*_{αβ} and *V* as a Taylor series fails. Therefore, the potential function of the LAM cannot be properly described, even by increasing the order of the polynomial expansion. In order to prevent such dramatic convergence issues, Hougen, Bunker, and Johns,^{75} in their seminal work on quasilinear triatomic molecules, suggested removing the bending motion from the vibrational part and incorporating this LAM into the rotational part. To this end, they (i) introduced a so-called nonrigid reference configuration $airef(\rho )$ instead of $aie$, where *ρ* is a coordinate describing the LAM, and (ii) treated the remaining 3*N* − 7 small amplitude vibrations using the normal mode coordinates on a numerical grid in *ρ*. By definition, the displacements of the coordinates describing the small amplitude motions are assumed to be zero along the reference configuration. The HBJ approach thus ensures that these 3*N* − 7 coordinates vary smoothly around $airef(\rho )$, making relevant a Taylor series expansion of both the KEO and *V* in terms of normal coordinates.

Though initially developed for quasilinear triatomic molecules, the HBJ approach can be finally applied to a wide class of nonrigid polyatomic molecules exhibiting different kinds of LAM (bending, inversion, or internal rotation). The major difference between the HBJ Hamiltonian^{75} and its rigid counterpart, namely the Eckart–Watson Hamiltonian,^{50} lies in the fact that all the coefficients involved in the KEO and in the Taylor expansion of the potential part will now be *ρ* dependent (see Secs. II B and III B).

**I**

^{ref}when all the displacements

**d**

_{i}are zero (see, e.g., the discussion of Flaud and Perrin in Ref. 91). The PAS allows us to remove the inertial cross-terms $I\alpha \beta ref$. In order to minimize the rotation-LAM coupling by removing the cross terms $I\alpha \rho ref$ (

*α*=

*x*,

*y*,

*z*) in

**I**

^{ref}, IAS is a good choice. Mathematically speaking, this latter condition implies that the angular-momentum-like vector will vanish due to

*ρ*,

_{1}, and

**U**(

*ρ*) is the

*ρ*-dependent rotation matrix to be determined. This choice is analogous to the original HBJ axis system for triatomic molecules,

^{75}followed soon after by a series of works where an auxiliary angle

^{79,81,92}

*ϵ*(

*ρ*) satisfying Eq. (1) was introduced instead of

**U**(

*ρ*). A recent method was proposed in Ref. 93 for the derivation of the rotation matrix

**U**(

*ρ*) and is briefly recalled in the Appendix. In the RAS, the angular momentum due to the LAM is not vanishing but is of constant magnitude along the

*z*axis.

^{6}

As an illustrative example, we have determined the rotation matrix **U**(*ρ*) in order to compute a new reference configuration satisfying (1) for the H_{2}O_{2} molecule. As a starting point, we have considered the reference configuration given in Ref. 94, called $airef,ini(\rho )$ such as *I*_{yρ} ≠ 0. By solving Eq. (A6), the new IAS configuration is given by $airef,IAS(\rho )=U(\rho )airef,ini(\rho )$. The *x* components of the reference configuration for the atoms H1 and O1, before and after transformation, are displayed in Fig. 1 as a function of *ρ*. It is worth mentioning that within the IAS treatment, the torsion angle *ρ* will be defined as half of the dihedral angle when studying internal rotation. As a direct consequence, the rotational and torsion angles, *χ* and *ρ*, will be double valued, making the use of irreps of extended MS groups necessary to classify the rotational and torsional energy levels.

Moreover, according to the relation between the laboratory and molecule-fixed axis systems, it turns out that the infinitesimal vibrational displacement vectors **d**_{i} must be subjected to seven constraints. There are three equations corresponding to the center of mass condition; three other ones have a form similar to the Eckart conditions; and the last one (the Sayvetz condition^{95}) minimizes the interaction between the small and large amplitude vibrations.

### B. Coordinate transformation and expansion of the potential energy surface

In Sec. II A, the LAM coordinate *ρ* as well as the notion of reference configuration attached to a molecule-fixed axis system have been introduced. As stated by Jensen and Szalay,^{81,96,97} *ρ* can be considered a kind of “effective” coordinate that is different from the geometrically defined curvilinear $\rho \u0304$ coordinate when the “rigid” vibrations are beyond their equilibrium geometry. By construction, the HBJ KEO depends explicitly on *ρ* and on 3*N* − 7 rectilinear, normal coordinates *Q*_{i}, while the *ab initio* potential function is built from a function $u\u0304(\rho \u0304)$, where $\rho \u0304$ describes a bond, torsion, or inversion angle, and from a set of curvilinear coordinates (i.e., valence, Morse, cosine, etc.) to describe the small amplitude vibrations. In this work, we assume that the potential function *V* is known and expressed in terms of 3*N* − 6 coordinates $Sk,\sigma k(\Gamma k)$ adapted to the symmetry of the molecule, where Γ_{k} is an irreducible representation of a molecular symmetry group and *σ*_{k} is a component when dim(Γ_{k}) > 1. Among these $Sk,\sigma k(\Gamma k)$ coordinates, one corresponds to $u\u0304(\rho \u0304)$ and needs to be transformed to a *ρ*-dependent function to be consistent with the KEO. To this end, it is more convenient to rewrite the PES in terms of the 3*N* − 7 coordinates $Sr,\sigma r(\Gamma r)$ for the rigid part and put the dependence in $\rho \u0304$ into the expansion coefficients. As a direct consequence, some linear terms like $\u2211rF\u0304r(\rho \u0304)Sr(\Gamma r)$ will now appear in the potential, while the quadratic part will be written as $\u2211rF\u0304r,r\u2032(\rho \u0304)Sr(\Gamma r)Sr\u2032(\Gamma r\u2032\u2032)$. Therefore, contrary to the “ordinary” treatment of semirigid molecules, vibrations with different symmetries can interact in the quadratic part of the potential for nonrigid molecules.

In order to solve the *GF* problem, we first need to express the PES as a function of the *ρ* coordinate instead of $\rho \u0304$. To this end, a general numerical recipe was proposed in Ref. 87 to find the relation $\rho \u0304=f(\rho ,Sr)$ for any polyatomic molecules so that the potential function transforms as $V(Sr,\rho \u0304)\u2192V(Sr,\rho )$. The procedure can be summarized as follows:

- First, we assume that all calculations are performed on a grid composed of
*M*_{ρ}points,where the choice of these points will be discussed in Sec. III C 4. As already mentioned, the main drawback of the HBJ approach lies in the fact that the kinematic matrix of dimension (3(2)$\rho m\u2208[\rho min,\rho max],m=1,\u2026,M\rho ,$*N*− 7) × (3*N*− 7) writes as**G**=**G**(*ρ*_{m}) due to its dependence on $ai\alpha ref(\rho )$, whereas the elements of the force constant matrix write $F\u0304rr\u2032=F\u0304rr\u2032(\rho \u0304)$. The aim is thus to find the transformation $F\u0304(\rho \u0304)\u2192F(\rho m)$ to properly solve the*GF*problem on the grid. To this end, we have shown in Ref. 87 that $u\u0304(\rho \u0304)$ must necessarily be of the formwhere the(3)$u\u0304(\rho \u0304)=u(\rho m)+\u2211rCr(\rho m)Sr(\Gamma r)+\u2211rr\u2032Crr\u2032(\rho m)Sr(\Gamma r)Sr\u2032(\Gamma r\u2032),$*C*_{r}and*C*_{rr′}coefficients are to be determined. In this approach, only the linear and quadratic terms are needed, unlike the MORBID approach, where higher-order expansion coefficients*C*_{rr′r″}, $Crr\u2032r\u2033r\u2032\u2032\u2032$, etc. were required.^{56} - In the next step, we expand the 3
*N*− 7 coordinates $Sr(\Gamma r)$ as well as the seven Eckart–Sayvetz relationsin terms of the 3(4)$T\alpha =1M\u2211imidi\alpha ,R\alpha =1I\alpha \alpha ref(\rho )\u2211imi\u03f5\alpha \beta \gamma ai\beta ref(\rho )di\gamma ,S=1I\rho \rho ref(\rho )\u2211i,\alpha midai\alpha ref(\rho )d\rho di\alpha ,$*N*small amplitude Cartesian displacements**d**_{i}. Here,*M*is the total mass of the system,*ϵ*corresponds to the Levi-Cività tensor, and $I\alpha \alpha ref$, $I\rho \rho ref$ are the diagonal elements of the 4 × 4 inertia matrix.^{75}In the standard treatment of semirigid molecules, only the linearized part**S**=**Bd**is required for solving the*GF*problem, but here we also need to compute the quadratic terms**d**_{i}**d**_{j}to determine the*C*_{r}and*C*_{rr′}coefficients in Eq. (3). To this end, we first construct the 3*N*linear and 3*N*(3*N*+ 1)/2 quadratic forms $Xslin$ and $Xs\u2033quad=XslinXs\u2032lin$ [*s*,*s*′ = 1, …, 3*N*,*s*≥*s*′;*s*″ = 1, …, 3*N*(3*N*+ 1)/2] with $Xlin=(S,T\alpha ,R\alpha ,S)$. Similarly, we define the vector**Y**^{lin}= (*d*_{1x}, …,*d*_{Nz}) and also construct the quadratic form $Ys\u2033quad=YsYs\u2032$. Finally, we follow closely the strategy established in Ref. 98 to make this problem linear by defining the transformation**B**_{quad}on the grid such thatThe Cartesian displacements contained in(5)$(Xlin,Xquad)t=Bquad(\rho m)(Ylin,Yquad)t.$**Y**can thus be deduced by inverting Eq. (5) and by the seven constraints $T\alpha =R\alpha =S=0$. We thus writewith now(6)$(d1,\u2026,dN)t=Aquadinv(\rho m)(Xlin,Xquad)t,$**X**^{lin}= (**S**, 0, 0, 0) and where $Aquadinv$ is a sub-matrix extracted from $Bquad\u22121$. In this scheme, the nuclear displacement vectors**d**_{i}are variables depending on both the symmetry and*ρ*coordinates. - In the last step, we express the function $u\u0304(\rho \u0304)$ in terms of Cartesian coordinates and make the substitution (6) to finally get the desired coefficients
*C*_{r}and*C*_{rr′}in Eq. (3). We are now able to find the eigensolutions of the**GF**matrix and convert the*ab initio*curvilinear potential function to a normal-mode polynomial expansion for each*ρ*_{m}using Eq. (3) as well as the 3(7)$V(Sr,\rho \u0304)\u2192V(Qr,\rho m)\u2261V\rho m(Qr),r=1,\u2026,3N\u22127,$*N*× (3*N*− 7) orthogonal transformation**l**(*ρ*_{m}) linking the mass-weighted Cartesian displacements to the normal coordinates.^{99}

In Ref. 87, we intended at this stage to fit all the matrix elements *l*_{iα,s}(*ρ*_{m}) as well as all the Coriolis coefficients $\zeta ss\u2032\alpha (\rho m)$ to a polynomial function in *u*(*ρ*) in order to expand both the KEO and potential parts in terms of the 3*N* − 6 coordinates (*Q*_{i}, *ρ*), like for semirigid molecules. In this work, the nonrigid part of the HBJ Hamiltonian will be treated numerically for more flexibility, which is the best choice if certain values of *ρ* are near a singularity.

## III. CONSTRUCTION OF A HYBRID HBJ-BASED MODEL

### A. Motivations

In the past few years, a general methodology based on the normal-mode Eckart–Watson formalism and on the extensive use of symmetry has been proposed for semirigid molecules. It also combines dimensionality reduction techniques to manage memory at each stage of calculation. All the procedures have been implemented in the TENSOR computer code and have led to the construction of accurate line lists for molecules up to seven atoms.^{13,21,51–53} The main motivations for the extension to nonrigid systems can be summarized as follows:

*Reduction-compression techniques*. The use of the HBJ theory is a quite natural choice because most of the reduction-compression tools previously developed for semirigid molecules can also be considered and adapted to tackle the problem of nonrigid polyatomic molecules.*A general formulation*. Like for the Watson KEO, the form of the HBJ KEO remains the same, whatever the number of atoms. Moreover, the use of the*ρ*coordinate allows an optimal separation between small and large amplitude nuclear motions, as achieved by the Sayvetz condition.*Variational calculation*. The HBJ model inspired a series of works and led to the development of both effective vibration-LAM-rotation Hamiltonians suitable for data fitting^{77–86}and nuclear-motion Hamiltonians, as those implemented in the MORBID^{56}and TROVE^{57}computer codes. In this work, we intend to use the HBJ Hamiltonian in its original formulation based on the 3*N*− 7 normal coordinates, including all terms in Eqs. (10)–(14) below, for solving variationally the vibration-LAM-rotation wave equation.*Extensive use of symmetry*. HBJ should also benefit from the use of the ITO formalism for a full account of symmetry properties. The use of ITOs is probably the most efficient way to deal with degenerate vibrations when the molecular symmetry group is non-Abelian. A suitable combination of ITOs with a numerical treatment of the LAM will lead to the introduction of a so-called*hybrid*Hamiltonian model for the spectra calculation of nonrigid molecules.*Toward a nonrigid effective model*. Effective Hamiltonian computer codes designed to calculate and fit energy levels and line transitions for molecules containing one or several internal rotors already exist and were developed over the past three decades (see, e.g., the XIAM,^{100}BELGI-C_{1},^{101}BELGI-C_{s},^{102}BELGI-C_{s}-2Tops,^{103}ERHAM,^{104}or RAM36^{85}codes). Recently, we have proposed a numerical approach for the construction of*ab initio*effective models,^{90}which is a clear alternative to the rather involved Van Vleck perturbation method. Undoubtedly, this approach could be of great help to extract an effective nonrigid Hamiltonian from our hybrid model and to refine parameters based on the experimental data, even when several interacting vibrational states are involved.

### B. Expansion of the HBJ kinetic energy operator

*ab initio*potential energy function

*V*in Eq. (7) can be expanded on a numerical grid in terms of 3

*N*− 7 normal coordinates

*Q*. The KEO can also be expressed as a sum of products in the same fashion. In that case, each integral is a product of integrals associated with each mode, and the computation of the matrix elements can benefit from the use of the Wigner–Eckart theorem (see Sec. III C 4). In this work, we propose to put each term of the expansion into a better suited normally ordered form of the type

*X*(

*Q*)

*Y*(

*P*)

*Z*(

*J*

_{α}), where

*X*and

*Y*are not necessarily commuting functions. Such a form makes it quite convenient to build the KEO in any order and in a systematic manner. Unless specified otherwise, the running indices

*j*,

*k*,

*l*,

*r*, and

*s*will take values 1, …, 3

*N*− 7, while the Greek letters

*α*and

*β*will equal

*x*,

*y*,

*z*, or

*ρ*. We can show that the HBJ KEO can be rewritten for each point

*ρ*

_{m}of the grid as

*P*

_{j}= −

*iℏ∂*/

*∂Q*

_{j}is the conjugate momentum of the normal coordinates

*Q*

_{j}, while

*J*

_{α}(

*α*=

*x*,

*y*,

*z*) and

*J*

_{ρ}= −

*iℏd*/

*dρ*are the components of the total angular momentum and LAM coordinate, respectively. If

*α*≠

*β*, the terms containing

*J*

_{α}

*J*

_{β}are sorted as

*J*

_{x}

*J*

_{y},

*J*

_{x}

*J*

_{z}, and

*J*

_{y}

*J*

_{z}. For the sake of simplicity, the symmetry labels have been omitted from all the coordinates and conjugate momenta that are simply denoted here by $Qk\sigma (\Gamma )\u2261Qj$, $J\rho (\Gamma )\u2261J\rho $, etc. To be rigorous, these symmetry labels should appear everywhere, and in that case, the

*k*index in $Qk\sigma (\Gamma )$ would run over the vibrational modes instead of the vibrational coordinates.

*Q*

_{i},

*P*

_{i},

*J*

_{ρ},

*J*

_{α}) and is computed for each

*ρ*

_{m}, like for the potential part. Using the well-known commutation rules and some algebra, it is straightforward to show that

*μ*is the determinant of the 4 × 4 reciprocal inertia matrix with the components

*μ*

_{αβ}computed from the inverse of $I\alpha \beta \u2032$ (see, e.g., Refs. 75 and 81 for all the definitions and for the expression of the Coriolis coefficients

*ζ*). The various quantities involved in Eqs. (10)–(14) are also evaluated for each

*ρ*

_{m}, that is,

*μ*

_{αβ}=

*μ*

_{αβ}(

**Q**,

*ρ*

_{m}),

*ζ*=

*ζ*(

*ρ*

_{m}), $\mu \alpha \beta \u2032=\mu \alpha \beta \u2032(Q,\rho m)$, or (ln

*μ*)′ = (ln

*μ*)′(

**Q**,

*ρ*

_{m}).

The expression (8) of the normally ordered HBJ KEO as well as the numerical procedure in Sec. II have been implemented in an updated version of our TENSOR computer code. In order to find the Taylor series of (7) and (8) in terms of the 3*N* − 7 normal coordinates, the crucial point is the computation of the successive derivatives.

The first and second derivatives in

*ρ*involved in (8) are evaluated numerically using the five-point central finite difference formula and quadruple-precision floating-point numbers.To circumvent the problem of round-off errors encountered in the finite difference method for high-order derivatives with respect to the 3

*N*− 7 normal coordinates, a modified version of the computer program COSY INFINITY,^{105}including the dynamic allocation and the quadruple precision, was implemented in the TENSOR code. It allows high-order multivariate automatic differentiation with almost no loss of precision. The Taylor series expansions of*V*,*μ*_{αβ}, and ln*μ*and their derivatives with respect to*Q*_{i}and*ρ*are thus all performed using the COSY algorithm. The mixed derivative $d(\u2202\u2061ln\u2061\mu \u2202Qs)/d\rho $ involved in (14) combines finite difference and COSY calculation. The convergence of both the KEO and potential parts will be discussed in Sec. IV in the case of 3- and 4-atomic molecules exhibiting a large bending, inversion, and internal rotation motion.

Finally, though the current version of the TENSOR code is dedicated to molecules exhibiting only one LAM, its extension to *N*_{ρ} large amplitude coordinates (e.g., the hydrazine or methylamine molecule^{63}) could be carried out with only a few modifications. As recalled by Bunker and Jensen^{74} in their Eqs. (15-14) and (15-15), the new inertia matrix will be of dimension (3 + *N*_{ρ}) × (3 + *N*_{ρ}), with some coupling terms between the different LAM coordinates *ρ*_{1}, *ρ*_{2}, etc. Therefore, the HBJ Hamiltonian could be expressed in terms of 3*N* − 6 − *N*_{ρ} normal-mode coordinates and of the 3 + *N*_{ρ} rotational components *J*_{x}, *J*_{y}, *J*_{z}, $J\rho 1$, $J\rho 2$, etc.

### C. Tensorial Hamiltonian and dipole moment operators on a numerical grid: The hybrid model

#### 1. Symmetry considerations

We just recall here the basic principles and main definitions that can be useful when dealing with symmetry, which plays a central role in nonrigid molecules. Indeed, the formidable development in experimental techniques makes it possible to resolve small tunneling splittings that are calculated quantum-mechanically, while group theory provides meaningful labels. It is not the purpose of the present work to give a detailed review of all the symmetry properties of a given molecule and to study how the vibration-LAM-rotation wave functions transform in the operations of a symmetry group. Extensive group-theoretical considerations can be found, for example, in the excellent reviews^{66,67,69,73,74,106–108} (and references therein).

By definition, a CNPI group is the direct product of permutation groups $Sm$ with the inversion group *ɛ* = {*E*, *E**}, where *E** is the laboratory-fixed inversion operation. For example, the CNPI group *G*_{n} consists of *n* nuclear permutations and nuclear permutation inversion operations. A molecular symmetry group is formed by all “feasible” operations of the CNPI group and can be denoted as *G*_{n′} (*n*′ ≤ *n*) or as $G(M)$, where $G$ is the corresponding point group. Typically, in the first version of the TENSOR computer code designed for semirigid molecules, both the Abelian (*C*_{2}, *C*_{i}, *C*_{s}, *C*_{2v}, *C*_{2h}, *D*_{2} and *D*_{2h}) and non-Abelian (*D*_{n}, *C*_{nv}, *D*_{nh} (*n* ≥ 3), *D*_{nd}, *T*_{d}, *O*_{h}) point groups were implemented.

*Non-tunneling molecules*: for molecules where no tunneling splitting is observed, the extension of TENSOR to the MS groups is direct because of their isomorphism with the point groups. However, the MS group is not necessarily the same as the CNPI group. For example, for ammonia, the CNPI group $S3\xd7\epsilon $ or *G*_{12} is isomorphic to the group *D*_{3h}(*M*) or to the point group *D*_{3h} of the planar configuration, while for phosphine, *G*_{12} is not isomorphic to the group *C*_{3v}(*M*) or *C*_{3v}.

*Tunneling molecules*: if at least one barrier in the description of the LAM is not “insuperable,” a tunneling splitting may be observed, and the analysis of the spectra of nonrigid molecules will require the use of MS groups that can be significantly different from the point groups. It may also arise that the MS group of a tunneling molecule is isomorphic to a point group, but most of the time, it is not really possible to classify the states of molecules with a LAM using a familiar point group alone. For molecules containing identical coaxial rotors such as ethane or hydrogen peroxide, it is necessary to use an *extended* MS group^{74} *G*_{m}(*EM*) = *G*_{m} × {*E*, *E*′}, where the operation *E*′ consists in increasing a rotational angle (*χ*_{a} or *χ*_{b}) by 2*π*. Such groups can also be isomorphic to a point group [e.g., *G*_{4}(*EM*) ∼ *D*_{2h} for hydrogen peroxide].

More generally, for symmetry groups consisting of many operations, we have to face the determination of the character tables and symmetry species, which can be tedious. In this context, Hougen^{109} presented a series of “recipes” aimed at helping the user handle permutation-inversion groups for practical spectroscopic applications. As stated by Woodman,^{67} the most common way to handle a large MS group is to find a structure with smaller subgroups. Indeed, many large groups can be expressed as a *group product*, making the determination of the symmetry species and generators easier. For example, for ethane-like molecules, we have the direct product structure^{72}^{,} $G36=C3v+\xd7C3v\u2212$, consisting of 36 = 6 × 6 feasible operations, where +/− is used to distinguish the two methyl groups. It was also shown in Ref. 67 that the symmetry groups of nonrigid molecules can be expressed as semi-direct products of the type $A=B\u2227C$ if only $B$ is an invariant subgroup of $A$. Alternatively, a method based on nonrigid group theory^{69,110} was also proposed to handle very large MS groups.

Once the character table and generators of a given MS group are known, the Clebsch–Gordon coupling coefficients as well as the recoupling 6*C* and 9*C* Wigner’s symbols involved in the calculation of matrix elements [see Eq. (30) below] can be thus deduced. For example, the Clebsch–Gordon coefficients adapted to an MS group can be computed by following the method presented in Ref. 111. A set of irreducible tensor operators for the small-amplitude vibrational, rotational, and large-amplitude vibrational motions, as well as symmetry-adapted vibration-LAM-rotation wavefunctions, can also be deduced from symmetry considerations.

We will assume hereafter that the character table, multiplication table, and symmetry species Γ_{i} of a given MS group are known in order to (i) compute and store the coupling, recoupling, and other related coefficients involved in the tensorial formalism and (ii) find the transformation laws of the wavefunctions in the operations of the MS group. The general theory presented in Sec. II as well as in Subsections III C 2 – III C 4 below is now implemented in the updated version of TENSOR.

#### 2. Hybrid Hamiltonian

For asymmetric top molecules with only one-dimensional irreps, accounting for symmetry is trivial because matrix and character coincide in group theory. Conversely, dealing with non-Abelian groups whose dimensions may be greater than 1 is somewhat more challenging and requires more attention. For example, two-fold or three-fold degenerate vibrations are involved for point groups characterized by the presence of a symmetry axis higher than two-fold or in MS (or double) groups like *G*_{6}, *G*_{12}, *G*_{24}, *G*_{36}, etc. Within this context, the use of ITOs is relevant to deal with degenerate vibrational modes like *ν*_{3} and *ν*_{4} of the NH_{3} molecule. Historically, and still nowadays, the ITOs were introduced to analyze vibration–rotation infrared spectra of methane-like *T*_{d} molecules using effective polyad models and the so-called tetrahedral formalism.^{112} They have been extended later on to other molecular species in various formalisms.^{113–116} More recently, we have shown that the ITOs could also be used for the construction of complete nuclear-motion Hamiltonians and *ab initio* dipole moment operators.^{88} They were thus employed with success in the construction of spectroscopic line lists of semirigid molecules from the first principles.^{13,21,51–53}

*ρ*

_{m}, the “initial” HBJ Hamiltonian $H\rho mHBJ=T\rho mHBJ+V\rho m$ writes as a sum of products and takes the simple form

*ρ*

_{m}in (2), the Hamiltonian parameters can be stored in a two-dimensional array, that is,

*h*

_{j}(

*ρ*

_{m}) ≡

*h*

_{jm}. Each row is a parameter associated with a given operator $Tj$, and each column is a point of the numerical grid (2). At this stage, we follow the strategy previously established for semirigid molecules and assume that the Hamiltonian (15) can be rearranged into a sum of vibration-LAM-rotation ITOs,

_{0}of the MS group, i.e., Γ

_{vρr}= Γ

_{0}. In Eq. (17),

*V*,

*R*, and

*L*stand for the vibrational, rotational, and LAM irreducible tensors of degree Ω

_{v}, Ω

_{r}, and Ω

_{ρ}in (

*Q*,

*P*), (

*J*

_{x},

*J*

_{y},

*J*

_{z}), and

*J*

_{ρ}, respectively.

*ɛ*is the parity in the conjugate momenta

*P*such that $\epsilon =(\u22121)\Omega r+\Omega \rho $ due to the time reversal invariance. More details about the definition and construction of

*V*and

*R*can be found elsewhere,

^{112,115–121}even for open-shell molecules in a degenerate electronic state,

^{122–126}but only in the frame of effective polyad Hamiltonians using creation-annihilation operators. Our vibrational ITOs can be built in various coordinate systems and account for all intra-polyad coupling terms.

*V*is built as the tensor product of

*N*

_{m}− 1 operators $Vk\Theta k(\Gamma k)$, where

*N*

_{m}is the total number of vibrational modes, while the single mode associated with the LAM is incorporated into

*L*. Each $Vk\Theta k(\Gamma k)$ is constructed from the symmetrized powers $Qk$ and $Pk$ of degree Ω

_{kq}and Ω

_{kp}in the small amplitude coordinates $Qs(\Gamma )$ and the conjugate momenta and $Ps(\Gamma )$ involved in the vibrational mode

*k*,

**α**

_{q}and

**α**

_{p}. Here, Ω

_{v}=

*∑*

_{k}Ω

_{qk}+ Ω

_{pk}and

*σ*

_{k}is a component when dim(Γ

_{k}) > 1. The rotational part

*R*is built recursively in terms of

*J*

_{x},

*J*

_{y}and

*J*

_{z},

*K*

_{r}is the rank of the tensor in SO(3), and

*α*

_{r}is a rotational labeling that will be different for spherical, symmetric, and asymmetric-top molecules (see Sec. III C 4).

_{ρ}in

*J*

_{ρ}and of symmetry Γ

_{ρ2}. To be consistent, the

*ρ*-dependent tensorial parameters

*s*

_{j}to be determined have been naturally included in

*L*. The identity operator

*I*introduced in (19) reflects the fact that the LAM function

*s*

_{j}(

*ρ*) possesses a certain symmetry Γ

_{ρ1}in the MS group (see Fig. 2). Usually, there are two types of symmetry species Γ

_{ρ1}for non-tunneling molecules or for molecules with one tunneling like ammonia for which $\Gamma \rho 1=A1\u2032$ or $A2\u2033$. When dealing with several tunnelings, more symmetry species Γ

_{ρ1}can be involved. For example, the allowed symmetry species Γ

_{ρ1}of the extended MS group

*G*

_{4}(

*EM*) will be

*A*

_{gs},

*A*

_{us},

*A*

_{gd}, and

*A*

_{ud}for the H

_{2}O

_{2}molecule. Note that H

_{2}O

_{2}is among these molecules for which the symmetry species of a point group (here,

*D*

_{2h}) can be used due to the one-to-one mapping of its elements with those of

*G*

_{4}(

*EM*).

*s*

_{j}(

*ρ*

_{m}). To this end, we expand the vibrational, rotational, and LAM ITOs in terms of elementary operators using the Clebsch–Gordan coefficients [e.g., see Eq. (21) of Ref. 98], sorted in such a way that the Hamiltonian (17) writes as

**Z**is a group symmetry transformation. By equating Eqs. (15) and (20) and bearing in mind that the two sums involved in these two equations do not necessarily contain the same number of terms if degenerate vibrations are involved, a set of parameters {

*s*

_{j}} is obtained by solving for each point

*ρ*

_{m}the following overdetermined system of equations:

*dgels*routine of the LAPACK (Linear Algebra PACKage) library was used. It generally takes few seconds per point

*ρ*

_{m}to solve (21) and obtain a set of tensorial parameters {

*s*

_{j}}. As an illustrative example, let us consider four vibration–rotation–torsion coupling terms between the

*ν*

_{4}(

*A*

_{1u}),

*ν*

_{5}(

*B*

_{1u}), and

*ν*

_{6}(

*B*

_{1u}) modes of H

_{2}O

_{2}written in the form (17) as

*D*

_{2h}are used here instead of those of

*G*

_{4}(

*EM*). The four parameters

*s*

_{j}are plotted in Fig. 2 for each value

*ρ*

_{m}from 0 to 4

*π*. In order to check the overall symmetry of the Hamiltonian (22), these four grids composed of 100 points have been fitted by analytical functions. When dealing with internal rotation, it is quite obvious that the Fourier series (cos(

*kρ*/2), sin(

*kρ*/2)) are well adapted to describe periodic functions, where the symmetry species will depend on the parity of

*k*, which is also related to the parity

*K*of a rotational wave function (

*K*is the projection of

*J*on the molecule axis

*z*). In our case, we obtain

*E**(

*ρ*/2) =

*π*−

*ρ*/2 and

*E*′(

*ρ*/2) =

*π*+

*ρ*/2, the functions

*s*

_{1}(

*ρ*),

*s*

_{2}(

*ρ*),

*s*

_{3}(

*ρ*), and

*s*

_{4}(

*ρ*) transform as

*A*

_{g},

*B*

_{2g},

*B*

_{2u}, and

*A*

_{u}, making the Hamiltonian (22) totally symmetric. We can thus see that the

*A*

_{g}/

*A*

_{u}functions are 2

*π*-periodic while the

*B*

_{2g}/

*B*

_{2u}functions are 4

*π*-periodic. Note that only the

*A*

_{g}/

*A*

_{u}species are involved in the vibration-LAM Hamiltonian.

#### 3. Hybrid dipole moment

*Z*component is required. In order to be consistent with the hybrid Hamiltonian model, (i) the Eckart–Sayvetz conditions are used for the orientation of the molecule-fixed frame (

*xyz*) to conveniently define the

*ab initio*components $\mu \alpha m$ and (ii) symmetry-adapted ITOs are used for the construction of the space-fixed frame

*Z*component as follows:

*H*. Moreover, we have the condition Ω

_{p}= Ω

_{ρ}= Ω

_{r}= 0. The tensor dipole moment parameters are determined by solving a system of equations like (21) for each

*ρ*

_{m}. In Eq. (24), $C$ is the direction cosine tensor.

^{112}

Finally, the quality of the line intensity calculation will be governed by two ingredients: the calculation of the eigenpair {*E*, Ψ} by solving the time-independent Schrödinger equation and an accurate dipole moment surface, whose construction is beyond the scope of this work.

#### 4. Hamiltonian matrix elements and choice of the LAM grid

*V*=

*R*=

*I*in (17),

*ρ*=

*ρ*

_{m}, where the

*s*

_{j}(

*ρ*) are not necessary known functions, though some of them could be easily related to the potential or to the tensor inertia matrix. For example,

*s*

_{1}(

*ρ*) is written as

*aρ*

^{4}−

*bρ*

^{2}+ ⋯ for an inversion potential, or

*s*

_{1}(

*ρ*) is expressed as a Fourier series for a torsional potential, while $s3(\rho )\u223c\mu \rho \rho ref(\rho )/2$. If we consider a set of LAM basis functions $\varphi v\rho (C\rho )$ of symmetry

*C*

_{ρ}, the matrix elements for each term of the Hamiltonian (26) can be computed numerically using the standard quadrature formula

*v*

_{ρ}is the principal LAM quantum number taking values 0, 1, …, and

*w*

_{m}are the quadrature weights. It is thus demonstrated after a brief inspection of Eqs. (25) and (27) that the computation of the matrix elements for the LAM problem is greatly simplified if a set of

*M*

_{ρ}quadrature points is used to define the grid (2). In that case, we just need to read the stored parameters

*s*

_{j}(

*ρ*

_{m}) of the LAM Hamiltonian (15) to compute the matrix elements (27). In this work, a Gauss–Legendre quadrature of ∼100 points at quadruple precision was considered to make a numerical integration between

*ρ*

_{min}and

*ρ*

_{max}. Moreover, the quadrature points generally go around singularity points; e.g., the singularity in

*μ*

_{αβ}at the linear configuration of a triatomic molecule is avoided because the first quadrature point takes a nonzero positive value when integrating from

*ρ*

_{min}= 0 to

*ρ*

_{max}.

Let us now focus on the calculation of the matrix elements for parts *V* and *R*. In order to compute the matrix elements of the full Hamiltonian (17), Eq. (26) can be generalized to the whole vibration-LAM-rotation problem, and the superscript *j* in Eq. (27) will run over all the vibration-LAM-rotation tensor operators. The key to the present approach is thus to compute and store in memory the matrix elements $Lvl\u2032vlj$ in a LAM basis set of all the *ρ*-dependent Hamiltonian parameters. This procedure takes a few seconds because the LAM basis contains generally less than 20 functions and the Hamiltonian has a few hundred or thousands of parameters.

_{v}, Φ

_{rot}, and $\Phi v\rho $ are symmetry-adapted vibrational, rotational, and LAM functions, respectively, described briefly as follows:

*Small amplitude basis functions*: The vibrational function Φ_{v} is a product of *N*_{m} − 1 harmonic oscillator functions. Their construction is trivial for non-degenerate modes. The treatment of two- and three-fold degenerate vibrations requires more attention. For example, a method was proposed in Ref. 115 in the case of *C*_{3v} and *T*_{d} molecules and extended later on to arbitrary point groups.^{127} The treatment of MS groups can be carried out in the same fashion.

*Rotational basis functions*: The construction of the symmetry-adapted Φ

_{rot}rotational functions will depend on which family the molecule belongs to. For spherical top molecules, a method was proposed in Ref. 128 to build such functions for both integer and half-integer

*J*values, that is, for both vectorial and spinorial representations.

^{129}In that case,

*α*in (28) plays the role of a multiplicity index. For symmetric and asymmetric tops, we can consider symmetrized Wang-type basis functions

*G*coefficients, including a consistent choice for the phase factors, have also been determined

^{127}for arbitrary point groups and implemented in TENSOR. Again, its extension to MS groups can be carried out quite easily.

*Large amplitude basis functions*: The choice of the functions Φ_{ρ} will depend on the nature of the LAM (Sec. IV).

^{130}for the small amplitude vibration and the rotation parts, combined with a numerical integration for the LAM part. The matrix elements for each term of the sum in Eq. (17) are given by

*V*‖⋯) and (⋯‖

*R*‖⋯) are vibrational and rotational reduced matrix elements,

^{112,115}and the 9

*C*are the recoupling Wigner’s symbols computed from the symmetry-adapted Clebsch–Gordan coefficients (see, e.g., Ref. 111 for their definition). The rotational reduced matrix elements involved the so-called isoscalar factors, which are computed from the

*G*coefficients in (29), the Clebsch–Gordon coefficients, as well as the 3

*j*Wigner symbols. [

*C*] is the dimension of the irrep

*C*, and $Mvj$ is a factor that contains a product of 9

*C*symbols characterizing the vibrational inner coupling scheme. $Lv\rho \u2032v\rho j$ are the matrix elements (27) of the LAM operators computed by numerical quadrature. From the selection rules, the use of symmetry thus allows us to diagonalize separately each block labeled by

*J*and

*C*

_{vρr}. The computation of the matrix elements for the dipole moment will be similar.

Finally, our hybrid model combines the rapidity of group-theoretical methods (via the use of ITOs and the Wigner–Eckart theorem) for the rigid part and the flexibility of numerical quadrature integration for the nonrigid part. All the elements (30) are generally computed from a few seconds up to some tens of minutes by using standard libraries for massive parallel programming, even for large basis sets composed of ∼100 000 functions. Here, all the matrices were diagonalized using direct eigensolvers.

#### 5. Reduced Hamiltonian and reduced vibrational eigenvectors

One of the main advantages of using the Watson or HBJ Hamiltonian lies in the fact that the form of their KEO remains unchanged, whatever the number of atoms. In turn, the number of terms in the Hamiltonian (15) may increase dramatically when performing a Taylor series expansion of the multivariate *μ*_{αβ} and *V* functions. A similar problem arises when managing the number of vibration-LAM functions in the direct-product basis (28) to compute variational solutions for *J* > 0. The number of rovibrational functions becomes rapidly huge, even by pruning the basis using Eq. (32).

In order to make vibration-LAM-rotation calculations as tractable as possible for the user, a strategy was proposed in Refs. 13, and 88–90 for semirigid molecules in order to drastically reduce the number of terms in the sum (15) as well as the number of basis functions, with a small loss of precision. We propose here to apply the same strategy for nonrigid molecules, with illustrative examples given in Sec. IV. Very briefly, as follows:

*Reduced Hamiltonian*: the

*p*-order polynomial expansion (15) expressed in terms of (

*Q*

_{i},

*P*

_{i},

*J*

_{ρ},

*J*

_{α}) is first transformed into a second-quantized form in terms of $(ai+,ai,J\rho ,J\alpha )$. Then, we define a new order

*p*′ <

*p*, and we discard during the second-quantized transformation all the terms of degree $>p\u2032$ in

*a*

^{+}and

*a*. In a final step, a backward transformation of the low-order second-quantized Hamiltonian allows us to define a so-called reduced Hamiltonian

*p*-order polynomial is transformed to a

*p*′-order polynomial. The advantage of this method is twofold: it allows reducing the number of terms by at least one order of magnitude, and it generally allows getting rid of artifacts and spurious minima inherent to the high-order Taylor expansion of multivariate functions.

*Reduced eigenvectors*: The *J* = 0 problem is first solved by using a suitable vibrational basis *F*(*m*) [see Eq. (32)] to properly converge the desired eigenvalues. Then, we select relevant vibrational eigenvectors; we define so-called reduced eigenvectors $\Psi red(m\u2192m\u2032)$, simply denoted as the *F*(*m* → *m*′) reduction with *m*′ < *m*, and we compute the rovibrational solutions by following the procedure described by Eqs. (7)–(11) of Ref. 90.

## IV. APPLICATIONS

In order to validate our HBJ-based hybrid model, we consider in this section four molecules with one “floppy” large amplitude vibration: the quasilinear CH_{2} methylene molecule, the planar CH_{3} methyl radical, the NH_{3} ammonia molecule with inversion, and the H_{2}O_{2} (hydrogen peroxide) molecule with internal rotation. It is worth mentioning that the proposed hybrid model was recently used to validate our new PESs,^{131–133} but it was described in a very sparse manner. In this work, the theory as well as the methodology illustrated with practical applications are presented in more detail.

Beyond the necessary validation step, this section also aims at carefully studying the impact of the Hamiltonian and basis set truncation on the energy level calculation. Recent papers^{89,134} studied the dependence of rovibrational calculations with respect to the truncation order of the Watson Hamiltonian. A similar convergence study is proposed here for nonrigid molecules and turns out to be very relevant for

controlling carefully the precision of the calculated eigenvalues and eigenvectors and,

constructing a tailor-made model suitable either for the analysis of high-resolution spectra or for the modeling of some planetary atmospheres requiring less accuracy.

*κ*

_{ρ}and

*κ*

_{i}are some weight coefficients. Combined with symmetry, this condition can drastically reduce the dimension of the Hamiltonian matrix to be diagonalized by several orders of magnitude. Note that in all the calculations below, the rigid part of the Hamiltonian (17) was rewritten in terms of dimensionless normal coordinates and conjugate momenta (

*q*,

*p*) using the standard relation

*q*=

*γ*

^{1/2}

*Q*.

### A. CH_{2} molecule

In the past four decades, published papers have been dedicated to the energy level calculation and line transition prediction of methylene in its ground electronic triplet state.^{136–140} Due to its large bending motion and low barrier to linearity around 2000 cm^{−1}, methylene turns out to be a good first candidate to validate our hybrid model. A recent accurate *ab initio* PES for CH_{2} has been calculated,^{133} including both relativistic and diagonal Born–Oppenheimer corrections as well as high-order electronic correlations. This PES will be naturally used in the present paper for various validation tests.

Our grid is composed of 100 quadrature Gauss–Legendre points, and the hybrid Hamiltonian (17) was computed from *ρ*_{min} = 0 to *ρ*_{max} = 3.1 radians. The construction of the symmetrized powers for the two stretching coordinates is trivial when dealing with one-dimensional irreps. Indeed, for a symmetry species Γ, the symmetrized powers will be just Γ_{0} (resp. Γ) for even (resp. odd) powers. The standard harmonic oscillator basis functions (HOBF) were used for the stretching modes, while the displaced HOBF $NvHv(\rho \u2032)e\u2212\rho \u20322/2$, with *ρ*′ = 5(*ρ* − 0.7), multiplied by a weight function $sin(\rho )$, was chosen for the bending mode to properly treat the overall range of *ρ*, even near *ρ* = 0. These functions are ortho-normalized using the standard Gram–Schmidt technique.

As a first test, the convergence of the KEO and potential developed in terms of *q*_{1} and *q*_{3} was studied. To this end, we have compared in Fig. 3 (left panel) the *J* = 0 energy levels obtained from the Hamiltonian (17) truncated at orders 10, 14, and 18 with respect to those obtained from the Hamiltonian truncated at order 20, taken as the benchmark calculation. The 10th, 14th, 18th, and 20th order Hamiltonians contain 378, 778, 1332, and 1648 ITOs, respectively. The *F*(27) basis set defined in (32) and composed of 2135 *A*_{1} functions and 1925 *B*_{2} functions was used to compute the vibrational energy levels, which are converged within 10^{−4} up to 10 000 cm^{−1} using the Hamiltonian truncated at order 18. In Fig. 3 (right panel), we show that our *J* = 0 and *J* = 10 levels are in very good agreement with those obtained from the DVR3D computer code.^{135} Indeed, using the 20th order Hamiltonian and the *F*(27) basis, the error on the *J* = 0 and *J* = 10 levels between DVR and our variational calculation is 10^{−5} up to 10 000 cm^{−1}. For *J* = 10, the basis contains about 22 000 functions but can be drastically reduced by choosing conveniently the weights *κ*_{i} in (32). For example, if *κ*_{1} = *κ*_{3} = 1.2, the number of basis functions is decreased by 35%, and the rms error becomes 0.0002 up to 10 000 cm^{−1}. By taking 50% fewer functions (*κ*_{1} = *κ*_{3} = 1.4), the convergence of the levels is not fully achieved (rms error of 0.001 cm^{−1}), but the results remain suitable for high-resolution spectroscopy.

### B. CH_{3} molecule

*D*

_{3h}(

*M*)) methyl radical CH

_{3}for which the out-of-plane

*ν*

_{2}mode exhibits a LAM behavior. Recently, new

*ab initio*CH

_{3}PES and DMS in their ground electronic states have been published

^{132}and are thus considered in this work. The hybrid Hamiltonian and dipole moment operator have been computed here on a numerical grid from

*ρ*

_{min}= 0.41 to

*ρ*

_{max}= 2.69 radians. The construction of the vibrational operators for the doubly degenerate modes

*ν*

_{3}and

*ν*

_{4}of symmetry

*E*′ in

*D*

_{3h}(

*M*) closely follows the method presented in Ref. 115 for the

*C*

_{3v}and

*T*

_{d}point groups. For twofold degenerate oscillators, it is common to introduce the auxiliary boson operators

^{141}

^{,}

*b*

_{1}and

*b*

_{2}and form the products $b1lb2m\u2261(lm)$, like in the

*u*(2) formalism.

^{121}The same holds for the construction of the creation part with $b1+$ and $b2+$. It is quite straightforward to see how these operators transform using the generators of

*D*

_{3h}(

*M*) and to show that the two quantities

*j*=

*l*−

*m*, we have $\Gamma 1=A1\u2032$ when

*l*=

*m*, $\Gamma 1=A1\u2032$, $\Gamma 2=A2\u2032$ for

*j*= 3

*p*(

*p*≠ 0), and (Γ

_{1},

*σ*

_{1}) = (

*E*′,

*a*), (Γ

_{2},

*σ*

_{2}) = (

*E*′,

*b*) for

*j*= 3

*p*+ 1 or 3

*p*+ 2, with

*p*an integer number. All the factors

*e*

^{iϕ(l,m)}equal 1, except for

*j*= 3

*p*+ 2, where it equals −1. For a mode of symmetry

*E*″, the symmetry rules will slightly change. In the ITO formalism, the vibrational reduced matrix elements involved in (30) take a quite simple form when using the creation-annihilation operators involved in effective Hamiltonians. For this work, we have shown

^{127}that the reduced matrix elements of symmetrized powers for doubly degenerate modes can also be put in a closed form using the (

*q*,

*p*) formalism. Concerning the rotational part, the symmetries Γ

_{r}involved in (29) will differ according to the values of

*k*; that is, we will have $\Gamma r=A1\u2032+A2\u2032$ for

*k*= 6

*p*(

*p*≠ 0), $\Gamma r=A1\u2033+A2\u2033$ for

*k*= 6

*p*+ 3, Γ

_{r}=

*E*′ for

*k*= 6

*p*+ 2 or

*k*= 6

*p*+ 4, and Γ

_{r}=

*E*″ for

*k*= 6

*p*+ 1 or

*k*= 6

*p*+ 5. Concerning the bending mode

*ν*

_{2}, simple HOBF $NvHv(\rho \u2032)e\u2212\rho \u20322/2$, with

*ρ*′ = 8.1

*ρ*, of symmetry $A1\u2032$ for

*v*even and $A2\u2033$ for

*v*odd, were used.

As a first test, the number of quadrature points *M*_{ρ} was gradually increased from 50 to 150 in order to see the impact on the energy level calculation. A variational calculation was performed for different *M*_{ρ} using the *F*(23) basis with *κ*_{1} = *κ*_{2} = 1.4 and *κ*_{3} = 1.3, resulting in 11 010, 9857, 20 849, 7869, 8861, and 16 719 basis functions of symmetry $A1\u2032$, $A2\u2032$, *E*′, $A1\u2033$, $A2\u2033$, and *E*″, respectively. The errors between the *J* = 0 levels computed from *M*_{ρ} = 50, 80, and 150 and those computed from *M*_{ρ} = 100, taken as the benchmark, are depicted in Fig. 4 (top left panel). We clearly see that *M*_{ρ} = 50 is not accurate enough for practical applications, while the results become very similar at *M*_{ρ} = 80.

In a second step, we studied the convergence of the Hamiltonian expansion. Like for CH_{2}, we have gradually increased the order of the expansion (17) and computed the vibrational levels for a given basis set [here, *F*(23)]. The errors between the energy levels computed from the Hamiltonian truncated at orders 6, 8, and 11 and those computed from the Hamiltonian truncated at order 12 are given in Fig. 4 (top right panel). We can see quite large fluctuations beyond 5000 cm^{−1} between orders 12 and 11, proving that this latter was not yet properly converged to study a highly excited (ro)vibrational state.

In a third step, the convergence of the variational calculation has been studied using different basis sets, namely *F*(12), *F*(14), *F*(17), *F*(20), and *F*(23). The number of vibrational basis functions of symmetry *E*′ increases to 525, 1071, 2556, 5536, and 11 010, respectively. It is seen in Fig. 4 (bottom left panel) that the convergence error of the vibrational levels up to 5000 cm^{−1} is about 10^{−3} cm^{−1} using the *F*(20) basis. The use of the *F*(23) basis allows for converging the levels within 10^{−4} up to 5000 cm^{−1} and within 10^{−3} up to 6000 cm^{−1}, which turns out to be enough for the modeling of high-resolution spectra. Finally, a room-temperature CH_{3} line list up to 6000 cm^{−1}, including both cold and hot bands (see Fig. 4, bottom right panel), was constructed using the DMS presented in Ref. 132. Comparisons with the results obtained in Ref. 142 can also be found in Ref. 132.

### C. NH_{3} molecule

The ammonia molecule is a good candidate to be tested using the hybrid model because its umbrella inversion is among the most prototypical cases of large amplitude motion. In the past two decades, many have been were devoted to the construction of refined PESs and published in the ExoMol group for the construction of molecular line lists^{143} combined with a MARVELization procedure.^{144} In this work, we have considered our recently published “pure” (i.e., in the sense without applying empirical corrections) *ab initio* PES,^{131} which turns out to be accurate up to 6600 cm^{−1} with a root-mean-square error of 0.4 cm^{−1} between the calculated and observed levels (see Table 3 of Ref. 131).

For this molecule, a grid *ρ* ∈ [0.37, 2.83] radians was chosen to ensure a consistent energy calculation up to at least 6000 cm^{−1}. The construction of ITOs adapted to both small amplitude vibrations and rotation is similar to that of the CH_{3} molecule. The major difference lies in the treatment of the LAM, in particular in the choice of the basis functions adapted to the symmetric double-well potential. As usual, the two wells correspond to the two equivalent pyramidal *C*_{3v}(*M*) equilibrium configurations. However, *C*_{3v}(*M*) is no longer suitable to describe tunneling because no operation in this group allows switching from one configuration to another. Therefore, for a proper classification of the overall vibration–inversion–rotation energy levels, we need to consider the larger permutation-inversion group *G*_{12}, which is isomorphic to the MS group *D*_{3h}(*M*) = *C*_{3v}(*M*) × {*E*, *E**} and consists of eight feasible operations. In the case of ammonia, the tunneling produces a splitting of the rigid-molecule *C*_{3v}(*M*) states into two components, which are deduced from the reverse correlation (or induction) table between *D*_{3h}(*M*) and *C*_{3v}(*M*). For example, the *A*_{1} states split into $A1\u2032$ and $A2\u2033$ components; the *A*_{2} states split into $A1\u2033$ and $A2\u2032$ components; and the *E* states split into *E*′ and *E*″ components. The question of how the energy levels of the infinite and low barrier limits are correlated was first raised by Watson.^{68}

*ν*

_{2}. In order to make the analogy with the field of quantum optics, the basis functions used here are nothing but a linear combination of displaced squeezed states (see, e.g., Ref. 145 and references within) of the type

*λ*

_{j}) and squeezing (

*κ*

_{j}) parameters play the role of variational parameters to be optimized. Here, the functions

*ϕ*

_{v}(

*ρ*) are simply HOBF, but Morse-like or Pöschl–Teller-like functions could also be used instead. The matrix elements

*D*

_{vv′}and

*S*

_{vv′}of the displaced and squeeze operators can be computed analytically.

^{146}

Like for CH_{2} and CH_{3}, the convergence of the Hamiltonian model is first studied and shown in Fig. 5 (left panel). To this end, we have computed the energy levels for the $A1\u2032$ and $A2\u2033$ blocks composed of 17 739 functions by using the *F*(17) basis and by gradually increasing the Hamiltonian truncation order from 6 to 12. For the inversion motion, the new coordinate *ρ*′ = *ρ* − *π*/2 was introduced for convenience so that the planar configuration is now defined at *ρ*′ = 0. In this new coordinate system, the parameters of the basis function (34) have been defined as follows: *κ*_{1} = *κ*_{2} = 9.45 and *λ*_{1} = −*λ*_{2} = 0.338. We can thus estimate the precision of the Hamiltonian truncated at order 12 below 0.01 up to 6000 cm^{−1}. Starting from this Hamiltonian model, we have then examined the convergence of the variational calculation by gradually increasing the number of basis functions. The results are displayed in Fig. 5 (right panel). We can see that the energy levels using the *F*(16) basis, which is composed of 13 185, 11 700, and 24 864 $A1\u2032/A2\u2033$, $A1\u2033/A2\u2032$, and *E*′/*E*″ functions, respectively, are converged within 0.001 cm^{−1}, with the only exception for the level near 5356.1 cm^{−1} corresponding to $5\nu 2(A2\u2033)$. In that case, the use of the *F*(17) basis composed of 17 739, 15 924, and 33 642 $A1\u2032/A2\u2033$, $A1\u2033/A2\u2032$, and *E*′/*E*″ functions is required to properly converge the five-quanta vibrational bands. Needless to say, converging levels above 6000 within 10^{−3} cm^{−1} will probably require more basis functions, in particular for the inversion mode.

For *J* > 0 calculations, the *reduced basis functions* *F*(17) → *F*(*p*) initially introduced for the rovibrational energy level calculation of semirigid molecules (see Sec. III C 5) have been tested for the NH_{3} molecule for *J* = 10. The four weights *κ*_{i} associated with the *F*(*p*) basis used here are 1.4, 1, 1.4, and 1.3. In order to achieve a good convergence of the energy levels for the block $(J=10,A1\u2032)$, different values of *p* have been considered, namely *p* = 10, 12, 13, and 14. The dimensions of the corresponding blocks are 7278, 16 474, 23 555, and 34 111, respectively, while those without using reduced functions are 354 222. Therefore, such functions allow us to reduce the dimension of the problem by at least one order of magnitude and thus use direct eigensolvers for the diagonalization. As usual, the biggest calculation corresponding to *p* = 14 will be taken as the benchmark.

The convergence of the levels with respect to *p* is depicted in Fig. 6 (left panel). We can see that the average error between *p* = 13 and *p* = 14 is 0.001 up to 6000 cm^{−1}, which is the typical resolution of Fourier transform infrared spectra. Even the smaller basis (*p* = 10) can produce qualitatively correct results, with errors around 0.01 cm^{−1} or less up to 4000 cm^{−1}. We can conclude that the reduced basis functions can also be used for nonrigid molecules to manage memory during high *J* calculations. In Fig. 6 (right panel), we have plotted the reduced rovibrational energy levels up to *J* = 18 using *p* = 12. The different colors represent the mixing coefficients in the eigenvector decomposition.

### D. H_{2}O_{2} molecule

Hydrogen peroxide is the simplest molecule that undergoes an internal rotation motion. It consists of a torsional motion of the O–H groups around the O–O bond, with a height of the *trans* barrier (*ρ* = *π*, symmetry point group *C*_{2h}) around 400 cm^{−1} and of the *cis* barriers (*ρ* = 0, symmetry point group *C*_{2v}) around 2600 cm^{−1}. Recently, a variational calculation of the energy levels of H_{2}O_{2} was presented in Ref. 147, based on the PES of Malyszek and Koput.^{148} It is common to choose the torsion angle between the two “halves” as 2*ρ*, that is, a rotation of +*ρ* for the first half and of −*ρ* for the second half (see, e.g., Refs. 94 and 149), whereas a torsional angle of *ρ*, with rotations of ±*ρ*/2, was considered in Ref. 150 and in this work. It is important to note that, contrary to the three previous molecules, H_{2}O_{2} will pass through *cis* and *trans* configurations with different symmetries as *ρ* increases from 0 to *π*, with the point group for intermediate configurations being *C*_{2}.

If only one barrier is not “insuperable,” for example, the *trans* barrier, then the permutation-inversion “four-group” *G*_{4} must be used to classify the overall vibration–rotation–torsion energy levels. This group consists of four feasible operations and is isomorphic to the point groups *C*_{2v} and *C*_{2h}. However, if the energy splitting due to both the *cis* and *trans* tunnelings can be resolved in the observed spectra, as in the case of H_{2}O_{2}, Hougen showed that this group must be extended to *G*_{4}(*EM*) in order to properly describe separately the torsional and rotational states. This group is isomorphic to the MS group *D*_{2h}(*M*) and to the point group *D*_{2h}, which is already implemented in the TENSOR code. Accordingly, the symmetry species of *D*_{2h} can be used to classify the energy levels. If we refer now to the reverse correlation table between *D*_{2h} and *C*_{2}, the rigid-molecule *A* states of *C*_{2} split into *A*_{1g}, *B*_{2g}, *B*_{2u}, and *A*_{1u} components. Alternatively, a quantum number *τ* taking the values 1, 2, 3, and 4, is generally used to label these four components.^{94,149}

For hydrogen peroxide, our energy levels can be directly compared to those computed from the WAVR4 code by Kozin *et al.*^{151} To this end, the *ab initio* PES parameters determined by Malyszek and Koput^{148} were used to generate a grid of points that has been fitted using a set of symmetry coordinates in order to be compatible with both TENSOR and WAVR4. Trial calculations were first made using WAVR4 to ensure that the newly fitted PES has no hole. The Fortran code as well as the new fitted PES parameters are provided in supplementary material. It is worth mentioning at this stage that the aim of the present work was not to reproduce the observed data but to validate the convergence of our Taylor-type expanded hybrid model with respect to the exact KEO model implemented in the WAVR4 computer code. As already shown in Ref. 147, the accuracy of the vibrational band origins using the *ab initio* PES of Malyszek and Koput^{148} is within 1 cm^{−1} but can be drastically improved using a slightly adjusted PES.

In order to achieve good convergence, we have studied the dependence of the vibrational levels for the first symmetry block *A*_{g} with respect to the truncation order of (i) the kinetic energy operator (8), (ii) the potential (7), (iii) the reduced Hamiltonian (31), and (iv) the basis set (32). For the convergence studies (i)–(iii), the basis *F*(18) was employed using the weight coefficients *κ*_{1} = *κ*_{5} = 1.2 and *κ*_{2} = *κ*_{6} = 1.1, *κ*_{3} = *κ*_{ρ} = 1.

*Convergence study* (i). The convergence of the vibrational levels computed with a KEO truncated at orders 6, 8, and 10 with respect to that truncated at order 11 is shown in Fig. 7 (top, left panel). This amounts to studying the dependence of the energy levels with respect to the *μ*_{αβ} tensor truncation order. For these four calculations, the potential expansion was truncated at order 16. We can see that the levels are converged within 0.001 cm^{−1} (and even below) up to 4000 cm^{−1} using a KEO truncated at order 10.

*Convergence study* (ii). In Fig. 7 (top, right panel), we give the convergence of the vibrational levels computed with the potential truncated at orders 10, 12, and 14 with respect to that truncated at order 16 and by using the 10th order KEO expansion. Similarly to the Watson Hamiltonian,^{89} we can conclude that the convergence of *μ*_{αβ} is faster than the convergence of *V* in the HBJ Hamiltonian.

*Convergence study* (iii). We have also tested the convergence of the reduced (31) and non-reduced (15) Hamiltonian expansions with respect to the non-reduced 16th order Hamiltonian *H*(16), consisting of a potential truncated at order 16 and a KEO truncated at order 10. We can see in Fig. 7 (bottom, left panel) that the reduced Hamiltonian *H*(16 → 8) is more accurate than the Taylor-based expansion *H*(10), using fewer parameters. The reduced Hamiltonian *H*(16 → 10) is converged within 0.001 up to 3000 cm^{−1} and within 0.01 up to 4000 cm^{−1} with respect to the full *H*(16) expansion.

*Convergence study* (iv). From the Hamiltonian *H*(16), the convergence of the variational calculation using the basis sets *F*(14), *F*(16), *F*(18), and *F*(19) with respect to the basis *F*(20) is displayed in Fig. 7 (bottom, right panel). These basis sets contain 10 871, 20 961, 37 820, 49 780, and 64 721 vibrational functions of symmetry *A*_{g}, respectively. We can see that the *F*(16) basis would be suitable up to 4000 cm^{−1} for the modeling of low and medium resolution infrared spectra, while the *F*(19) basis allows converging vibrational levels within 0.001 up to 3500 cm^{−1}.

Finally, we give in Tables I and II the differences between the vibrational energy levels computed from the TENSOR and WAVR4^{151} computer codes for the symmetry blocks (*A*_{g}, *B*_{1u}) and (*A*_{u}, *B*_{1g}), respectively. We can see that the results obtained from WAVR4 are quite sensitive to the input parameters (*j*_{max}, *n*, *icut*) (see Ref. 151 for more explanations), in particular for converging the overtones of the torsional mode *ν*_{4}. Several days of calculation were necessary to compute the *J* = 0 levels using WAVR4 for (*j*_{max}, *n*, *icut*) = (25, 14, 300) against a few hours using TENSOR. Finally, we can note the very good agreement between both calculations, with errors ∼0.001 up to 3000 cm^{−1}, which is much better than the agreement between TROVE and WAVR4 in Ref. 147, with deviations up to 0.5 cm^{−1}.

Band . | Sym . | (i) . | (ii) . | (iii) . | (iv) . | (v) . | (vi) . | (vii) . |
---|---|---|---|---|---|---|---|---|

ν_{4} | A_{g} | 255.7882 | 0.0025 | 0.0019 | 0.0018 | 0.0008 | 0.0002 | 0.0001 |

2ν_{4} | A_{g} | 570.4892 | 0.0069 | 0.0043 | 0.0041 | 0.0025 | −0.0002 | −0.0004 |

ν_{3} | A_{g} | 865.9407 | 0.0056 | 0.0044 | 0.0044 | 0.0013 | 0.0002 | 0.0001 |

3ν_{4} | A_{g} | 1001.7398 | 0.0087 | 0.0047 | 0.0045 | 0.0034 | −0.0007 | −0.0008 |

ν_{3} + ν_{4} | A_{g} | 1120.4606 | 0.0075 | 0.0058 | 0.0051 | 0.0026 | 0.0008 | 0.0002 |

ν_{6} | B_{1u} | 1264.3401 | 0.1387 | 0.1266 | 0.1263 | 0.0132 | 0.0003 | 0.0001 |

ν_{2} | A_{g} | 1394.5277 | 0.1271 | 0.1154 | 0.1152 | 0.0128 | 0.0003 | 0.0001 |

ν_{3} + 2ν_{4} | A_{g} | 1431.7821 | 0.0199 | 0.0148 | 0.0140 | 0.0057 | 0.0003 | −0.0005 |

4ν_{4} | A_{g} | 1478.6183 | 0.0320 | 0.0241 | 0.0240 | 0.0076 | −0.0010 | −0.0011 |

ν_{4} + ν_{6} | B_{1u} | 1505.3746 | 0.1540 | 0.1354 | 0.1351 | 0.0203 | 0.0003 | 0.0001 |

ν_{2} + ν_{4} | A_{g} | 1683.5634 | 0.1592 | 0.1349 | 0.1345 | 0.0269 | 0.0006 | 0.0002 |

2ν_{3} | A_{g} | 1714.5097 | −0.0157 | −0.0182 | −0.0187 | 0.0029 | 0.0005 | 0.0001 |

2ν_{4} + ν_{6} | B_{1u} | 1854.0556 | 0.2180 | 0.1825 | 0.1821 | 0.0364 | 0.0002 | −0.0003 |

ν_{3} + 3ν_{4} | A_{g} | 1861.6873 | 0.0304 | 0.0214 | 0.0203 | 0.0093 | −0.0002 | −0.0011 |

5ν_{4} | A_{g} | 1945.9637 | 0.0574 | 0.0299 | 0.0296 | 0.0282 | −0.0011 | −0.0014 |

2ν_{3} + ν_{4} | A_{g} | 1967.5418 | −0.0269 | −0.0317 | −0.0329 | 0.0055 | 0.0009 | −0.0002 |

ν_{2} + 2ν_{4} | A_{g} | 1976.1079 | 0.2358 | 0.1855 | 0.1848 | 0.0561 | 0.0007 | −0.0001 |

ν_{3} + ν_{6} | B_{1u} | 2113.7491 | 0.2389 | 0.2156 | 0.2148 | 0.0263 | 0.0019 | 0.0011 |

ν_{2} + ν_{3} | A_{g} | 2241.9684 | 0.2514 | 0.2236 | 0.2225 | 0.0317 | 0.0022 | 0.0009 |

2ν_{3} + 2ν_{4} | A_{g} | 2277.3907 | 0.0012 | −0.0143 | −0.0157 | 0.0150 | −0.0005 | −0.0018 |

3ν_{4} + ν_{6} | B_{1u} | 2306.8829 | 0.3843 | 0.3169 | 0.3139 | 0.0733 | 0.0028 | −0.0007 |

6ν_{4} | A_{g} | 2316.7995 | 0.0758 | 0.0315 | 0.0302 | 0.0463 | −0.0003 | −0.0017 |

ν_{3} + ν_{4} + ν_{6} | B_{1u} | 2354.2088 | 0.2863 | 0.2430 | 0.2421 | 0.0465 | 0.0016 | 0.0007 |

ν_{3} + 4ν_{4} | A_{g} | 2355.9264 | 0.0609 | 0.0266 | 0.0254 | 0.0354 | −0.0006 | −0.0018 |

ν_{2} + 3ν_{4} | A_{g} | 2392.7082 | 0.3737 | 0.2860 | 0.2817 | 0.0958 | 0.0043 | −0.0006 |

2ν_{6} | A_{g} | 2505.3506 | 1.0302 | 0.9121 | 0.9043 | 0.1386 | 0.0119 | 0.0030 |

ν_{2} + ν_{3} + ν_{4} | A_{g} | 2530.0509 | 0.3083 | 0.2195 | 0.2176 | 0.0989 | 0.0028 | 0.0009 |

3ν_{3} | A_{g} | 2545.5883 | −0.6603 | −0.6702 | −0.6718 | 0.0030 | −0.0021 | −0.0027 |

ν_{2} + ν_{6} | B_{1u} | 2648.8057 | 1.8124 | 1.6103 | 1.5977 | 0.2427 | 0.0202 | 0.0061 |

ν_{3} + 2ν_{4} + ν_{6} | B_{1u} | 2701.2782 | 0.3994 | 0.2994 | 0.2979 | 0.1012 | 0.0008 | −0.0008 |

2ν_{3} + 3ν_{4} | A_{g} | 2705.2429 | 0.0315 | −0.0012 | −0.0030 | 0.0301 | −0.0037 | −0.0054 |

7ν_{4} | A_{g} | 2720.3876 | 0.1485 | 0.0579 | 0.0549 | 0.0964 | −0.0011 | −0.0043 |

ν_{4} + 2ν_{6} | A_{g} | 2743.9942 | 1.4453 | 1.1785 | 1.1665 | 0.3147 | 0.0169 | 0.0040 |

2ν_{2} | A_{g} | 2765.5520 | 1.0040 | 0.7899 | 0.7788 | 0.2444 | 0.0153 | 0.0028 |

3ν_{3} + ν_{4} | A_{g} | 2796.6891 | −0.6625 | −0.6821 | −0.6856 | 0.0051 | −0.0066 | −0.0082 |

4ν_{4} + ν_{6} | B_{1u} | 2799.9239 | 0.5421 | 0.2867 | 0.2796 | 0.3039 | 0.0054 | −0.0022 |

ν_{3} + 5ν_{4} | A_{g} | 2809.5394 | 0.2537 | 0.1122 | 0.1106 | 0.1563 | −0.0011 | −0.0029 |

ν_{2} + ν_{3} + 2ν_{4} | A_{g} | 2824.6633 | 0.3854 | 0.1918 | 0.1894 | 0.2082 | 0.0012 | −0.0016 |

ν_{2} + 4ν_{4} | A_{g} | 2861.3768 | 0.6177 | 0.2763 | 0.2684 | 0.3969 | 0.0066 | −0.0017 |

ν_{2} + ν_{4} + ν_{6} | B_{1u} | 2914.7957 | 2.6087 | 1.8111 | 1.7898 | 0.8872 | 0.0279 | 0.0072 |

2ν_{3} + ν_{6} | B_{1u} | 2945.3719 | 0.2881 | 0.2364 | 0.2341 | 0.0554 | 0.0018 | 0.0000 |

rms error | 0.6452 | 0.5064 | 0.5018 | 0.1851 | 0.0071 | 0.0026 |

Band . | Sym . | (i) . | (ii) . | (iii) . | (iv) . | (v) . | (vi) . | (vii) . |
---|---|---|---|---|---|---|---|---|

ν_{4} | A_{g} | 255.7882 | 0.0025 | 0.0019 | 0.0018 | 0.0008 | 0.0002 | 0.0001 |

2ν_{4} | A_{g} | 570.4892 | 0.0069 | 0.0043 | 0.0041 | 0.0025 | −0.0002 | −0.0004 |

ν_{3} | A_{g} | 865.9407 | 0.0056 | 0.0044 | 0.0044 | 0.0013 | 0.0002 | 0.0001 |

3ν_{4} | A_{g} | 1001.7398 | 0.0087 | 0.0047 | 0.0045 | 0.0034 | −0.0007 | −0.0008 |

ν_{3} + ν_{4} | A_{g} | 1120.4606 | 0.0075 | 0.0058 | 0.0051 | 0.0026 | 0.0008 | 0.0002 |

ν_{6} | B_{1u} | 1264.3401 | 0.1387 | 0.1266 | 0.1263 | 0.0132 | 0.0003 | 0.0001 |

ν_{2} | A_{g} | 1394.5277 | 0.1271 | 0.1154 | 0.1152 | 0.0128 | 0.0003 | 0.0001 |

ν_{3} + 2ν_{4} | A_{g} | 1431.7821 | 0.0199 | 0.0148 | 0.0140 | 0.0057 | 0.0003 | −0.0005 |

4ν_{4} | A_{g} | 1478.6183 | 0.0320 | 0.0241 | 0.0240 | 0.0076 | −0.0010 | −0.0011 |

ν_{4} + ν_{6} | B_{1u} | 1505.3746 | 0.1540 | 0.1354 | 0.1351 | 0.0203 | 0.0003 | 0.0001 |

ν_{2} + ν_{4} | A_{g} | 1683.5634 | 0.1592 | 0.1349 | 0.1345 | 0.0269 | 0.0006 | 0.0002 |

2ν_{3} | A_{g} | 1714.5097 | −0.0157 | −0.0182 | −0.0187 | 0.0029 | 0.0005 | 0.0001 |

2ν_{4} + ν_{6} | B_{1u} | 1854.0556 | 0.2180 | 0.1825 | 0.1821 | 0.0364 | 0.0002 | −0.0003 |

ν_{3} + 3ν_{4} | A_{g} | 1861.6873 | 0.0304 | 0.0214 | 0.0203 | 0.0093 | −0.0002 | −0.0011 |

5ν_{4} | A_{g} | 1945.9637 | 0.0574 | 0.0299 | 0.0296 | 0.0282 | −0.0011 | −0.0014 |

2ν_{3} + ν_{4} | A_{g} | 1967.5418 | −0.0269 | −0.0317 | −0.0329 | 0.0055 | 0.0009 | −0.0002 |

ν_{2} + 2ν_{4} | A_{g} | 1976.1079 | 0.2358 | 0.1855 | 0.1848 | 0.0561 | 0.0007 | −0.0001 |

ν_{3} + ν_{6} | B_{1u} | 2113.7491 | 0.2389 | 0.2156 | 0.2148 | 0.0263 | 0.0019 | 0.0011 |

ν_{2} + ν_{3} | A_{g} | 2241.9684 | 0.2514 | 0.2236 | 0.2225 | 0.0317 | 0.0022 | 0.0009 |

2ν_{3} + 2ν_{4} | A_{g} | 2277.3907 | 0.0012 | −0.0143 | −0.0157 | 0.0150 | −0.0005 | −0.0018 |

3ν_{4} + ν_{6} | B_{1u} | 2306.8829 | 0.3843 | 0.3169 | 0.3139 | 0.0733 | 0.0028 | −0.0007 |

6ν_{4} | A_{g} | 2316.7995 | 0.0758 | 0.0315 | 0.0302 | 0.0463 | −0.0003 | −0.0017 |

ν_{3} + ν_{4} + ν_{6} | B_{1u} | 2354.2088 | 0.2863 | 0.2430 | 0.2421 | 0.0465 | 0.0016 | 0.0007 |

ν_{3} + 4ν_{4} | A_{g} | 2355.9264 | 0.0609 | 0.0266 | 0.0254 | 0.0354 | −0.0006 | −0.0018 |

ν_{2} + 3ν_{4} | A_{g} | 2392.7082 | 0.3737 | 0.2860 | 0.2817 | 0.0958 | 0.0043 | −0.0006 |

2ν_{6} | A_{g} | 2505.3506 | 1.0302 | 0.9121 | 0.9043 | 0.1386 | 0.0119 | 0.0030 |

ν_{2} + ν_{3} + ν_{4} | A_{g} | 2530.0509 | 0.3083 | 0.2195 | 0.2176 | 0.0989 | 0.0028 | 0.0009 |

3ν_{3} | A_{g} | 2545.5883 | −0.6603 | −0.6702 | −0.6718 | 0.0030 | −0.0021 | −0.0027 |

ν_{2} + ν_{6} | B_{1u} | 2648.8057 | 1.8124 | 1.6103 | 1.5977 | 0.2427 | 0.0202 | 0.0061 |

ν_{3} + 2ν_{4} + ν_{6} | B_{1u} | 2701.2782 | 0.3994 | 0.2994 | 0.2979 | 0.1012 | 0.0008 | −0.0008 |

2ν_{3} + 3ν_{4} | A_{g} | 2705.2429 | 0.0315 | −0.0012 | −0.0030 | 0.0301 | −0.0037 | −0.0054 |

7ν_{4} | A_{g} | 2720.3876 | 0.1485 | 0.0579 | 0.0549 | 0.0964 | −0.0011 | −0.0043 |

ν_{4} + 2ν_{6} | A_{g} | 2743.9942 | 1.4453 | 1.1785 | 1.1665 | 0.3147 | 0.0169 | 0.0040 |

2ν_{2} | A_{g} | 2765.5520 | 1.0040 | 0.7899 | 0.7788 | 0.2444 | 0.0153 | 0.0028 |

3ν_{3} + ν_{4} | A_{g} | 2796.6891 | −0.6625 | −0.6821 | −0.6856 | 0.0051 | −0.0066 | −0.0082 |

4ν_{4} + ν_{6} | B_{1u} | 2799.9239 | 0.5421 | 0.2867 | 0.2796 | 0.3039 | 0.0054 | −0.0022 |

ν_{3} + 5ν_{4} | A_{g} | 2809.5394 | 0.2537 | 0.1122 | 0.1106 | 0.1563 | −0.0011 | −0.0029 |

ν_{2} + ν_{3} + 2ν_{4} | A_{g} | 2824.6633 | 0.3854 | 0.1918 | 0.1894 | 0.2082 | 0.0012 | −0.0016 |

ν_{2} + 4ν_{4} | A_{g} | 2861.3768 | 0.6177 | 0.2763 | 0.2684 | 0.3969 | 0.0066 | −0.0017 |

ν_{2} + ν_{4} + ν_{6} | B_{1u} | 2914.7957 | 2.6087 | 1.8111 | 1.7898 | 0.8872 | 0.0279 | 0.0072 |

2ν_{3} + ν_{6} | B_{1u} | 2945.3719 | 0.2881 | 0.2364 | 0.2341 | 0.0554 | 0.0018 | 0.0000 |

rms error | 0.6452 | 0.5064 | 0.5018 | 0.1851 | 0.0071 | 0.0026 |

Band . | Sym . | (i) . | (ii) . | (iii) . | (iv) . | (v) . | (vi) . | (vii) . |
---|---|---|---|---|---|---|---|---|

GS | A_{u} | 11.2122 | −0.0034 | −0.0034 | 0.0033 | −0.0003 | 0.0000 | −0.0001 |

ν_{4} | A_{u} | 371.3261 | −0.0014 | −0.0015 | 0.0015 | −0.0001 | −0.0001 | −0.0002 |

2ν_{4} | A_{u} | 776.8274 | 0.0119 | 0.0103 | −0.0101 | 0.0014 | −0.0018 | −0.0006 |

ν_{3} | A_{u} | 877.6885 | 0.0087 | 0.0077 | −0.0078 | 0.0010 | −0.0011 | 0.0000 |

ν_{3} + ν_{4} | A_{u} | 1227.0708 | 0.0134 | 0.0112 | −0.0109 | 0.0020 | −0.0022 | −0.0005 |

3ν_{4} | A_{u} | 1245.4558 | 0.0155 | 0.0131 | −0.0128 | 0.0023 | −0.0025 | −0.0006 |

ν_{6} | B_{1g} | 1284.6829 | 0.1640 | 0.1548 | −0.1546 | 0.0099 | −0.0097 | 0.0000 |

ν_{2} | A_{u} | 1400.7098 | 0.1674 | 0.1574 | −0.1571 | 0.0107 | −0.0104 | 0.0000 |

ν_{3} + 2ν_{4} | A_{u} | 1638.1640 | 0.0202 | 0.0170 | −0.0162 | 0.0039 | −0.0038 | −0.0007 |

ν_{4} + ν_{6} | B_{1g} | 1648.4799 | 0.2196 | 0.2066 | −0.2064 | 0.0138 | −0.0136 | 0.0001 |

4ν_{4} | A_{u} | 1720.1736 | 0.0330 | 0.0261 | −0.0260 | 0.0066 | −0.0077 | −0.0011 |

3ν_{3} | A_{u} | 1726.8764 | −0.0076 | −0.0101 | 0.0102 | 0.0025 | −0.0025 | 0.0000 |

ν_{2} + ν_{3} | A_{u} | 1772.8478 | 0.2064 | 0.1894 | −0.1891 | 0.0181 | −0.0178 | 0.0000 |

2ν_{4} + ν_{6} | B_{1g} | 2072.6916 | 0.2268 | 0.2033 | −0.2030 | 0.0269 | −0.0269 | −0.0004 |

2ν_{3} + ν_{4} | A_{u} | 2076.4227 | 0.0016 | −0.0028 | 0.0034 | 0.0044 | −0.0043 | −0.0005 |

ν_{3} + 3ν_{4} | A_{u} | 2103.3759 | 0.0231 | 0.0160 | −0.0152 | 0.0072 | −0.0075 | −0.0010 |

ν_{3} + ν_{6} | B_{1g} | 2134.6189 | 0.2206 | 0.2036 | −0.2025 | 0.0193 | −0.0178 | 0.0005 |

ν_{2} + 2ν_{4} | A_{u} | 2170.7709 | 0.2092 | 0.1735 | −0.1730 | 0.0400 | −0.0398 | −0.0004 |

5ν_{4} | A_{u} | 2205.8643 | 0.0567 | 0.0305 | −0.0303 | 0.0267 | −0.0280 | −0.0014 |

ν_{2} + ν_{3} | A_{u} | 2248.5455 | 0.2091 | 0.1873 | −0.1864 | 0.0246 | −0.0230 | 0.0008 |

2ν_{3} + 2ν_{4} | A_{u} | 2483.5870 | −0.0086 | −0.0154 | 0.0168 | 0.0074 | −0.0076 | −0.0015 |

ν_{3} + ν_{4} + ν_{6} | B_{1g} | 2496.0545 | 0.2672 | 0.2396 | −0.2389 | 0.0297 | −0.0282 | 0.0008 |

2ν_{6} | A_{u} | 2538.4081 | 0.9592 | 0.8801 | −0.8778 | 0.0921 | −0.0864 | 0.0027 |

3ν_{4} + ν_{6} | B_{1g} | 2551.8234 | 0.3660 | 0.3021 | −0.2999 | 0.0689 | −0.0670 | −0.0006 |

3ν_{3} | A_{u} | 2558.3483 | −0.6483 | −0.6569 | 0.6582 | 0.0035 | −0.0048 | −0.0015 |

ν_{3} + 4ν_{4} | A_{u} | 2579.2230 | 0.0632 | 0.0440 | −0.0435 | 0.0193 | −0.0201 | −0.0012 |

ν_{2} + ν_{3} + ν_{4} | A_{u} | 2614.6564 | 0.3396 | 0.2791 | −0.2773 | 0.0658 | −0.0633 | 0.0005 |

ν_{2} + 3ν_{4} | A_{u} | 2628.8437 | 0.3410 | 0.2729 | −0.2703 | 0.0741 | −0.0713 | −0.0001 |

ν_{2} + ν_{6} | B_{1g} | 2660.3104 | 1.4349 | 1.3261 | −1.3199 | 0.1318 | −0.1203 | 0.0045 |

6ν_{4} | A_{u} | 2705.4012 | 0.1297 | 0.0716 | −0.0697 | 0.0603 | −0.0604 | −0.0023 |

2ν_{2} | A_{u} | 2768.9855 | 0.8934 | 0.7426 | −0.7364 | 0.1715 | −0.1623 | 0.0022 |

3ν_{2} + ν_{4} | A_{u} | 2908.1590 | −0.5381 | −0.5535 | 0.5561 | 0.0049 | −0.0089 | −0.0050 |

ν_{4} + 2ν_{6} | A_{u} | 2910.8426 | 1.1266 | 0.9862 | −0.9787 | 0.1694 | −0.1578 | 0.0027 |

ν_{3} + 2ν_{4} + ν_{6} | B_{1g} | 2918.9996 | 0.4153 | 0.3473 | −0.3461 | 0.0781 | −0.0770 | −0.0001 |

2ν_{3} + 3ν_{4} | A_{u} | 2945.9496 | −0.0863 | −0.1090 | 0.1109 | 0.0182 | −0.0222 | −0.0053 |

2ν_{3} + ν_{6} | B_{1g} | 2966.8810 | 0.2657 | 0.2261 | −0.2239 | 0.0418 | −0.0398 | 0.0001 |

rms error | 0.4347 | 0.3932 | 0.3915 | 0.0569 | 0.0539 | 0.0018 |

Band . | Sym . | (i) . | (ii) . | (iii) . | (iv) . | (v) . | (vi) . | (vii) . |
---|---|---|---|---|---|---|---|---|

GS | A_{u} | 11.2122 | −0.0034 | −0.0034 | 0.0033 | −0.0003 | 0.0000 | −0.0001 |

ν_{4} | A_{u} | 371.3261 | −0.0014 | −0.0015 | 0.0015 | −0.0001 | −0.0001 | −0.0002 |

2ν_{4} | A_{u} | 776.8274 | 0.0119 | 0.0103 | −0.0101 | 0.0014 | −0.0018 | −0.0006 |

ν_{3} | A_{u} | 877.6885 | 0.0087 | 0.0077 | −0.0078 | 0.0010 | −0.0011 | 0.0000 |

ν_{3} + ν_{4} | A_{u} | 1227.0708 | 0.0134 | 0.0112 | −0.0109 | 0.0020 | −0.0022 | −0.0005 |

3ν_{4} | A_{u} | 1245.4558 | 0.0155 | 0.0131 | −0.0128 | 0.0023 | −0.0025 | −0.0006 |

ν_{6} | B_{1g} | 1284.6829 | 0.1640 | 0.1548 | −0.1546 | 0.0099 | −0.0097 | 0.0000 |

ν_{2} | A_{u} | 1400.7098 | 0.1674 | 0.1574 | −0.1571 | 0.0107 | −0.0104 | 0.0000 |

ν_{3} + 2ν_{4} | A_{u} | 1638.1640 | 0.0202 | 0.0170 | −0.0162 | 0.0039 | −0.0038 | −0.0007 |

ν_{4} + ν_{6} | B_{1g} | 1648.4799 | 0.2196 | 0.2066 | −0.2064 | 0.0138 | −0.0136 | 0.0001 |

4ν_{4} | A_{u} | 1720.1736 | 0.0330 | 0.0261 | −0.0260 | 0.0066 | −0.0077 | −0.0011 |

3ν_{3} | A_{u} | 1726.8764 | −0.0076 | −0.0101 | 0.0102 | 0.0025 | −0.0025 | 0.0000 |

ν_{2} + ν_{3} | A_{u} | 1772.8478 | 0.2064 | 0.1894 | −0.1891 | 0.0181 | −0.0178 | 0.0000 |

2ν_{4} + ν_{6} | B_{1g} | 2072.6916 | 0.2268 | 0.2033 | −0.2030 | 0.0269 | −0.0269 | −0.0004 |

2ν_{3} + ν_{4} | A_{u} | 2076.4227 | 0.0016 | −0.0028 | 0.0034 | 0.0044 | −0.0043 | −0.0005 |

ν_{3} + 3ν_{4} | A_{u} | 2103.3759 | 0.0231 | 0.0160 | −0.0152 | 0.0072 | −0.0075 | −0.0010 |

ν_{3} + ν_{6} | B_{1g} | 2134.6189 | 0.2206 | 0.2036 | −0.2025 | 0.0193 | −0.0178 | 0.0005 |

ν_{2} + 2ν_{4} | A_{u} | 2170.7709 | 0.2092 | 0.1735 | −0.1730 | 0.0400 | −0.0398 | −0.0004 |

5ν_{4} | A_{u} | 2205.8643 | 0.0567 | 0.0305 | −0.0303 | 0.0267 | −0.0280 | −0.0014 |

ν_{2} + ν_{3} | A_{u} | 2248.5455 | 0.2091 | 0.1873 | −0.1864 | 0.0246 | −0.0230 | 0.0008 |

2ν_{3} + 2ν_{4} | A_{u} | 2483.5870 | −0.0086 | −0.0154 | 0.0168 | 0.0074 | −0.0076 | −0.0015 |

ν_{3} + ν_{4} + ν_{6} | B_{1g} | 2496.0545 | 0.2672 | 0.2396 | −0.2389 | 0.0297 | −0.0282 | 0.0008 |

2ν_{6} | A_{u} | 2538.4081 | 0.9592 | 0.8801 | −0.8778 | 0.0921 | −0.0864 | 0.0027 |

3ν_{4} + ν_{6} | B_{1g} | 2551.8234 | 0.3660 | 0.3021 | −0.2999 | 0.0689 | −0.0670 | −0.0006 |

3ν_{3} | A_{u} | 2558.3483 | −0.6483 | −0.6569 | 0.6582 | 0.0035 | −0.0048 | −0.0015 |

ν_{3} + 4ν_{4} | A_{u} | 2579.2230 | 0.0632 | 0.0440 | −0.0435 | 0.0193 | −0.0201 | −0.0012 |

ν_{2} + ν_{3} + ν_{4} | A_{u} | 2614.6564 | 0.3396 | 0.2791 | −0.2773 | 0.0658 | −0.0633 | 0.0005 |

ν_{2} + 3ν_{4} | A_{u} | 2628.8437 | 0.3410 | 0.2729 | −0.2703 | 0.0741 | −0.0713 | −0.0001 |

ν_{2} + ν_{6} | B_{1g} | 2660.3104 | 1.4349 | 1.3261 | −1.3199 | 0.1318 | −0.1203 | 0.0045 |

6ν_{4} | A_{u} | 2705.4012 | 0.1297 | 0.0716 | −0.0697 | 0.0603 | −0.0604 | −0.0023 |

2ν_{2} | A_{u} | 2768.9855 | 0.8934 | 0.7426 | −0.7364 | 0.1715 | −0.1623 | 0.0022 |

3ν_{2} + ν_{4} | A_{u} | 2908.1590 | −0.5381 | −0.5535 | 0.5561 | 0.0049 | −0.0089 | −0.0050 |

ν_{4} + 2ν_{6} | A_{u} | 2910.8426 | 1.1266 | 0.9862 | −0.9787 | 0.1694 | −0.1578 | 0.0027 |

ν_{3} + 2ν_{4} + ν_{6} | B_{1g} | 2918.9996 | 0.4153 | 0.3473 | −0.3461 | 0.0781 | −0.0770 | −0.0001 |

2ν_{3} + 3ν_{4} | A_{u} | 2945.9496 | −0.0863 | −0.1090 | 0.1109 | 0.0182 | −0.0222 | −0.0053 |

2ν_{3} + ν_{6} | B_{1g} | 2966.8810 | 0.2657 | 0.2261 | −0.2239 | 0.0418 | −0.0398 | 0.0001 |

rms error | 0.4347 | 0.3932 | 0.3915 | 0.0569 | 0.0539 | 0.0018 |

## V. CONCLUSION

In this paper, we have proposed a hybrid nuclear-motion Hamiltonian based on the HBJ formalism but written in terms of ITOs adapted to a given molecular symmetry group. The treatment of the nonrigid coordinates is performed numerically on a grid of points suitably chosen in order to accelerate the computation of the energy levels. For each point of the grid, a set of vibrational–rotational ITOs is constructed. In the introduction part, we asked a question about the pertinence of using the tools that were developed for semirigid molecules. We can answer positively to this question. Indeed, we have shown that the ITO formalism initially introduced in the construction of the effective Hamiltonian for semirigid spherical top molecules, as well as the Hamiltonian and basis-set reduction procedure introduced to accelerate calculations, are also both designed for nonrigid molecules. The hybrid model has been validated on small nonrigid systems (CH_{2}, CH_{3}, NH_{3}, H_{2}O_{2}) and will pave the way for the study of more complex molecules, with possibly several large amplitude motions, for the analysis of high-resolution spectra, and for the construction of molecular line lists.

Undoubtedly, the construction of nonrigid effective models, as proposed in Ref. 90, would be of great help in refining parameters quite easily on observed data. Moreover, an upcoming version of TENSOR is in preparation and will benefit from the use of iterative eigensolvers and nested contracted basis functions in conjunction with group-theoretical considerations in order to consider molecules with more than seven atoms as ethane-like molecules.

## SUPPLEMENTARY MATERIAL

Fortran code including the H2O2 PES in symmetry coordinates used in this work and generated from the ab initio PES by Malyszek and Koput.148

## ACKNOWLEDGMENTS

The authors acknowledge support from the Romeo Computer Center in Reims-Champagne-Ardenne, from the program IRP “SAMIA2,” and from the French ANR TEMMEX project (Grant No. 21-CE30-0053-01). This work was supported by the Russian Scientific Foundation (RSF, Grant No. 22-42-09022). This work was dedicated to Professor Per Jensen, who is gratefully acknowledged for the fruitful discussions he shared with the authors during his visit to Reims.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Michaël Rey**: Conceptualization (lead); Funding acquisition (equal); Methodology (lead); Software (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). **Dominika Viglaska**: Software (equal); Writing – review & editing (equal). **Oleg Egorov**: Software (equal); Validation (equal); Writing – review & editing (equal). **Andrei V. Nikitin**: Funding acquisition (equal); Software (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: DERIVATION OF THE ROTATION MATRIX **U**(*ρ*) IN EQ. (1)

^{152}already derived the rotation matrix

**U**(

*ρ*) involved in Eq. (1) using the Floquet theory.

^{153}More recently, an alternative method was proposed in Ref. 93 and consists of integrating the evolution equation

**Ω**(

*ρ*) is a skew-symmetric matrix involving the LAM dependent angular velocity vector components (

*ω*

_{x},

*ω*

_{y},

*ω*

_{z}) defined by $\omega (\rho )=\u2212(I\u2009Rot-Rot)\u22121I\u2009Rot-LAM$, where

**I**

_{Rot-Rot}(

*ρ*) is the rotational tensor of inertia. Introducing the three generators of the so(3) Lie algebra

*α*=

*x*,

*y*, and

*z*). Here, the three components $px(\rho )$, $py(\rho )$, and $pz(\rho )$ are to be determined. Similarly, we can write

**Ω**(

*ρ*) =

*∑ω*

_{α}(

*ρ*)

*S*

_{α}. After some algebra, we can show that

*B*

_{k}is the Bernoulli number. Bearing in mind that the only non-vanishing term

*B*

_{2r−1}is given for

*r*= 1 and using the expansion in the Laurent series of cot(

*z*) around

*z*= 0, the elements

*p*

_{α}are determined by solving a system of three ordinary differential equations

*p*

_{α}’s as

^{71}we could also introduce the Euler angles as

## REFERENCES

*Higher-Order Effects in the Vibration-Rotation Spectra of Semirigid Molecules*

*Molecular Vibrational-Rotational Spectra*

*Molecular Symmetry and Spectroscopy*

*Vibration-Rotational Spectroscopy and Molecular Dynamic–Advances in Quantum Chemical and Spectroscopical Studies of Molecular Structures and Dynamics*

*Advanced Series in Physical Chemistry*

*Fundamentals of Molecular Symmetry*

*Principles of Symmetry, Dynamics, and Spectroscopy*

*Spherical Top Spectra*

*Group Theory and Its Applications to the Quantum Mechanics of Atomic Spectra*