A numerical-tensorial “hybrid” nuclear motion Hamiltonian and dipole moment operator for spectra calculation of polyatomic nonrigid molecules

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⁵⁰ 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,⁷⁴ 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⁵⁶ and TROVE⁵⁷ 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₂, CH₃, NH₃, and H₂O₂.

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,⁷⁵ 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(ρ)$ instead of $aie$⁠, where ρ is a coordinate describing the LAM, and (ii) treated the remaining 3N − 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 3N − 7 coordinates vary smoothly around $airef(ρ)$⁠, 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⁷⁵ and its rigid counterpart, namely the Eckart–Watson Hamiltonian,⁵⁰ 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).

As an illustrative example, we have determined the rotation matrix U(ρ) in order to compute a new reference configuration satisfying (1) for the H₂O₂ molecule. As a starting point, we have considered the reference configuration given in Ref. 94, called $airef,ini(ρ)$ such as I_yρ ≠ 0. By solving Eq. (A6), the new IAS configuration is given by $airef,IAS(ρ)=U(ρ)airef,ini(ρ)$⁠. 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⁹⁵) minimizes the interaction between the small and large amplitude vibrations.

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 $ρ̄$ coordinate when the “rigid” vibrations are beyond their equilibrium geometry. By construction, the HBJ KEO depends explicitly on ρ and on 3N − 7 rectilinear, normal coordinates Q_i, while the ab initio potential function is built from a function $ū(ρ̄)$⁠, where $ρ̄$ 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 3N − 6 coordinates $Sk,σk(Γ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,σk(Γk)$ coordinates, one corresponds to $ū(ρ̄)$ 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 3N − 7 coordinates $Sr,σr(Γr)$ for the rigid part and put the dependence in $ρ̄$ into the expansion coefficients. As a direct consequence, some linear terms like $∑rF̄r(ρ̄)Sr(Γr)$ will now appear in the potential, while the quadratic part will be written as $∑rF̄r,r′(ρ̄)Sr(Γr)Sr′(Γr′′)$⁠. 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 $ρ̄$⁠. To this end, a general numerical recipe was proposed in Ref. 87 to find the relation $ρ̄=f(ρ,Sr)$ for any polyatomic molecules so that the potential function transforms as $V(Sr,ρ̄)→V(Sr,ρ)$⁠. The procedure can be summarized as follows:

• First, we assume that all calculations are performed on a grid composed of M_ρ points,

  $ρm∈[ρmin,ρmax],m=1,…,Mρ,$

  (2)

  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 (3N − 7) × (3N − 7) writes as G = G(ρ_m) due to its dependence on $aiαref(ρ)$⁠, whereas the elements of the force constant matrix write $F̄rr′=F̄rr′(ρ̄)$⁠. The aim is thus to find the transformation $F̄(ρ̄)→F(ρm)$ to properly solve the GF problem on the grid. To this end, we have shown in Ref. 87 that $ū(ρ̄)$ must necessarily be of the form

  $ū(ρ̄)=u(ρm)+∑rCr(ρm)Sr(Γr)+∑rr′Crr′(ρm)Sr(Γr)Sr′(Γr′),$

  (3)

  where the 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′r″r′′′$⁠, etc. were required.⁵⁶ 
• In the next step, we expand the 3N − 7 coordinates $Sr(Γr)$ as well as the seven Eckart–Sayvetz relations

  $Tα=1M∑imidiα,Rα=1Iααref(ρ)∑imiϵαβγaiβref(ρ)diγ,S=1Iρρref(ρ)∑i,αmidaiαref(ρ)dρdiα,$

  (4)

  in terms of the 3N small amplitude Cartesian displacements d_i. Here, M is the total mass of the system, ϵ corresponds to the Levi-Cività tensor, and $Iααref$⁠, $Iρρref$ are the diagonal elements of the 4 × 4 inertia matrix.⁷⁵ 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_id_j to determine the C_r and C_rr′ coefficients in Eq. (3). To this end, we first construct the 3N linear and 3N(3N + 1)/2 quadratic forms $Xslin$ and $Xs″quad=XslinXs′lin$ [s, s′ = 1, …, 3N, s ≥ s′; s″ = 1, …, 3N(3N + 1)/2] with $Xlin=(S,Tα,Rα,S)$⁠. Similarly, we define the vector Y^lin = (d_1x, …, d_Nz) and also construct the quadratic form $Ys″quad=YsYs′$⁠. 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 that

  $(Xlin,Xquad)t=Bquad(ρm)(Ylin,Yquad)t.$

  (5)

  The Cartesian displacements contained in Y can thus be deduced by inverting Eq. (5) and by the seven constraints $Tα=Rα=S=0$⁠. We thus write

  $(d1,…,dN)t=Aquadinv(ρm)(Xlin,Xquad)t,$

  (6)

  with now X^lin = (S, 0, 0, 0) and where $Aquadinv$ is a sub-matrix extracted from $Bquad−1$⁠. 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 $ū(ρ̄)$ 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

  $V(Sr,ρ̄)→V(Qr,ρm)≡Vρm(Qr),r=1,…,3N−7,$

  (7)

  using Eq. (3) as well as the 3N × (3N − 7) orthogonal transformation l(ρ_m) linking the mass-weighted Cartesian displacements to the normal coordinates.⁹⁹ 

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 $ζss′α(ρm)$ to a polynomial function in u(ρ) in order to expand both the KEO and potential parts in terms of the 3N − 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.

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:

1. 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.

2. 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.

3. 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⁵⁶ and TROVE⁵⁷ computer codes. In this work, we intend to use the HBJ Hamiltonian in its original formulation based on the 3N − 7 normal coordinates, including all terms in Eqs. (10)–(14) below, for solving variationally the vibration-LAM-rotation wave equation.

4. 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.

5. 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,¹⁰⁰ BELGI-C₁,¹⁰¹ BELGI-C_s,¹⁰² BELGI-C_s-2Tops,¹⁰³ ERHAM,¹⁰⁴ or RAM36⁸⁵ codes). Recently, we have proposed a numerical approach for the construction of ab initio effective models,⁹⁰ 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.

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 3N − 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 3N − 7 normal coordinates, a modified version of the computer program COSY INFINITY,¹⁰⁵ 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(∂⁡ln⁡μ∂Qs)/dρ$ 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⁶³) could be carried out with only a few modifications. As recalled by Bunker and Jensen⁷⁴ 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 ρ₁, ρ₂, etc. Therefore, the HBJ Hamiltonian could be expressed in terms of 3N − 6 − N_ρ normal-mode coordinates and of the 3 + N_ρ rotational components J_x, J_y, J_z, $Jρ1$⁠, $Jρ2$⁠, etc.

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₂, C_i, C_s, C_2v, C_2h, D₂ 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×ε$ or G₁₂ is isomorphic to the group D_3h(M) or to the point group D_3h of the planar configuration, while for phosphine, G₁₂ 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⁷⁴  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₄(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¹⁰⁹ presented a series of “recipes” aimed at helping the user handle permutation-inversion groups for practical spectroscopic applications. As stated by Woodman,⁶⁷ 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+×C3v−$⁠, 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∧C$ 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 6C and 9C 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.

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₆, G₁₂, G₂₄, G₃₆, etc. Within this context, the use of ITOs is relevant to deal with degenerate vibrational modes like ν₃ and ν₄ of the NH₃ 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.¹¹² 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.⁸⁸ They were thus employed with success in the construction of spectroscopic line lists of semirigid molecules from the first principles.^{13,21,51–53}

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.

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′vlj$ 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.

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.¹²⁷ The treatment of MS groups can be carried out in the same fashion.

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

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.

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 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 $Ψred(m→m′)$⁠, 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.

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₂ methylene molecule, the planar CH₃ methyl radical, the NH₃ ammonia molecule with inversion, and the H₂O₂ (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.

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⁻¹, methylene turns out to be a good first candidate to validate our hybrid model. A recent accurate ab initio PES for CH₂ has been calculated,¹³³ 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 Γ₀ (resp. Γ) for even (resp. odd) powers. The standard harmonic oscillator basis functions (HOBF) were used for the stretching modes, while the displaced HOBF $NvHv(ρ′)e−ρ′2/2$⁠, with ρ′ = 5(ρ − 0.7), multiplied by a weight function $sin(ρ)$⁠, 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₁ and q₃ 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₁ functions and 1925 B₂ functions was used to compute the vibrational energy levels, which are converged within 10⁻⁴ up to 10 000 cm⁻¹ 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.¹³⁵ 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⁻⁵ up to 10 000 cm⁻¹. 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.2, the number of basis functions is decreased by 35%, and the rms error becomes 0.0002 up to 10 000 cm⁻¹. By taking 50% fewer functions (κ₁ = κ₃ = 1.4), the convergence of the levels is not fully achieved (rms error of 0.001 cm⁻¹), but the results remain suitable for high-resolution spectroscopy.