Romain Garnier

Marc Odunlami
email: marc.odunlami@univ-pau.fr

Vincent Le Bris

Didier Bégué

Isabelle Baraille

Olivier Coulaud

Adaptive vibrational configuration interaction (A-VCI): a posteriori error estimation to efficiently compute anharmonic IR spectra

Keywords: numbers: 31, 15, xt, 33, 20, Tp

A new variational algorithm called adaptive vibrational configuration interaction (A-VCI) intended for the resolution of the vibrational Schrödinger equation was developed. The main advantage of this approach is to efficiently reduce the dimension of the active space generated into the configuration interaction (CI) process. Here, we assume that the Hamiltonian writes as a sum of products of operators. This adaptive algorithm was developed with the use of three correlated conditions i.e. a suitable starting space ; a criterion for convergence, and a procedure to expand the approximate space. The velocity of the algorithm was increased with the use of a posteriori error estimator (residue) to select the most relevant direction to increase the space. Two examples have been selected for benchmark. In the case of H 2 CO, we mainly study the performance of A-VCI algorithm: comparison with the variationperturbation method, choice of the initial space, residual contributions. For CH 3 CN, we compare the A-VCI results with a computed reference spectrum using the same potential energy surface and for an active space reduced by about 90 %.

I. INTRODUCTION

Vibrational spectroscopy, a widely used experimental technique for the identification and characterization of molecules in complex chemical environments, has been applied successfully in fields as diverse as biochemistry, agribusiness, interstellar chemistry, as well as the chemistry of materials. With the joint progresses of modeling techniques and computer technology, quantum mechanical calculations of vibrational spectra have become an essential and powerful contributor for the identification and discrimination of diverse molecules in complex chemical systems 1,2 . The support of theoretical methods is used to fully assign the bands of the vibrational spectra and identify the impact of anharmonicity, presence of non-fundamental transitions, solvent or aggregation effects. Nevertheless, accurate calculation of vibrational properties for large molecular systems remains a very challenging problem for computational chemistry. The common starting point for these calculations is the time-independent Schrödinger equation and the standard separation of electronic and nuclear motions within the framework of the Born-Oppenheimer approximation. As a consequence, the wave function of the molecule is written as a product of a nuclear and an electronic function, the nuclear motion being determined by the potential electronic energy, usually called potential energy surface (PES). The nuclear part describes both the rotation and the vibration of the molecular system. The most general and efficient way of calculating the vibrational properties of a molecule is to compute the eigenvalues of the vibrational Hamiltonian in a discrete space spanned by basis functions built as products of functions of a single variable. The exact solutions of the harmonic equation are the quantum harmonic oscillators written as functions of reduced normal coordinates. They are often considered as a natural basis set in current methods admitting that the discretization of the anharmonic Hamiltonian should give a good approximation of the vibrational spectra. The computational cost of the calculation and the accuracy of the computed vibrational frequencies depend essentially on two limiting steps:

(i) The PES is generally expressed as a polynomial determined by adjustment of a set of molecular structures and depends on the accuracy of the calculations of the electronic wave functions. While it is now common to determine the PES accurately for small molecules (3-5 atoms), this becomes almost unfeasible when solving anharmonic approximations for systems of larger size. The truncation of the analytic shape of the PES (degree of the potential) is unavoidable and the series is often truncated at the fourth order.

(ii) For the resolution of the vibrational Schrödinger equation, the anharmonic vibrational wave function is expanded as a linear combination of products of 1-D functions. This point raises two problems: the choice of the 1-D functions and the dimension of the basis directly connected to the treatment of the discretized Hamiltonian.

Unlike the state-specific methods like vibrational coupled cluster (VCC) 3,4 , the variational methods such as vibrational configuration interaction (VCI) [5][6][7][8] , parallel vibrational multiple window configuration interaction (P-VMWCI) 9 , vibrational self-consistent field (VSCF 10,11 , VMCSCF 12 ,• • • ) are based on the diagonalization of the discretized Hamiltonian matrix in the most extended basis set of products. These approaches commonly use Lanczos 13,14 method or Davidson 14,15 algorithms to diagonalize the Hamiltonian matrix in the full approximate space. As the number of degrees of freedom exponentially increases with the number of atoms in the molecule, these methods quickly encounter the curse of dimensionality for large molecules (more than 10 atoms). In addition, the perturbation approach [16][17][18] (second order perturbation theory VMP2, for example) is easier to implement than variational methods but fails to reproduce the strong anharmonic couplings such as Fermi and Coriolis resonances. Another way to compute eigenvalues of the vibrational Hamiltonian is to combine the VCI method with basis selection techniques as in the so-called variationperturbation methods (VP) [19][20][21][22] . Then, the challenging key point is to find a suitable basis set to generate an active space of minimal dimension containing all the vibrational states of interest and the most pertinent interacting ones as shown in [23][24][25][26] .

The analysis of the vibrational model proposed in this paper focuses on the choice of the discretized harmonic space for the computation of the eigenvalues of the Hamiltonian matrix. In this paper we consider that the vibrational Hamiltonian is a sum of products (SOP) of unidimensional operators.

We develop an adaptive approach of the variational formulation of the problem to efficiently reduce the dimension of the active space. In the proposed algorithm, we completely use the properties of the orthonormal basis of the harmonic operator to construct a sequence of nested spaces in the full approximation space. In these spaces we search the eigenvalues of the Hamiltonian and we present a criterion that ensures the convergence of our method.

By the way, we never construct the complete matrix. This last point removes the memory bottleneck for medium-sized molecules. Section II is devoted to the theoretical background used in this paper and to the related works. Developing an adaptive algorithm to find both a minimal approximation space and some eigenpairs of the Hamiltonian requires at least three conditions: 1) a suitable starting space ; 2) a criterion for convergence and 3) a procedure to expand the approximate space. In section III, we present a posteriori error estimator to select the most relevant direction to enlarge the search space. Based on this estimator we introduce the A-VCI algorithm (namely adaptive vibrational configuration interaction). In section IV, we describe the results obtained for the H 2 CO and CH 3 CN infrared spectra with the A-VCI method. In the case of the H 2 CO molecule (6 vibrational coordinates) we present the influence of the A-VCI parameters with respect to the convergence of the method, and we compare our algorithm with the standard variation-perturbation method. Finally, we show for the CH 3 CN molecule (12 vibrational coordinates) that the A-VCI results obtained in a reduced basis set agree with the reference spectrum computed by Avila et. al 27 using the same PES.

II. THEORETICAL BACKGROUND, NOTATIONS, DEFINITIONS

Let consider the dimensionless normal coordinates q = (q 1 , q 2 , q 3 , . . . q N ), where N is the number of vibrational degrees of freedom of the system. In the Born-Oppenheimer approximation, the vibrational Hamiltonian H(q) is the sum of the harmonic operator H 0 and the anharmonic part V and it writes as

H(q) = H 0 + V(q) = N i=1 ω i 2 - ∂ 2 ∂q 2 i + q 2 i + S s 1 =3 K s N i=1 q s i i , (1)
where ω i is the harmonic frequency (in cm -1 ) associated with the q i coordinate, s 1 is the sum of elements of multi-index s = (s 1 , s 2 , . . . , s N ) and S is the maximal degree of the PES.

The eigenfunctions of H 0 noted φ 0 n (q) are products of the one-dimensional normalized Hermite functions ψ n i (q i ):

φ 0 n (q) = N i=1 ψ n i (q i ), for all n = (n 1 , . . . , n N ), (2)
and the corresponding eigenvalues (in cm -1 ) are

E 0 n = N i=1 (n i + 1 2 )ω i ,
n i being the quantum number associated with the ω i frequency.

Let define d i the maximal degree of the Hermite function associated with the i-coordinate in ( 2) and Π d the product space of size M of all quantum numbers defined by

Π d = N i=1 [0, d i ].
We introduce θ a mapping between Π d and I = {0, . . . , M -1}, and for i ∈ I there exists a

unique n ∈ Π d such that i = θ(n).
We seek the eigenfunctions φ of H as linear combinations of the eigenfunctions of H 0 as

φ(q) = n∈Π d x n φ 0 n (q).
By using the variational method on the eigenvalue problem for the operator defined in (1), we obtain the matrix form of the problem

H f ull X = EX (3)
where H f ull is a M × M symmetric matrix and (E, X) is an eigenpair of H f ull . An element of H f ull is calculated as

H θ(m),θ(n) = φ 0 m |Hφ 0 n = E 0 m δ m,n + φ 0 m |Vφ 0 n . (4)
Thanks to the properties of Hermite functions and the PES (SOP form) the matrix coefficients ( 4) are analytic and easily calculated 28 .

A. Related works

Solving Eq. ( 3) by iterative methods needs to compute fast matrix-vector products. Even if the Hamiltonian matrix is very sparse, the size of the vectors, M , grows as d N if all d i = d in each direction. As shown in Ref. 29 , if N = 12 and d = 10, then for a single vector using 64 bits for each entry one needs 8 TB of memory. To overcome the curse of dimensionality, several approaches have been developed, which can be classified in two categories following how the basis is constructed: (i) fixed basis set (full VCI) ; (ii) adaptive basis set construction (variation-perturbation methods).

In the first approach, the goal is to choose the smallest subset of Π d giving a good approximation of the eigenpairs. The first idea is to consider a complete active space B d including all the harmonic vibrational states φ 0 n such that

B d = {φ 0 n / |n| = N i=1 n i ≤ d}. (5)
The total number of elements of this binomial basis is

(d + N )! N !d! .
In Ref. When an adaptive basis set construction is used, the eigenvalue computation of matrix H is usually performed by a variational method combined with an iterative process. The matrix H is diagonalized in embedded basis sets B (j) ⊂ B (j+1) ⊂ Π d , where B (j) is the basis set at iteration j. In the classical variation-perturbation algorithm 19,20,22 configurations φ 0 m of B (j+1) are selected by imposing the following second order perturbation criterion on the energy:

φ 0 m |Vφ 0 n 2 E 0 n -E 0 m > ε VP , E 0 n -E 0 m = 0,
where φ 0 n ∈ B (j) and ε VP is a given energetic threshold depending on the accuracy we want to reach. Usually the initial subspace B (0) includes the ground state and its monoexcitations to reproduce the fundamental transitions, which mainly govern the shape of the vibrational spectrum.

The convergence of the eigenvalue E (j) at iteration j is classically checked by evaluating |E (j) -E (j-1) |. This constitutes a misleading criterion 33,34 , some plateaus and oscillations can be observed between two iterations. To converge to highly excited energy levels the authors proposed an approach based on a modified Davidson algorithm. The convergence of the iterative process is evaluated from the norm of a residual vector. They empirically concluded that residue values of 200, 100 and 20 cm -1 corresponded to absolute errors of 10, 0.1 and 0.01 cm -1 .

B. Mathematical background

The Bauer-Fike theorem 35,36 is a classical tool in spectral perturbation theory to localize eigenvalues. This theorem is useful to characterize the error on a perturbation of the eigenvalue of a diagonalizable matrix. There are several forms of the theorem, and we consider here the two following one.

Theorem II.1 (Bauer-Fike) Let ( Ẽ, X) be an approximated eigenpair of the real symmetric matrix H with X = 1, and let R = H X -Ẽ X the residual vector. H is diagonalizable satisfying H = P DP T , where D is a diagonal matrix and P an unitary matrix then

∃E ∈ σ(H) such that | Ẽ -E| ≤ R
where σ(H) is the spectrum of H and . the Euclidean norm.

This formulation justifies the use of the residue to check the convergence of the eigenvalues with an iterative solver like Davidson schemes 33,34 .

In the context of the perturbation theory, we prefer the following alternative formulation of the theorem. If ∆H is a symmetric perturbation of a symmetric matrix H then for any eigenpairs ( Ẽ, X) of H + ∆H with X = 1, there is an eigenvalue E of H such that

| Ẽ -E| ≤ ∆H X ≤ ∆H . ( 6
)
If the Frobenius norm of ∆H is small, then the theorem says that each eigenvalue of H +∆H is close to an eigenvalue of H. However, the theorem says nothing about the correspondence between an eigenvalue of the perturbed matrix and the eigenvalue of the unperturbed matrix.

III. A FAST RESIDUE-BASED ADAPTIVE VCI ALGORITHM A. A posteriori error estimator

Classically an eigenvalue algorithm checks the convergence of a pair (E, X) using the scaled norm of the residual vector R = HX -EX. In such algorithm the size of the approximate space (i.e., full space) is fixed and could be very large. Now, consider an adaptive process in which both the size of each vector and the approximate space dimension increase during the iterations. Then which residue should be considered? Do we need to extend it in the full space? How do we efficiently compute it? These are the main issues of the adaptive procedure we will address here.

Let B a subset of the full approximate space. A function φ of B writes as a linear combination of 1-D Hermite function products. As the operator H is a polynomial operator the function Hφ can be expressed as the sum of two linear combinations of functions:

Hφ = φ 0 n ∈B α n φ 0 n + φ 0 n / ∈B β n φ 0 n .
Let H(B) the image space of B by H. This space is spanned by all harmonic functions φ 0 m such that φ 0 m |Hφ = 0 for all φ in B. As operator H 0 is diagonal in the 1-D Hermite function products, the image space is decomposed into two orthogonal spaces as

H(B) = B ⊕ B R . (7)
We denote by m the number of elements in B and by m R the number of elements in B R .

As φ 0 n is the product of 1-D Hermite functions and the degree of the PES, S, is small it is easy to see that H(B) is a small subset of the full approximate space.

Let introduce the Rayleigh matrix H of operator H in the space H(B). Thanks to the above space decomposition this matrix is decomposed as follows

H =         H H T R H R H C         (8)
where H = B T HB is a m × m Rayleigh matrix that approximates the Hamiltonian operator

H in the orthonormal basis B, H R = B T R HB a m R ×m matrix and H C = B T R HB R a m R ×m R matrix.
Due to the space decomposition (7) we only need to extend a vector X of space B in space H(B). This is easily done by adding m R zeros to X. The extended vector X writes (X, 0 m R ) T .

Consider (E, X) an eigenpair of the Rayleigh matrix H, the issue is to know if its corresponding extended eigenpair is a good approximation of an eigenpair of H. More precisely, we would like to know how far is (E, X) from an eigenpair of matrix H. The residual vector R of the extended eigenpair in H(B) writes

R = H X -E X = (0 m , H R X) T . ( 9
)
and then its norm is just the norm of the vector H R X. By using the alternative formulation of the Bauer-Fike theorem (Eq. 6) we obtain the estimation between the computed eigenvalue

E and Ẽ the target eigenvalue of H |E -Ẽ| ≤ H R X . (10)
Inequality ( 10) has to be interpreted as a posteriori error estimator to measure if the computed eigenpair is a good approximation of the target eigenvalue in a larger space. If H R X is small enough then (E, X) is a good approximation of an eigenpair of H. In an iterative process based on nested basis B (0) ⊂ B (1) ⊂ . . . ⊂ B (j) , the evaluation of H R X is a good convergence criterion to find the first F eigenpairs of H f ull . Furthermore, we will use it to select the appropriate functions in B R to enlarge the search space and decrease the residue norm.

To use in practice the a posteriori error estimation (10), we have to compute H R X at a reasonable cost. Different strategies are considered depending on the dimension, m R , of the space B R . First, due to the properties of the basis functions, the H R is a tall and skinny sparse matrix (i.e. m m R ). For a given basis element, the number of terms on its associated line in the full matrix depends on the number of terms of the PES and its global degree. In the block decomposition (8) H C is a huge m R × m R sparse matrix but in our approach we never store it. Therefore there is no memory bottleneck and we can easily store all the data in memory. We construct the matrix H R and store it in the CSR-format to efficiently perform matrix-vector products. On the other hand, for large molecules if we cannot store the full H R matrix on one machine, we can either parallelize on several nodes to distribute the matrix or perform the matrix-vector product on the fly, i.e., without constructing the matrix. This matrix-free approach increases the CPU-time of a matrixvector product but classical optimization (blocking, vectorization, parallelization, ...) can be performed to reduce the overhead 37 .

B. The adaptive vibrational configuration interaction

At each iteration j of the algorithm the elements of the block decomposition ( 8) are indexed by the superscript (j) as follows

H(j) =         H (j) H (j) R T H (j) R H (j) C         , (11)
with the corresponding decomposition space

H(B (j) ) = B (j) ⊕B (j) R . We denote by m (j) (resp. m (j) R ) the number of elements in B (j) (resp. B (j) R ). Then, H (j) = B (j) T HB (j) is a m (j) × m (j) matrix and H (j) R = B (j) R T HB (j) a m (j) R × m (j) matrix.
We introduce an adaptive process where the selection of the basis functions is controlled by the residual contributions defined in (9). The algorithm builds a suitable basis set, B (j) , to compute the smallest F eigenvalues and eigenvectors of H (j) at the desired accuracy ε. We denote this adaptive scheme with enrichment by residual contributions A-VCI, for adaptive vibrational configuration interaction. for j ≥ 0 do 1 Build the Rayleigh matrix H (j) = B (j) T HB (j) 2 Compute the first F eigenpairs of H (j) denoted (E

(j) , X (j) ) F -1 =0 . 3 Construct the first F residual vectors R (j) = H (j)
R X // Check the convergence of the algorithm.

if max ∈[0,F -1] R (j) /E (j) ≤ ε then //
The method has converged exit else 4 Compute the new basis elements by selecting an orthonormal subset

A (j) of B (j)
R such that some components of R (j) will vanish at next step.

B (j+1) = B (j) ⊕ A (j)
The procedure of A-VCI is written in Algorithm 1. Lines ( 1) to (3) correspond to the classical Rayleigh-Ritz procedure 14 . First, we construct the Rayleigh matrix of the operator j) in line (1). At line (2) the eigenvalue problem can be solved by your favorite solver, here we use the Krylov-Schur procedure 38,39 to compute the F smallest eigenpairs. Then, we extend these eigenvectors in B (j) R to compute the residues (line ( 3)). Unlike Davidson or Lanczos procedures, the residues are not computed in the full space but only in the space H(B (j) ) . The size of the residual vectors is only m

H in B (
(j)
R which is really small compared to the size of Π d . The convergence criterion is based on the relative residue which is standard for iterative procedures. In line (4) the new selected directions are incorporated into the space B (j) . Here we select some basis elements (orthonormal vectors) of B (j) R that avoid a Gram-Schmidt process for orthonormalising the selected directions. This is another important difference with the Davidson procedure. We point out that only the last rows and columns of H (j) and H (j) R have to be computed at iteration j. The efficiency of our adaptive algorithm depends on the way we select the new basis functions (line 4) to add to B (j) in order to decrease the residue. If we add too many functions in B (j) at once then the memory and CPU requirements to perform the next eigensolver iteration will be prohibitive. On the other hand, if we add too few functions then we increase the number of iterations in the adaptive solver and the final cost will also be large. The way to construct the space is a tradeoff between the cost of the eigenvalue solvers at each iteration and the number of iterations of the A-VCI method to reach the convergence.

Component-wise procedure to increase the approximation space

The classical strategy is to consider for the next approximate space the image by H of the current one. This leads to the following choice

A (j) = B (j) R = H(B (j) ) \ B (j) . ( 12
)
With this approach there is no selection among the interacting configurations of the previous iteration to build the current basis. All the functions of B (j)

R are selected even if they do not contribute to the residues. Then B (j+1) = H(B (j) ) is used to construct the Hamiltonian at the next iteration. Here, the basis grows quickly and leads to prohibitive memory and CPU costs both in matrix construction and eigenvalue solver. In the other hand, we can only add one orthogonal direction to the approximate space like in the Lanczos procedure. In this case, the direction that maximizes the error on the residues is chosen.

By this way at the convergence we should obtain the minimal approximate space we can find. However, the number of iterations will be too large and it corresponds to the size of the final approximate space minus the size of B (0) . Our goal is both to reduce the size of the basis without altering the accuracy of the results and to reach the convergence in a small number of iterations. Eq. (9) shows that the residue of an extended eigenpair is orthogonal to the space B (j) . The main idea to build A (j) is to select the most relevant directions in order to improve the approximation of the eigenvector at next step. Many procedures can be imagined to enlarge the approximate space depending on how we choose the directions and how many functions are kept at each step.

Here we propose a component-wise approach based on the residual vector as described in Algorithm 2 to select the functions we add to the space B (j) . (j) , R (j) ) at iteration j Output: The set A (j) of functions to incorporate in B (j) begin

1 forall the such that R (j) /E (j) > ε do Determine the nnz non-zero entries in R (j) K = 0 2 foreach i in [m, m + m R -1] such that |(R (j) ) i | > ε √ nnz E (j) do K = K + 1 ; C K = |(R (j) ) i |/E (j) 3 Compute the average η = 1 K K k=1 C k . 4 For all k in [1, K] such that C k > η find its associated multi-index k and add φ 0 k in A (j)
First of all, we select relative residues greater than the threshold ε (line 1) i.e., the eigenpairs that have not yet converged. When we add a direction (i.e., a component) of the residue in the search space its contribution will be zero at the next iteration. So, the main idea is to select all the components of the residue such that the Euclidean norm of the remaining components will be less than the prescribed threshold. For that purpose, we introduce in line 2 the new threshold ε/ √ nnz with nnz the non-zero element number of the current residue to take the highest residue components in magnitude. In order to limit the number of added functions, we build the average of all previous contributions (line 3), and we only select the components corresponding to the values over this average (line 4). At the end, the space A (j) contains all directions that have the most important contributions on the residues.

Construction of the initial subspace

The last point is to explain how the initial space is chosen. Generally, we are interested by all the eigenvalues lower than E max . As we assume that the perturbation between H 0 and H is small, the eigenpairs of H 0 , are relatively close to the eigenpairs of H. Then it is natural to consider the space spanned by the harmonic functions with their associated energy lower than E max and writes

B (0) = span{φ 0 m /E 0 m ≤ E max }. Each frequency ν 0 m corresponding to E 0 m verifies E 0 m = E 0 0 + ν 0 m
, where E 0 0 represents the ground state energy. Introducing the frequency upper limit ν max such that E max = E 0 0 +ν max , we provide the following alternate definition of the initial subspace:

B (0) = span{φ 0 m /ν 0 m ≤ ν max }. ( 13
)

IV. RESULTS AND DISCUSSION

A. Formaldehyde molecule, H 2 CO

The proposed method is tested on the benchmark molecule H 2 CO to which numerous studies have been realized 34,40,41 . In this first example, we focus on the numerical efficiency of the A-VCI algorithm. The results will be compared to a reference calculation in a large basis set (296 010 functions), using the same PES developed by Le Bris et al. 42 . This potential energy surface is a polynomial quartic function as in Seidler et al. 41 , and has been determined at the CCSD(T)/aug-cc-pVTZ level. The electronic energies have been obtained with the Gaussian software 43 . All the force constants are provided in the supplemental material 44 .

The corresponding fundamental harmonic frequencies are:

ω 1 = 1 181.2, ω 2 = 1 261.1, ω 3 = 1 529.0, ω 4 = 1 765.4
, ω 5 = 2 932.5, ω 6 = 3 000.5 (cm -1 ).

VCI reference calculation

We perform a vibrational configuration interaction (VCI) calculation on H 2 CO in the harmonic binomial basis of maximal degree d = 21, defined in (5) and noted B 21 . The VCI space B 21 contains 296 010 configurations, and we consider it as a reference basis. Anharmonic frequencies, eigenvector contributions and intensities are presented in Table I along with the corresponding harmonic frequencies in the limit of ν max = 4 000 cm -1 . In this example, it is quite straightforward to associate an anharmonic frequency with its harmonic equivalent despite the discrepancies between the two approximations. The anharmonic frequencies are in agreement with those previously reported in 41 . All the fundamental transitions have an IR activity. The inclusion of the effects from electric anharmonicity leeds to active overtones (2ω 1 , 2ω 3 and 2ω 4 ) and combination bands (ω 2 + ω 4 ). We observe a strong mixing of the states ω 6 and ω 2 + ω 3 , with important IR activities for the two transitions.

Choice of the initial space

One key feature of the A-VCI approach is the choice of the initial subspace B (0) of the iterative process. In the classical variation-perturbation (VP) technique, the ground state and its monoexitations generate the initial space. According to Table I, it contains 7 basis functions, each monoexitation corresponding to an harmonic frequency ω i , i = 1 • • • 6. So the initial space includes states 1-4, 10, 11 and the ground state, with an energetic hole of 1170 cm -1 between ω 4 and ω 5 that includes five harmonic states. In the A-VCI approach, the starting subspace has no energetic hole: it contains all the functions with harmonic frequencies between 0 and the maximal boundary of the domain of interest. Therefore, 23 basis functions form B (0) to study the 0-4 000 cm -1 frequency range with A-VCI.

Here we compare the results given by the first iterations of VP and A-VCI for the 23 lowest eigenvalues of the vibrational problem. Our goal is to evaluate how fast the convergence is according to the growth of the basis. The absolute errors on eigenvalues are computed with respect to the reference values of Table I and absolute errors on eigenvectors are estimated by

1 -|X ref • X|, (14)
where X ref is the reference eigenvector computed in B 21 and X the one obtained with VP or A-VCI algorithm. The relative residue R (j) /E (j) of the eigenvalue E (j) of H (j) at iteration j is computed. A comparison with VP is presented on Fig. 1 for the 8 frequencies having an IR intensity greater than 2 km/mol. At iteration 1 VP needs 455 functions compared to the 1 218 of A-VCI and there is no significant gap between the two methods in the calculation of the fundamental transitions represented by frequencies 1-4, 10 and 11. On the contrary A-VCI gives a much better approximation of frequencies 9 and 13 than VP with a difference of one order of magnitude on the different error indicators. The particularity of these two frequencies is that their corresponding harmonic states are outside the starting space of VP while contained in the one of A-VCI. The 9th harmonic frequency belongs to the energetic FIG. 1. Absolute error on (a) H 2 CO active frequencies (cm -1 ) with an IR intensity higher than 2 km/mol, (b) associated eigenvectors and (c) relative residue for VP and A-VCI, on a logarithmic scale. At each iteration all interacting configurations are added to the basis: the size of the resulting subspace is given in parentheses. The initial space contains 7 or 23 vibrational configurations for VP (the fundamental and its monoexitations) and A-VCI (all the harmonic functions φ 0 n for which ν 0 n ≤ ν max ) respectively. The reference calculation for absolute errors is presented in Table I.

hole between between ω 4 and ω 5 , and the 13th one is greater than ω 6 . Therefore, the VP calculation of frequencies ranging from 0 to 4000cm -1 requires additional iterations to compensate for the lack of states 5-9 and 12-23 in its initial subspace. Whereas VP needs two iterations and 5 002 basis elements to get an error less than 1 cm -1 on the 23 targeted eigenvalues, a maximal error of 1.6 cm -1 is already obtained with A-VCI from iteration 1, with 1 218 basis functions. We observe the same behavior on eigenvectors and relative residues. Compared to VP, A-VCI provides a good accuracy with smaller VCI subspaces.

The fast convergence of eigenvectors impacts directly the convergence of intensities. Since the A-VCI starting space has no energetic hole, the minimal needed information is already available at the beginning of the calculation, and then the convergence can be achieved faster with less basis functions than with VP. Finally, we show here that residues, providing error information on eigenvalues and eigenvectors, are good indicators to control the convergence.

Influence of the number of elements in the initial subspace

As shown on Fig. 1, iteration 1 of A-VCI using B (1) = H(B (0) ) gives relative residues slightly higher than 10 -2 with 1 218 basis functions. What happens if we increase the size m (0) of this initial subspace? This is equivalent to enlarge the frequency range of study [0, ν max ]. To answer this last question, we use in A-VCI the strategy described in ( 12).

We compute frequencies in H(B (0) ) starting from B (0) subspaces corresponding to different values of m (0) . Fig. 2 shows how for the frequencies 2 and 10 the residue decreases with respect to the dimension of the initial basis size. The convergence of ν 2 with a relative residue less than 10 -2 is achieved with 2 162 basis functions which corresponds to a starting space of size m (0) = 63 (ν max = 5 500 cm -1 ). The highest frequencies are the most difficult to converge. To lower the residue under 10 -2 for ν 10 we need 9 085 basis functions with m (0) = 533 (ν max = 10 000 cm -1 ). Even without energetic holes, increasing the size of the initial subspace is obviously not the best way to perform a VCI process efficiently.

Increase of the subspace based on the residual contributions

In this section we consider the A-VCI Algorithm 1 with Algorithm 2 to increase the subspace. Setting the following parameters: ν max = 4 000 cm -1 , F = 23, ε = 10 -2 , we aim to compute the 22 frequencies listed in Table I with a relative residue less than 10 -2 . We remind that maximal frequency ν max of 4 000 cm -1 corresponds to a starting space B (0) of m (0) = 23 elements.

On Fig. 3 we observe the residue is below the requested precision at iteration 3, using 2 516 (0) . The frequency ν max (cm -1 ) corresponding to m (0) is reported on the upper axis. The residue is represented for ν 2 and ν 10 computed in H(B (0) ).

basis functions, for the 8 active frequencies with an IR intensity greater than 2 km/mol and the highest one ν 22 . It should be noted that the 2 243 basis functions of iteration 2 already ensure the desired convergence for all frequencies except for the last one ν 22 . To verify the convergence of the algorithm we plot, in Fig. 4, the gaps with the frequencies computed in the reference basis B 21 . Using algorithm 2 to enrich the basis the A-VCI method allows to get eigenvalues with absolute errors below 1 cm -1 with only 2 516 basis functions. To reach such a precision with the strategy ( 12) the 1 218 elements of iteration 1 are not enough and another iteration is needed at the cost of 10 962 basis functions. Compared to the ). The y-axis represents the relative residue and the x-axis the frequency number. For each iteration j the size of the subspace B (j) is given in parentheses.

classical variation-perturbation method, these results confirm the efficiency of the A-VCI approach, where the B (j) basis is enriched by residue-based selected functions. Our new selection method allows to compute frequencies with a good approximation by using less basis functions.

As shown on Fig. 5, using A-VCI or VP without any criterion to select relevant functions dramatically increases the size of the basis over iterations. On the other hand, small residues A-VCI: ε =10 -2 , it. 3 (2516) A-VCI: B (j) =H(B (j-1) ), it. 1 (1218) A-VCI: B (j) =H(B (j-1) ), it. 2 (10962)

FIG. 4. Absolute error (logarithmic scale in cm -1 ) of the A-VCI frequencies computed for H 2 CO with respect to the calculation in the reference basis (see Table I). The size of the subspace B (j) at iteration j is given in parentheses. We compare results obtained with the enrichment strategy (12) (B (j) = H(B (j-1) )) at iteration 1 and 2 and with the Algorithm 1 (ε = 10 -2 ) at iteration 3.

are obtained with few basis elements with a A-VCI process control by a residue-based approach. In addition to a better control of the convergence, the residue is clearly a powerful tool to efficiently grow the successive subspaces in an adaptive VCI process. As it can be seen on Fig. 4, the largest errors are related to the highest frequencies whatever the size of the basis. A relevant description of high frequencies (ν 22 in this example) with the least expanded basis set is a crucial stake in theoretical spectroscopy. The adaptive approach of A-VCI based on the residue is therefore a suitable way to address these challenging problems. Basis size VP (B (j) =H(B (j-1) )) for ν 10 A-VCI (ε =10 -2 ) for ν 10 FIG. 5. Size of the basis for VP (B (j) = H(B (j-1) )) and A-VCI (ε = 10 -2 ) in the case of H 2 CO, on a logarithmic scale and as a function of the relative residue of ν 10 .

B. Acetonitrile molecule, CH 3 CN

In the previous section, we have shown how the A-VCI can accurately compute frequencies of a 6-degree-of-freedom molecule. With these promising results one question arises: does the method also work well on larger systems? One can wonder how the method behaves when the spectrum becomes denser and denser as the molecule grows. In this section we answer this matter by applying the A-VCI to the 12-degree-of-freedom acetonitrile molecule CH 3 CN. The relative residue threshold of the method is fixed to ε = 10 -2 in Algorithm 1 and Algorithm 2. We use the normal coordinate quartic potential of CH 3 CN used in Ref. 27,29,32 , initially introduced in Ref. 30 with a method mixing coupled-cluster and density functional calculations. This potential counts 311 terms (12 quadratic, 108 cubic and 191 quartic) and its quadratic part is inferred from the harmonic frequencies calculated at the CCSD(T)/cc-pVTZ level: 3 065,ω 2 = 2 297,ω 3 = 1 413,ω 4 = 920,ω 5 = ω 6 = 3 149,ω 7 = ω 8 = 1 487, ω 9 = ω 10 = 1 061, ω 11 = ω 12 = 361 (cm -1 ).

ω 1 =

Influence of the frequency domain of interest

In the A-VCI method the frequency range of study [0, ν max ] is directly related to the size m (0) of the initial space to compute the F = m (0) smallest eigenpairs of H. Table II shows for different values of F how the frequency ν max and the sizes of the final spaces change. As there is almost a ratio of 10 between the lowest and the largest harmonic frequencies, even for F = 121 the maximal reached frequency is under 2 600 cm -1 . The two spaces B (j) and TABLE II. Number of iterations and sizes of the final spaces at the convergence of A-VCI (ε = 0.01) for different numbers of computed eigenvalues (F ). Here F = m (0) is also the size of the initial subspace.

F = m (0) ν max (cm -1 ) Number of iterations Size of B (j) Size of B (j) R 3
B (j)
R grow quickly when F increases but the number of iterations to reach the convergence remain small. Even if at the convergence the size of B (j) is large it is smaller than the size of the binomial basis. For instance, B 19 is the smallest binomial basis containing all the functions of B (j) for F = 121 at the last iteration. The total size of B (j) and B Moreover, the choice of 1-D Hermite function products and a PES in SOP form leads to sparse matrices as shown in Table III. For F = 121 all non-zero elements of H

R only take 372 MB at the convergence. The matrix storage is only mandatory for H (j) in order to compute its smallest F eigenvalues. But the memory size of H (j) is small compared to the memory needed for H (j) R . Table III shows that there is one order of magnitude between the number of non-zero elements of H (j) and H (j) R . As the PES has a small number of terms we get very sparse Hamiltonian matrices. Therefore the memory cost is not an issue in the A-VCI procedure. The storage increases with the number of iterations but according to Table II this number does not seem to blow up when requiring more and more frequencies.

Comparison with a reference calculation

In this part we compare the CH 3 CN frequencies obtained with A-VCI to a spectrum computed by Avila et al. 27 The reference frequencies result from calculations in a VCI space of size 743 103, which can be decomposed into two independent spaces of 374 251 and 368 852 basis functions thanks to the molecule symmetry. The initial subspace B (0) of the A-VCI is given by the first m (0) = 121 harmonic functions of CH 3 CN with associated frequencies ranging from 0 to 2 527 cm -1 in order to compute the first 120 frequencies of CH 3 CN.

The largest absolute error and residue of all computed frequencies are presented in Table IV for each iteration. The growth rate of the basis is important at the beginning of the algorithm, and decreases from iteration 3 where 92 frequencies have converged over the 120 targeted ones. When the main part of the spectrum has converged few contributions are added to the basis in the last iterations to decrease the residue under ε. A total number of 78 859 states gives at iteration 4 of the method a maximal relative residue lower than 10 -2 with a corresponding absolute error below 1 cm -1 .

The absolute error and the residue at iteration 4 of A-VCI are reported in

V. CONCLUSION

In this work, we have proposed an efficient VCI algorithm with an iterative basis selection technique based on the computation of the residual contribution between two nested spaces.

Compared to the classical approach using the image space we control the basis size to avoid a too quick growth of the space. Moreover, we have shown that the choice of the initial space including all harmonic configurations corresponding to an energetic interval without powerful strategies should be tested to build a smaller space that ensures the convergence at the targeted accuracy and with an affordable computational cost.

The AVCI-algorithm is based on two key ideas: 1) for a given space B, we can decompose its image by the Hamiltonian in a direct sum between the two spaces B and B R ; 2) we can construct the residue in B R to select some unit-vectors to increase the active space. This algorithm could be easily adapted to other sets of coordinates, Hamiltonian structure, basis functions ... as soon as the two above ingredients are valid. Typically, our approach works if we are able to store or to compute the sparse matrix representing the discretization of the Hamiltonian in the full approximation space.

Algorithm 1 :

 1 Adaptive vibrational configuration interaction (A-VCI) algorithm beginInitialization: construct the initial orthonormal basis B(0)

Algorithm 2 :

 2 Procedure to enlarge the search space at each iteration Input: The smallest F eigenpairs (E

FIG. 2 .

 2 FIG.2. Relative residue (blue curves) and corresponding H(B (0) ) basis size (red curve) as function of the size m (0) of B(0) . The frequency ν max (cm -1 ) corresponding to m (0) is reported on the upper

  4 % of the 141 120 525 functions of B 19 .

  any energetic hole is a crucial point. Another key ingredient is the selection technique implemented to extend iteratively the basis of the harmonic configurations. The selection is based on a posteriori error estimator that allows to add only the most relevant directions in the search space. The first tests on the vibrational spectra of H 2 CO (6-D space) and CH 3 CN (12-D space) confirm the good performance of the A-VCI algorithm. The quality of the results is the same as those obtained with the same PES in reference calculations performed in large basis sets. The reduced dimension of the A-VCI final basis significantly lowers the computational cost of the IR frequencies. The results for CH 3 CN show in high dimension that our selection criterion adds too many directions to the VCI space. More

  27 , Avila et al. introduce the following energetic criterion

			N
			α i n i ≤ d, with
	α i =	ω i ω min	i=1 and ω min the lowest harmonic frequency, to drastically reduce the degree
	of Hermite functions in high energy directions. The last trick to decrease the size of the
	problem is to take advantage of the molecular symmetry when possible 27,30 . For larger
	molecules the tensor decomposition was introduced by 31 in the context of the VCC theory
	and some numerical problems related to the canonical decomposition have been pointed
	out. Recently, Leclerc et al. 29 use the tensor structure of the operator (1) and the basis
	functions (2) to propose a promising approach based on tensor decomposition to solve the
	problem in Π d : the reduced-rank block power method (RRBPM). Although RRBPM has
	difficulties to converge for highest eigenvalues, significant improvements are obtained by
	using nested contraction and hierarchical tensor format 32 .

TABLE I .

 I Computed harmonic, anharmonic (ν i ) frequencies (in cm -1 ), and anharmonic infrared intensities (I ν i ) (in km/mol) for H 2 CO in the reference basis B 21 (296 010 harmonic configurations).Expansion coefficients on the harmonic configurations ω k (k = 1 . . . 6) are reported to support the assignments of the transitions.

	Frequency number i Harmonic frequency	ν i	I ν i	Description
	1	1181.2	1145.8	4.20	ω 1 (0.99)
	2	1261.1	1234.8	11.44	ω 2 (0.99)
	3	1529.0	1492.3	7.50	ω 3 (0.98)
	4	1765.4	1740.6	65.45	ω 4 (0.97)
	5	2362.4	2289.3	1.40	2ω 1 (0.98)
	6	2442.3	2385.2	0.00	ω 1 + ω 2 (0.99)
	7	2522.2	2468.6	1.14	2ω 2 (0.98)
	8	2710.2	2635.9	0.01	ω 1 + ω 3 (0.99)
	9	2790.1	2700.7	21.44	ω 2 + ω 3 (0.87), ω 6 (0.44)
		2932.5	2787.9	51.30	ω 5 (0.91)
		3000.5	2828.4	72.62	ω 6 (0.81), ω 2 + ω 3 (0.46)
		2946.6	2880.8	0.00	ω 1 + ω 4 (0.97)
		3026.5	2977.4	6.09	ω 2 + ω 4 (0.87), 2ω 3 (0.37)
		3058.0	2983.3	0.93	2ω 3 (0.91), ω 2 + ω 4 (0.35)
		3294.4	3226.1	0.30	ω 3 + ω 4 (0.96)
		3543.6	3422.2	0.00	3ω 1 (0.95)
		3530.8	3469.8	0.31	2ω 4 (0.93)
		3623.5	3529.7	0.00	2ω 1 + ω 2 (0.97)
		3703.4	3621.1	0.00	ω 1 + 2ω 2 (0.97)
		3783.3	3694.5	0.01	3ω 2 (0.93)
		3891.4	3774.3	0.00	2ω 1 + ω 3 (0.96)
		3971.3	3822.2	0.00 ω 1 + ω 2 + ω 3 (0.71), ω 1 + ω 6 (0.65)

TABLE III .

 III Number of non-zero elements (NNZ) in H and H R matrices and memory used to store them. In parenthesis, we give the percentage of non-zero elements in the matrix. Results are given at the convergence for different values of F .

		H (j) at the convergence	H R at the convergence (j)
	F	NNZ	Memory (MB)	NNZ	Memory(MB)
	3	18 080 (6%)	0.2	280 743 (0.8%)	3
	7	75 614 (3%)	0.6	884 620 (0.2%)	8
	27	835 579 (0.8%)	7	5 826 320 (0.06%)	46
	51	2 417 778 (0.4%)	19	13 599 134 (0.03%)	109
	121	11 383 967 (0.2%)	91	46 447 775 (0.01%)	372

TABLE IV .

 IV Convergence of A-VCI frequencies. The absolute error is evaluated with respect to Ref. 27 .

	Iteration number	B (j)	Largest abs.	Largest rel.	Number of
	j	size	error (cm -1 )	residue	converged frequencies
	0	121	363.37	1.13 × 10 -1	0
	1	7 694	234.12	7.77 × 10 -2	7
	2	52 007	176.69	7.09 × 10 -2	92
	3	70 080	20.33	3.68 × 10 -2	114
	4	78 859	0.91	8.85 × 10 -3	120

  Table V for the 120 targeted frequencies. While 743 103 basis functions are needed in Ref.27 , A-VCI only use 78 859 basis functions to provide computational results of the same accuracy: the absolute error between the frequencies calculated by the two methods does not exceed 1 cm -1 and the associated residue is always under 10 -2 . It should also be noted that the symmetry properties of CH 3 CN taken into account in Ref.27 are not for A-VCI. This comparison with a reference calculation on CH 3 CN illustrates the ability of A-VCI to study large vibrational molecular systems.

TABLE V :

 V Absolute error on frequencies and relative residue at iteration 4 of A-VCI (ε = 10 -2 ) applied to CH 3 CN. Ref. 27 10 -3 ω 10 + ω 11 + 2ω 12 (0.61), ω 9 + 2ω 11 + ω 12 (0.60) 10 -3 ω 9 + ω 11 + 2ω 12 (0.61), ω 10 + 2ω 11 + ω 12 (0.55) 10 -3 ω 10 + 2ω 11 + ω 12 (0.57), ω 9 + ω 11 + 2ω 12 (0.52) 10 -3 ω 4 + ω 11 + 3ω 12 (0.67), ω 4 + 3ω 11 + ω 12 (0.62) + 2ω 11 + 2ω 12 (0.67), ω 4 + 4ω 11 (0.45) 10 -3 ω 4 + 3ω 11 + ω 12 (0.63), ω 4 + ω 11 + 3ω 12 (0.63) + 4ω 11 (0.63), ω 4 + 4ω 12 (0.63) 94 2346.59 (A 1 ) 0.38 6.85 × 10 -3 ω 4 + 2ω 11 + 2ω 12 (0.61), ω 4 + 4ω 12 (0.46) + ω 10 + ω 11 (0.67), 2ω 10 + ω 12 (0.49) 102 2425.48 (A 2 ) 0.44 6.67 × 10 -3 ω 9 + ω 10 + ω 12 (0.60), 2ω 9 + ω 11 (0.52) + ω 10 + ω 11 (0.57), ω 3 + ω 10 (0.55) a + ω 10 + ω 12 (0.66), ω 3 + ω 9 (0.52) a 10 -3 ω 3 + 2ω 11 + ω 12 (0.65), ω 9 + 3ω 11 + ω 12 (0.38) 10 -3 ω 10 + 2ω 11 + 2ω 12 (0.58), ω 9 + ω 11 + 3ω 12 (0.48) 10 -3 ω 9 + 2ω 11 + 2ω 12 (0.56), ω 10 + ω 11 + 3ω 12 (0.55) 111 2485.09 (A 1 ) 0.03 3.92 × 10 -3 ω 9 + ω 11 + 3ω 12 (0.59), ω 10 + 4ω 11 (0.56) 112 2485.09 (A 2 ) 0.75 5.91 × 10 -3 ω 10 + 4ω 11 (0.49), ω 10 + 3ω 11 + ω 12 (0.46) 113 2485.87 (E) 0.03 5.77 × 10 -3 ω 10 + 4ω 11 (0.52), ω 10 + 3ω 11 + ω 12 (0.41) 10 -3 ω 3 + 2ω 11 + ω 12 (0.45), ω 10 + 2ω 11 + 2ω 12 (0.43) 10 -3 ω 9 + 2ω 11 + 2ω 12 (0.54), ω 3 + ω 11 + 3ω 12 (0.45) 10 -3 2ω 4 + ω 11 + ω 12 (0.85), 3ω 4 + ω 11 + ω 12 (0.34) + 2ω 11 (0.62), 2ω 4 + 2ω 12 (0.59)

	21 51 81	1451.10 (E) 1981.85 (A 1 ) 2206.78 (E)	0.29 0.06 0.49	6.19 × 10 -3 2.70 × 10 -3 4.70 × 10 -3	2ω 11 + 2ω 12 (0.78), 4ω 11 (0.44) ω 4 + 2ω 11 + ω 12 (0.81), ω 4 + 3ω 12 (0.44) ω 7 + 2ω 11 (0.76), ω 8 + ω 11 + ω 12 (0.49)
	22 52 82	1452.83 (E) 1981.85 (A 2 ) 2207.56 (E)	0.08 0.09 0.03	2.15 × 10 -3 4.24 × 10 -3 3.85 × 10 -3	4ω 12 (0.68), 4ω 11 (0.63) ω 4 + ω 11 + 2ω 12 (0.77), ω 4 + 3ω 11 (0.51) ω 8 + 2ω 12 (0.77), ω 8 + 2ω 11 (0.52)
	23 53 83	is used to compute the absolute error. The reference uses 1452.83 (E) 0.07 2.37 × 10 -3 3ω 11 + ω 12 (0.66), ω 11 + 3ω 12 (0.65) 1982.86 (E) 0.03 3.15 × 10 -3 ω 4 + 3ω 12 (0.80), ω 4 + 2ω 11 + ω 12 (0.44) 2207.56 (E) 0.62 8.85 × 10 -3 ω 7 + 2ω 12 (0.9), ω 7 + 2ω 11 (0.34)
	743 103 basis functions and takes advantage of the symme-try properties of the molecule, written in parentheses. At iteration 4, A-VCI uses 78 859 functions, with no symmetry assumptions. The last column reports the normal mode de-scription. A-VCI (ε = 10 -2 ), iteration 4 Absolute error (cm -1 ) Rel. residue Description 360.99 (E) Ref. 27 (cm -1 ) 0.07 2.72 × 10 -3 ω 12 (0.97) 360.99 (E) 0.06 3.01 × 10 -3 ω 11 (0.99) 723.18 (E) 0.04 3.12 × 10 -3 2ω 11 (0.69), 2ω 12 (0.67) 723.18 (E) 0.06 5.24 × 10 -3 ω 11 + ω 12 (0.97) 723.83 (A 1 ) 0.05 3.14 × 10 -3 2ω 12 (0.69), 2ω 11 (0.67) 900.66 (A 1 ) 0.08 2.24 × 10 -3 ω 4 (0.95) 1034.13 (E) 0.09 1.81 × 10 -3 ω 9 (0.97) 1034.13 (E) 0.08 2.11 × 10 -3 ω 10 (0.97) 1086.55 (A 1 ) 0.05 2.76 × 10 -3 2ω 11 + ω 12 (0.83), 3ω 12 (0.49) 1086.55 (A 2 ) 0.20 6.04 × 10 -3 ω 11 + 2ω 12 (0.80), 3ω 11 (0.54) 1087.78 (E) 0.08 2.21 × 10 -3 3ω 12 (0.83), 2ω 11 + ω 12 (0.49) 1087.78 (E) 0.11 5.19 × 10 -3 3ω 11 (0.80), ω 11 + 2ω 12 (0.54) 1259.82 (E) 0.08 2.13 × 10 -3 ω 4 + ω 12 (0.94) 1259.82 (E) 0.08 2.18 × 10 -3 ω 4 + ω 11 (0.94) 1388.97 (A 1 ) 0.08 1.99 × 10 -3 ω 3 (0.73), ω 9 + ω 11 (0.45) 1394.69 (E) 0.07 2.23 × 10 -3 ω 9 + ω 11 (0.68), ω 10 + ω 12 (0.68) 1394.69 (E) 0.03 4.03 × 10 -3 ω 9 + ω 12 (0.83), ω 10 + ω 11 (0.49) 1394.91 (A 2 ) 0.06 4.93 × 10 -3 ω 10 + ω 11 (0.83), ω 9 + ω 12 (0.49) 1397.69 (A 1 ) 0.08 2.02 × 10 -3 ω 3 (0.63), ω 10 + ω 12 (0.51) 1451.10 (E) 0.07 2.36 × 10 -3 ω 11 + 3ω 12 (0.69), 3ω 11 + ω 12 (0.68) 1453.40 (A 1 ) 0.08 4.73 × 10 -3 2ω 11 + 2ω 12 (0.56), 4ω 12 (0.55) 1483.23 (E) 0.06 2.44 × 10 -3 ω 8 (0.97) 1483.23 (E) 0.12 4.99 × 10 -3 ω 7 (0.97) 1620.22 (E) 0.08 2.24 × 10 -3 ω 4 + ω 11 + ω 12 (0.93) 1620.22 (E) 0.07 2.43 × 10 -3 ω 4 + 2ω 11 (0.67), ω 4 + 2ω 12 (0.65) 1620.77 (A 1 ) 0.07 2.51 × 10 -3 ω 4 + 2ω 12 (0.66), ω 4 + 2ω 11 (0.65) 1749.53 (E) 0.03 3.06 × 10 -3 ω 3 + ω 12 (0.65), ω 10 + 2ω 12 (0.44) 1749.53 (E) 0.01 3.30 × 10 -3 ω 3 + ω 11 (0.65), ω 9 + 2ω 11 (0.44) 1756.43 (A 1 ) 0.09 4.83 × 10 -3 ω 9 + ω 11 + ω 12 (0.70), ω 10 + 2ω 12 (0.51) 1756.43 (A 2 ) 0.31 6.25 × 10 -3 ω 10 + ω 11 + ω 12 (0.71), ω 9 + 2ω 11 (0.54) 1757.13 (E) 0.19 6.01 × 10 -3 ω 10 + 2ω 11 (0.88), ω 9 + ω 10 + ω 11 (0.31) 1757.13 (E) 0.61 8.03 × 10 -3 ω 9 + 2ω 12 (0.91) 1759.77 (E) 0.03 3.15 × 10 -3 ω 3 + ω 12 (0.62), ω 10 + 2ω 12 (0.59) a 1759.77 (E) 0.04 4.09 × 10 -3 ω 3 + ω 11 (0.61), ω 9 + 2ω 11 (0.57) a 1785.21 (A 1 ) 0.12 2.86 × 10 -3 2ω 4 (0.89), 3ω 4 (0.31) 1816.80 (E) 0.16 5.05 × 10 -3 2ω 11 + 3ω 12 (0.78), 4ω 11 + ω 12 (0.49) 1816.80 (E) 0.38 6.57 × 10 -3 3ω 11 + 2ω 12 (0.71), ω 11 + 4ω 12 (0.59) 1818.95 (A 1 ) 0.13 4.60 × 10 -3 5ω 11 (0.61), ω 11 + 4ω 12 (0.60) 1818.95 (A 2 ) 0.20 5.21 × 10 -3 4ω 11 + ω 12 (0.63), 5ω 12 (0.60) 1820.03 (E) 0.07 4.53 × 10 -3 5ω 12 (0.65), 2ω 11 + 3ω 12 (0.51) 1820.03 (E) 0.18 5.50 × 10 -3 5ω 11 (0.63), 3ω 11 + 2ω 12 (0.58) 1844.26 (A 2 ) 0.05 2.10 × 10 -3 ω 7 + ω 12 (0.96) 1844.69 (A 1 ) 0.43 2.29 × 10 -3 ω 7 + ω 11 (0.70) , ω 8 + ω 12 (0.66) 1844.33 (E) 0.15 5.57 × 10 -3 ω 8 + ω 11 (0.96) 1844.33 (E) 0.29 2.31 × 10 -3 ω 8 + ω 12 (0.70), ω 7 + ω 11 (0.66) 1931.55 (E) 0.01 3.25 × 10 -3 ω 4 + ω 9 (0.94) 1931.55 (E) 0.03 3.65 × 10 -3 ω 4 + ω 10 (0.94) 1982.86 (E) 0.04 3.90 × 10 -3 ω 4 + 3ω 11 (0.76), ω 4 + ω 11 + 2ω 12 (0.50) 2057.07 (A 1 ) 0.02 3.04 × 10 -3 2ω 9 (0.68), 2ω 10 (0.69) 2065.29 (E) 0.04 2.78 × 10 -3 2ω 10 (0.68), 2ω 9 (0.68) 2065.29 (E) 0.57 7.54 × 10 -3 ω 9 + ω 10 (0.96) 2111.38 (E) 0.05 3.65 × 10 -3 ω 3 + 2ω 12 (0.49), ω 3 + 2ω 11 (0.49) 2111.38 (E) 0.09 4.26 × 10 -3 ω 3 + ω 11 + ω 12 (0.69), ω 10 + ω 11 + 2ω 12 (0.40) 2112.30 (A 1 ) 0.02 3.30 × 10 -3 ω 3 + 2ω 11 (0.55), ω 3 + 2ω 12 (0.52) 2119.33 (E) 0.15 4.72 × 62 Frequency number 1 2 3 4 5 6 7 8 9 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 54 55 56 57 58 59 60 61 2119.33 (E) 0.41 6.33 × 63 2120.54 (E) 0.30 5.60 × 10 -3 ω 10 + 3ω 11 (0.68), ω 9 + 3ω 12 (0.67) 64 2120.54 (E) 0.38 6.12 × 65 2120.91 (A 2 ) 0.25 5.38 × 10 -3 ω 9 + 3ω 12 (0.61), ω 10 + 3ω 11 (0.59) 66 2122.83 (E) 0.08 4.15 × 10 -3 ω 3 + 2ω 12 (0.44), ω 10 + 3ω 12 (0.43) a 67 2122.83 (E) 0.10 4.38 × 10 -3 ω 3 + ω 11 + ω 12 (0.61), ω 9 + 2ω 11 + ω 12 (0.54) a 68 2123.30 (A 1 ) 0.07 4.05 × 10 -3 ω 9 + 3ω 11 (0.44), ω 3 + 2ω 11 (0.43) 69 2142.61 (E) 0.15 3.90 × 10 -3 2ω 4 + ω 12 (0.87) 70 2142.61 (E) 0.13 4.05 × 10 -3 2ω 4 + ω 11 (0.87) 71 2183.64 (E) 0.15 4.40 × 10 -3 3ω 11 + 3ω 12 (0.75), ω 11 + 5ω 12 (0.44) 72 2183.64 (E) 0.78 8.13 × 10 -3 4ω 11 + 2ω 12 (0.66), 2ω 11 + 4ω 12 (0.64) 73 2186.14 (E) 0.14 4.35 × 10 -3 ω 11 + 5ω 12 (0.69), 5ω 11 + ω 12 (0.65) 74 2186.14 (E) 0.53 6.02 × 10 -3 6ω 11 (0.6), 6ω 12 (0.6) 75 2187.64 (E) 0.10 4.17 × 10 -3 6ω 11 (0.53), 6ω 12 (0.53) 76 2187.64 (E) 0.18 4.61 × 10 -3 3ω 11 + 3ω 12 (0.57), 5ω 11 + ω 12 (0.55) 77 2188.14 (A 1 ) 0.56 7.41 × 10 -3 2ω 11 + 4ω 12 (0.54), 4ω 11 + 2ω 12 (0.52) 78 2206.63 (A 1 ) 0.01 3.29 × 10 -3 ω 7 + ω 11 + ω 12 (0.83), ω 8 + 2ω 12 (0.44) 79 2206.63 (A 2 ) 0.13 4.75 × 10 -3 ω 8 + ω 11 + ω 12 (0.83), ω 7 + 2ω 11 (0.49) 80 2206.78 (E) 0.12 4.85 × 10 -3 ω 8 + 2ω 11 (0.78), ω 7 + ω 11 + ω 12 (0.42) 84 2250.75 (A 1 ) 0.04 3.17 × 10 -3 ω 2 (0.90) 85 2287.94 (A 1 ) 0.02 3.93 × 10 -3 ω 4 + ω 9 + ω 11 (0.61), ω 4 + ω 10 + ω 12 (0.59) 86 2290.22 (E) 0.01 3.31 × 10 -3 ω 4 + ω 9 + ω 12 (0.91) 87 2290.22 (E) 0.02 4.04 × 10 -3 ω 4 + ω 10 + ω 12 (0.66), ω 4 + ω 9 + ω 11 (0.65) 88 2290.38 (A 2 ) 0.24 6.14 × 10 -3 ω 4 + ω 10 + ω 11 (0.91) 89 2297.19 (A 1 ) 0.03 3.22 × 10 -3 ω 3 + ω 4 (0.82) 90 2344.75 (E) 0.11 2.78 × 91 2344.75 (E) 0.60 6.96 × 10 -3 ω 4 92 2346.13 (E) 0.10 3.04 × 93 2346.13 (E) 0.07 4.45 × 10 -3 ω 4 95 2384.87 (E) 0.00 3.61 × 10 -3 ω 4 + ω 8 (0.95) 96 2384.87 (E) 0.02 3.84 × 10 -3 ω 4 + ω 7 (0.95) 97 2415.02 (E) 0.01 3.66 × 10 -3 2ω 10 + ω 12 (0.68), ω 3 + ω 10 (0.51) 98 2415.02 (E) 0.08 4.26 × 10 -3 2ω 9 + ω 11 (0.68), ω 3 + ω 9 (0.52) 99 2420.14 (E) 0.01 3.96 × 10 -3 2ω 9 + ω 12 (0.65), ω 3 + ω 10 (0.60) 100 2420.14 (E) 0.24 5.60 × 10 -3 2ω 10 + ω 11 (0.65), ω 3 + ω 9 (0.61) 101 2425.44 (A 1 ) 0.07 4.33 × 10 -3 ω 9 103 2427.96 (E) 0.04 4.10 × 10 -3 ω 9 104 2427.96 (E) 0.38 6.50 × 10 -3 ω 9 105 2474.51 (A 1 ) 0.04 2.82 × 106 2474.51 (A 2 ) 0.26 5.16 × 10 -3 ω 3 + ω 11 + 2ω 12 (0.63), ω 3 + 3ω 11 (0.42) 107 2476.17 (E) 0.02 2.82 × 10 -3 ω 3 + 3ω 12 (0.65), ω 10 + 4ω 12 (0.45) 108 2476.17 (E) 0.15 4.35 × 10 -3 ω 3 + 3ω 11 (0.64), ω 9 + 4ω 11 (0.44) 109 2483.41 (E) 0.02 3.81 × 110 2483.41 (E) 0.76 2485.87 (E) 0.91 7.67 × 10 -3 ω 9 + 4ω 12 (0.77) 115 2486.96 (A 1 ) 0.11 4.10 × 116 2486.96 (A 2 ) 0.42 6.77 × 117 2487.79 (E) 0.01 3.24 × 10 -3 ω 10 + 4ω 12 (0.55), ω 3 + 3ω 12 (0.49) 118 2487.79 (E) 0.28 5.74 × 10 -3 ω 9 + 4ω 11 (0.49), ω 3 + 3ω 11 (0.46) 7.55 × 114 119 2501.26 (E) 0.32 3.32 × 120 2501.26 (E) 0.00 5.20 × 10 -3 2ω 4

a The main assignation slightly differs from Ref.

27

.

ACKNOWLEDGMENTS

We thank G. Avila and T. Carrington for providing us the potential energy surface of the CH 3 CN molecule that was used for the comparison tests. We also thank L. Giraud and the two anonymous reviewers for their constructive suggestions on preliminary versions of this manuscript. A part of this work was funded by the grand CRA-HPC of the Conseil Régional d'Aquitaine. Experiments presented in this paper were carried out using the PlaFRIM experimental computing platform, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI, IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux, CNRS and ANR in accordance to the programme Investissements d'Avenir (see www.plafrim.fr). Some numerical computations were also carried out using the MCIA platform (Mésocentre de Calcul Intensif Aquitain, see www.mcia.univ-bordeaux.fr). Finally we acknowledge the Direction du Numérique of the Université de Pau et des Pays de l'Adour for the computing facilities it provided us.