Yann Garniron

Emmanuel Giner

Jean-Paul Malrieu

Anthony Scemama

Alternative definition of excitation amplitudes in multi-reference state-specific coupled cluster

A central difficulty of state-specific Multi-Reference Coupled Cluster (MR-CC) in the multiexponential Jeziorski-Monkhorst formalism concerns the definition of the amplitudes of the single and double excitation operators appearing in the exponential wave operators. If the reference space is a complete active space (CAS), the number of these amplitudes is larger than the number of singly and doubly excited determinants on which one may project the eigenequation, and one must impose additional conditions. The present work first defines a state-specific reference-independent operator T m which acting on the CAS component of the wave function |Ψ m 0 maximizes the overlap between (1 + T m )|Ψ m 0 and the eigenvector of the CAS-SD (Singles and Doubles) Configuration Interaction (CI) matrix |Ψ m CAS-SD . This operator may be used to generate approximate coefficients of the triples and quadruples, and a dressing of the CAS-SD CI matrix, according to the intermediate Hamiltonian formalism. The process may be iterated to convergence. As a refinement towards a strict coupled cluster formalism, one may exploit reference-independent amplitudes provided by (1 + T m )|Ψ m 0 to define a reference-dependent operator T m by fitting the eigenvector of the (dressed) CAS-SD CI matrix. The two variants, which are internally uncontracted, give rather similar results. The new MR-CC version has been tested on the ground state potential energy curves of 6 molecules (up to triple-bond breaking) and two excited states. The non-parallelism error with respect to the full-CI curves is of the order of 1 mE h .

I. INTRODUCTION

The single-reference Coupled Cluster (CC) formalism [1][2][3][4] is the standard technique in the study of the ground state of closed-shell molecules, i.e., those for which a mean-field treatment provides a reasonable zeroth-order single-determinant wave-function Φ 0 . This method incorporates the leading contributions to the correlation energy in a given basis set; it is based on the linked-cluster theorem 5 and is size-extensive since it is free from unlinked contributions. The method generates an approximate wave function under the action of a wave operator Ω acting on the single-determinant reference Φ 0 , and assumes an exponential character to the wave operator

Ψ = ΩΦ 0 = e T Φ 0 . (1)
The most popular version only introduces single and double excitation operators in T , and is known as the Coupled Cluster Singles and Doubles (CCSD) approximation. It incorporates the fourth-order correction of the quadruply excited determinants. The lacking fourth-order contribution concerns the triply excited determinants, which may be added in a perturbative manner. The CC equations, obtained by projecting the eigenequation on each of the Singles and Doubles (SD), lead to coupled quartic equations. In practice, guess values of a) Electronic mail: scemama@irsamc.ups-tlse.fr the amplitudes of the T0→i operators appearing in the T operator may be taken as the coefficients of the singles and doubles |i in the intermediate normalization of the SD Configuration Interaction (CI) vector. The solution of the CC equations may be obtained by treating the effect of the triples and quadruples as an iterative dressing of the SD CI matrix, 6 according to the Intermediate Effective Hamiltonian (IEH) theory. 7,8 The field of application of this method, which satisfies formal requirements and is numerically efficient, is however limited to the systems and the situations where a single-determinant zerothorder description is relevant. This is no longer the case when chemical bonds are broken, creating open shells, as occurs in most of the chemical reactions. The magnetic systems generally present several open shells, and the low spin-multiplicity states are inherently of multiple-determinant character. Due to near degeneracies, most of the excited states are not only of multi-determinantal but of multi-configurational character. The conception of a multi-reference (MR) counterpart of the CCSD formalism is highly desirable, and has been the subject of intense research. The most comprehensive review has been given by Bartlett and his colleagues. 9 For formal reasons and in particular to treat correctly the breaking of bonds, the reference space, or model space, is usually taken as a Complete Active Space (CAS), i.e., the Full-CI (FCI) of a well-defined number of electrons (the active electrons) in a well-defined set of orbitals (the active MOs). The other MOs are called inactive. Let us label |I , |J , . . . the reference determinants. The determinants |i , |j , . . . which interact with the reference space are obtained under purely inactive or semi-active single and double excitations; they generate the CAS-SD CI space, the diagonalization of which provides a size-inconsistent energy E m CAS-SD and the corresponding eigenvector,

|Ψ m CAS-SD = |Ψ m 0 + |Ψ m SD = I ∈CAS C m I |I + i CAS c m i |i (2)
with Ψ m CAS-SD |Ψ m CAS-SD = 1. One strategy, which is not very aesthetic since it breaks the symmetry between degenerate reference determinants, but which has given rather satisfactory results, consists in selecting (eventually in an arbitrary manner) a specific single reference and in introducing in the wave operator the multiple excitations which generate the other references (the other determinants of the model space). 10 A similar procedure was proposed by Li and Paldus which uses specific three and four body amplitudes issued from a MR-CISD function. 11 The other strategies consider all the references on an equal footing, and are really multi-reference. Let us call N the number of references, and n the number of SD determinants. If the treatment pretended to provide N eigenvectors simultaneously, one might define the N × n amplitudes sending from the references to the outerspace determinants, in a unique manner but this state-universal approach is not practicable when the model space is a CAS.

Most of the proposed formalisms are state-specific. In this case one faces the famous multi-parentage problem. This problem is recalled in Section II A. Sufficiency conditions have to be imposed. 12 One solution was proposed by Mukherjee and co-workers, and has been widely tested. [13][14][15] Another one had been proposed earlier by one of us (JPM) and co-workers. 16 It consists, for a given outer-space determinant, in scaling the amplitudes of the various excitation operators TI→i on the interaction between the outer-space determinant and its parents. A recent work has implemented this second solution of the state-specific Multi-Reference Coupled Cluster (MR-CC) problem and has tested its accuracy and robustness on a series of molecular benchmarks, comparing its results to the full-CI energies. 17 In the text, we will refer to this method as λ-MR-CCSD. The present work proposes an alternative process to define the amplitudes of the excitation operators, and this new method will be called µ-MR-CCSD.

The state-specific MR-CC formalisms are usually based on the Jeziorski-Monkhorst 18

We shall leave in a first time this assumption and define in Section II B a reference-independent operator T which acting on the component of the desired state in the model space, |Ψ m 0 , provides a vector as close as possible to the CAS-SD eigenvector. This solution, defining reference-independent amplitudes of the excitations, may be exploited directly to generate approximate values of the coefficients of the triply and quadruply excited determinants, according to the exponential structure of the wave operator. From these coefficients, one may dress the CAS-SD CI matrix, redefine amplitudes, and iterate the process to convergence. This solution, presented in Section II C, is not an MR-CC technique; one may call it an exponential dressing of the CAS-SD CI matrix. Section II D redefines reference-dependent excitation amplitudes from the reference-independent amplitudes by a fitting of the previous amplitudes on the coefficients of the singles and doubles of the (dressed) CAS-SD CI eigenvector. This represents an alternative solution to multi-parentage problem and opens the way to a strict MR-CC formalism. Section III presents a series of numerical tests on the bond breaking of single, double, and triple bonds in ground states of molecules as well as a few tests on excited states. The results are compared to our previous proposal and with Full Configuration Interaction (FCI) results.

II. FORMALISMS

In this section, all the presented formalisms are statespecific. To simplify the notations we will consider that the state superscript m is implicit for the wave functions (Ψ m → Ψ) and for the excitation operators ( T m → T ).

A. The multi-parentage problem in the Jeziorski-Monkhorst approach

Since one wants to produce a MR-CCSD method, one may start from a preliminary CAS-SD CI calculation which will help to fix guess values of the amplitudes of the excitation operators. Let us call |I , |J , . . . the determinants of the CAS, i.e., the so-called reference vectors, and |i , |j , . . . the singles and doubles which do not belong to the CAS and interact with them. The resulting approximate wave function of the targeted state |Ψ is written as

|Ψ CAS-SD = I C I |I + i c i |i . ( 4
)
Although this function is not size consistent, one may note that the coefficients on the CAS determinants are no longer those of the CAS-CI: they incorporate the effect of the dynamical correlation on the composition of the CAS component of the wave function.

In CC formalisms the wave operator Ω is assumed to take an exponential form Ω = exp( T ) (5) and in our previous MR-CC formalism 17 the Jeziorski-Monkhorst multi-exponential structure of the wave operator was adopted, introducing reference-specific wave operators acting specifically on each reference vector (Eq. ( 3)). One may exploit the knowledge of the CAS-SD CI eigenvector to determine guess operators TI defined in such a manner that

|Ψ CAS-SD = I C I TI |I . ( 6
)
Each of the TI operators is a sum of single and double excitations TI→i possible on |I , multiplied by an amplitude t I→i , TI

= i t I→i TI→i . ( 7
)
In the single-reference CC, the amplitudes of the excitation operators are obtained by projecting the eigenequation on the singly and doubly excited determinants; the number of unknowns is equal to the number of equations. This is no longer the case in the MR context: projecting the eigenequation on each of the singly or doubly excited vectors |i is not sufficient to define the amplitudes t I→i since for many classes of excitation, an outer-space determinant interacts with several references, |i = TI→i |I = TJ→i |J . The condition

C i = I t I→i C I ( 8
)
is not sufficient to define the amplitudes, even if one restricts the excitation operators to single and double excitations. Additional constraints have to be introduced to fix the amplitudes, and this is the famous multi-parentage problem. The number of amplitudes is larger than the number of outer-space determinants so that one cannot determine directly the guess values of the amplitudes from Eq. ( 6). Different additional constraints have been proposed. One of them consists in scaling the amplitudes on the Hamiltonian interactions between the references and the outer space determinants,

t I→i t J→i = i| Ĥ |I i| Ĥ |J . ( 9
)
This constraint is expressed as

t I→i = λ i i| Ĥ |I , (10)
where

λ i = c i i| Ĥ |Ψ 0 . ( 11
)
This solution has been recently implemented 17 and shown to provide excellent agreements with full-CI results on a series of molecular problems. From now on, we will refer to this method as the λ-MR-CCSD.

When the term i| Ĥ |Ψ 0 is small, the λ-MR-CCSD presents minor stability problems which may introduce some jitter in the potential energy surfaces. A more important problem of the λ-MRCC is illustrated by considering the case of a non-interacting A • • • B system with localized MOs on A and B. The 2-hole 2-particle inactive double excitations of the type Ti A j B →r A s B have zero amplitude in the λ-MRCC formalism since the integral i A j B |r A s B = 0. The coefficient of the determinant Ti A j B →r A s B |I is not zero but it is equal to the product t i A →r A t j B →s B c I . In Sec. II B, we propose an alternative solution to the multi-parentage problem to define amplitudes which do not suffer from this pathological behavior.

B. Introduction of reference-independent amplitudes

The present method differs from the λ-MR-CCSD in the definition of the amplitudes, introduced in this section. The formalism will leave in the first step the Jeziorski-Monkhorst formulation of the wave operator and will consider the possibility to define a unique state-specific reference-independent operator T , written as a sum of single and double excitation operators,

T = mnpq t mn→pq a † p a † q a n a m + mp t m→p a † p a m , (12a)
= mnpq t mn→pq Tmn→pq + mp t m→p Tm→p , (12b)
where the indices p and q run on the virtual and active MOs and the indices m and n run on the inactive occupied and active MOs, excluding the possible occurrence of 4 active MOs. An operator of this kind (but keeping only the linearly independent combinations of the elementary operators) is used in the internally contracted MR-CC method (ic-MRCC) by Evangelista and Gauss, 19 and by Hanauer and Köhn. 20 A similar and more compact formulation was already suggested by Mahapatra et al. [START_REF] Mahapatra | Advances in Quantum Chemistry[END_REF] Our formalism differs by both the determination of the amplitudes and by the way we use them, as will appear later. The ic-MRCC method determines the amplitudes of the excitations by solving the projected coupled cluster equations, where the amplitudes appear up to quartic terms. Hereafter we exploit the knowledge of the CAS-SD CI eigenvector to determine the guess values of the reference-independent amplitudes. These excitation amplitudes will be used later on to estimate the coefficients of the triples and quadruples, and perform an iterative dressing of the CAS-SD CI matrix introducing the coupling between the singles and doubles with the triples and quadruples.

We propose a criterion to fix the amplitudes t = {t mn→pq , t m→p }. Given the fact that we have at our disposal the CAS-SD wave function, a natural way to solve this overdetermined problem is to minimize the distance between the CAS-SD vector and the vector obtained by applying the (1+ T ) operator on the CAS wave function

arg min t (1 + T )|Ψ 0 -|Ψ CAS-SD = arg min t T |Ψ 0 -|Ψ SD , (13)
T |Ψ 0 being normalized such that T |Ψ 0 = |Ψ SD .

To perform the minimization, we build the N SD × N t transformation matrix A i,mn→pq = i| Tmn→pq |Ψ 0 which maps from the outer space of determinants {|i } to the space of excited wave functions { Tmn→pq |Ψ 0 }, and we search for the vector of amplitudes t which minimizes A.tc by solving the normal equations,

(A † A)t = A † c. (14)
Note that in the single-reference case, A is a permutation matrix and the CAS-SD wave function is exactly recovered. The matrix A is usually so large that the use of standard singular value decomposition (SVD) routines to obtain the least squares solution is prohibitive.

Let us first consider the most numerous 2-hole-2-particle inactive double excitations Tjk→rs . These excitations consist in creating two holes in the doubly occupied orbitals and two particles in the unoccupied orbitals. For each excitation of this kind, as all the involved orbitals are outside of the active space, the number of determinants originating from this process is equal to the number of determinants in the reference. Moreover, each one of these excited determinants is doubly excited with respect to only one determinant |I of the reference, and the excitation degree with respect to all other reference determinants is necessarily higher than two. Therefore, all excited determinants created by such a 2-hole-2-particle process have only one parent in the reference, and the corresponding rows of A contain only one non-zero element located in the jk → rs column with value A i,jk→rs = C I . The condition fixing the amplitude t jk→rs is given by arg min t jk→rs Tjk→rs |Ψ 0 t jk→rs -|Ψ SD (15) which is obtained by minimizing min

t jk→rs I C I t jk→rs - i c i i| Tjk→rs |I 2 ( 16
)
using Eq. ( 14), and this condition turns out to be satisfied using only one non-zero coefficient c i with

t jk→rs = I C I c i I C 2 I . ( 17
)
One may notice that this is the weighted average of the ratios between the coefficients of the doubly excited determinants |i and the coefficient of their unique reference generator,

t jk→rs = 1 I C 2 I I C 2 I c i C I . ( 18
)
The maximum number of non-zero elements per row of A is equal to the number of reference determinants since each excitation operator applied on a reference produces no more than one excited determinant. Hence, for all the remaining active excitations, A remains sparse and we solve Eq. ( 14) using Richardson's iterative procedure 22

     t 0 = A † c, t n+1 = A † c + I -A † A t n , (19)
which may be implemented very efficiently using sparse matrix products.

There are cases where multiple amplitudes applied to different references lead to same determinant: Tjk→rs |I = Tlm→tv |J = |i . The linear system is underdetermined, so there are infinitely many possible amplitudes verifying the equations. Among the infinity of possibilities, the SVD picks one particular solution given by A + c, where A + is the pseudoinverse of A. As this solution minimizes the norm of the amplitude vector, 23 the arbitrariness brought by the null space of A is minimized and one obtains the most sensible solution.

C. Evaluation of the coefficients of triples and quadruples and iterative dressing of the CAS-SD CI matrix

This section recalls the procedure described in our previous work. 17 The so-determined excitation operator T may be used to generate the approximate values of the coefficients of the triples and quadruples as obtained by the action of 1 2 T 2 . Actually one may assume, in the spirit of the internally contracted MR-CC methods, that the wave operator Ω generating the correlated wave function Ψ from Ψ 0 ,

Ψ = ΩΨ 0 , (20)
has an exponential structure,

Ω = exp( T ). ( 21
)
But this form will be simply used to estimate the coefficients of the triply and quadruply excited determinants {|α }, leaving the internally contracted structure of the outer-space. The coefficients of these determinants are estimated as

c α = 1 2 α| T 2 |Ψ 0 . ( 22
)
All the determinants {|α } are generated by applying all the single and double substitutions on the singles and doubles, and filtering out the determinants which are already in the wave function. For each |α one searches for the reference determinants {|I } α which differ by no more than 4 orbital substitutions from |α (its grand-parents). One then identifies the set of all possible complementary excitations as the products of excitations Tp and Tq , which generate |α from every member |I of the set {|I } α , i.e.,

S α = (p, q, I) : ∀|I ∈ {|I } α , Tp Tq |I = |α . ( 23
)
It now straightforward to find the set of singles and doubles with which |α interacts through the matrix elements i| Ĥ |α , namely,

{|i } α = |i : ∀(p, •, •)∈ S α , |i = T † p |α . ( 24
)
For each |i , in the eigenequation

i| Ĥ |i -E c i + J i| Ĥ |J C J + j i i| Ĥ |j c j + α i| Ĥ |α c α = 0, (25)
the coefficient c α is given by the genealogy of |α ,

c α = (p,q,I)∈S α (-1) n(I→α) t p t q C I , (26)
n(I → α) being given by the number of permutations needed to go from |I to |α . One may replace the sum over the α by a dressing of the matrix elements between the determinant |i and the references which are grand-parents of |α ,

i| ∆|I = α i| Ĥ |α (p,q,J)∈S α :(J=I) (-1) n(I→α) t p t q (27) since I i| ∆|I C I = α i| Ĥ |α c α . ( 28
)
The effect of the triples and quadruples is thus incorporated as a change of the columns of the CAS-SD CI matrix corresponding to the interaction between the references and the singles and doubles,

i| Ĥ |i -E c i + J i|( Ĥ + ∆)|J C J + j i| Ĥ |j c j . ( 29
)
This type of dressing was already employed in our previous MR-CC implementation. 17 One will find in the same reference the practical procedure to make the dressed matrix Hermitian without any loss of information. Of course the whole process may be iterated. The diagonalization of the dressed CAS-SD CI matrix provides new values of the coefficients, not only of the singles and doubles which no longer suffer from the truncation but also those of the references: the method is fully decontracted. From the new wave function new amplitudes are obtained, a new dressing is defined and the process is repeated till convergence, which is usually rapidly obtained (3-4 iterations).

As opposed to the λ-MR-CCSD method which uses reference-specific amplitudes, the amplitudes introduced in Sec. II B are reference-independent. As a consequence, the formalism is not a strict MR-CC method since we exploit the CAS-SD CI function which slightly differs from the vector resulting from the action of T on the vector. Although the distance between these two vectors has been minimized, they are not identical, (1 + T )|Ψ 0 |Ψ CAS-SD . Once the T operator has been obtained, one might imagine a contracted exponential formalism calculating T 2 |Ψ 0 and the interaction between T |Ψ 0 and T 2 |Ψ 0 , but this calculation requires to return to the determinants. This formalism would remain internally contracted and would be less accurate than the decontracted procedure we propose. Actually in this version, the deviations of the approximate reference-independent amplitudes from optimal ones, those which would generate the exact coefficients of the singles and doubles, only affect the evaluation of the coefficients of the triples and quadruples, and these deviations represent a minor source of error in the correction restoring the size extensivity. This reliability will be illustrated in the numerical tests.

D. State-specific MR-CC variant

In order to return to MR-CC formalism, one may simply exploit the reference-independent amplitudes as an initial guess to define reference-dependent amplitudes. Currently the determinant |i belonging to the singles and doubles has a coefficient ci in T |Ψ 0 , ci = i| T |Ψ 0 = {(p,I):

( Tp |I = |i )} t p C I , (30)
which differs from the coefficient c i in |Ψ SD . One can define a parameter µ i , specific of the determinant |i ,

µ i = c i ci , (31)
which multiplying with ci will produce the exact coefficient c i of |i in the (dressed) CAS-SD CI eigenvector. So the previous reference-independent amplitudes have now become reference-dependent. The excitation Tp which excites |I to |i (|i = Tp |I ) receives a reference-dependent amplitude

t I→i = t p,I = µ i t p . ( 32
)
The same excitation will receive a somewhat different amplitude when it acts on another reference t p,J t p,I . This version returns to the Jeziorski-Monkhorst formalism as the wave operator again is a sum of reference-specific operators. The so-obtained amplitudes may be exploited to generate the coefficients of the triples and quadruples, and one may follow the same strategy as in our previous formalism, with an iterative column dressing of the interactions between the singles and doubles and the references. In what follows, we will refer to this method as µ-MR-CCSD as it involves the µ i (Eq. ( 31)).

As the overlap between (1 + T )|Ψ 0 and |Ψ CAS-SD has been maximized, the coefficients ci and c i are expected to be very close in particular if c i is large, and the parameter µ i should be close to 1, at least for the determinants which contribute significantly to the wave function. In practice we observe this tendency, but the smallest coefficients are sacrificed during the maximization of the overlap and their µ i can be very far from 1. This introduces some instabilities in the iterations, so we chose to limit the values of µ i in the [-µ max i , µ max i ] range, with

µ max i = 2 + 100 × exp -20 |c i | max j |c j | . ( 33
)
In this way, when c i is large µ i is constrained in the [2,2] range, and when c i is small, µ i is constrained in [102,102].

The effect on the stability of the iterations is significant, and the effect on the energy differences is not noticeable, as seen in Sec. III. Our procedure makes use of an (non-compulsory but convenient) approximation, namely, the fact that we have not subtracted the product of the single excitations (the T 2 1 contributions) from coefficients of the doubles to fix the T2 amplitudes. From a perturbation expansion, one sees that this neglect only introduces fifth-order errors on the energy, which are responsible for small deviations from strict additivity of the energies. Notice that a correct treatment of the T n 1 operations, although tedious, is perfectly conceivable in our formalism and would insure a perfect MR-CC character.

III. NUMERICAL TESTS

In this section, we first numerically evaluate the errors made by the different approximations. Then, we compare the here-proposed dressed CAS-SD and µ-MR-CCSD to the λ-MR-CCSD presented in Ref. 17 on standard benchmark systems. 13,15,19,20,[24][25][26][27][28][29][30][31][32][33] The basis set used is Dunning's cc-pVDZ, 34 and the molecular orbitals were obtained using the CAS-SCF code present in GAMESS. 35 All the following calculations were made using the Quantum Package, [START_REF] Scemama | Quantum package v1.0[END_REF] an open-source program developed in our group. Full-CI energies were obtained using the CIPSI algorithm. [START_REF] Huron | [END_REF][38][39] In all the calculations (full-CI, CAS-SD, and MR-CC), only the valence electrons are correlated (frozen core approximation).

A. Approximations

T 2 1 × T2 , T 3 1 , and T 4 1

To estimate the errors due to the approximate treatment of the T 2 1 × T2 , T 3 1 , and T 4 1 operators, we chose a single-reference example in which the single excitations are important at the CISD level. In the single reference case, all the excitations are of the 2 hole-2 particle type, so the normal equations (Eq. ( 14)) are solved exactly and all the values of µ are equal to 1. The only difference to standard CCSD is the approximate treatment of the T 3 1 and T 4 1 , as explained in Section II D. We have calculated the energy of the FH molecule at a distance of 1.2 Å with the single-reference CCSD programs of GAMESS 40 and Gaussian 09, [START_REF] Frisch | gaussian 09[END_REF] and our µ-MR-CCSD implementation using the Hartree-Fock determinant as reference. The results are presented in Table I. The Hartree-Fock and CISD energies agree up to 10 8 a.u., but our implementation differs from the CCSD by almost ∼-10 -6 a.u. We attribute this difference to the approximation in the T 2 1 × T2 , T 3 1 , and T 4 1 operators, and it represents a relative error of 4.4 × 10 6 on the correlation energy.

To measure the effect of this approximation on the size consistency, we have calculated the energy of the CH 3 radical, with C-H bond lengths of 1.103 Å and H-C-H bond angles of 107.69 • in the 6-31G basis set. We have also calculated the energy of the dimer with an intermolecular distance of 100 Å. The active space of the monomer contains only the singly occupied orbital, and the dimer is an open shell singlet with a CAS(2,2) wave function. For the µ-MR-CCSD calculation, we have constrained the µ i as in Eq. ( 33) or we have let it unconstrained. The results are given in Table II and show that the deviation to additivity of the energy is reduced by an order of magnitude going from the CAS-SD to the dressed CAS-SD, and by two orders of magnitude when including the µ i factors. The constraint on the µ i introduces a small error which is below 10 4 a.u.

B. Bond breaking

For all the applications we compare the dressed CAS-SD and µ-MR-CCSD with the λ-MR-CCSD and the CAS-SD values. Results are also given using the reference-independent dressing of the CAS-SD CI matrix. All the applications are presented as energy differences with respect to the full-CI energy estimated by a CIPSI calculation with a second-order perturbative correction. The smallest and largest values of the CIPSI perturbative corrections along the curves are given in Table III. We empirically estimate the error to the FCI energy to be in the order of 10% of the largest contribution. For ethane and twisted ethylene, which have the largest perturbative corrections, we have performed a larger CIPSI calculation at the point with the largest PT2 contribution. For ethylene the PT2 contribution was reduced to 5.6 mE h , but the total energy changed by only 0.08 mE h . In this case, we can consider that the CIPSI energy of ethylene is converged. In the case of ethylene, the For the full series of compounds, Figure 1 shows the energy difference with respect to the full-CI along the reaction coordinate. Table IV summarizes the non-parallelism errors (NPEs) and the maximum of the error obtained along the curve. The MR-CC treatment reduces the average and maximum error of the CAS-SD with respect to the full-CI by a factor close to 4. The correction is larger when the system involves an important number of inactive electrons (F 2 , C 2 H 6 ) than when this number is small (BeH 2 , N 2 ). One actually knows that the size-consistency error of the CAS-SD treatment increases with the number of inactive electrons; this error disappears in the MRCC treatment, which essentially misses some fourth-order connected effects of the triples.

Single-bond breaking

We present here the single bond breaking of the σ bonds of the C 2 H 6 and F 2 molecules and of the π bond of ethylene. The active spaces were chosen with two electrons in two MOs, the minimal wavefunctions to describe properly the dissociation of the molecules. In the case of ethane, the NPE of the CAS-SD is 5.1 mE h , and is reduced to 3.5 mE h with the µ-MR-CCSD. The curve of the dressed CAS-SD has the lowest NPE (1.3 mE h ). The curves obtained by both MR-CCSD methods give equivalent results, with NPEs of 3.5 and 3.6 mE h . In the case of F 2 the NPE of the dressed CAS-SD is 0.9 mE h and the NPE of the µ-MR-CCSD is 1.5 mE h ; both are better than the NPE of the λ-MR-CCSD which is 3.1 mE h . Also, one can remark here some numerical instabilities in the λ-MR-CCSD where the curve is not smooth.

In the next example, the π bond of ethylene is broken by the rotation of the CH 2 fragments. The CAS-SD has an NPE of 1.5 mE h , and using the dressed CAS-SD reduces the NPE to 0.7 mE h . The µ-MR-CCSD gives an NPE of 0.5 mE h , and the NPE obtained with the λ-MR-CCSD is slightly better with a value of 0.3 mE h .

Insertion of Be in H 2

We present the results obtained by the insertion of a beryllium atom into the H 2 molecule, which is a popular benchmark for MR-CC methods. The reference is still a CAS(2,2) for comparison with the literature, even though this choice of reference is not the most appropriate for a correct description of the reaction. The geometries are given by the relation

z = 2.54 -0.46x (a.u.), (34)
where the beryllium atom is at the origin and the hydrogen atoms are at the coordinates (x, 0, ±z). In this particular case, the µ-MR-CCSD gives an NPE of 1.7 mE h which is larger than the NPE of 1.3 mE h obtained by the λ-MR-CCSD. This is due to only one point of the curve, the maximum which is higher by 0.3 mE h ; all the other points being very close by less than 0.1 mE h . Here, the dressed CAS-SD and the µ-MR-CCSD are equivalent.

Two bond breaking

For breaking two bonds we have used CAS(4,4) wave functions as the reference space. The first example is the simultaneous breaking of the two O-H bonds of the water molecule by stretching. Here, the CAS-SD exhibits an NPE of 1.9 mE h which is significantly improved to 0.2 mE h with the dressed CAS-SD. The µ-MR-CCSD, with an NPE of 0.5 mE h , is slightly more parallel to the full-CI curve than the λ-MR-CCSD which has an NPE of 0.7 mE h .

The second example is the double-bond breaking of ethylene by stretching. One should first clarify that the energy differences in the figure do not match those of the torsion along the bond because in the former example the reference was a CAS (2,2), and here it is a CAS (4,4). Dressing the CAS-SD reduces the NPE from 2.7 mE h to 1.7 mE h . One can remark a discontinuity in the curve at large distances. The µ-MR-CCSD and λ-MR-CCSD slightly improve the NPE to values of 1.7 mE h and 1.8 mE h .

Triple-bond breaking

N 2 is the typical benchmark for breaking a triple bond. Here, we have used a CAS (6,6) reference wave function. At the CAS-SD level, the NPE is 1.7 mE h , and the dressed CAS-SD does not reduce the NPE. Here, it is necessary to use referencedependant amplitudes to recover a low NPE: 0.9 mE h with the λ-MR-CCSD, and 0.3 mE h with the µ-MR-CCSD.

C. Excited states

Triplet state of F 2

We report here calculations on the triplet state 3 Σ + u of F 2 . The reference wave function was prepared in two different ways, both using restricted open-shell Hartree-Fock molecular orbitals. The first reference wave function labeled m s = 1 is a single open-shell determinant, and the second wave function is the triplet m s = 0, made of two determinants 1/ √ 2(αβ -βα). To ensure that the CAS-SD is a strict eigenfunction of the Ŝ2 operator, we have included in Ψ SD all the determinants with the same space part as the singles and doubles with respect to the CAS. These determinants are treated in the same way as singles and doubles and are treated variationally in the diagonalizations. Of course, those which are triples or quadruples with respect to Ψ ref are excluded from the set of the {α} and have no effect in the dressing.

To reduce the computational cost, the triples and quadruples were not augmented with all the determinants with the same space part. The absence of some determinants gives rise to a slight deviation (<10 -6 a.u.) of Ŝ2 from the desired eigenvalue, and it is expected to have some impact on the iterative dressing. It is worth checking the effect of this deviation from the exact spin multiplicity. The first test concerns the comparison of the m s = 0 and m s = 1 components of a triplet state.

According to Figure 2, in all the cases, the NPE of the CAS-SD (2.6 mE h ) is improved to a value of 1.1-1.3 mE h . As expected the dressed CASSD and the µ-MR-CCSD are strictly equivalent for m s = 1. Indeed, for both variants, the usual single-reference amplitudes c i /c 0 are recovered. The amplitudes of the λ-MR-CCSD lower the curve by 1 mE h when going from m s = 1 to m s = 0. The dressed CAS-SD gives a slightly higher energy by only 0.3 mE h , and introducing the reference-dependence via the µ i reduces the difference to 0.2 mE h .

If one considers the error on the singlet-triplet gap with respect to the full-CI reference, it appears clearly that the µ-MR-CCSD with m s = 0 gives the most accurate results, with errors lying between 0 and 0.9 mE h along the curve.

Avoided crossing in FH and LiF

We have calculated the potential energy surfaces of the two lowest 1 Σ + states of FH, using as reference wave function the CAS(2,2) with state-averaged CAS-SCF molecular orbitals in the aug-cc-pVDZ basis set. Figure 3 shows the NPEs of the ground and excited states. In the ground state, the NPE is 1.8 mE h for both MR-CCSD variants, but the λ-MR-CCSD shows some numerical instabilities, as opposed to the µ-MR-CCSD which gives a very smooth curve.

In the excited state, the situation is different: surprisingly the best NPE is obtained by the CAS-SD. The reason is the following. At large distances, the CAS-SD description is correct, but the size-consistency error raises the energy. At short distances, the CAS-SD description of the excited state is not as accurate as for large distances since there are determinants with a large coefficient which are not in the active space: the CAS contributes by 0.85 to the norm of the FCI wave function, but in the case of the CAS-SD its contribution is 0.91. This bad description of the CAS-SD raises also the energy at short distances, and the errors compensate along the curve. All the other methods correct the size-consistency error, so the long range errors are corrected but not the error due to the incompleteness of the CAS at short range. This explains why the deviations to the FCI decrease with the distance, and why the NPEs are so large. The two variants of the MR-CCSD agree at short and long distances, but they differ significantly between 1.5 and 2.5 Å; the region of the avoided crossing. To understand these differences, we have plotted in Figure 4 the two eigenvalues of the two state-specific Hamiltonians-one dressed for the ground state and one dressed for the excited state. It appears that the state of interest is very well described, but the dressing for other root has a much lower quality. This strong state-specific character is due to the fitting procedure which is implicitly weighted by the state of interest. The large coefficients have a higher quality in the amplitudes, but the FIG. 4. Potential energy surfaces of the two lowest 1 Σ + states of FH with the µ-MR-CCSD method. The energy of the state corresponding to the dressing is plotted in plain curves, and the energy of the other state is plotted in dashed curves.

important determinants of the second root are usually not the same as in the state of interest. The λ-MR-CCSD has amplitudes which depend less on the wave function, so the quality is comparable on both states, and the choice of these amplitudes is better suited for calculating excited states within the same symmetry, as will be confirmed by the next example.

In Figure 5 we have represented the avoided crossing of LiF, also calculated with the aug-cc-pVDZ basis set. The physical situation is similar to FH, but the energy difference between the ground and the excited states is much smaller. A striking result is that the λ-MR-CCSD, although being statespecific, is able to reproduce very well the whole potential energy surfaces of both states. The position of the avoided crossing is very well reproduced by the three methods: the CAS-SD crosses at 6.3 Å, the full-CI crosses at 6.8 Å, and the dressed CAS-SD and the two MR-CCSD variants cross at 6.9 Å. The µ-MR-CCSD and λ-MR-CCSD coincide in the short-range (≤5 Å) and in the long range (≥7.2 Å), but when the two states become very close in energy in the region of the crossing the dressed CAS-SD and the µ-MR-CCSD are unable to give reasonable values. This disappointing result motivates a future work on a multi-state µ-MR-CCSD.

D. Sensitivity to the choice molecular orbitals

The µ-MR-CCSD algorithm we propose is in the Jeziorski-Monkhorst framework, so it is not invariant with respect to the choice of molecular orbitals. In this section, we checked its sensitivity to the choice of the MO set by FIG. 5. Potential energy surfaces of the two lowest 1 Σ + states of LiF. FIG. 6. Comparison between pseudo-canonical (dashed curves) and localized (plain curves) MOs in F 2 and N 2 . Difference with respect to the full-CI energy.

comparing results obtained with pseudo-canonical CAS-SCF orbitals and with localized MOs in the F 2 and N 2 molecules (Figure 6).

In the F 2 molecule, using localized MOs is a better choice than the pseudo-canonical MOs. The best NPE is obtained by the µ-MR-CCSD method with a value of 0.9 mE h . In the case of N 2 , the situation is different: the NPE of the λ-MR-CCSD goes from 0.9 mE h to 3.6 mE h , and the NPE of the µ-MR-CCSD is relatively stable around 0.3-0.4 mE h . On the other hand, the dressed CAS-SD gives a better NPE with local orbitals, going from 1.7 mE h to 1.1 mE h .

The fact that the µ-MR-CCSD is less sensitive to the MO set than the λ-MR-CCSD can be understood. By changing the MO set, a single excitation rotates into a combination of single and double excitations. In the λ-MR-CCSD method, the amplitudes are calculated by taking into account the matrix elements of the Hamiltonian, which are of different nature depending on the degree of excitation, so the amplitudes are expected to change significantly. In the µ-MR-CCSD variant, the amplitudes are adjusted in such a way that they fit the CAS-SD wave function, which is invariant by rotation of the MOs. Therefore, it is expected to be more robust with respect to the change of MO set.

IV. CONCLUSIONS

We have proposed a method to determine referenceindependent amplitudes by fitting the CAS-SD CI vector. These amplitudes may be used to perform a state-specific iterative dressing of the CAS-SD Hamiltonian in order to take into account the effect of the triples and quadruples in the spirit of the coupled cluster formalism. Alternatively, these amplitudes may be rescaled to reproduce the exact coefficients of the singles and doubles to introduce a reference-dependent character.

In that case, the CAS-SD CI vector is recovered by the application of (1 + T ) on the reference wave function, so we reach here the Jeziorski-Monkhorst coupled cluster formalism.

The CAS-SD dressed with reference-independent amplitudes gives excellent results for single-bond breaking (F 2 and ethane) and the simultaneous breaking of the two O-H bonds of water, with a non-parallelism error lower than the milli-Hartree. When the active space becomes larger, it is necessary to go to the reference-dependent MR-CCSD introducing the µ factors in Eq. (31). In the case of ethylene and N 2 , this keeps the NPE to a value close to the milli-Hartree.

We have shown numerically that the here-proposed amplitudes are not very sensitive to the value of m s for open-shell systems, and to the choice of the molecular orbitals. This is clearly an improvement compared the amplitudes proposed earlier. 17 But we have also shown that the former amplitudes are a better choice when computing excited states of the same symmetry because the here-proposed amplitudes have a much more pronounced state-specific character which may be disadvantageous if the states are too close in energy. This problem can be cured by leaving the state-specific formalism for a multi-state formalism, [START_REF] Malrieu | [END_REF] and this will be the object of a future work.

  splitting of the wave operator into a sum of operators acting individually on the various references T m = I T m I |I I |.

FIG. 2

 2 FIG. 2. F 23 Σ + u . Difference with respect to the full-CI energy for the m s = 0 and m s = 1 wave functions (top), and error on the singlet-triplet gap∆E = E( 3 Σ + u ) -E( 1 Σ + g ) (bottom). On both graphics, the two curves of the dressed CAS-SD coincide.

TABLE I .

 I Comparison of the single-reference energies obtained with Gaussian 09, GAMESS, and the Quantum Package for FH at 1.2 Å, cc-pVDZ. All energies are converged below 10 10 a.u. 526 065 100.170 216 086 100.178 425 629 Quantum package 99.959 526 065 100.170 216 097 100.178 426 538

		Hartree-Fock	CISD	CCSD
	Gaussian 09	99.959 526 039 100.170 216 059 100.178 425 609
	GAMESS	99.959

TABLE III .

 III Second-order perturbative correction in the CIPSI calculations. Minimum and maximum values among all the points of the potential energy curve.largest calculation dropped the PT2 value to 6.7 mE h and the total energies differed by 1.4 mE h . So we estimate that the CIPSI curve of ethane has an accuracy less than 2 mE h in total energy, and all the other curves have an accuracy below the mE h . The non-parallelism errors (NPEs) of all the CIPSI curves are estimated with an error below the mE h .

			E PT2 (E h )
			Smallest	Largest
	C 2 H 6	13.4 × 10 3	17.8 × 10 3
	C 2 H 4 twisted	3.37 × 10 3	9.86 × 10 3
	C 2 H 4	2.41 × 10 3	6.85 × 10 3
	F 2 F 2	3 Σ + u	0.31 × 10 3 0.13 × 10 3	1.42 × 10 3 0.47 × 10 3
	N 2		61.1 × 10 6	0.41 × 10 3
	BeH 2	12.3 × 10 6	35.4 × 10 6
	H 2 O	1.59 × 10 6	69.1 × 10 6
	FH	0.23 × 10 6	55.1 × 10 6
	LiF	0.17 × 10 6	12.3 × 10 6

TABLE II .

 II Evaluation of the size-consistency error in the dressed CAS-SD and the µ-MR-CCSD, 6-31G basis set.

						FIG. 1. Dissociation curves. Differ-
						ence with respect to the full-CI energy
						using the MR-CCSD method presented
						in Ref. 17 and with the MR-CCSD
						method proposed in this work, as well as
						the CAS-SD and the dressed CAS-SD.
				Dressed	µ-MR-CCSD	µ-MR-CCSD
		CAS-SCF	CAS-SD	CAS-SD	µ ∈ Eq. (33)	unconstrained µ
	•CH 3	39.528 586 5	39.622 437 9	39.625 570 2	39.625 570 2	39.625 570 2
	•CH 3 ×2	79.057 173 0	79.244 875 8	79.251 140 3	79.251 140 3	79.251 140 3
	H 3 C• • • CH 3	79.057 173 0	79.237 098 2	79.250 695 7	79.251 039 3	79.251 107 3
	Error		7.78 × 10 3	4.44 × 10 4	1.01 × 10 4	3.29 × 10 5

TABLE IV .

 IV Non-parallelism errors (NPEs) and maximum errors with respect to the full-CI potential energy surface (mE h ).

				CAS-SD	λ-MR-CCSD	Dressed CAS-SD	µ-MR-CCSD
			NPE	Max error	NPE	Max error	NPE	Max error	NPE	Max error
	C 2 H 6	5.1	35.5	3.6	8.4	1.3	8.3	3.5	8.3
	F 2		3.8	19.8	3.1	4.0	0.9	4.2	1.5	3.8
	C 2 H 4 twist	1.5	27.7	0.3	6.7	0.7	7.1	0.5	6.5
	BeH 2	2.9	4.1	1.3	1.8	2.0	2.4	1.7	2.1
	H 2 O	1.9	4.8	0.7	1.2	0.2	1.3	0.5	1.3
	C 2 H 4 stretch	2.7	20.0	1.6	5.2	1.7	6.2	1.8	6.0
	N 2		1.7	9.0	0.9	2.2	1.7	3.8	0.3	2.3
	F 2	3 Σ + u (m s = 1)	2.6	18.6	1.3	3.3	1.2	3.5	1.2	3.3
	F 2	3 Σ + u (m s = 0)	2.6	18.6	1.2	1.8	1.3	3.5	1.1	3.3
	FH (ground state)	2.6	14.6	1.8	4.0	2.1	4.5	1.8	4.0
	FH (excited state)	3.3	20.9	8.8	8.5	10.5	10.1	7.1	8.3
	F 2 (local)	3.8	19.8	1.2	3.2	1.5	3.5	0.9	3.0
	N 2 (local)	1.7	9.0	3.6	5.0	1.1	3.5	0.4	1.8

ACKNOWLEDGMENTS

This work has been made through generous computational support from CALMIP (Toulouse) under the Allocation No. 2015-0510, and GENCI under the Allocation No. x2015081738. We are also very grateful to the two anonymous reviewers for their comments, which helped us improve the manuscript significantly.


Anthony Scemama

Yann Garniron

Michel Caffarel

Pierre-François Loos

Deterministic construction of nodal surfaces within quantum Monte Carlo: the case of FeS

Keywords: quantum Monte Carlo, diffusion Monte Carlo, full configuration interaction, multireference trial wave function

In diffusion Monte Carlo (DMC) methods, the nodes (or zeroes) of the trial wave function dictate the magnitude of the fixed-node (FN) error. Within standard DMC implementations, they emanate from short multideterminant expansions, stochastically optimized in the presence of a Jastrow factor. Here, following a recent proposal, we follow an alternative route by considering the nodes of selected configuration interaction (sCI) expansions built with the CIPSI (Configuration Interaction using a Perturbative Selection made Iteratively) algorithm. In contrast to standard implementations, these nodes can be systematically and deterministically improved by increasing the size of the sCI expansion. The present methodology is used to investigate the properties of the transition metal sulfide molecule FeS. This apparently simple molecule has been shown to be particularly challenging for electronic structure theory methods due to the proximity of two low-energy quintet electronic states of different spatial symmetry. In particular, we show that, at the triple-zeta basis set level, all sCI results -including those extrapolated at the full CI (FCI) limit -disagree with experiment, yielding an electronic ground state of 5 Σ + symmetry. Performing FN-DMC simulation with sCI nodes, we show that the correct 5 ∆ ground state is obtained if sufficiently large expansions are used. Moreover, we show that one can systematically get accurate potential energy surfaces and reproduce the experimental dissociation energy as well as other spectroscopic constants.

I. INTRODUCTION

From an experimental point of view, transition metal sulfides have proven to be useful in a variety of fields including biological chemistry, 1 catalysis, 2 and electrochemistry. 3 From the computational side, the apparently simple FeS diatomic molecule has been giving nightmares to computational chemists. The challenging features of the electronic structure of FeS originate from the energetic proximity of two electronic states with same multiplicity that compete for being the ground state. To make things worse, the equilibrium bond lengths associated with these two states are extremely close to each other.

Experimentally, the ground state of FeS is assigned to be 5 ∆, 4,5 with an equilibrium bond length r e = 2.017 Å, 5 and a dissociation energy D 0 = 3.31 (15) eV. 6 For this state, the harmonic frequency ω e has been estimated to 518 ± 5 cm -1 . 7 Very recently, a much more accurate value of the dissociation energy D 0 = 3.240(3) eV has been obtained by Matthew et al. using predissociation threshold technique. 8 FeS has been extensively studied by density-functional theory (DFT) and post-Hartree-Fock methods. In short, a) Corresponding author: scemama@irsamc.ups-tlse.fr b) Corresponding author: loos@irsamc.ups-tlse.fr most (but not all) DFT functionals correctly predict a 5 ∆ ground state, [9][10][11][12][13] while CAS-based multireference methods such as CASSCF/ACPF, 14 CASPT2, 15 or CASSCF/ICACPF 16 systematically predict5 Σ + lower than 5 ∆.

Here, we investigate this problem using quantum Monte Carlo (QMC). In recent years, QMC has been applied with great success to a large variety of main group compounds (see e.g. 17-20 for recent applications). Transition metal systems are more challenging but a number of successful studies have also been reported. When multireference effects are weak, QMC is seen as a very accurate method providing benchmark results of a quality similar or superior to the gold-standard CCSD(T). However, when multireference effects are dominantas it is usually the case for metallic compounds with partially-filled d shells -the situation is more complicated. Indeed, the results may depend significantly on the trial wave function Ψ T used to guide the walkers through configuration space. In theory, QMC results should be independent of the choice of Ψ T . However, it is not true in practice because of the fixed-node approximation which imposes the Schr ödinger equation to be solved with the additional constraint that the solution vanishes at the zeroes (nodes) of the trial wave function. Using an approximate Ψ T leads to approximate nodes and, thus, to an approximate energy, known as the FN energy. The FN energy being an upper bound of the exact energy, this gives us a convenient variational principle for characterizing the nodal quality: "the lower the FN energy, the better the nodes". In situations where multireference effects are strong, getting accurate nodes may be difficult. As we shall see, this is the main challenge we are facing in the present work.

Most QMC studies for transition metal-containing systems have been performed with pseudopotentials. In this case, an additional source of error, the so-called localization error, is introduced. This error, specific to QMC, adds up to the standard error associated with the approximate nature of pseudopotentials. Similarly to the FN error, the localization error depends on Ψ T and vanishes only for the exact wave function. Therefore, to get accurate and reliable QMC results, both sources of error have to be understood and controlled.

In 2011, Petz and L üchow reported a FN diffusion Monte Carlo (FN-DMC) study of the energetics of diatomic transition metal sulfides from ScS to FeS using pseudopotentials and single-determinant trial wave functions. 34 The pseudopotential dependence was carefully investigated, and comparisons with both DFT and CCSD(T) as well as experimental data were performed. In short, it was found that FN-DMC shows a higher overall accuracy than both B3LYP and CCSD(T) for all diatomics except for CrS and FeS that appeared to be particularly challenging.

Very recently, Haghighi-Mood and L üchow had a second look at the difficult case of FeS. 41 In particular, they explored the impact of the level of optimization on the parameters of multideterminant trial wave functions (partial or full optimization of the Jastrow, determinant coefficients and molecular orbitals) on both the FN and localization errors. Their main conclusions can be summarized as follows. Using a single-determinant trial wave function made of B3LYP orbitals or fully-optimized orbitals in the presence of a Jastrow factor is sufficient to yield the correct state ordering. However, in both cases, the dissociation energy is far from the experimental value and thus multideterminant trial wave functions must be employed. Although a natural choice would be to take into account the missing static correlation via a CASSCF-based trial wave function, they showed that it is insufficient and that a full optimization is essential to get both the correct electronic ground state and reasonable estimates of the spectroscopic constants.

In the present study, we revisit this problem within the original QMC protocol developed in our group these last few years. 37,[42][43][44][45][START_REF] Caffarel | Using cipsi nodes in diffusion monte carlo[END_REF][START_REF] Caffarel | [END_REF] In the conventional protocol, prevailing in the QMC community and employed by Haghighi-Mood and L üchow, the nodes of the Slater-Jastrow (SJ) trial wave function

Ψ SJ T = exp(J) Ψ det (1)
are obtained by partially-or fully-optimizing the Jastrow factor J and the multiderminant expansion Ψ det (containing typically a few hundreds or thousands of determinants). This step is performed in a preliminary variational Monte Carlo calculation by minimizing the energy, the variance of the local energy (or a combination of both) employing one of the optimization methods developed within the QMC context. [START_REF] Umrigar | [END_REF][49][50][51] We note that, in practice, the optimization must be carefully monitored because of the large number of parameters (several hundreds or thousands), the nonlinear nature of most parameters (several minima may appear) and the inherent presence of noise in the function to be minimized. Within our protocol, we rely on configuration interaction (CI) expansions in order to get accurate nodal surfaces, without resorting to the stochastic optimization step. Our fundamental motivation is to take advantage of all the machinery and experience developed these last decades in the field of wave function methods. In contrast to the standard protocol described above, the CI nodes can be improved deterministically and systematically by increasing the size of the CI expansion. In the present work, we do not introduce any Jastrow factor, essentially to avoid the expensive numerical quadrature involved in the calculation of the pseudopotential, and to facilitate the control of the localization error. To keep the size of the CI expansion reasonable and retain only the most important determinants, we propose to use selected CI (sCI) algorithms, such as CIPSI (Configuration Interaction using a Perturbative Selection made Iteratively). 42 Using a recently-proposed algorithm to handle large numbers of determinants in FN-DMC 44 we are able to consider up to a few million determinants in our simulations.

Over the last few years, we have witnessed a re-birth of sCI methods. 36,37,[42][43][44][45] Although these various approaches appear under diverse acronyms, most of them rely on the very same idea of selecting determinants iteratively according to their contribution to the wave function or energy, an idea that goes back to 1969 in the pioneering works of Bender and Davidson, 52 and Whitten and Hackmeyer. 53 Importantly, we note that any sCI variants can be employed here.

The price to pay for using sCI expansions instead of optimized SJ trial wave functions is the need of much larger multideterminant expansions as well as the presence of larger statistical and systematic errors (such as time-step and basis set incompleteness errors). However, these disadvantages are compensated by the appealing features of sCI nodes: i) they are built in a fully-automated way; ii) they are unique and reproducible; iii) they can be systematically improved by increasing the level of selection and/or the basis set (with the possibility of complete basis set extrapolation 45 ), and iv) they easily produce smooth potential energy surfaces. 43 Unless otherwise stated, atomic units are used throughout.

II. COMPUTATIONAL DETAILS

All trial wave functions have been generated with the electronic structure software QUANTUM PACKAGE, [START_REF] Scemama | Quantum package v1.0[END_REF] while the QMC calculations have been performed with the QMC=CHEM suite of programs. [START_REF] Scemama | Qmc=chem[END_REF][START_REF] Scemama | [END_REF] Both softwares are developed in our laboratory and are freely available. For all calculations, we have used the triple-zeta basis sets of Burkatzki et al. 79,80 (VTZ-ANO-BFD for Fe and VTZ-BFD for S) in conjunction with the corresponding Burkatzki-Filippi-Dolg (BFD) small-core pseudopotentials including scalar relativistic effects. For more details about our implementation of pseudopotentials within QMC, we refer the interested readers to Ref. 46. As pointed out by Hammond and coworkers, 81 thanks to the absence of Jastrow factor in our trial wave functions, the non-local pseudopotential can be localized analytically and the usual numerical quadrature over the angular part of the non-local pseudopotential can be eschewed.

In practice, the calculation of the localized part of the pseudopotential represents only a small overhead (about 15%) with respect to a calculation without pseudopotential (and the same number of electrons).

In order to compare our results for the dissociation energy of FeS with the experimental value of Matthew et al. 8 and the (theoretical) benchmark value of Haghighi-Mood and L üchow, 41 we have taken into account the zero-point energy (ZPE) correction, the spin-orbit effects as well as the core-valence correlation contribution the same way as Ref. 41. For the 5 ∆ state, this corresponds to an increase of the dissociation energy by 0.06 eV, and a 0.02 eV stabilization of 5 ∆ compared to 5 Σ + .

A. Jastrow-free trial wave functions

Within the spin-free formalism used in QMC, a CIbased trial wave function is written as

Ψ T (R) = N det ∑ I=1 c I D I (R) = N det ∑ I=1 c I D ↑ I (R ↑ )D ↓ I (R ↓ ), (2)
where R = (r 1 , ..., r N ) denotes the full set of electronic spatial coordinates, R ↑ and R ↓ are the two subsets of spin-up (↑) and spin-down (↓) electronic coordinates, and

D σ I (R σ ) (σ = ↑ or ↓) are spin-specific determinants.
In practice, the various products D ↑ I D ↓ I contain many identical spin-specific determinants. For computational efficiency, it is then advantageous to group them and compute only once their contribution to the wave function and its derivatives. 44 Therefore, the Jastrow-free CI trial wave functions employed in the present study are rewritten in a "spin-resolved" form

Ψ T (R) = N ↑ det ∑ i=1 N ↓ det ∑ j=1 c ij D ↑ i (R ↑ )D ↓ j (R ↓ ), (3)
where D σ i i=1,...,N σ det denotes the set of all distinct spinspecific determinant appearing in Eq. ( 2).

B. Quantum Monte Carlo calculations

To avoid handling too many determinants in Ψ T , a truncation scheme has to be introduced. In most CI and/or QMC calculations, the expansion is truncated by either introducing a cutoff on the CI coefficients or on the norm of the wave function. Here, we use an alternative truncation scheme knowing that most of the computational effort lies in the calculation of the spin-specific determinants and their derivatives. Removing a product of determinants whose spin-specific determinants are already present in other products does not change significantly the computational cost. Accordingly, a natural choice is then to truncate the wave function by removing independently spin-up and spin-down determinants. To do so, we decompose the norm of the wave function as

N = N ↑ det ∑ i=1 N ↓ det ∑ j=1 c ij 2 = N ↑ det ∑ i=1 N ↑ i = N ↓ det ∑ j=1 N ↓ j . (4)
A determinant

D ↑ i is retained in Ψ T if N ↑ i = N ↓ det ∑ j=1 c ij 2 > , (5)
where is a user-defined threshold. A similar formula is used for D ↓ j . When = 0, the entire set of determinants is retained in the QMC simulation.

In order to treat the two electronic states ( 5 Σ + and 5 ∆) on equal footing, a common set of spin-specific determinants D σ i i=1,...,N σ det is used for both states. In addition, a common set of molecular orbitals issued from a preliminary state-averaged CASSCF calculation is employed. These CASSCF calculations have been performed with the GAMESS package 82 while, for the atoms, we have performed ROHF calculations. The active space contains 12 electrons and 9 orbitals (3d and 4s orbitals of Fe and 3p orbitals of S). The multideterminant expansion (2) has been constructed using the sCI algorithm CIPSI, 54,55 which uses a second-order perturbative criterion to select the energetically-important determinants D I in the FCI space. 36,37,[42][43][44][45]65 A n s -state truncated sCI expansion (here n s = 2) is obtained via a natural generalization of the state-specific criterion introduced in Eq. ( 5): a determinant

D ↑ i is retained in Ψ T if N ↑ i = 1 n s n s ∑ k=1 N ↓ det ∑ j=1 c (k) ij 2 > , (6)
with a similar formula for D ↓ j . The characteristics of the various trial wave functions considered here (and their acronyms) at r FeS = 2.0 Å are presented in Table I. For other r FeS values, the numbers of determinants are slightly different. Our largest sCI trial wave function contains 8 388 608 determinants and is labeled sCI(∞). The sCI(n) wave functions with n = 4, 5, and 6 are obtained by truncation of the sCI(∞) expansion setting = 10 -n . They contain respectively 15 723, 269 393, and 1 127 071 determinants. At this stage, we are not able to use the entire 8 388 608 determinants of the sCI(∞) wave function within our FN-DMC simulations. In comparison, Haghighi-Mood and L üchow's CASSCFbased trial wave function (labelled as HML in Table II) only contains 630 and 500 determinants for the 5 Σ + and 5 ∆ states, respectively. 41 However, as discussed in the introduction, fully-optimized SJ trial wave functions require much smaller multireference expansions.

Based on these trial wave functions, we performed FN-DMC calculations with the stochastic reconfiguration algorithm developed by Assaraf et al. 83 In order to remove the time-step error, all our FN-DMC results have been extrapolated to zero time-step using a two-point linear extrapolation with τ = 2 × 10 -4 and 4 × 10 -4 . 84 Note that, because the variance of the local energy is larger than in SJ calculations, time-step errors are enhanced and shorter time steps are required.

III. RESULTS AND DISCUSSION

In Table II, we report FN-DMC energies at equilibrium geometry as well as other quantities of interest such as the dissociation energy D 0 , the equilibrium distance r e and harmonic frequency ω e obtained with various trial wave functions. These values are obtained via the standard four-parameter Morse potential representation of the numerical values gathered in Tables III andIV. 85 For comparison purposes, Haghighi-Mood and L üchow's results are also reported based on their best trial wave function. 41 When available, the experimental result is also reported. [5][6][7][8] The value of D 0 is always calculated with respect to the 5 ∆ state adding the corresponding corrections for ZPE, spin-orbit effects and core-valence correlation, as described above (see Sec. II A). The dissociation profile of FeS obtained with FN-DMC is depicted in Fig. 1 for various trial wave functions.

For the variational results gathered in Tables III and IV, the FCI limit has been reached by the method recently proposed by Holmes, Umrigar and Sharma 74 in the context of the (selected) heat-bath CI method. 70,73,74 In order to obtain FCI results, they proposed to linearly extrapolate the sCI energy E sCI as a function of the second-order Epstein-Nesbet energy

E PT2 = ∑ α α| Ĥ|Ψ sCI 2 E sCI -α| Ĥ|α , ( 7
)
which is an estimate of the truncation error in the sCI algorithm, i.e E PT2 ≈ E FCI -E sCI . 54 In Eq. ( 7), the sum runs over all external determinants |α (i.e. not belonging to the sCI expansion) connected via Ĥ to the sCI wave function Ψ sCI , i.e. α| Ĥ|Ψ sCI = 0. When E PT2 = 0, the FCI limit has effectively been reached. In our case, E PT2 is efficiently evaluated thanks to our recently-proposed hybrid stochastic-deterministic algorithm, 65 which explains the error bar on E PT2 in Tables III andIV. The extrapolated FCI results are labeled exFCI from hereon.

To obtain the FN-DMC curve with an effective FCI trial wave function, we have generalized the extrapolation procedure described above, and we have performed a three-point linear extrapolation of the FN-DMC energy as a function of E exFCI -E sCI using the sCI(4), sCI(5) and sCI( 6) results (see Fig. 2).

The first observation we would like to make is that, at the variational level, the 5 ∆ state is never found lower in energy than the 5 Σ + state, even after performing the extrapolation to the FCI limit (see Table IV). Because all post-Hartree-Fock methods are indeed an approximation of FCI, they are expected to predict a 5 Σ + ground state for this particular basis set. This observation is in agreement with the CASPT2 results previously published in the literature. [14][15][16] Thus, one can attribute the wrong state ordering to basis set incompleteness, the only remaining approximation.

At the FN-DMC level, one must include at least a few hundred thousand determinants in order to find the proper ground state. For larger values (10 -4 and 10 -5 ), D 0 is underestimated due to the unbalanced treatment of the isolated atoms compared to the dimer at equilibrium geometry. Indeed, for a given number of determinants, the energy of the atomic species is much closer to the FCI limit than the energy of FeS.

For = 10 -6 , our approach correctly predicts a 5 ∆ ground state. However, although our FN-DMC energies are much lower than those obtained with the HML trial wave function, our estimate of the dissociation energy (D 0 = 3.062(39) eV) is still below the experimental value. This underestimation of D 0 can be ultimately tracked down to the lack of size-consistency of the truncated CI wave function. With more than 10 6 determinants in the variational space, the wave function is still 150 mE h higher than the exFCI wave function, while the atoms are much better described by the sCI wave function. To remove the size-consistency error, we then extrapolate the FN-DMC energies to the (size-consistent) FCI limit of the trial wave function, as described above.

In that case, using the extrapolated FN-DMC energies of the molecule and isolated atoms reported in Table II, we obtain a value of D 0 = 3.267 (49) eV, which nestles nicely between the experimental values of Matthew et al. 8 (3.240(3) eV) and Drowart et al. 6 (3.31(15) eV). As a final remark, we note that other spectroscopic constants, such as the equilibrium geometry and the harmonic frequency, are fairly well reproduced by our approach.

IV. CONCLUSION

In this article, the potential energy curves of two electronic states -5 ∆ and 5 Σ + -of the FeS molecule have been calculated using the sCI algorithm CIPSI and the stochastic FN-DMC method. In all our sCI calculations, 5 Σ + is found to be the ground state, in disagreement with experiment. It is not only true for all CIPSI expansions with up to 8 million determinants but also for the estimated FCI limit obtained using the extrapolation procedure recently proposed by Holmes et al. 74 This conclusion agrees with other high-level ab initio wave function calculations which all wrongly predict a ground state of 5 Σ + symmetry. FN-DMC calculations have been performed using CIPSI expansions including up to 1 127 071 selected determinants as trial wave functions. Contrary to standard QMC calculations, we do not introduce any Jastrow factor: the CI expansions have been used as they are (no optimization). It is found that, when the number of determinants in the trial wave function reaches few hundred thousands, the FN-DMC ground state switches from the 5 Σ + state to the correct 5 ∆ state, as predicted experimentally.

Generalizing the extrapolation procedure of Holmes et al., 74 an estimate of the FN-DMC potential energy curves corresponding to the FCI nodes can be obtained. The resulting dissociation energy is found to be 3.267 (49) eV, in agreement with the recent experimental value of Matthew et al. (3.240(3) eV). 8 As already observed in previous applications, the FN-DMC energy obtained with CIPSI nodes is found to systematically decrease as a function of the number of selected determinants. 36,37,42,43,45,[START_REF] Caffarel | Using cipsi nodes in diffusion monte carlo[END_REF] For the largest expansion, our fixed-node energies are lower than the values recently reported by Haghighi-Mood and L üchow 41 using a fully-optimized SJ trial wave function. This important result illustrates that "pure" sCI nodes is a realistic alternative to stochasticallyoptimized SJ trial wave functions, even for a challenging system such as FeS. A similar conclusion had already been drawn in our recent study of the water molecule. 45

FIG. 1 .FIG. 2 .

 12 FIG.1. E DMC (in Hartrees) for the5 Σ + (solid) and 5 ∆ (dashed) states of FeS as a function of r FeS (in Å) for various trial wave functions. The error bar corresponds to one standard error.

TABLE I .

 I Characteristics of the various sCI expansions at r FeS = 2.0 Å for various levels of truncation. The characteristics of the extrapolated FCI (exFCI) expansion are also reported.

	Method sCI	10 -4	N det 15 723	N det ↑ 191	N det ↓ 188	acronym sCI(4)
		10 -5	269 393	986	1 191	sCI(5)
		10 -6	1 127 071	3 883	4 623	sCI(6)
		0	8 388 608 364 365 308 072	sCI(∞)
	exFCI	-	∼ 10 27	∼ 10 16	∼ 10 11	exFCI

TABLE II .

 II FN-DMC energies E DMC (in Hartrees) at equilibrium geometry, dissociation energy D 0 (in eV), equilibrium distance r e (in Å), harmonic frequency ω e (in cm -1 ) for the 5 Σ + and 5 ∆ of FeS obtained with various trial wave functions Ψ T . The error bar corresponding to one standard error is reported in parenthesis.

	Ψ T		FeS ( 5 Σ + )			FeS ( 5 ∆)		Fe ( 5 D)	S ( 3 P)	D 0	Ref.
		E DMC	r e	ω e	E DMC	r e	ω e	E DMC	E DMC
	HML -134.057 1(4) 2.00(1) 518(7) -134.057 9(4) 2.031(7) 499(11) -123.812 6(4) -10.131 4(1) 3.159(15)	41
	sCI(4) -134.010 1(8) 1.994(7) 532(20) -134.004 0(7) 2.029(7) 502(15) -123.802 8(9) -10.127 9(2) 2.055(20) This work
	sCI(5) -134.047 9(10) 1.992(8) 551(24) -134.040 2(10) 2.048(11) 489(21) -123.823 4(10) -10.131 2(2) 2.389(28) This work
	sCI(6) -134.061 7(14) 1.994(12) 497(35) -134.067 1(14) 2.004(11) 550(32) -123.830 0(12) -10.133 4(3) 3.062(39) This work
	exFCI -134.086 3(15) 1.990(12) 523(37) -134.088 5(18) 2.016(14) 525(40) -123.837 2(12) -10.133 6(3) 3.267(49) This work
	Exp.	-	-	-	-	2.017	518(5)	-	-	3.240(3) 5, 7, and 8

TABLE III .

 III sCI energy E sCI , second-order perturbation correction E PT2 and FN-DMC energy E DMC (in Hartrees) for the Fe ( 5 D state) and S ( 3 P state) atoms obtained with various methods. The error bar corresponding to one standard error is reported in parenthesis.

	Method		Fe ( 5 D)			S ( 3 P)
		E sCI	E PT2	E DMC	E sCI	E PT2	E DMC
	sCI(4)	-123.418 124	-0.316 44(3)	-123.802 8(9)	-10.093 850	-0.023 348(2)	-10.127 9(2)
	sCI(5)	-123.607 608	-0.127 40(1)	-123.823 4(10)	-10.108 576	-0.007 941(1)	-10.131 2(2)
	sCI(6)	-123.673 435	-0.063 698(6)	-123.830 0(12)	-10.113 926	-0.002 179(0)	-10.133 4(3)
	sCI(∞)	-123.720 629	-0.016 987(2)	-	-10.115 844	-0.000 188(0)	-
	exFCI	-123.738 264	0	-123.837 2(12)	-10.115 996	0	-10.133 4(3)

TABLE IV .

 IV sCI energy E sCI , second-order perturbation correction E PT2 and FN-DMC energy E DMC (in Hartrees) for the 5 Σ + and 5 ∆ of FeS obtained with various methods. The error bar corresponding to one standard error is reported in parenthesis.

	Method	r FeS		FeS ( 5 Σ + )			FeS ( 5 ∆)
			E sCI	E PT2	E DMC	E sCI	E PT2	E DMC
	sCI(4)	1.8	-133.322 421	-0.602 63(6)	-133.989 0(12)	-133.287 542	-0.633 28(6)	-133.978 1(12)
		1.9	-133.338 106	-0.603 60(6)	-134.007 4(12)	-133.304 261	-0.635 26(6)	-133.998 2(12)
		2.0	-133.344 260	-0.600 73(6)	-134.009 2(12)	-133.316 924	-0.624 38(6)	-134.001 8(12)
		2.1	-133.343 984	-0.597 12(6)	-134.005 5(12)	-133.318 747	-0.617 18(6)	-134.002 9(12)
		2.2	-133.337 856	-0.595 60(6)	-133.998 5(12)	-133.321 678	-0.603 51(6)	-133.997 4(12)
		2.3	-133.312 474	-0.606 32(6)	-133.987 5(12)	-133.334 006	-0.587 80(6)	-133.985 3(12)
		2.4	-133.309 214	-0.601 37(6)	-133.977 0(12)	-133.323 583	-0.587 32(6)	-133.975 2(12)
		2.5	-133.295 921	-0.607 48(6)	-133.963 3(13)	-133.305 592	-0.595 18(6)	-133.963 5(12)
	sCI(5)	1.8	-133.607 003	-0.290 22(3)	-134.026 4(15)	-133.574 416	-0.309 14(3)	-134.014 1(16)
		1.9	-133.623 662	-0.290 03(3)	-134.043 1(15)	-133.595 588	-0.305 97(3)	-134.034 6(16)
		2.0	-133.631 686	-0.286 36(3)	-134.047 9(16)	-133.605 430	-0.302 44(3)	-134.038 9(17)
		2.1	-133.631 076	-0.283 06(3)	-134.045 2(16)	-133.608 180	-0.297 83(3)	-134.038 7(17)
		2.2	-133.622 872	-0.283 00(3)	-134.034 0(15)	-133.600 172	-0.299 37(3)	-134.029 5(13)
		2.3	-133.600 838	-0.289 70(3)	-134.023 1(17)	-133.618 393	-0.276 15(3)	-134.027 1(15)
		2.4	-133.591 037	-0.289 75(3)	-134.012 1(15)	-133.602 716	-0.279 50(3)	-134.010 3(16)
		2.5	-133.583 356	-0.286 68(3)	-133.998 9(15)	-133.590 900	-0.278 56(3)	-134.992 0(15)
	sCI(6)	1.8	-133.742 093	-0.159 64(2)	-134.044 8(21)	-133.725 217	-0.162 79(2)	-134.043 7(21)
		1.9	-133.759 798	-0.157 96(2)	-134.057 4(21)	-133.746 303	-0.160 67(2)	-134.061 7(22)
		2.0	-133.766 027	-0.155 98(2)	-134.062 6(22)	-133.756 008	-0.157 98(2)	-134.068 5(22)
		2.1	-133.764 192	-0.154 31(2)	-134.057 1(23)	-133.758 610	-0.154 21(2)	-134.062 5(22)
		2.2	-133.759 178	-0.151 39(2)	-134.052 9(22)	-133.755 183	-0.151 49(2)	-134.053 7(22)
		2.3	-133.749 068	-0.148 46(1)	-134.040 2(21)	-133.750 644	-0.148 52(1)	-134.043 5(22)
		2.4	-133.741 085	-0.146 21(1)	-134.028 6(23)	-133.740 633	-0.146 76(1)	-134.030 0(22)
		2.5	-133.731 347	-0.145 00(1)	-134.018 9(23)	-133.729 703	-0.145 52(1)	-134.013 3(22)
	sCI(∞)	1.8	-133.850 093	-0.071 305(7)	-	-133.836 804	-0.073 473(7)	-
		1.9	-133.868 551	-0.069 523(7)	-	-133.857 909	-0.071 490(7)	-
		2.0	-133.874 845	-0.067 992(7)	-	-133.866 607	-0.069 933(7)	-
		2.1	-133.873 424	-0.066 144(7)	-	-133.867 130	-0.068 250(7)	-
		2.2	-133.868 534	-0.063 444(6)	-	-133.864 246	-0.065 508(7)	-
		2.3	-133.856 234	-0.063 467(6)	-	-133.858 730	-0.061 578(6)	-
		2.4	-133.849 753	-0.058 745(6)	-	-133.849 494	-0.059 879(6)	-
		2.5	-133.839 274	-0.056 324(6)	-	-133.840 801	-0.057 100(6)	-
	exFCI	1.8	-133.924 967	0	-134.068 7(24)	-133.913 076	0	-134.065 9(25)
		1.9	-133.941 175	0	-134.080 6(24)	-133.932 524	0	-134.083 6(25)
		2.0	-133.945 859	0	-134.087 5(24)	-133.939 854	0	-134.090 8(26)
		2.1	-133.942 722	0	-134.083 3(25)	-133.938 821	0	-134.084 0(25)
		2.2	-133.934 789	0	-134.074 4(24)	-133.932 709	0	-134.072 3(23)
		2.3	-133.922 799	0	-134.061 8(24)	-133.923 119	0	-134.069 5(24)
		2.4	-133.910 763	0	-134.050 2(24)	-133.911 892	0	-134.051 0(25)
		2.5	-133.897 709	0	-134.039 1(24)	-133.900 371	0	-134.028 7(24)

∆: σ 2 π 4 σ 2 δ 3 σ 1 π 2 , 5 Σ + : σ 2 π 4 σ 2 δ 2 σ 2 π 2 ,

ACKNOWLEDGMENTS

The authors would like to thank Arne L üchow for numerous stimulating discussions. This work was performed using HPC resources from CALMIP (Toulouse) under allocation 2016-0510 and from GENCI-TGCC (Grant 2016-08s015).