Michaël Rey
email: michael.rey@univ-reims.fr

Iana S Chizhmakova

Andrei V Nikitin

Vladimir G Tyuterev

Towards a complete elucidation of the ro-vibrational band structure in the SF 6 infrared spectrum from full quantum-mechanical calculations

First accurate and complete theoretical room-temperature rotationally-resolved spectra in the range 300-3000 cm -1 are reported for the three most abundant isotopologues ( 32 SF 6 , 33 SF 6 and 34 SF 6 ) of the sulfur hexafluoride molecule. The SF 6 was widely used in the literature as prototype molecule for studying the multi-photon excitation processes with powerful lasers in the infrared range. On the other hand, SF 6 is an important greenhouse molecule with a very long lifetime in the atmosphere. Because of relatively low vibrational frequencies, the hot bands of this molecule contribute significantly to the absorption infrared spectra even at room temperature. This makes the calculation of complete rovibrational line lists required for fully converged opacity modeling extremely demanding. In order to reduce the computational costs, symmetry was exploited at all stages of the first global variational nuclear motion calculations by means of irreducible tensor operators. More than 2600 new vibrational band centers were predicted using our empirically-refined ab initio potential energy surface. Highly-excited rotational states were calculated up to J = 121, resulting in 6 billion transitions computed from an ab initio dipole moment surface and distributed over more than 500 cold and hot bands. The final line lists are made available through the The-oReTS information system (http://theorets.univ-reims.fr, http://theorets.tsu.ru). For the first time, the major (ro)vibrational band structures in the wavenumber range corresponding to the strongest absorption in the infra-red are completely elucidated for a seven-atomic molecule, providing excellent agreement with observed spectral patterns. It is shown that the obtained results are more complete than all available line lists, permitting reliable modelling of the temperature dependence of the molecular opacity.

Introduction

SF 6 is a very important greenhouse gas with a global warming potential ∼ 29000 times larger than that of CO 2 and its atmospheric lifetime can last up to 3200 years 1,2 since it is chemically inert. SF 6 is an important prototype molecule for infrared laser chemistry and laser isotope separation and for ultrahigh resolution spectroscopy. In particular, various problems related to the laser excitation of SF 6 in the quasi-continuum band formed by congested spectral patterns have been discussed [3][START_REF] Alimpiev | Excitation Spectrum of SF6 Irradiated by an Intense IR Laser Field[END_REF][START_REF] Lyman | Single-Infrared-Frequency Studies of Multiple-Photon Excitation and Dissociation of Polyatomic Molecules[END_REF] . As was reviewed in Ref. [START_REF] Lyman | Single-Infrared-Frequency Studies of Multiple-Photon Excitation and Dissociation of Polyatomic Molecules[END_REF] , a study of the literature indicated that, in many cases, the experimental data related with multiple-photon excitation differed by factors of 2 or 3 in absolute value. This makes important a more profound insight in the ro-vibrational patterns of levels and transitions of SF 6 in the infrared range. This molecule is mainly used as a gaseous dielectric insulating medium in the electrical industry 6 as well as as an atmospheric 7,8 and ocean [START_REF] Fine | CFCs and SF 6 as tracers of ocean processes In: Encyclopedia of Ocean Sciences[END_REF] tracer while its concentration is rapidly increasing in the atmosphere. Vertical profiles have been retrieved by Rinsland et al. [START_REF] Rinsland | [END_REF][11][12] from a series of spectroscopic measurements. Though detected as a trace species, this gas thus plays an important role in the climate changeas recognized by the Kyoto protocolmaking the need to have consistent theoretical models to monitor emissions. SF 6 belongs to the O h point group at the equilibrium structure and possesses 15 vibrational degrees of freedom. Due to its high octahedral symmetry, the molecule has doubly and triply degenerate vibrational modes of E and F-types. Among the 6 fundamental vibrations: three stretching modes (ν 1 , ν 2 , ν 3 ) transform according the A 1g , E g and F 1u irreducible representations and three bending modes (ν 4 , ν 5 , ν 6 ) according to F 1u , F 2g and F 2u . The g type modes are Raman active while ν 3 and ν 4 are infrared (IR) active; only ν 6 being inactive. A brief inspection of the harmonic vibrational modes (see Table 1) shows no simple relations between the ω i . Consequently, the SF 6 spectrum does not exhibit a clear polyad structure, contrary to many other spherical top molecules like methane [START_REF] Champion | Spherical Top Spectra[END_REF] . Moreover, the presence of relatively low degenerate vibrational modes, combined with a small rotational constant (B e ≈ 0.09 cm -1 ), make the rotationally-resolved IR spectrum very dense and crowded, even at room-temperature (RT). This particular feature gives rise to a tremendous number of degenerate and quasi-degenerate states which strongly complicate analysis of high-resolution spectra. Indeed, it can be shown that the density of lines in SF 6 spectra at T =296 K is similar to that of small hydrocarbons like CH 4 , C 2 H 2 , C 2 H 4 14-16 at much higher temperatures (T ≈ 700 -1500 K). The strongest absorption region corresponds to the stretching fundamental IR band ν 3 (F 1u ) located in the atmospheric windows around 10.5 µm is detectable by IR remote-sensing techniques 7,[START_REF] Rinsland | [END_REF] . However, the presence of low-lying vibrational energy levels makes the Boltzmann population of the excited states very important at RT leading to a substantial number of hot band (HB) transitions that significantly contribute to the overall opacity. For example, the HB transitions ν k + ∑ n i ν i -∑ n i ν i (k =3 or 4, i = k) falling in the two most congested ν 3 and ν 4 regions represent 66 and 70%, respectively, of the integrated intensity at RT. As expected, the "forbidden" mode ν 6 of lowest frequency is involved in most of the HB transitions (n 6 = 1, 2, 3, • • • ) allowed by symmetry. This paper is devoted to the construction of a first comprehensive, accurate and complete theoretical spectrum in the range 300-3000 cm -1 from variational nuclear motion calculations. To achieve this end, our recent full-dimensional ab initio potential energy (PES) and dipole moment (DMS) surfaces 17 were used. The quadratic force constants the PES were empirically optimized for this work. The final aim is to fill the gap between what is available in the current spectroscopic databanks and what is really needed to properly model spectra, in particular in the 10 µm windows of atmospheric interest. The paper is structured as follows. The state-of-the art of the SF 6 infrared spectroscopy as well as the main motivation for using global variational calculations are presented in Section 2. Section 3 is devoted to the construction of the Hamiltonian model, the choice of basis sets and the computation of ro-vibrational energy levels.

The full calculated line lists for the two most abundant species of the sulfur hexafluoride molecule, namely 32 SF 6 and 34 SF 6 , as well as validation on raw experimental spectra are presented in Section 4.

Effective versus global variational approach

Undoubtedly, the biggest challenge concerns the modelling of all HBs, in particular those falling in spectral windows of atmospheric interest. To this end, highly-excited ro-vibrational states have to be properly modelled to predict all necessary transitions belonging to cold and hot bands, even at RT conditions. The past four decades, a series of experimental and theoretical studies was devoted to this task to improve analysis and spectra predictions (see e.g. Refs. [18][19][20][21] ). Intensities of the infrared-active fundamental vibrations have been also studied 22,23 . Recently, much effort was made to analyse high-resolution SF 6 spectra using phenomenological effective Hamiltonians whose parameters were empirically-adjusted to experimental data [24][25][26][27][28][29] . These data have been included in line list compilations like HITRAN 30 or GEISA 31 . The most extensive empirical model to date has been constructed for the main isotopologue 32 SF 6 by the Dijon group for modelling the band system 28

(ν 1 , ν 2 , ν 3 , ν 1 + ν 3 , ν 2 + ν 3 , 2ν 3 , ν 1 + ν 2 + ν 3 , 3ν 3 , ν 3 -ν 1 , ν 3 -ν 2 )
, recently updated with the set 32 (ν 1 + ν 4 , ν 2 + ν 6 , ν 5 + ν 6 ). All these combination bands and overtones are of key importance for computing a part of HBs falling in the strongly absorbing ν 3 and ν 4 regions. However, a purely empirical approach based on effective Hamiltonians faces several conceptual difficulties with the increasing number of atoms. First, a new analysis has to be conducted for each isotopologue, unlike global variational calculations for which all species can be predicted almost simultaneously (see e.g. Ref. 33 ). Second, effective spectroscopic models are based on a line-by-line assignment procedure that suffers from the completeness issues in spectra analyses of heavy atmospheric molecules like SF 6 because of the presence of very few isolated lines in the highly congested spectra. The lack of assignments as well as ambiguities in the fits of the effective models could often lead to a poor determination of some fitted resonance coupling and empirical dipole transition moment parameters. In such cases, the ro-vibrational effective wavefunctions are not reliably determined and do not permit a correct description of the intensity borrowing effects induced by accidental resonances. In order to partly get around the problem, the most recent analyses have been realized in a pragmatic way: all experimental spectra have been cooled down (∼ 150 K). This permitted diminishing the density of line by removing high-J transitions as well as many HBs. In turn, a large amount of information concerning HBs is clearly missing at 296 K. For this reason, the recent line list in SHeCaSDa SF 6 database 32 as well as the closely related HITRAN 30 and GEISA 31 compilations do not contain enough line parameters to reproduce and explain the vibrational band structure at 296 K, in particular around ν 3 . Completeness of rovibrational bands will be thus the key requirement for a successful modelling of the radiative transfer.

To analyse rotationally resolved SF 6 absorption spectra in the infrared, it is imperative to elucidate the vibrational structure first. As explained above, all the previous studies have been done essentially by using the empirical approach with ad hoc effective models which clearly showed their limits. An important drawback of a purely empirical approach was that no systematic line intensity measurements have been carried out so far for this molecule because of the difficulty to extract non-blended lines among the highly congested spectral patterns. In this context, the present paper reports for the very first time full-dimensional, accurate and complete variational calculations in the range 300-3000 cm -1both for line positions and line intensitiesfor the three most abundant species of the sulfur hexafluoride molecule, namely 32 SF 6 (∼ 95%), 34 SF 6 (∼ 4%) and 33 SF 6 (∼ 0.7%). Though they may suffer from convergence issues concerning line positions, variational global calculations proved to be very efficient to accurately predict line intensities for both cold and hot bands [START_REF] Rey | [END_REF]16,[34][35][36][37][38] . An importance of anharmonicity in vibrational spectra of polyatomic molecules have been studied in many previous works (see e.g. Refs. [39][40][41][42][43][44] ). The focus of our study is different. Here, the aim was to produce complete line lists by (i) including line-by-line intensities for a variationally converged set of spectral transitions at medium or high spectral resolution, Table 1 Ab initio (CCSD(T)/CVQZ(DK)) 17 and optimized (this work, TW) symmetrized force constants

F i j = ∂ 2 V /∂ S (Γ) iσ ∂ S (Γ)
jσ (in mdyn, 1dyn=10 -5 N) and equilibrium geometry r e (in Å, 1 Å= 100 pm) for SF 6 . The harmonic wavenumbers ω i (in cm -1 ) are given for 32 SF 6 . The values in the column "best" correspond to a final empirical optimization (see text) and was used to compute all energy levels. Here ω opt and ω best are the same. ab initio 17 (ii) including all hot bands and (iii) providing the temperature dependence of the global opacity. As discussed in a recent review on theoretical rotationally resolved spectra predictions 37 , only very few groups [START_REF] Tennyson | Proc. Int. Astron. Union[END_REF][START_REF] Huang | [END_REF][47] (and see references therein) have recently developed specific tool for such applications, but for smaller molecules. In this paper, we focus on the main results permitting to explain the complicated features in observed SF 6 spectra. Technical details about the extension of the methods previously developed for smaller molecules 34,[48][49][50][51] are beyond the scope of this work.

Quantum-mechanical nuclear motion calculations

Refined potential energy surface

As a starting point, the full-dimensional ab initio potential and dipole moment surfaces have been computed at the CCSD(T)/CVQZ level of the theory (including Douglas-Kroll relativistic corrections), as reported in Ref. 17 . However, the harmonic frequencies computed from the pure ab initio force constant matrix F ab of the PES, together with the geometry-based G Wilson matrix, were off by about ∼1 to 10 cm -1 . The ab initio equilibrium geometry r ab e was not accurate enough to reproduce the rotational structure. In order to better match experimental fundamental band centers, the approach proposed in Ref. 34 was adopted to refine the initial solutions {F ab ,r ab e }.

A set of optimized quadratic force constants F opt for each symmetry block as well as a refined value r opt e were thus obtained (see Table 1) while the high-order ab initio force constants were kept unchanged and fixed to the values of Ref. 17 . The axis convention, the atom numbering as well as the definition of the symmetry coordinates can be found in Refs. 17,52 .

Hamiltonian model and symmetry

Since SF 6 is a semirigid molecule without large amplitude vibrations, the normal mode coordinates q

(Γ k )
kσ k were used to build the Hamiltonian in the Watson-Eckart formalism 53 implemented in the TENSOR computer code. Here, Γ k denotes an irreducible representation and σ k a label to distinguish the degenerate components of a vibrational mode k. Simple analytically solvable potential representations for independent vibrational modes do not permit approaching experimental accuracy in high-resolution spectra. The aim of this work was to account as fully as possible for the high-order anharmonicity of the PES and for the coupling among ro-vibrational states. To this end, we do not use perturbation theory, as often employed in the literature for polyatomic molecules, but the symmetry-adapted variational approach in 18 dimensions: 15 vibrational and 3 rotational degrees of freedom. For a full account of symmetry and to build linearly independent polynomials in q at a given order of the expansion of the complete nuclear motion Hamiltonian, the use of irreducible tensor operators (ITOs) is an optimal choice. The ITO formalism in the version proposed by Nikitin et al. 54 turned out to be very appropriate to describe the coupling of degenerate vibrations. Note that the ITO formalism, initially designed for effective models [START_REF] Champion | Spherical Top Spectra[END_REF] , is extended here for the first time for global variational nuclear motion calculations of XY 6 -type molecules. Concerning the construction of rotational basis functions adapted to the non-Abelian O h point group for very high J values, the symmetrization procedure proposed by Rey et al. 55 was employed. In summary, all the (3+15) rotation-vibrational degrees of freedom were simultaneously taken into account to build the sum-of-product Hamiltonian in the ITO representation-including all intra-molecular coupling terms and anharmonicities at a given order-as well as the direct-product symmetry-adapted basis functions.

At the next step, a "brute-force" 18-dimensional variational calculation of complete rotationally resolved spectra using so many terms in the Hamiltonian up to J = 120 with standard tools and direct eigensolvers was not feasible. To overcome the first obstacle related to the large number of Hamiltonian terms, the reduction technique described in Refs. 48,50,51 was applied in a slightly modified form: instead of expressing the reduced tensor Hamiltonian in terms of harmonic oscillator creation and annihilation operators, this latter was expressed directly in q (Γ) i and in conjugate momenta p

(Γ) i = -id/dq (Γ)
i . This allowed us to drastically reduce the number of terms and thus to accelerate the calculation of all matrix elements via the Wigner-Eckart theorem. Briefly, the initial m-th order polynomial expansion of the Eckart-Watson Hamiltonian in (q

(Γ) i , p (Γ)
i ) was reduced to a transformed polynomial a smaller order m < m without loss of accuracy for the eigenvalues. It was then converted to the ITO representation 34 . For this work, both the normal-mode, refined PES and kinetic energy part have been expanded at order m = 9 and reduced at order m = 6. A final ITO Hamiltonian H (9→6) composed of ∼ 5600 non-zero ro-vibrational parameters was built. In the case of SF 6 , such normal mode representation proved to be very efficient and allowed a fast convergence of the Hamiltonian expansion. This type of calculation inherently and systematically accounts for all intra-molecular coupling terms and for anharmonicities through a power series expansion in normal coordinates q in both the kinetic and potential parts.

Computation of vibrational energy levels

At the next step, we have solved the eigenvalue/eigenfunction problem for the full anharmonic potential with J = 0. The vibrational energies were obtained by the exact numerical diagonalization of the Hamiltonian matrix built in a direct-product "pruned" basis defined by a set of primitive harmonic oscillator vibrational functions such that 34

F κ (n) ⇔ 6 ∑ i=1 κ i v i ≤ n, (1)
with

v i = 0, • • • , n.
Here, κ i (≥ 1) are weight coefficients chosen to select an appropriate number of stretching and bending basis functions.

The maximum number of quanta v max i for each mode is given by the integer part n/κ i . For this work, n = 10 and κ = {1.05, }. This compressed anharmonic vibrational basis will be used in the next step for rovibrational calculations. A comparison of our variationally computed energies with the empirical levels derived from effective models is given in Table 2 for 32 SF 6 . The error in the fundamental bands is below 10 -3 cm -1 . A very good agreement proves the high accuracy of our refined PES, on the one hand, and validates all the theoretical "machinery" for solving the nuclear Schrödinger equation, on the other hand. However, there remain relatively large discrepancies (e.g. -0.4 cm -1 for 2ν 1 + ν 3 ) that can be partly explained by the fact that only the quadratic force constants of the PES have been empirically refined using fundamental band origins. By doing this, we assumed that a large part of the errors was due to uncertainty in the quadratic terms ("harmonic" contribution to the PES). Higher-order terms in the PES also play a substantial role to properly model highly-excited states but they were fixed to their ab initio values in this study.

Consequently, though they were improved some overtones and combination bands did not reach the same accuracy as the fundamentals bands.

In order to check the basis set convergence of J = 0 energies, we have also computed the Γ = A 1g levels by using a larger set of primitive functions with the cut-off n = 11 in Eq. ( 1). In comparison with n = 10, the convergence errors for 1-, 2-, 3-, 4-, 5-and 6-quanta levels turned out to be 0.0001, 0.0003, 0.006, 0.15, 0.9 cm -1 in average. Clearly, though their contributions to the absorption spectrum are almost negligible, very highly-excited vibrational levels are not fully converged as expected. Unfortunately, the use of the F κ (11) basis setcomposed of about twice more basis functions than F κ (10)goes beyond our computer resources to treat all symmetry blocks. Finally, the most extensive list composed of 2686 (in comparison with ∼25 ones available in the literature) selected vibrational energy (sub)levels up to 3000 cm -1 is given in Supplementary Material S 1 . Selected vibrational band centers for 33 SF 6 and 34 SF 6 as well as the isotopic shifts ∆(32 → 33) and ∆(32 → 34) are also given in Table 2 . As already noted in Refs. 21 , "anomalous" negative frequency shifts are observed. Such a feature is generally possible when some harmonic frequencies ω i remain unchanged during isotopic substitutions. For SF 6 , anomalous shifts are thus expected for ν 1 , ν 2 , ν 5 and ν 6 . As stated in Ref. 21 , "the change in the mixed-mode anharmonicity constants can lower the ground-state energy more than the excited-state energy, resulting in an increase in ν".

The negative shift in ν 1 between 32 SF 6 and 34 SF 6 was predicted to -0.044 cm -1 in Ref. 21 while Volkov et al. 57 have observed a shift of -0.0574 cm -1 . In this work, the shift is given by -0.0514 cm -1 which is close to observation. Similar shifts were predicted for ν 2 , ν 5 and ν 6 , increasing the band origins from 0.025 to 0.04 cm -1 . In this work, the band origins increase by 0.0269 to 0.0412 cm -1 . Note Fig. 1 Dimension of the Hamiltonian matrices to be diagonalized for the symmetry block F 1u as a function of J when using the full and reduced vibrational eigenfunctions Ψ [START_REF] Rinsland | [END_REF] and

Ψ (10→7) v
, respectively.

that a similar situation was reported in Ref. 58 between 12 CH 4 and 13 CH 4 for nν 2 because ω 2 is unchanged during the substitution 12 C → 13 C. Note that recently, J = 0 energy levels were reported 59 for the UF 6 molecule that belongs to the same O h point group, by using a rectangular collocation approach and product of harmonic basis functions but without exploiting the full symmetry.

Rovibrational energy levels for high-J values

For J > 0, "brute-force" calculations from the direct-product basis set (1) coupled with symmetric-top wave functions are clearly proscribed using available direct eigensolvers. Instead of considering iterative eigensolvers, as it is often recommended 60 , we have followed the method developed in our previous work for smaller molecules 49 . To this end, a set of reduced eigenfunctions Ψ Γ(n→m) v associated with a basis F κ (m) was obtained using the projection technique onto an appropriately chosen smaller subspace with m < n for which J > 0 calculations are feasible. For this work, rovibrational energy levels were computed up to J max = 121 using m = 7 and in a last step a final compression of the basis set was applied by selecting p% of eigenvectors where p = 100 up to J = 61 and then gradually decreases to p = 33 at J = 121. This strategy allows to keep approximately the same memory requirement as J increases, bearing in mind that accurate calculations are generally not necessary for J > 80 since the corresponding transitions are very weak at room temperature. Note that the J max value was chosen to ensure converging the partition function Q(T ) at T = 296 K. This gives Q(296) = 1607253.7 for 32 SF 6 , within 0.05 % of error. To better illustrate the use of the reduction procedure, the maximum dimension of the eigenvector matrix for J = 121 is 106784×35311; without reduction it would have been (23 × 23) millions. Figure 1 displays the size of the symmetry blocks F 1u as a function of J when using or not the reduced vibrational eigenfunctions Ψ (10→7) v

. In this figure (see the red part), we clearly see that the number of rovibrational basis functions decreases beyond J = 61 because of reduction of the vibrational subspace which vary from 100% to 33%.

In order to validate the optimized value r opt e

for the equilibrium structure (see Table 1), let us first focus on the accuracy of our variationally predicted rotational energy levels E varia. rot . To this end, we did a direct comparison with energies E emp.

rot obtained from the was optimized only for low-J values by using a small basis set during our refinement procedure. A way to improve the accuracy was to introduce an empirical correction to E varia. rot by the term δ rot =1.8×

10 -6 J(J + 1) cm -1 . This correction is equivalent to a shift in the equilibrium geometry by δ r e = -1.5×10 -5 Å to give the best estimated value r best e = r opt e + δ r e = 1.555938 Å. As a direct consequence, F opt had to be slightly modified (≡ F best ) to preserve the accuracy of the vibrational levels. Finally, both the Hamiltonian model and energy levels were recalculated using the newly optimized couple {F best , r best e } given in Table 1. This results in the improved set of rotational levels (see Figure 2(d)) whose uncertainty does not exceed 0.0035 cm -1 at J = 80. Above this value, the extrapolation error of the effective model beyond the range of observed data also contributes to ∆ rot , which increases up to 0.1 cm -1 for J = 120.

Selected sets of computed rovibrational energy levels are displayed in Figure 2(b,c) in the region of ν 2 + ν 4 , 2ν 2 , ν 3 + ν 6 , ν 1 + ν 5 , ν 4 + 2ν 6 and ν 2 + ν 6 , listed in increasing energies. The sequence of quasi-degenerate groups of quantum levels called cluster is well known for SF 6 and was explained by semiclassical analysis 61 . The different colours give information on how the energy levels are distributed and on how the states are mixed altogether. Except strong local resonances (e.g. between ν 3 + ν 6 and ν 1 + ν 5 at J = 7, see figure 2(c)), the inter-band couplings appear to be relatively weak for SF 6 . As an indicator of the degree of the coupling, the average value of the absolute maximum contribution in the wavefunction decomposition for the first 100 vibrational levels is 0.922, close to 1. This may partly explain why the convergence of the variational calculation for SF 6 is rather fast compared to other molecules like methane or ethylene.

Construction of the line listsvalidation

Methodology

The ultimate aim of this work was the construction of the room-temperature ab initio based line lists for 32 SF 6 , 33 SF 6 and 34 SF 6 in the range 300-3000 cm -1 . Here, we briefly recall the methodology obtained by the extension of the approach previously developed in our previous works 34,[48][49][50][51] for smaller molecules.

-Dipole moment: the line intensities were computed using the home-made TENSOR code from the fourth order normal-mode DMS reported in Ref. 17 . In order to be consistent with the Hamiltonian, the dipole moment component µ z (q) which transforms as the irreducible representation F 1u in the molecular-fixed frame was also converted to ITOs. A set of 213 ab initio dipole moment parameters in the ITO representation was thus deduced algebraically and used to compute line intensities. The explicit use of symmetry considerations when building the dipole moment turns out to be of great help to deduce all the symmetry-allowed transitions (J ,C ) ← (J ,C ) via the selection rules ∆J = 0, ±1 and C ×C ⊃ A 1u . For intensity calculations, the nuclear spin statistical weights for 32 SF 6 and 34 SF 6 are 2,

species A 1τ , A 2τ , E τ , F 1τ and F 2τ (τ = g, u), respectively while those of 33 SF 6 have to be multiplied by a factor 4.

-Completeness: the construction of complete line lists is essential for atmospheric retrievals and opacity calculations. Completeness is controlled by two cut-off criteria, namely the lower-state energy and the intensity cut-offs. A prior in-depth study allowed us to fix the maximum lower state energy E max low at [2600 + 0.09J(J + 1)] cm -1 to control all necessary HBs. Figure 3 (left) displays all lower state energies involved in rovibrational transitions around the 10 µm by taking an intensity cut-off of 10 -25 cm/molecule at 296 K. In this figure are depicted all cold band transitions (E low ≡ E rot ) denoted as P i -P 0 as well as all hot band transitions starting from a lower state with n quanta of excitation, denoted as P i -P n (n > 0). To illustrate the gap with the literature, a comparison of the lower state energies involved in our calculations and in the currently available lists HITRAN and SHeCaSDa are shown in Figure 3 (right).

In order to properly converge the integrated intensities using an optimal number of lines, we chose to split the whole wavenumber range W in three intervals. For each interval of 32 SF 6 , a specific intensity cut-off I max cut was defined as follows: for 300 < W ≤ 910 cm -1 , Range 0-300 cm -1 Range 300-700 cm -1 Range 700-1000 cm -1 Range 1000-2000 cm -1 (X 5000) (X 60) (X 20) -Strong and weak lines: at this stage, the major limitation when building the SF 6 line lists is the huge number of transitions that becomes hardly manageable for practical applications. The full line list for 32 SF 6 contains 6 billion lines when using the previously defined cut-offs. To make the problem tractable for the user, we have adopted the same strategy as for CF 4 34 and hot methane [START_REF] Rey | [END_REF] .

I max cut = 5 ×
This consists in extracting from the full list a smaller subset composed of strong and medium lines which are responsible for all sharp absorption features in the spectrum. In our calculations, this set of "strong" lines is composed of few millions of transitions that contribute to the largest absorption in spectra. They are selected by applying the condition I ≥ αI max cut where α is a cut-off factor. For this work, the value α = 1000 was chosen for W ≤ 300 cm -1 and α = 5000 otherwise. The remaining set is formed by an accumulation of several millions (and even billions in our case) of very weak overlapping lines, resulting in a smooth quasi-continuum (QC). Though the contribution of each individual weak line could be very small, the whole set plays a substantial role in the modelling of the overall opacity. 32 is also given in the range 925-960 cm -1 .

Validation and distribution

Figures 567gives examples of comparison between theoretical (including the CB, HB, and total CB+HB contributions shown in different colors) and experimental PNNL raw spectra [START_REF] Johnson | Proceedings of SPIE -The International Society for Optical Engineering[END_REF] which have been recently included in HITRAN as absorption cross-sections without lineby-line assignments. All simulations were made at a resolution of 0.1 cm -1 and a pressure of 760 Torr using a path length of 20 cm and a Voigt profile. The main absorption bands are indicated in these figures and the integrated intensities for the strongest absorption bands in the ν 4 and ν 3 regions of 32 SF 6 are given in Table 3 . Though a good agreement can be seen in Figs. 567, plotting residuals between the PNNL and the TheoReTS data is quite a tricky task. On the one hand, the line positions, in particular for hot bands, may be shifted and on the other hand, the PNNL line profile was difficult to reproduce, making thus difficult a proper comparison. Instead, accuracy of the different patterns can be established by a direct comparison of the absorption cross-sections of SF 6 in the regions displayed in Figs. 567. The results are summarized in the table 4. The overall agreement between theory and the observed PNNL spectra in the strongest absorption region around 947 cm -1 is confirmed by an error of about 4%. Assuming the value of the integrated intensity for ν 4 consistent in SHeCaSDa (see Tab. 3), half of the 12% of error between PNNL and TheoReTS around 615 cm -1 can be possibly The QC line lists was compressed to 1.81 millions of lines using the super-line approach.

explained by a deviation (about 2.9%) of our dipole moment first derivative ∂ µ/∂ q 4 and the other half by line strength errors for some hot bands. Simulation from SHeCaSDa is also displayed in the region of ν 4 (around 600 cm -1 ) and ν 3 (around 950 cm -1 ). This latter one does not reproduce the overall absorption shape of the room-temperature observed spectra due to the lack of HBs, as indicated by the tables 3 and 4 and corroborated by the table 5. Among the 500 cold and hot bands contained in the "strong" line list of 32 SF 6 , the integrated intensities for the 260 strongest ones can be found in Supplementary Material S 2 . For example, the most important HB ν 3 + ν 6ν 6 that represents half of the ν 3 contribution to the absorption is clearly lacking in SHeCaSDa and HITRAN (see table 3).

A detailed description of the final line lists for 32 SF 6 is given in Table 5. The number of transitions in the "strong lines subset" represents 0.25% of the full line list file for 32 SF 6 , providing 86% of the opacity. The remaining 5.999 billion lines of the QC were compressed using our super-line approach [START_REF] Rey | [END_REF] , permitting to reduce the number of lines by 3 orders of magnitude with a minimal loss of accuracy. The final compressed QC line lists is composed of 1.81 million super-lines. For theoretical simulation of absorption or emission cross-sections at a given temperature, the super-lines can be treated in the same manner as the true lines contained in the "strong line subset". A comparison with the recent SHeCaSDa database is also given in this table. As already seen in the figure 3 and confirmed the table 4, the lack of opacity in the ν 4 and ν 3 regions in SHeCaSDa is due to 80% of missing HBs. Both the strong and QC line lists are freely downloadable from the TheoReTS information system 47 (http://theorets.univ-reims.fr, http://theorets.tsu.ru).

Conclusion

To summarize, we present in this work the first global theoretical predictions of infrared rotationally resolved spectra for the sulfur hexafluoride molecule. The rotation-vibration energy levels are determined from variational calculations using our refined PES obtained by the optimisation of the ab initio one of Ref. 17 . Most importantly, absolute line intensities are calculated and extensive theoretical

Fig. 2

 2 Fig. 2 (a) Cut of the SF 6 PES 17 along the S-F 1 and S-F 2 bonds . (b) Reduced energy levels for ν 2 + ν 4 , 2ν 2 , ν 3 + ν 6 , ν 1 + ν 5 , ν 4 + 2ν 6 and ν 2 + ν 6 , listed by increasing energy order. (c) Portion of reduced energy levels for ν 3 + ν 6 and ν 1 + ν 5 . (d) Rotational errors using our optimized (r opt e ) and "best" (r best e , see text) equilibrium geometry values.

Fig. 3

 3 Fig.3Comparison of the lower state energies involved in the TheoReTS database47 (left panel, this work) and in HITRAN30 /SHeCaSDa 32 (right panel). The intensity cut-off was fixed to 10 -25 cm/molecule. For the definition of P i -P j , see text.

  empirical effective Hamiltonian24 . Figure2(d) displays the rotational errors ∆ rot = E emp. rot -E varia. rot up to J = 80. We clearly see in this figure that the discrepancy increases with J because r opt e

  10 -29 cm/molecule; for 910 < W ≤ 965 cm -1 , I max cut = 7 × 10 -28 cm/molecule and I max cut = 1 × 10 -29 cm/molecule otherwise. It is shown in Fig. 4 that the integrated intensities are nearly converged at J max = 121 using these cut-offs. Concerning 33 SF 6 and 34 SF 6 , since only strong lines are present in spectra with natural abundances, a constant intensity cut-off of 10 -26 cm/molecule was applied.

Fig. 4

 4 Fig. 4 Integrated intensities (in cm/molecule) for different spectral ranges at different values of the rotational angular momentum J.

  619.6474) and (945.5481, 618.2594) cm -1 for 33 SF 6 and 34 SF 6 , respectively. The other ω's remain unchanged during isotopic substitutions x SF 6 → y SF 6 .

				Opt. (TW) "Best" (TW)
	F 11	0.224671	0.220313	0.220314
	F 22	0.151637	0.151728	0.151731
	F 33	0.172397	0.169213	0.169213
	F 44	0.306367	0.295676	0.295672
	F 34	0.174623	0.163428	0.163423
	F 55	0.226574	0.217858	0.217858
	F 66	0.200283	0.192309	0.192034
	r e		1.557027	1.555953	1.555938
	ω 1		794.8578	787.1117	/
	ω 2 ω 3 ω 4	a a	653.0081 964.7491 633.3353	653.2106 963.7003 621.0293	/ / /
	ω 5		539.9829	529.5244	/
	ω 6		358.9892	351.5388	/
	a (ω 3 , ω opt 4 ) = (954.3495, opt

Table 2

 2 Selected vibrational band centers (in cm -1 ) for 32 SF 6 computed from variational calculations using our refined PES (this work, TW) and from empirical effective models20,[26][27][28][29]56 . Calculated band centers for 33 SF 6 and 34 SF 6 as well as the isotopic shifts ∆ (in cm -1 ) with respect to 32 SF 6 are also given (see text for explanations).

			32 SF 6		33 SF 6		34 SF 6
	Band	Empirical	Variational (TW) Emp.-Var. Variational (TW) ∆(32 → 33) Variational (TW) ∆(32 → 34)
	ν 6 (F 2u )	347.7367	347.7366	0.0001	347.7504	-0.0137	347.7635	-0.0269
	ν 5 (F 2g )	523.4490	523.4494	-0.0004	523.4705	-0.0215	523.4906	-0.0412
	ν 4 (F 1u )	614.9819	614.9821	-0.0002	613.6757	1.3062	612.3611	2.6210
	ν 2 (E g )	643.3737	643.3734	0.0003	643.3908	-0.0171	643.4073	-0.0339
	ν 1 (A 1g )	774.5453	774.5454	-0.0001	774.5721	-0.0268	774.5968	-0.0514
	ν 3 (F 1u )	948.1025	948.1024	0.0001	939.0272	9.0753	930.4646	17.6378
	ν 2 + ν 6 (F 2u )	989.1190	989.1653	-0.0463	989.2029	-0.0839	989.2324	-0.0671
	ν 2 + ν 6 (F 1u )	991.2740	991.3071	-0.0331	991.2641	0.0099	991.2401	0.0670
	ν 2 + ν 4 (F 1u )	1257.0849 1257.0190	0.0659	1255.7262	1.3587	1254.4199	2.5991
	ν 1 + ν 4 (F 1u )	1388.3820 1388.3192	0.0628	1387.0563	1.3257	1385.7712	2.5480
	ν 2 + ν 3 (F 1u )	1587.7173 1587.9089	-0.1916	1578.9063	8.8110	1570.3867	17.5222
	ν 2 + ν 3 (F 2u )	1591.1923 1591.0909	0.1014	1582.1352	9.0571	1573.6517	17.4392
	ν 1 + ν 3 (F 1u )	1719.7390 1719.8226	-0.0836	1710.8106	8.9284	1702.2811	17.5415
	2ν 3 (A 1g )	1889.0500 1889.3553	-0.3053	1871.3404	17.7096	1854.3259	35.0294
	2ν 3 (E g )	1891.6000 1892.0052	-0.4052	1873.9120	17.6880	1856.8083	35.1969
	2ν 3 (F 2g )	1896.5300 1896.5375	-0.0075	1878.4443	18.0857	1861.3746	35.1629
	ν 1 + ν 2 + ν 4 (F 1u ) 2028.1549 2027.9459	0.2090	2026.7878	1.3671	2025.5056	2.4403
	2ν 1 + ν 4 (F 1u )	2160.0257 2159.9263	0.0994	2158.3789	1.6468	2158.2626	1.6637
	2ν 1 + ν 3 (F 1u )	2489.4682 2489.8765	-0.4083	2480.5044	8.9638	2473.5020	16.3745

  1τ , A 2τ , E τ , F 1τ , F 2τ (τ = g and u), respectively. For a given block of symmetry Γ, we have thus included in the calculations the set of eigensolutions up to {E

	Γ(n=10) v	, Ψ Γ(n=10) v

1, 1.05, 1, 1, 1} were chosen, resulting in 32579, 31263, 63824, 91728, 93029, 29960, 30640, 60582, 94470 and 93775 functions for the symmetry blocks A

  Comparison of theoretical absorption spectra of 32 SF 6 + 34 SF 6 with experimental medium resolution PNNL records[START_REF] Johnson | Proceedings of SPIE -The International Society for Optical Engineering[END_REF] at room temperature the range 590-895 cm -1 . Contribution of the cold (CB) and hot (HB) bands are plotted and the strongest bands are indicated (the complete list can be found in Supplementary Material S 2 ). A comparison with the SHeCaSDa database32 is also given in the range 590-640 cm -1 . Comparison of theoretical absorption spectra of 32 SF 6 + 33 SF 6 + 34 SF 6 with experimental medium resolution PNNL records[START_REF] Johnson | Proceedings of SPIE -The International Society for Optical Engineering[END_REF] at room temperature the range 925-1005 cm -1 . Contribution of the cold (CB) and hot (HB) bands are plotted and the strongest bands are indicated (the complete list can be found in Supplementary Material S 2 ). A comparison with the SHeCaSDa database

	-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 ×10 -17 ×10 -18 5 PNNL spectra TheoReTS 32 SF 6 (CB+HB) TheoReTS 32 SF 6 (CB) TheoReTS 32 SF 6 (HB) TheoReTS 34 SF 6 (CB+HB) SHeCaSDa 32 SF 6 (CB+HB) 34 SF 6 (ν 4 ) SHeCaSDa 32 SF 6 (CB+HB) -4 PNNL cross-sections (cm 2 / molecule) -3 -2 -1 0 1 2 3 4 ν 3 +3ν 6 -3ν 6 ν 3 +2ν 6 -2ν 6 ν 3 +ν 5 +ν 6 -ν 5 +ν 6 ν 1 +ν 5 +ν 6 -2ν 6 ν 3 +ν 5 -ν 5 ν 3 +ν 6 -ν 6 ν 2 +ν 3 -ν 2 2ν 3 -ν 3 PNNL spectra TheoReTS 32 SF 6 + 33 SF 6 (CB+HB) TheoReTS 32 SF 6 (CB) TheoReTS 32 SF 6 (HB) TheoReTS 34 SF 6 (CB+HB) PNNL cross-sections (cm / molecule) 34 SF 6 33 SF 6	ν 4 +ν 5 -ν 5 ν 3 +ν 4 +ν 6 -ν 4 +ν 6	ν 4 2ν 4 -ν 4 +ν ν 6 -ν 6 4 4 +ν ν 5 +ν 6 -ν 6 5 -ν 2ν 4 +ν 6 -ν 4 -ν 6 6 2ν 2 -ν 6 ν 3 +ν 4 -ν 4 ν 2 +ν 3 +ν 6 -ν 2 +ν	ν 4 +2ν 6 -2ν 6	ν 3
	590 925 -5 -8 -6 -4 -2 0 2 4 6 8 PNNL cross-sections (cm 2 / molecule) ×10 -20 ×10 -19 2 PNNL spectra 595 600 930 PNNL spectra TheoReTS 32 SF 6 (CB+HB) 605 935 ν 5 +3ν 6 -2ν 6 2ν 5 +2ν Wavenumber (cm 610 615 620 Wavenumber (cm 940 -1 ) -1 ) ν 5 +ν 6 ν 5 +2ν 6 -ν 6 2ν 5 +ν 6 -ν 5 ν 4 +ν 5 +ν 6 -ν 4 6 -ν 6 5 -ν TheoReTS 32 SF 6 (CB) TheoReTS 32 SF 6 (HB) TheoReTS 34 SF 6 (CB+HB) 34 SF 6 -1 0 1 PNNL cross-sections (cm 2 / molecule) ν 2 +ν 5 +ν 6 -ν 5 ν 2 +2ν 6 -ν 6 ν 6 TheoReTS 32 SF 6 (CB+HB) 2 +ν 2ν 2 +ν 6 -ν 2 ν 2 +ν 4 +ν 6 -ν 4 2ν 2 +2ν 6 -ν 2 -ν 6 ν ν 2 +ν 2 +4ν 6 -3ν 6 4 +2ν 6 -ν 6 TheoReTS 32 SF 6 (CB) 4 -ν TheoReTS 32 SF 6 (HB)	625 945	630	950	635	640
	970 -2	850	855 975	860 980	865 985	870	875 990	880 995	885 1000	890	1005
					Wavenumber (cm Wavenumber (cm -1 -1 ) )
	Fig. 5 1-23 | 15

Fig.

6

Table 3

 3 Integrated intensities ∑ S i j (×10 17 cm/molecule) for the strongest absorption bands in the ν 4 and ν 3 regions of 32 SF 6 and comparison with the SHeCaSDa database 32 . + ν 4 + ν 6ν 2 + ν 6 0.017 / 948.393 ν 2 + ν 3 + ν 6ν 2 + ν 6 0.188 / 614.355 2ν 4 + ν 5ν 4 + ν 5

	ν 4 region

Table 4

 4 Sum of PNNL cross-sections (cm 2 /molecule) for SF 6 and comparison with TheoReTS (this work) and SHeCaSDa 32 .

	Region (cm -1 ) TheoReTS (TW)	PNNL	Error (%)
	560-650	1.96×10 -16	2.23×10 -16	12.1
	820-910	1.46×10 -17	1.30×10 -17	-11.0
	910-965	3.06×10 -15	3.20×10 -15	4.4
	965-1005	4.07×10 -17	4.01×10 -17	-1.5
	1005-1500	7.39×10 -18	7.65×10 -18	3.4
	1550-1730	3.16×10 -17	3.47×10 -17	8.9
	Region (cm -1 )	SHeCaSDa 32	PNNL	Error (%)
	560-650	0.86×10 -16	2.23×10 -16	61.4
	910-965	1.02×10 -15	3.20×10 -15	68.4

Table 5

 5 Summary of the final line lists for 32 SF 6 and comparison with the SHeCaSDa database32 . The number of strong and quasi-continuum (QC) lines (#lines) as well as the integrated intensities (∑lines) are indicated. Here, CB=cold bands, HB=hot bands.

	Strong lines		QC	SHeCaSDa	32	Error a
	Range / Type #lines	∑lines	#lines	∑lines	#lines	∑lines
	×10 6	×10 -16	×10 9	×10 -16	×10 5	×10 -16
	300-910 CB 0.0426 3.828×10 -2	0.0005	4.273×10 -5 1.3165 3.942×10 -2	-2.8%
	HB 2.7963 7.337×10 -2	1.6862	1.718×10 -2 0.3851 2.201×10 -2	75.6%
	910-965 CB 0.0231 0.615	0.0002	3.170×10 -4 0.5310 0.553	10.1%
	HB 2.9370 1.008	0.9543	0.253	1.6568 0.175	86.1%
	965-3000 CB 0.6444 1.443×10 -2	0.0217	2.650×10 -4 /	/	/
	HB 7.8830 2.254×10 -2	2.9239	1.213×10 -2 /	/	/
	Total (CB+HB) 15.470 1.771	5.9991 b 0.283	3.3584 0.789	/
	a ∑ (Strong+QC-SHeCaSDa) / ∑ (Strong+QC) ×100.

b

Acknowledgements

The authors acknowledge support from the Romeo computer center of Reims-Champagne-Ardenne and of French-Russian collaboration program LIA "SAMIA". The support of RFBR grant (Grant No. 20-35-90075) is also acknowledged.

molecular line lists are constructed for the three most abundant species 32 SF 6 , 33 SF 6 and 34 SF 6 using ab initio dipole moment functions.

Room-temperature calculated spectra match observed ones with unprecedented accuracy for ro-vibrational patterns of seven-atomic molecules. This is, in particular, the case of the most important range in the atmospheric windows around 10 µm where the error between the sum of absorption cross-sections of PNNL and TheoReTS is about 4%, against 68% between PNNL and the most recent analyses. This work permitted to elucidate for the first time all complicated cold and hot band structures in the rotationally resolved room-temperature SF 6 spectra. For applications requiring accurate opacity calculations and its temperature dependence, our complete theoretical line lists should be preferred to partial lists included in the HITRAN or SHeCaSDa databases.


Oleg Egorov
email: oleg.egorov@iao.ru

Michaël Rey

Andrei V Nikitin

Dominika Viglaska

New ab initio potential energy surfaces for NH3 constructed from explicitly correlated coupled cluster methods

Keywords: NH3, ab initio, potential energy surface, CCSD(T)-F12, CCSDT, CCSDTQ, CCSDTQP, aug-cc-pCV6Z methods

Accurate calculation of ab initio potential energy surfaces (PES) for the NH3 molecule is a difficult task because of the poor convergence of the standard CCSD(T) method with respect to the basis size. Currently, the best available NH3 PESs contain empirically refined parameters. In this paper, we show that CCSD(T)/aug-cc-pCV6Z calculations are not sufficient to properly describe the PES over a large range of nuclear configurations. However, the PES obtained by the extrapolation of the CCSD(T)/aug-cc-pCVXZ (X=T, Q, 5, and 6) energies to the complete basis set limit is closer to that based on the explicitly correlated CCSD(T)-F12a method using the cc-pCVQZ-F12 orbital basis set. All the ab initio PESs constructed in this work include the following corrections: one electron relativistic effects, diagonal Born-Oppenheimer correction and high-order electronic correlations (CCSDT, CCSDTQ and CCSDTQP). Finally, the root mean square deviation between the predicted band centers obtained from our final "pure" ab initio PES and the experimental ones in the spectral region [0-7000] cm -1 is divided by two compared to the most accurate ab initio PES available in the literature.

INTRODUCTION

Although the use of empirically-refined potential energy surfaces (PESs) allows computing some energy levels close to the spectroscopic accuracy, there are at least two reasons for which the development of pure ab initio PESs is necessary: a) ab initio PESs generally better extrapolate energy levels beyond the range of observed data, making more consistent the prediction of strong interacting bands or of hot bands for high temperature calculations; b) testing and searching for accurate ab initio methods to build new PESs is very stimulating and allows checking consistency those PESs already published.

Currently, the single-and double-excitation coupled-cluster method with a perturbative treatment of triple excitations [CCSD(T)] is widely used for constructing the PESs of molecules below their dissociation energies. In early works, this method was combined with the augmented correlation-consistent polarized basis sets of triple (aug-cc-pVTZ) or quadruple (aug-cc-pVQZ) zeta quality. However, the errors on the energy levels may reach 20 cm -1 , and even higher, with respect to the NH3 experimental band origins. [START_REF] Rajamäki | Vibrational energy levels for symmetric and asymmetric isotopomers of ammonia with an exact kinetic energy operator and new potential energy surfaces[END_REF] More advanced calculations were based on larger basis sets, including the core-valence correlation effects and extrapolation to the Complete Basis Set (CBS) limit. Huang et al. [START_REF] Huang | An accurate global potential energy surface, dipole moment surface, and rovibrational frequencies for NH3[END_REF] obtained CCSD(T) energies extrapolated to the CBS limit using the aug-cc-pCVXZ (X=T, Q, and 5) basis sets. Although diagonal Born-Oppenheimer breakdown corrections (DBOC), relativistic and high-order correlation effects were included in ref 2, the deviations from the experimental band origins were 2.5 cm -1 in average. Such deviations turn out (i) smaller than those obtained from the previous ab initio PES of Yurchenko et al. [START_REF] Yurchenko | Potential-energy surface for the electronic ground state of NH3 up to 20000 cm -1 above equilibrium[END_REF] where only the corrections for the corevalence correlation and relativistic effects were considered, but (ii) bigger than the work of Rajamӓki et al. [START_REF] Rajamäki | High excitations in coupled-cluster series: vibrational energy levels of ammonia[END_REF] In the work by Polyansky et al., [START_REF] Polyansky | Calculation of rotation-vibration energy levels of the ammonia molecule based on an ab initio potential energy surface[END_REF] the multireference configuration interaction (MRCI) method was considered with the core-valence aug-cc-pCVQZ and aug-cc-pCV5Z basis sets and basis set extrapolation to improve the calculation of highly-excited vibrational states of NH3. More than ten possible ab initio contributions to the dissociation energy were analyzed and the mean errors on the band origins was of 1 cm -1 for energies up to 7000 cm -1 . Marquardt et al. [START_REF] Marquardt | Global Analytical Potential Energy Surface for the Electronic Ground State of NH3 from High Level ab Initio Calculations[END_REF] have constructed a global analytical PES of NH3 by combining both the CCSD(T) and MRCI methods: the CCSD(T) energies were extrapolated to the CBS limit using the aug-cc-pVXZ (X=T, Q, and 5) basis sets as in ref 3 and the MRCI energies were calculated using the aug-cc-pVQZ basis set. The final PES was used to compute the band origins of four isotopologues but the Obs.-Calc. deviation up to 5000 cm -1 was of 11 cm -1 on the average for the main isotopologue.

In ref 4, authors showed for NH3 that the explicitly correlated CCSD(T)-R12 method [START_REF] Kutzelnigg | Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. I. General theory[END_REF][START_REF] Noga | Coupled cluster theory that takes care of the correlation cusp by inclusion of linear terms in the interelectronic coordinates[END_REF] correctly describes the basis set incompleteness error. According to ref 4, the symmetric part of the NH3 PES including the radial and inversion modes converges very slowly with respect to the basis set size. For this reason, the CCSD(T)-R12 method could provide more reliable results than those obtained from the standard CCSD(T) method using a large basis set size. Another conclusion of ref 4 concerns the electronic correlation contribution which is important for levels up to five quanta of excitation.

In this work, we consider only "pure" ab initio PESs of NH3 i.e. without any empirical corrections. To achieve this end and unlike ref 4, all the corrections beyond the CCSD(T) level were systematically included for the six vibrational degrees of freedom. Here, the more advanced ab initio method for the DBOC calculations have been performed on a grid composed of ~40000 reference ab initio points. The Root Mean Square (RMS) deviation between the calculated and observed band origins is two times less than that obtained in refs 4 and 5. Section 2 gives all details about the construction of the final PES where various ab initio methods and basis sets, the analytical form of the PES as well as the grid of the reference ab initio points are discussed. The validation of our final PES is given in Section 3 before concluding in Section 4.

Concerning the nuclear motion calculations, our variational computer code TENSOR [START_REF] Rey | Ab initio ro-vibrational Hamiltonian in irreducible tensor formalism: a method for computing energy levels from potential energy surfaces for symmetric-top molecules[END_REF][START_REF] Rey | Complete nuclear motion Hamiltonian in the irreducible normal mode tensor operator formalism for the methane molecule[END_REF][START_REF] Rey | Group-theoretical formulation of an Eckartframe kinetic energy operator in curvilinear coordinates for polyatomic molecules[END_REF] has been adapted to molecules exhibiting one large amplitude motion. [START_REF] Viglaska | Derivation of ρ-dependent coordinate transformations for nonrigid molecules in the Hougen-Bunker-Johns formalism[END_REF] This code was initially designed to the treatment of rigid molecules which was applied with success to the construction of accurate spectroscopic line lists (see e.g. CH4, [START_REF] Rey | Ab initio variational predictions for understanding highly congested spectra: rovibrational assignment of 108 new methane subbands in the icosad range (6280-7800 cm -1 )[END_REF][START_REF] Rey | TheoReTS -An information system for theoretical spectra based on variational predictions from molecular potential energy and dipole moment surfaces[END_REF][START_REF] Rey | Accurate Theoretical Methane Line Lists in the Infrared up to 3000 K and Quasi-continuum Absorption/Emission Modeling for Astrophysical Applications[END_REF] CF4, [START_REF] Rey | Understanding global infrared opacity and hot bands of greenhouse molecules with low vibrational modes from first-principles calculations: the case of CF4[END_REF] PH3, 17 NF3, [START_REF] Egorov | Global modeling of NF3 line positions and intensities from far to mid-infrared up to 2200 cm -1[END_REF] C2H4, [START_REF] Viglaska-Aflalo | A global view of isotopic effects on rovibrational spectra of six-atomic molecules: a case study of eleven ethylene species[END_REF] and SF6 [START_REF] Rey | Towards a complete elucidation of the ro-vibrational band structure in the SF6 infrared spectrum from full quantum-mechanical calculations[END_REF] ). Concerning the ab initio calculations, the MOLPRO, [START_REF] Werner | The Molpro quantum chemistry package[END_REF][START_REF] Werner | version 2019.2, a package of ab initio programs[END_REF] MRCC, [START_REF] Kállay | The MRCC program system: Accurate quantum chemistry from water to proteins[END_REF][START_REF] Kállay | a quantum chemical program suite[END_REF] and CFOUR [START_REF] Matthews | Coupled-Cluster Techniques for Computational Chemistry: the CFOUR Program Package[END_REF][START_REF] Stanton | a quantum chemical program package[END_REF] packages have been considered in this work.

CONSTRUCTION OF THE POTENTIAL ENERGY SURFACE

Choice of the ab initio methods and basis sets. Since the symmetric part of the NH3

PES is very sensitive to the choice of the basis set, one way to study the convergence is to analyze the height of the inversion barrier. In the following, the results of ref 4 computed with the CCSD(T)-R12 explicitly correlated method will be considered as the reference value. The R12 [START_REF] Kutzelnigg | Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. I. General theory[END_REF][START_REF] Noga | Coupled cluster theory that takes care of the correlation cusp by inclusion of linear terms in the interelectronic coordinates[END_REF] and F12 [START_REF] Knizia | Simplified CCSD(T)-F12 methods: Theory and benchmarks[END_REF] families of explicitly correlated methods include explicitly the inter-electronic coordinate in the wave function, allowing to reach the complete basis set values without extrapolations. Two of the most popular simplified versions, called F12a and F12b, [START_REF] Knizia | Simplified CCSD(T)-F12 methods: Theory and benchmarks[END_REF] are implemented in the MOLPRO package. In order to make consistent the comparison with ref 4, the correlation of all electrons was necessary. For this reason, the CCSD(T) method combined with the aug-cc-pCVXZ (X=T, Q, 5, and 6) basis sets (hereafter, ACVTZ, ACVQZ, ACV5Z, and ACV6Z) are considered for this work. Concerning the CCSD(T)-F12q (q=a, b) methods, the cc-pCVQZ-F12 basis set (hereafter, CVQZ-F12) was used. It is currently the bigger basis implemented in MOLPRO for the explicitly correlated methods for which the core-valence correlation is taken into account.

To extrapolate the energies to the CBS limit, we used the following functional

  exp X CBS E E A C X      , (1)
which is available in MOLPRO, where X E and CBS E are the total energies (Hartree-Fock(reference)

+ correlation) corresponding to the ACVXZ basis set and CBS limit. The choice (1) differs from the inverse powers functional used in refs 2 and 3 but is in a better agreement with that of Klopper et al. [START_REF] Klopper | Equilibrium inversion barrier of NH3 from extrapolated coupled-cluster pair energies[END_REF] According to our tests for NH3, the converged values obtained from the inverse powers functional are below those provided by the explicitly correlated methods. This is probably one of the reasons why the ab initio band origins in refs 2 and 3 are systematically higher than the observed values (see in Table 3). 4). The result for CCSD(T)-R12 was taken from ref 4. Note, here the "CBS" results come from the extrapolation of the total CCSD(T)/ACVXZ energies (given in parenthesis in cm -1 ). The results based on the extrapolation of the correlation energy are given in supplementary material I. The NH3 PES has two stationary points: the first one corresponds to the global minimum at the pyramidal (C3v) geometry while the second one corresponds to the saddle point at the planar (D3h)

geometry. The height of the inversion barrier can be determined from these two geometries. Although they can be optimized in MOLPRO, we preferred optimizing the geometries using the constructed PESs.

According to Figure 1 and Table 1, the height of the inversion barrier computed at the CCSD(T)/ACV6Z level of the theory is about 4 cm -1 higher than the reference value provided by the explicitly correlated CCSD(T)-R12 method. A good agreement between the F12b method with the CVQZ-F12 basis set and CCSD(T)/ACV6Z is confirmed by the comprehensive study conducted in ref 29. On the contrary, The F12a method provides a value in excellent agreement with the CBS limit taken from ref 4.

The three points extrapolation scheme, which involves the CCSD(T) energies computed with the ACVXZ (X=T, Q, and 5) basis sets, allows approaching the reference value of the barrier within 1.7 cm -1 using the functional (1). Another scheme, which additionally includes the CCSD(T)/ACV6Z energy, turns out to be very promising since the agreement is within 0.5 cm -1 in this case. An exclusion of the CCSD(T)/ACVTZ energies does not allow improving the results: the CBS(Q56Z) extrapolation gives a value that is 1.3 cm -1 higher than the reference one. As it can be seen in Section 3.2 of the present paper, the convergence test presented in this section mainly explains the final agreement between the calculated and observed band origins of NH3.

Note that by default MOLPRO extrapolates the correlation energy. From Table S1 of Supporting File I, the reader, however, can see that the agreement with the experimental data is slightly better in terms of RMS when the total energies are extrapolated instead of the correlation ones (RMSr+c=0.480 cm -1 and RMSc = 0.522 cm -1 for the CCSD(T)/CBS(TQ56Z) PES). Although a difference of 0.042 cm -1 is not big and clearly indicates that the Hartree-Fock energies are well converged at the CCSD(T)/ACV6Z level of the theory, we will consider below the results based on the extrapolation of the total energies.

Analytical form of the PES.

The NH3 molecule has six vibrational degrees of freedom where there are two symmetric and two doubly degenerate asymmetric vibrations. For a full account of symmetry, a set of suitable coordinates has to be introduced. Actually, the symmetry coordinates are preferred to the internal ones because each coordinate transforms as an irreducible representation of the group. [START_REF] Papoušek | Molecular Vibrational-Rotational Spectra. Theory and Applications of High Resolution Infrared[END_REF] For NH3, such a set of coordinates can be built with the help of the projection technique. [START_REF] Bunker | Molecular Symmetry and Spectroscopy[END_REF] For stretching and bending vibrations, the symmetry coordinates were defined as the following linear combinations:

         
                            , ( 2
)
where the Morse-type functions for stretching vibrations

  1 exp k k e y b r r        , (3)
and cosine functions for bending ones

    cos cos k k e      , (4)
with k = 1, 2, and 3; b = 1.9 Å -1 , were used for a correct asymptotic behavior. Concerning the inversion coordinate  , the sine function is a good choice because it equals to one at the planar configuration

( o 90  )     1 2 sin sin A e S      , (5)
where

    1 2 3 2 sin sin / 6 3          . (6)
Finally, the analytical form of the PES can be expressed as a linear combination of irreducible

tensor operators p j R   max 11 1 2 3 3 4 4 0 , , , , , p AAE E E E p a b a b j j pj V S S S S S S C        R (7)
where j C are the expansion coefficients to be determined. Since the NH3 PES must be totally symmetric with respect to the symmetry operations of the D3h group, each irreducible tensor has to transform as 1 A and it can be written as a power series expansion in the symmetry coordinates:

ii i i i i n pp p p p p AA p E E E E j i a b a b i T S S S S S S                                       R . ( 8
)
Here, the i T 's are the D3h Clebsch-Gordan coefficients and each term in Eq. ( 8) is of order p. In this work, the polynomial (7) was expanded up to max 8 p  , resulting in 590 irreducible tensor operators. However, only 305 tensors were included in the final fits of the PESs. It is worth mentioning that due to the use of the symmetry coordinates, the equilibrium geometry only depends on the two-dimensional symmetric part of the PES. points with a maximum energy of 14000 cm -1 for the coupled coordinates.

Construction of the

To properly fit the ab initio points, we have considered a weight function similar to that previously used for CH4 [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF] where the weights decrease asymptotically from 1.0 to 0.9 for energies in the range of 0-8000 cm -1 , and then lower to 0.5 and 0.1 at 10000 and 15000 cm -1 , respectively. For all energies above 20000 cm -1 , a fixed weight of 0.001 was applied.

Since the ab initio calculations using the ACV6Z basis set are quite demanding (up to 2 (resp.

7) hours for one point using the Cs (resp. C1) symmetry), our strategy was to build another fulldimensional but smaller grid composed of 2235 points (called Grid_2) and used to compute the energy difference between the CCSD(T) calculations made with the ACV6Z and ACV5Z basis sets.

It turns out that the behavior of CCSD(T)/[ACV6Z-ACV5Z] is quite smooth with an energy difference going from -412 cm -1 to -485 cm -1 and a mean value of -434 cm -1 . Thus, only 40 irreducible tensors in the expansion (7) were necessary to fit this difference with an RMS error less than 0.08 cm -1 .

Finally, the correction factor CCSD(T)/[ACV6Z-ACV5Z] was applied to the whole Grid_1 leading to an extrapolated CCSD(T)/ACV6Z PES. The 305 irreducible tensors of the expansion (7) reproduced the Grid_1 of the reference ab initio points with the RMS error less than 0.5 cm -1 .

Corrections beyond the CCSD(T) method.

With the increase of the computer capabilities, it is now common to take into account corrections beyond CCSD(T) for computing the vibrational band origins within 1 cm -1 of accuracy. Generally, there are three types of corrections:

one electron relativistic effects (Relativistic), DBOC, and high-order electronic correlation effects (in this work CCSDT, CCSDTQ, and CCSDTQP). In the case of NH3, all these corrections were analyzed in details in ref 4. Therefore, in the following we shall mainly focus on the differences between this work and that of Rajamäki et al. [START_REF] Rajamäki | High excitations in coupled-cluster series: vibrational energy levels of ammonia[END_REF] Unlike ref 4, our aim was to apply these three corrections to the six degrees of freedom of NH3. In particular, the energies at the CCSDT/cc-pVTZ and CCSDT(Q)/cc-pVTZ levels of the theory were calculated for the whole Grid_1 with the help of the MRCC code. The calculations of the Relativistic and DBOC contributions were also performed for Grid_1. In this work, the CFOUR code was used to compute DBOC from the CCSD method and the cc-pVTZ basis set while only HF/cc-pVDZ was considered in ref 4. Concerning the relativistic corrections in MOLPRO, the second-order Douglas-Kroll-Hess Hamiltonian was used in conjunction with the ACVTZ-DK basis. It turns out that the "DK" extension in the name of the basis set, which is important to optimize the inner-core orbitals, [START_REF] Werner | The Molpro quantum chemistry package[END_REF] was omitted in ref 4, resulting to a difference of +1 cm -1 on the inversion barrier height: 19 cm -1 versus 18 cm -1 (see in Section 3.3).

The smaller grid (Grid_2) was considered for the two most time demanding corrections, namely CCSDTQ and CCSDTQP. The ab initio calculations at the CCSDTQ/cc-pVTZ level took 3 (resp. 13) hours for one point in the Cs (resp. C1) symmetry. Since the convergence of high-order coupled-cluster methods is quite fast with the respect to the basis set size, the CCSDTQP method was combined here with the cc-pVDZ basis set in MRCC, as in ref 4. At this stage, it is worth mentioning that bigger basis sets require much more computational effort when considering higher order coupled cluster methods. For example, the cost of calculation using CCSDT(Q)/cc-pVQZ is about 11 times more demanding than that using CCSDT(Q)/cc-pVTZ employed in this work for the same nuclear configuration. For one point at the CCSDTQ(P)/cc-pVTZ level, up to 42 hours are necessary (about 2 days), against 10 minutes using CCSDTQ(P)/cc-pVDZ. A good compromise would be to decrease drastically the number of points in the grid but this may lead to a wrong estimation of the contribution of the correction because of the coupling between the vibrational modes. Undoubtedly, the complete construction of the PES using high-order electronic correlation methods will become possible in a near future.

The scaling factor from CCSD(T) to CCSDT has been evaluated by computing the [CCSDT-CCSD(T)]/cc-pVTZ energy difference which was equal to -46 cm -1 on the average. It turns out that the CCSDT(Q) correction gives the biggest contribution to the inversion barrier height with a shift of +3 cm -1 . The [CCSDT(Q)-CCSDT]/cc-pVTZ energy difference was equal to -76 cm -1 on the average while the mean difference [CCSDTQ-CCSDT(Q)]/cc-pVTZ was found to be +3 cm -1 on the average with a smaller impact on the barrier height (see in Section 3.3). Concerning the CCSDTQ→CCSDTQP scaling, the energies were calculated using the cc-pVDZ basis set. The energy differences going from -2 cm -1 to -11 cm -1 have been successfully fitted with an RMS error of 0.04 cm -1 using only 50 tensor parameters up to 4 th order. There are less experimental data concerning the cold bands falling in the spectral region 4000-6000 cm -1 . In particular, the combination bands v1+v2, v2+v3, and v2+2v4 were studied in refs 72 and 73 from Fourier transform spectra. Brown et al. [START_REF] Brown | Empirical line parameters of NH3 from 4791 to 5294 cm -1[END_REF] analyzed the v1+v4 and v3+v4 regions that are important for remote sensing of planetary atmospheres. Čermák et al. [START_REF] Čermák | Spectroscopy of 14 NH3 and 15 NH3 in the 2.3 μm spectral range with a new VECSEL laser source[END_REF] reported new experimental spectra of the v1+v2/v2+v3 bands based on new laser sources. However, the very congested spectra in the high-energy region of NH3 make analysis very difficult using effective Hamiltonians, as discussed in ref 76.

Around 6600 cm -1 , the spectra contain many blended lines that are difficult to resolve, even by Doppler-limited spectroscopy. Lundsberg-Nielsen et al. [START_REF] Lundsberg-Nielsen | Analysis of the High-Resolution Spectrum of Ammonia ( 14 NH3) in the Near-Infrared Region, 6400-6900 cm -1[END_REF] assigned the experimental transitions of the v1+v3 and 2v3 bands from the room-temperature Fourier spectra. Some lines of these bands were also measured in ref 78 using the external cavity tunable diode laser to provide the higher sensitivity.

Berden et al. [START_REF] Berden | Cavity-enhanced absorption spectroscopy of the 1.5 μm band system of jet-cooled ammonia[END_REF] applied the cavity-enhanced absorption technique combined with the slit-nozzle expansion for rotationally cooling the spectra. As a result, some lines of the v1+2v4 band were additionally detected. Li et al. [START_REF] Li | External cavity tunable diode laser spectra of the ν1 + 2ν4 stretch-bend combination bands of 14 NH3 and 15 NH3[END_REF] In a more recent work, [START_REF] Földes | Low-temperature high-resolution absorption spectrum of 14 NH3 in the ν1+ν3 band region (1.51 μm)[END_REF] the v1+v3 band was studied from the jet-cooled cavity ring-down (CRDS) and cavity enhanced absorption (CEAS) spectroscopy. Twagirayezu et al. [START_REF] Twagirayezu | Quadrupole splittings in the near-infrared spectrum of 14 NH3[END_REF] recorded sub-Doppler, saturation dip, spectra of the v1+v3, v1+2v4, and v3+2v4 bands with partially resolved hyperfine splittings due to the nitrogen quadrupole coupling.

A large number of experimental transitions of NH3 has been collected via the MARVEL (Measured Active Rotational-Vibrational Energy Levels) spectroscopic network protocol. In the most recent release, Furtenbacher et al. [START_REF] Furtenbacher | Empirical rovibrational energy levels of ammonia up to 7500 cm -1[END_REF] applied MARVEL with 28 new experimental data sources below 11 7500 cm -1 , resulting in the most reliable list of experimental band origins for NH3. This list is used below in Tables 2 and3 to validate our vibrational band centers calculated from our ab initio PESs. The original ab initio PESs (before adding the corrections) were obtained with the CCSD(T) and CCSD(T)-F12q (q = a, b) methods and with the basis sets presented in the Table. According to Table 2, the largest deviations from the CCSD(T)/ACV6Z PES correspond to the fundamental and overtone frequencies of the inversion mode. The CCSD(T)-F12b/CVQZ-F12 PES produces similar results which are confirmed by our convergence test in Figure 1 and Table 1.

These two methods give an inversion barrier height that is about 4 cm -1 higher than the values provided by the CCSD(T)-R12, 4 CCSD(T)-F12a/CVQZ-F12, and CCSD(T)/CBS(TQ56Z) methods.

Although the deviations of the fundamental frequencies for the other modes are within 1 cm -1 , our two PESs based on CCSD(T)/CBS(TQ56Z) and CCSD(T)-F12a/CVQZ-F12 are about 2-3 times more accurate than the best available ones.

When comparing the results obtained from the CCSD(T)/CBS(TQ56Z) and CCSD(T)-F12a/CVQZ-F12 PESs, one can see the levels where the inversion overtone is involved are less accurate when using the first PES. Conversely, the energy levels near 6600 cm -1 are in better agreement with observation than for the CCSD(T)-F12a/CVQZ-F12 PES. Finally, the RMS error using these two PESs is rather small (within 0.1 cm -1 ). Since the calculations of the ab initio points using the CCSD(T)/ACVXZ (X=T, Q, 5, and 6) methods are about 300 times more demanding, the use of CCSD(T)-F12a/CVQZ-F12 is a good compromise and turns out more appropriate for the construction of the full NH3 PES with much less effort within 1 cm -1 of accuracy. The CCSD(T)/CBS(Q56Z) PES provides an intermediate result in terms of RMS. From Table 2, it is seen that the frequencies of the inversion mode become more accurate using the CBS(Q56Z) extrapolation in comparison with the CCSD(T)/ACV6Z PES but they are less accurate compared with CCSD(T)/CBS(TQ56Z) calculations.

According to the benchmarks, [START_REF] Knizia | Simplified CCSD(T)-F12 methods: Theory and benchmarks[END_REF] the F12a method gives slightly overestimated values of the correlation energy whereas the F12b methodthe underestimated ones. However, there is no reliable information on using the F12a method in a combination with the core-valence correlation CVQZ-F12 basis set. Results of this work indicate that the "overestimation property" of the F12a method is suitable for NH3 because the standard CCSD(T) method has poor convergence with respect to the basis set size for this molecule.

There are only five bands in Table 2 for which the band centers calculated with the CCSD(T)-F12a/CVQZ-F12 PES have a deviation bigger than 1.0 cm -1 . Among these, four bands are associated with the third overtone of the inversion mode, namely [4v2( 1 A , 2 A ) and 4v2+v4 ( E , E )] while the last one is associated with the inversion-bending combination band 2v2+3v4( E ). According ref 85, all these bands, with the exception of 4v2( 1 A ), are related with only 1 or 2 rovibrational levels in the MARVEL database, making their accuracy quite poor compared to v2+v3 ( E , E ), for example, whose band origins were confirmed by more than 100 rovibrational observed levels.

As it can be seen in Table 3, are less than their critical values (of 0.02 and 0.05) that confirms the reliability of coupled cluster method. For further details about the T1 and D1 diagnostics, the original works [START_REF] Lee | A diagnostic for determining the quality of single-reference electron correlation methods[END_REF][START_REF] Janssen | New diagnostics for coupled-cluster and Møller-Plesset perturbation theory[END_REF] can be suggested.

Height of the inversion barrier: contribution of the corrections. The inversion barrier

separates two wells in the potential of NH3 and for a correct calculation of the inversion splitting, one needs to accurately estimate its height during the construction of the PES. The influence of different corrections on the equilibrium and saddle point geometries as well as on the inversion barrier height is given in Table 4. Although our analysis is based on the CCSD(T)-F12a/CVQZ-F12 PES, such corrections give similar contributions for the other PESs that can be seen in Tables S2-S7 of Supporting File I. Details of the calculation for the three corrections are given in Section 2.4.

According to Table 4, the relativistic correction shifts the barrier height by +19.39 cm -1 while its contribution was estimated at +18.1 cm -1 in ref 4. We have checked that the difference of about +1 cm -1 was caused by omitting the "DK" extension in the basis set name in ref 4. Our relativistic corrections using the ACVTZ basis set gave a contribution of +17.7 cm -1 to the barrier height. Since the "DK" extension is of key importance for optimizing the orbitals when using the relativistic Douglas-Kroll-Hess Hamiltonian, our test calculations of the band origins have demonstrated the less good agreement with the MARVEL database.

The DBOC corrections lower the barrier height by 10 cm -1 which is quite consistent with refs In this work, we have calculated the CCSDT(Q) energies for the whole Grid_1 composed of 40148 nuclear configurations. As it can be seen in the last column of Table 4, the difference between the CCSDT(Q) and CCSDTQ corrections is less than +0.2 cm -1 compared to the contribution of +1.4

cm -1 when scaling from CCSD(T) to CCSDT. The final contribution of the high-order electronic correlations to the barrier height is of +5.03 cm -1 . Note that the excitations of the quintuple order

[CCSDTQP] shift the barrier by +0.40 cm -1 against +1.0 cm -1 as estimated in ref 4. This difference could be explained by a larger number of points in our grid, the basis set size being similar in both cases.

Effect of the corrections on the band origins.

As expected, the three corrections considered in this work have played a major role in the prediction of the band origins associated with the radial and inversion vibrations (see Table 5). The contribution of the relativistic and DBOC corrections is up to 9.31 cm -1 and -4.80 cm -1 in the 2v2( 1 A ) band while that due to the electronic correlations has a bigger impact on the radial vibrations (up to -4.66 cm -1 and -4.9 cm -1 in the v1( 1 A

) and v3( E ) bands) and becomes larger for combination bands (up to -5.92 cm -1 and -5.58 cm -1 in the v1+v4( E ) and v3+2v4( E ) bands). Given values are the differences in the band origins (in cm -1 units) when the corresponding correction is taken into account. For example, the "CCSDT" column means the difference in the band origins calculated at the [CCSD(T)- when scaling to the CCSDTQP level from CCSDT(Q)→CCSDTQ. On the other hand, if the electronic correlations of the quadruple order would have been calculated directly using the CCSDTQ method without their estimations by the perturbation theory [CCSDT(Q)] then the CCSDTQ→CCSDTQP scaling would have lowered the RMS by 0.2 cm -1 , which is considerable.

SUMMARY AND CONCLUSION

In this work, five new ab initio potential energy surfaces (PESs) for NH3 have been constructed and tested. To reach an accuracy of 1 cm -1 in the prediction of the band origins of NH3 located in the 0-7000 cm -1 region from the reliable CCSD(T) method, the relativistic, DBOC and high-order electronic correlation effects have been added. Despite all these corrections, the largest aug-cc-pCV6Z basis set was not enough to obtain a good agreement with the experimental data tabulated in the MARVEL database. This was partly demonstrated by the convergence test shown in Figure 1. Note that the results using the explicitly correlated CCSD(T)-F12b method with the CVQZ-F12 orbital basis set were close to those obtained from CCSD(T)/aug-cc-pCV6Z. This last statement will be valuable when developing the PESs of large polyatomic molecules, for which the use of the aug-cc-pCV6Z basis set is not feasible.

It was also shown that for a correct extrapolation to the complete basis set (CBS) limit, the energies calculated with the aug-cc-pCV6Z basis set have to be included through the extrapolation functional. This was confirmed by our CCSD(T)/CBS(TQ56Z) PESs from which the RMS deviation with the MARVEL levels is divided by 2 in comparison with the best available ab initio PES. The The CCSD(T)-F12a/CVQZ-F12 PES gave an inversion height of 1788.58 cm -1 which is about 2 cm -1 larger than that already published. This is mainly explained by more accurate calculation of the relativistic and DBOC corrections in this work. As to the electronic correlations, the quadruple excitations gave the largest contribution on both the inversion barrier height and band origins. The RMS deviation with the experimental band origins of 1 cm -1 can be already achieved after adding the electronic correlations of the perturbative quadruple level of the excitations [CCSDT(Q)]. We hope this work will be of great help for future ro-vibrational calculations and spectra analysis of NH3. The vibrational energy levels with their assignments can be found in Supplementary Information as well as the program code to compute potential energies from the constructed PESs.

Supporting Information

Supporting File I includes the following: Table S1

Figure 1 .

 1 Figure 1. Height (in cm -1 ) of the inversion barrier for NH3 computed from the CCSD(T) and CCSD(T)-F12q (q=a, b) methods as a function of the cardinal number for the basis set (the effect of other contributions is given in Table4). The result for CCSD(T)-R12 was taken from ref 4. Note, here the "CBS" results come

a, b

 b Geometries and Einv were obtained: ain MOLPRO; bfrom the fitted PES. c Ref 4.

  grid and fit of the PES. The grid of the nuclear configurations was first generated by fixing e    and varying simultaneously the remaining five symmetry coordinates. The resulting grid was then prolonged to e       , where   took the following values: ±5, ±7, ±14, and ±22. Since the use of the sine function in Eq. (5) may be problematic at the two limit angles (  0° and  180°), the range of the 1D cut for the inversion angle was from -1 ) by step of o 1 . The final grid (called Grid_1) thus contains 40148

1 .

 1 Variational calculations. The Hamiltonian model was built using the Hougen-Bunker-Johns formalism[START_REF] Hougen | The vibration-rotation problem in triatomic molecules allowing for a large-amplitude bending vibration[END_REF] where the relation between the curvilinear bar{ρ} and effective coordinates ρ was explicitly established in ref 12. For this work, a hybrid model was implemented in our homemade computer code TENSOR which consists in Taylor-expanding the kinetic energy and potential parts in terms of normal mode coordinates associated with the ν1, ν3, and ν4 modes for each point of a numerical grid describing the nonrigid coordinate ρ.[START_REF] Viglaska | Derivation of ρ-dependent coordinate transformations for nonrigid molecules in the Hougen-Bunker-Johns formalism[END_REF] Here, a 11 th order Hamiltonian composed of 8648 D3h irreducible tensor operators was built on a grid of ρ ϵ [0.5, 2.7] rad by step of 0.02. In order to properly converge the vibrational energy levels (with errors < 0.01 cm -1 ) in the range of observed values, 6945, 5985 and 12915 basis functions have been used in the variational calculation for the symmetry blocks ( 1 A , 2 A ), ( 2 A , 1 A ) and ( E , E ), respectively. In this work, the atomic masses were employed in the kinetic energy operator while the nuclear masses were used in the DBOC calculation, as it is implemented in the CFOUR code (more details about the role of the masses with/without non-adiabatic contributions can be found in refs 34 and 35).

3. 2 .

 2 Comparison with experimental data. Investigation of the high-resolution spectra of NH3 has a long-standing history. Below, we give a non-exhaustive list of the main experimental studies. Since this molecule is of great interest in astrophysical observations, its inversion-rotation spectra within ground and excited vibrational states were the subject of many works (see e.g. refs 36-43, and references therein). Among the early results on the study of infrared spectra, let us cite the works by Urban et al.[START_REF] Urban | Coriolis and l-type interactions in the ν2, 2ν2, and ν4 states of 14 NH3[END_REF][START_REF] Urban | The ν2 band of 14 NH3: A calibration standard with better than 1×10 -4 cm -1 precision[END_REF][START_REF] Urban | High-J transitions in the ν2 bands of 14 NH3 and 15 NH3[END_REF] where the low-lying fundamental band v2 was studied using an effective vibration-inversion-rotation Hamiltonian. The hyperfine structure observed in the ground state and the v2 bands due to the nitrogen quadrupole splitting was then analyzed in refs 47-50. In addition, we can also mention some other studies where the transitions of the v2 band were measuredrefs 40, 51-57, and references therein; while refs 44, 58-61, and references therein were devoted to the overtone of the inversion mode as well as to the v4 band. The 3v2/v2+v4 system was studied in ref 62 while the interacting band system v1/v3/2v4 was recorded and analyzed in ref 63 together with the 4v2 band. A re-analysis of 2v2/v4 was made in ref 64 and by Pearson et al. 43 who also re-analyzed the 3v2/v2+v4 band system in ref 65. Some efforts were made in refs 42, 65-71, and references therein to assign hot band transitions. Previous works on v1/v3/2v4 can be found in the references cited in this paragraph.

  assigned new transitions of the v1+2v4 band with the help of an external cavity tunable diode laser spectrometer. The spectral region around v3+2v4 was analyzed in ref 81 by isotopic shift labeling while Sung et al.[START_REF] Sung | Extended line positions, intensities, empirical lower state energies and quantum assignments of NH3 from 6300 to 7000 cm -1[END_REF] conducted an analysis in the whole 6300-7000 cm -1 region. The line positions and intensities of all previously measured bands (v1+v3, 2v3, v1+2v4, v3+2v4, and 2v2+3v4) were reanalyzed from the room and low temperature Fourier spectra. According to ref 82, among the 4800 resolved features in the recorded spectra, only 1096 lines were assigned.

b

  The experimental values were taken from the MARVEL database; squared brackets denote the values with low spectroscopic accuracy (see details in ref 85).c The CCSD(T)/ACV6Z PES was obtained from the CCSD(T)/ACV5Z PES by the extrapolation of the CCSD(T)/[ACV6Z-ACV5Z] energy difference computed on a grid of the 2235 points (see details in Section 2.3).dThe CCSD(T)/CBS PESs were obtained by the extrapolation of the CCSD(T)/ACVXZ (X=Q, 5, and 6; X=T, Q, 5, and 6) energies to the complete basis set limit (see Section 2.1).The marked by bold-italic font band origins with the deviations bigger than 1.1 cm -1 for the CCSD(T)-F12a/CVQZ-F12 PES are related with only 1 or 2 experimental rovibrational levels (see details in ref 85).

  the ab initio PESs developed in this work allow decreasing the RMS error by a factor 2 compared to refs 4 and 5. The ab initio results from refs 2 and 3, which are based on the extrapolations of the CCSD(T)/ACVXZ (X = T, Q, 5) and CCSD(T)/AVXZ (X = T, Q, 5) energies respectively, are systematically bigger than the experimental values. According to our convergence test (see Figure 1), the NH3 energies computed with the ACV6Z orbital basis set size should be taken into account when extrapolating to CBS. It is confirmed by our CCSD(T)/CBS(TQ56Z) PES which is as accurate as the PES based on the explicitly correlated coupled cluster method. Although the multireference methods were used in ref 5, the ACV6Z basis set was also not considered. Finally, the band origins computed from the ab initio PES in ref 5 are in less good agreement with MARVEL than those given in ref 4 and they are systematically overestimated, like in refs 2 and 3.Concerning the importance of the multireference effects for NH3, we put FigureS1with the T1 and D1 diagnostic numbers in Supporting File I. The T1 and D1 diagnostic numbers computed and printed by the MOLPRO package are used to estimate the multireference behavior of the coupled cluster method. It is seen that for the nuclear geometries considered in this work for NH3, T1 and D1

2 and 28 .

 28 The decrease of 11 cm -1 obtained in ref 4 can be explained by the use of the Hartree-Fock method (see the end of Section 2.1 in ref 4).

F12a/

  CVQZ-F12+Relativistic+DBOC+CCSDT] and [CCSD(T)-F12a/CVQZ-F12+Relativistic+DBOC] levels of the theory We can see that the CCSDTQ and CCSDTQP contributions are very close in absolute values but with the opposite signs. The positive signs in the "CCSDTQ" column indicate that the CCSDT(Q) method overestimates the correlation energy from the determinants of the quadruple order of the excitation. For this reason, the RMS deviation with the MARVEL database is lowered by 0.05 cm-1

  CCSD(T)/CBS(Q56Z) PES provides the result that lies between the CCSD(T)/ACV6Z and CCSD(T)/CBS(TQ56Z) PESs in terms of RMS. Another PES based on the explicitly correlated CCSD(T)-F12a/CVQZ-F12 method gave quite similar to the CCSD(T)/CBS(TQ56Z) PES results but was about 300 less demanding in terms of CPU time. Both these PESs can be recommended for future theoretical calculations. Note also that two previously published PESs of NH3 constructed by the extrapolation of the CCSD(T) energies calculated with the aug-cc-pVXZ and aug-cc-pCVXZ (X=3, 4, and 5) basis sets gave systematically bigger values of the band centers.

  FigureS1(The T1 and D1 diagnostic numbers in dependence on the energy). The C++ code for computing potential energies from the constructed PESs can be found in Supporting File II.

Table 1 . Equilibrium and saddle point geometries as well as height of the inversion barrier (Einv, in cm -1 ) using different methods and basis sets

 1

	Method/Basis	re (Angs.)	ρe (deg.)	αe (deg.)	rplanar (Angs.)	Einv
	CCSD(T)/ACVTZ a	1.013429	67.76788	106.57474	0.996436	1879.07
	CCSD(T)/ACVQZ a	1.011504	67.89177	106.71031	0.995047	1806.97
	CCSD(T)/ACV5Z a	1.010867	67.94277	106.76596	0.994563	1783.45
	CCSD(T)/ACV6Z b	1.010755	67.96003	106.78478	0.994484	1778.11
	CCSD(T)/CBS(TQ5Z) b	1.010555	67.97033	106.79600	0.994310	1772.34
	CCSD(T)/CBS(Q56Z) b	1.010667	67.96433	106.78946	0.994413	1775.32
	CCSD(T)/CBS(TQ56Z) b	1.010654	67.96949	106.79509	0.994403	1774.50
	CCSD(T)-F12a/CVQZ-F12 b	1.010705	67.98181	106.80851	0.994464	1774.18
	CCSD(T)-F12b/CVQZ-F12 b	1.010680	67.97524	106.80135	0.994422	1778.04
	CCSD(T)-R12 c	1.01069	67.974	106.80	0.99445	1774.0
	For extrapolations, the total energies (reference + correlation) were used with the functional (1) (see details in ref 22).

Table 2 . Differences between experimental and theoretical band origins of NH3 computed from the ab initio PESs including Relativistic, DBOC, CCSDT, CCSDTQ, and CCSDTQP corrections a

 2

	Band					Exp. b -Calc.
	V1 V2 V3 V4 Γvib	Exp. b	CCSD(T)/ ACV6Z c	CCSD(T)-F12b/ CVQZ-F12	CCSD(T)/ CBS(Q56Z) d	CCSD(T)/ CBS(TQ56Z) d	CCSD(T)-F12a/ CVQZ-F12
	0 0 0 0	2 A 0.7934030(10)	0.016	0.009	0.010	0.007	0.001
	0 1 0 0	1 A	[932.4]	-1.184	-0.663	-0.693	-0.509	0.202
	0 1 0 0	2 A 968.121930(28) -0.567	-0.182	-0.328	-0.204	0.397
	0 2 0 0	1 A	[1597.5]	-2.482	-2.276	-1.356	-1.047	-0.424
	0 0 0 1 E 1626.27472(29) 0.181	-0.518	-0.001	-0.056	0.092
	0 0 0 1 E 1627.37238(20) 0.203	-0.506	0.012	-0.044	0.094
	0 2 0 0	2 A 1882.177445(30) -0.561	-0.524	-0.301	-0.162	0.399
	0 3 0 0	1 A	[2384.2]	-0.421	-0.767	-0.223	-0.063	0.326
	0 1 0 1 E 2540.5243(27) -1.150	-1.221	-0.812	-0.715	0.283
	0 1 0 1 E 2586.12855(50) -0.426	-0.677	-0.388	-0.350	0.480
	0 3 0 0	2 A 2895.52192(42) 0.387	-0.088	0.172	0.309	0.754
	0 2 0 E	[3189.4]	-2.346	-2.708	-1.561	-1.334	-0.384
	0 0 0 2	1 A	[3215.9]	0.281	-0.939	-0.098	-0.171	0.109
	0 0 0 2	2 A 3217.57919(50) 0.356	-0.874	-0.033	-0.107	0.161
	0 0 0 2 E 3240.16297(53) 0.370	-0.978	0.015	-0.089	0.189
	0 0 0 2 E 3241.59828(65) 0.402	-0.960	0.034	-0.072	0.193
	1 0 0 0	1 A	[3336.1]	0.510	-0.471	0.002	0.183	-0.191
	1 0 0 0	2 A 3337.09694(53) 0.528	-0.451	0.021	0.198	-0.177
	0 0 1 0 E 3443.62788(54) 0.577	-0.485	0.023	0.184	-0.471
	0 0 1 0 E 3443.98765(53) 0.585	-0.482	0.028	0.187	-0.471
	0 4 0 0	1 A	[3462.4]	0.909	0.266	0.184	0.560	1.064
	0 2 0 1 E	[3502.6]	-0.336	-0.902	-0.360	-0.279	0.529
	0 4 0 0	2 A	[4061.6]	1.517	0.786	0.319	0.881	1.557
	0 1 0 2	1 A	[4115.9]	-0.941	-1.510	-0.825	-0.742	0.559
	0 1 0 2 E	[4135.8]	-1.243	-1.874	-1.100	-1.035	0.248
	0 1 0 2	2 A	[4173.1]	-0.250	-1.017	-0.449	-0.444	0.613
	0 1 0 2 E	[4192.9]	-0.423	-1.279	-0.599	-0.618	0.422
	1 1 0 0	1 A	[4294.5]	-0.336	-0.681	-0.545	-0.056	0.301
	1 1 0 0	2 A 4320.030806(90) 0.088	-0.364	-0.271	0.149	0.439
	0 1 1 0 E 4416.915053(90) -0.302	-0.768	-0.571	-0.194	-0.015
	0 1 1 0 E 4435.446476(90) 0.040	-0.521	-0.363	-0.028	0.072
	0 2 0 2 E	[4773.8]	-2.277	-3.208	-2.040	-1.679	-0.442
	0 0 0 3 E	[4799.2]	0.739	-1.023	0.182	0.078	0.530
	0 0 0 3 E	[4801.4]	0.921	-0.862	0.344	0.234	0.666
	0 0 0 3	1 A	[4841.6]	0.462	-1.492	-0.076	-0.214	0.193
	1 0 0 1 E 4955.75609(50) 0.705	-0.996	-0.010	0.138	-0.095
	1 0 0 1 E	[4956.9]	0.709	-0.994	-0.002	0.137	-0.101
	1 2 0 0	1 A	[5000.2]	-2.016	-2.474	-1.580	-0.826	-0.366
	0 0 1 1	2 A	[5052.0]	0.710	-1.071	-0.049	0.093	-0.405
	0 0 1 1	1 A	[5052.7]	0.808	-0.944	0.045	0.195	-0.320
	0 0 1 1 E	[5052.6]	0.706	-1.082	-0.056	0.083	-0.447
	0 0 1 1 E	[5053.2]	0.716	-1.044	-0.052	0.095	-0.453

Table 3 . Comparison of calculated band origins of NH3 using different ab initio PESs with experiment.

 3 The experimental values were taken from the MARVEL database; squared brackets denote the values with low spectroscopic accuracy (see details in ref 85). I -Rajamӓki et al.[START_REF] Rajamäki | High excitations in coupled-cluster series: vibrational energy levels of ammonia[END_REF] ; II -Yurchenko et al.[START_REF] Yurchenko | Potential-energy surface for the electronic ground state of NH3 up to 20000 cm -1 above equilibrium[END_REF] ; III -Huang et al.[START_REF] Huang | An accurate global potential energy surface, dipole moment surface, and rovibrational frequencies for NH3[END_REF] ; IV -Marquardt et al.[START_REF] Marquardt | Global Analytical Potential Energy Surface for the Electronic Ground State of NH3 from High Level ab Initio Calculations[END_REF] ; V -Polyansky et al.[START_REF] Polyansky | Calculation of rotation-vibration energy levels of the ammonia molecule based on an ab initio potential energy surface[END_REF]

	Band						Exp. -Calc.
	V1 V2 V3 V4 Γvib	Exp. a	I	II	III	IV	V	CCSD(T)/	CCSD(T)-F12a/
									CBS(TQ56Z)	CVQZ-F12
	0 0 0 0	2 A 0.7934030(10) 0.00 0.00 -0.02 0.04 0.01	0.01	0.00
	0 1 0 0	1 A	[932.4]	0.39 -1.85 -0.48 0.13 -0.64	-0.51	0.20
	0 1 0 0	2 A 968.121930(28) 0.49 -1.61 -1.11 1.90 -0.38	-0.20	0.40
	0 2 0 0	1 A	[1597.5]	0.22 -3.81 0.48 -3.96 -1.12	-1.05	-0.42
	0 0 0 1 E 1626.27472(29) 0.65 -2.57 -2.05 -2.18 0.09	-0.06	0.09
	0 0 0 1 E 1627.37238(20) 0.67 -2.58 -2.07 -1.96 0.11	-0.04	0.09
	0 2 0 0	2 A 1882.177445(30) 0.21 -3.33 -2.06 3.19 -0.44	-0.16	0.40
	0 3 0 0	1 A	[2384.2]	0.01 -4.35 -2.58 3.17 -0.41	-0.06	0.33
	0 1 0 1 E 2540.5243(27) 1.31 -4.25 -1.89 -19.94 -1.13	-0.72	0.28
	0 1 0 1 E 2586.12855(50) 1.45 -4.08 -2.74 -13.35 -0.62	-0.35	0.48
	0 3 0 0	2 A 2895.52192(42) 0.00 -4.97 -4.07 6.71 -0.21	0.31	0.75
	0 0 0 2	1 A	[3215.9]	0.92 -5.25 -3.24 -4.97 -0.14	-0.17	0.11
	0 0 0 2	2 A 3217.57919(50) -0.75 -5.19 -3.23 -4.51 -0.07	-0.11	0.16
	0 0 0 2 E 3240.16297(53) 0.87 -5.25 -3.68 -6.71		-0.09	0.19
	0 0 0 2 E 3241.59828(65) 1.40 -5.28 -3.70 -7.13 -0.17	-0.07	0.19
	1 0 0 0	1 A	[3336.1]	-0.02 -6.14 -0.98 1.37 -1.52	0.18	-0.19
	1 0 0 0	2 A 3337.09694(53) -0.09 -6.12 -0.98 1.49 -1.50	0.20	-0.18
	0 0 1 0 E 3443.62788(54) -1.65 -5.72 -0.06 -2.22 -0.50	0.18	-0.47
	0 0 1 0 E 3443.98765(53) -1.55 -5.72 -0.06 -2.00 -0.50	0.19	-0.47
	0 4 0 0	1 A	[3462.4]	-0.18 -5.93 -5.21 9.32 -0.12	0.56	1.06
	0 1 0 2	1 A	[4115.9]	1.69 -7.11 -2.76 -18.51 -1.62	-0.74	0.56
	0 1 0 2 E	[4135.8]	1.53 -7.20 -3.24 -31.76 -2.05	-1.03	0.25
	0 1 0 2	2 A	[4173.1]	1.05 -6.96 -3.96 -11.30 -0.99	-0.44	0.61
	0 1 0 2 E	[4192.9]	1.52 -7.14 -4.30 -21.41 -1.28	-0.62	0.42
	1 1 0 0	1 A	[4294.5]	0.24 -7.96 -1.23 20.43 -2.16	-0.06	0.30
	1 1 0 0	2 A 4320.030806(90) 0.29 -7.91 -1.81 17.16 -1.99	0.15	0.44
	0 1 1 0 E 4416.915053(90) -1.22 -7.49 -0.80 7.94 -0.92	-0.19	-0.01
	0 1 1 0 E 4435.446476(90) -1.20 -7.48 -1.20 7.54 -0.82	-0.03	0.07
	1 0 0 1 E 4955.75609(50) 0.41 -9.27 -2.90 6.17 -1.45	0.14	-0.10
	1 0 0 1 E	[4956.9]	0.82 -9.34 1.14 6.38 -1.44	0.14	-0.10
	0 0 1 2 E 6677.43169(50) -0.54 -1.03 -2.22 -4.71 -1.11	-0.06	-0.86
	Total RMS		0.938 5.744 2.538 11.042 1.048	0.401	0.412

a

Table 4 . Effect of the relativistic and DBOC corrections as well as high-order electronic correlations on the equilibrium and saddle point geometries and on the height of the inversion barrier (Einv in cm -1 ) for NH3

 4

	Ab initio level	re (Angs.)	ρe (deg.)	αe (deg.)	rplanar (Angs.)	Einv
	CCSD(T)-F12a 1.010705	67.98181	106.80851	0.994464	1774.18, 1774.0 a
	+Relativistic	1.010722	67.92817	106.75003	0.994346	1793.57 (+19.39), 1792.1 a
	+DBOC	1.010750	67.95004	106.77389	0.994416	1783.56 (-10.01), 1781.1 a
	+CCSDT	1.010767	67.94875	106.77248	0.994426	1784.96 (+1.40), 1782.7 a
	+CCSDT(Q)	1.010981	67.92761	106.74942	0.994581	1788.02 (+3.06)
	+CCSDTQ	1.010969	67.92791	106.74975	0.994571	1788.18 (+0.16), 1785.7 a
	+CCSDTQP	1.010990	67.92572	106.74736	0.994584	1788.58 (+0.40), 1786.7 a
	a Rajamӓki et al. 4

Table 5 . Contributions of the corrections beyond the CCSD(T)-F12a/CVQZ-F12 level to the band origins of NH3

 5

	Band V1 V2 V3 V4 Γvib	Relativistic	DBOC	CCSDT	CCSDT(Q)	CCSDTQ	CCSDTQP
	0 0 0 0	2 A	-0.05	0.02	0.00	-0.02	0.00	0.00
	0 1 0 0	1 A	4.68	-2.43	0.51	0.35	0.06	0.09
	0 1 0 0	2 A	2.90	-1.63	0.41	-0.25	0.07	0.02
	0 2 0 0	1 A	9.31	-4.80	0.91	1.02	0.08	0.21
	0 0 0 1 E	0.55	-0.27	0.25	-1.42	0.11	-0.10
	0 0 0 1 E	0.48	-0.24	0.25	-1.44	0.12	-0.10
	0 2 0 0	2 A	3.72	-2.32	0.67	-1.09	0.15	-0.04
	0 3 0 0	1 A	3.81	-2.53	0.79	-1.67	0.20	-0.09
	0 1 0 1 E	5.62	-2.85	0.79	-0.91	0.16	0.01
	0 1 0 1 E	3.42	-1.87	0.66	-1.67	0.18	-0.08
	0 3 0 0	2 A	1.79	-1.75	0.81	-2.94	0.27	-0.22
	0 0 0 2	1 A	0.79	-0.47	0.39	-2.86	0.23	-0.22

TOC Graphic

ACKNOWLEDGEMENTS

The supports from the French-Russian laboratory IRP "SAMIA2" and from the ROMEO computer center of Reims-Champagne-Ardenne are acknowledged. O.E. and A.V.N. acknowledge the supports from the Ministry of Science and Higher Education of the Russian Federation (V.E. Zuev Institute of