Alexandrina Stoyanova

Andrew M Teale

Julien Toulouse

Trygve Helgaker

Emmanuel Fromager

Alternative separation of exchange and correlation energies in multi-configuration range-separated density-functional theory

The alternative separation of exchange and correlation energies proposed by Toulouse et al. [Theor. Chem. Acc. 114, 305 (2005)] is explored in the context of multiconfiguration range-separated density-functional theory. The new decomposition of the short-range exchange-correlation energy relies on the auxiliary long-range interacting wavefunction rather than the Kohn-Sham (KS) determinant. The advantage, relative to the traditional KS decomposition, is that the wavefunction part of the energy is now computed with the regular (fully-interacting) Hamiltonian. One potential drawback is that, because of double counting, the wavefunction used to compute the energy cannot be obtained by minimizing the energy expression with respect to the wavefunction parameters. The problem is overcome by using short-range optimized effective potentials (OEPs). The resulting combination of OEP techniques with wavefunction theory has been investigated in this work, at the Hartree-Fock (HF) and multi-configuration self-consistent-field (MCSCF) levels. In the HF case, an analytical expression for the energy gradient has been derived and implemented. Calculations have been performed within the short-range local density approximation on H 2 , N 2 , Li 2 and H 2 O. Significant improvements in binding energies are obtained with the new decomposition of the short-range energy. The importance of optimizing the short-range OEP at the MCSCF level when static correlation becomes significant has also been demonstrated for H 2 , using a finite-difference gradient. The implementation of the analytical gradient for MCSCF wavefunctions is currently in

I. INTRODUCTION

The simultaneous description of dynamical and non-dynamical (static or strong) electron correlation in atomic and molecular systems, using low-cost methodologies, remains a challenge for electronic-structure theory. In particular, standard Kohn-Sham density-functional theory 1 (KS-DFT) approximations have enjoyed success in treating phenomena where the description of short-range dynamical correlation is paramount but they have been unable to provide reliable results whenever static correlation is important-including the description of bond breaking, transition-metal compounds, conjugated polymers, and magnetic materials.

Over the years, a great deal of effort has been devoted to understand the shortcomings of density-functional approximations (DFAs) for the treatment of static correlation and to develop new approximations to address this situation. Within the framework of KS-DFT, the pragmatic unrestricted approach can of course be used for describing bond breaking for example but, in this case, space and spin symmetries are broken which is fundamentally not satisfactory. Recent progress towards DFAs capable of treating static correlation has been made, for example, by Malet and Gori-Giorgi 2 and by Becke 3 . Other authors have focussed on ensemble-DFT (E-DFT). The utility of E-DFT has been analyzed by Schipper et al. 4 and recent implementations of E-DFT variants have been made by Filatov et al. 5 , by Chai 6 , and by Nygaard and Olsen 7 .

Beyond the framework of standard KS-DFT, a number of groups have pursued the idea of hybridizing density-functional and multi-configuration self-consistent-field (MCSCF) approaches. The goal of these approaches is to treat static correlation using the flexibility of the MCSCF expansion, whilst treating dynamical correlation using density functionals [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] . A number of MCSCF-DFT hybrid approaches have been proposed, including the complete-activespace-DFT (CAS-DFT) schemes of Gräffenstein et al. 12 , Gusarov et al. 18 , and Miehlich et al. 11 A multi-configuration extension of KS-DFT using optimized effective potential techniques has also been proposed by Weimer et al. 23 A challenge for these methodologies is to avoid double counting of correlation effects, as the density functionals utilized in these approaches depend on the MCSCF expansions. Tackling the double counting problem is a difficult task as illustrated by Kurzweil et al. 24 in their analysis on the mapping of interacting onto partially interacting system.

A key step in avoiding the doubling counting problems of the MCSCF-DFT hybrid approaches was the proposal of Savin et al. [START_REF] Savin | Recent Developments and Applications of Modern Density Functional Theory[END_REF][START_REF] Stoll | Density functional methods in physics[END_REF] to divide the Coulomb interaction into long-range (lr) and short-range (sr) contributions. The introduction of this range-separated approach has led to a wide variety of hybrid wave-function/DFT approaches; in the present context, we note the multi-configuration short-range DFT (MC-srDFT) approach of Refs. 27-29. Range separation has proven to be of great utility for the treatment of dispersion interactions, where the simple physical intuition that interactions are long-ranged can be leveraged to provide a clean division of labour between the density-functional and wavefunction contributions. [30][31][32][33][34][35][36] However, for the description of static correlation, this approach is less effective since static correlation may not be interpreted as predominantly long-ranged.

Thus, even for systems such as stretched H 2 , where static correlation is expected to be dominant, the corresponding short-ranged density-functional still plays a significant role 37 and errors associated with these approximations can lead to significant errors in binding energies, albeit with an improved shape for the potential energy curve. [START_REF] Fromager | [END_REF] These considerations are reflected in the fact that a recently proposed MC-DFT approach based on a simple linear (rather than range-dependent) decomposition of the Coulomb interaction delivers similar accuracy in practice. 38 To overcome some of the difficulties associated with the MC-srDFT models for the description of static correlation, Toulouse et al. 39 proposed an alternative separation of the short-range exchange and correlation energies, which may be more natural in the context of hybrid methodologies that incorporate multi-configurational components. The practical performance of methods using this partitioning is investigated in the present work.

The paper is organized as follows. In Sec. II, exact and approximate formulations of rangeseparated DFT are presented. The differences between the traditional and the alternative short-range exchange-correlation energy decompositions are highlighted. For the latter, long-range HF and MCSCF approximations are combined with short-range OEPs in order to overcome double counting problems. Details of the various range-separated schemes that have been implemented are given in Sec. III, while the results obtained for H 2 , N 2 , H 2 O, and Li 2 are discussed in Sec. IV. Finally, conclusions are given in Sec. V.

II. THEORY

In this section, we present the theory underlying the various multi-configuration rangeseparated DFT models assessed in Sec. IV. We first introduce in Sec. II A the exact multi-determinantal extension of KS-DFT based on range separation. The conventional KS decomposition of the short-range exchange-correlation density functional, as well as an alternative one that relies on the long-range interacting wavefunction (rather than the KS determinant), are next discussed in Secs. II B and II C. The combination of OEPs with multi-determinant wavefunctions and density-functionals is then considered in Sec. II D 1.

Models based on HF and MC wavefunctions are discussed in Sec. II D 2. The derivation of the analytical energy gradient for the OEP optimization at the HF level is finally presented in Sec. II E. A summary is given in Sec. II F.

A. Multi-determinant range-separated DFT

In multi-determinant range-separated DFT, [START_REF] Savin | Recent Developments and Applications of Modern Density Functional Theory[END_REF][START_REF] Stoll | Density functional methods in physics[END_REF][START_REF] Savin | Recent advances in density functional methods[END_REF] referred to as short-range DFT (srDFT) in the following, the regular two-electron Coulomb interaction is split into long-and shortrange parts, w ee (r 12 ) = 1/r 12 = w lr,µ ee (r 12 ) + w sr,µ ee (r 12 ),

where µ is a parameter that controls the range separation. In this work, the long-range interaction is based on the error function

w lr,µ ee (r 12 ) = erf(µr 12 ) r 12 . (2)
The universal Levy-Lieb functional 41,42

F [n] = min Ψ→n Ψ| T + Ŵee |Ψ , ( 3
)
where T is the kinetic energy operator and Ŵee the two-electron repulsion operator, can then be rewritten as

F [n] = F lr,µ [n] + E sr,µ Hxc [n], (4)
with the universal long-range functional defined as

F lr,µ [n] = min Ψ→n Ψ| T + Ŵ lr,µ ee |Ψ = Ψ µ [n]| T + Ŵ lr,µ ee |Ψ µ [n] . (5)
The minimizing wavefunction Ψ µ [n] in Eq. ( 5) corresponds to the ground state of the auxiliary long-range interacting system with density n. When connecting the auxiliary and physical systems by a generalized adiabatic-connection path, 37,43,44 the complementary short-range Hartree-exchange-correlation (srHxc) energy can be expressed as

E sr,µ Hxc [n] = 1 µ/(1+µ) W ν Hxc [n] dν, (6)
where the srHxc integrand is given by

W ν Hxc [n] = ∂ Ŵ lr,ν/(1-ν) ee ∂ν Ψ ν/(1-ν) [n] . (7)
Note that the expression in Eq. ( 6) relies on the density constraint

n Ψ ν/(1-ν) [n] (r) = n(r), 0 ≤ ν ≤ 1. ( 8
)
According to the variational principle 45 and Eq. ( 4), the exact expression for the ground-state energy of an electronic system becomes

E = min n F [n] + dr v ne (r) n(r) = min n Ψ µ [n]| T + Ŵ lr,µ ee + Vne |Ψ µ [n] + E sr,µ Hxc [n] , (9)
where Vne = dr v ne (r) n(r) is the nuclear potential operator. The exact ground-state energy can also be obtained by a minimization over wavefunctions

E = min Ψ Ψ| T + Ŵ lr,µ ee + Vne |Ψ + E sr,µ Hxc [n Ψ ] = Ψ µ | T + Ŵ lr,µ ee + Vne |Ψ µ + E sr,µ Hxc [n Ψ µ ], (10)
where the minimizing wavefunction Ψ µ fulfills the self-consistent equation

( T + Ŵ lr,µ ee + V µ [n Ψ µ ])|Ψ µ = E µ |Ψ µ , (11) with
V µ [n] = dr v ne (r) + δE sr,µ Hxc δn(r) [n] n(r). ( 12
)
While Ŵ lr,µ ee vanishes and Eq. ( 11) reduces to the conventional KS equation at µ = 0, the full Schrödinger equation is recovered in the µ → +∞ limit, as Ŵ lr,µ ee reduces to the regular two-electron repulsion and the short-range interaction vanishes. For intermediate µ values, 0 < µ < +∞, a hybrid wave-function/DFT description is obtained. Then, in contrast to traditional KS-DFT, the exact auxiliary wavefunction Ψ µ is generally multi-determinantal owing to the explicit description of long-range interactions.

The srDFT approximation obtained by restricting the minimization in Eq. ( 10) to single determinants is in the following referred to as HF-srDFT; this approximation was referred to as range-separated hybrid (RSH) theory in Ref. 30. To describe multi-configurational electronic systems, a long-range MCSCF description has also been proposed, [START_REF] Pedersen | [END_REF][START_REF] Fromager | [END_REF] leading to the MC-srDFT model. In this work, we consider both schemes.

B. KS decomposition of the short-range energy

The conventional decomposition of the srHxc energy is analogous to that in standard KS-DFT 46 ,

E sr,µ Hxc [n] = E sr,µ H [n] + E sr,µ x [n] + E sr,µ c [n], (13)
where the short-range Hartree energy is defined as

E sr,µ H [n] = 1 2 drdr ′ n(r)n(r ′ )w sr,µ ee (|r -r ′ |) , (14)
and the exact short-range exchange energy is calculated from the non-interacting KS determinant Φ KS [n] with density n. The Hartree-exchange integrand is then obtained from Eq. ( 7) by replacing the long-range interacting wavefunction by Φ KS [n]

W ν Hx [n] = ∂ Ŵ lr,ν/(1-ν) ee ∂ν Φ KS [n] , (15)
which defines the short-range exchange energy, according to Eq. ( 6), as

E sr,µ x [n] = 1 µ/(1+µ) W ν Hx [n] dν -E sr,µ H [n] (16) = Φ KS [n]| Ŵ sr,µ ee |Φ KS [n] -E sr,µ H [n].
The corresponding correlation integrand is then given by

W ν c [n] = W ν Hxc [n] -W ν Hx [n] = ∂ Ŵ lr,ν/(1-ν) ee ∂ν Ψ ν/(1-ν) [n] - ∂ Ŵ lr,ν/(1-ν) ee ∂ν Φ KS [n] , (17)
leading to the following expression for the exact complementary short-range correlation energy

E sr,µ c [n] = 1 µ/(1+µ) W ν c [n] dν. ( 18
)
Note that, according to Eqs. ( 4), ( 5) and ( 16), this energy can also be expressed as

E sr,µ c [n] = E c [n] + Φ KS [n]| T + Ŵ lr,µ ee |Φ KS [n] -Ψ µ [n]| T + Ŵ lr,µ ee |Ψ µ [n] , (19)
where E c [n] is the regular correlation density-functional energy, recovered in the µ = 0 limit. It is clear from Eq. ( 19) that the complementary short-range correlation density functional contains the purely short-range correlation effects as well as their coupling with long-range correlation. 47 Various short-range functionals have been developed for practical srDFT calculations at the local density (srLDA), [47][48][49] generalized gradient, 31,46,[50][51][52][53] and meta-generalized gradient 54 levels of approximation. These functionals have been successfully

employed with a number of post-HF-srDFT and post-MC-srDFT long-range correlation treatments to describe dispersion. [30][31][32][33][34][35] However, for systems with significant static correlation, they are usually not accurate enough. [START_REF] Fromager | [END_REF]29 To improve upon the description of the short-range energy, we consider in the following an alternative separation of exchange and correlation energies.

C. Alternative decomposition of the short-range energy

As pointed out by Toulouse, Gori-Giorgi, and Savin, 39,55 it is more natural, in the context of srDFT, to define the short-range exchange energy in terms of the multi-determinantal wavefunction Ψ µ [n] introduced in Eq. ( 5). This observation leads to the following decomposition of the srHxc density-functional energy

E sr,µ Hxc [n] = E sr,µ H [n] + E sr,µ x,md [n] + E sr,µ c,md [n], (20)
where, what is referred to as the short-range multideterminantal (MD) exchange functional in Ref. 55, is defined as

E sr,µ x,md [n] = Ψ µ [n]| Ŵ sr,µ ee |Ψ µ [n] -E sr,µ H [n]. (21)
This expression arises naturally from Eqs. ( 6) and ( 7) when replacing the ν-dependent wavefunction in the integrand by the one obtained with the lower integration limit ν = µ/(1 + µ), namely Ψ µ [n]. We thus define the MD Hartree-exchange integrand as

W µ,ν Hx,md [n] = ∂ Ŵ lr,ν/(1-ν) ee ∂ν Ψ µ [n] . (22)
The short-range MD exchange energy is then obtained as follows, according to Eq. ( 6),

E sr,µ x,md [n] = 1 µ/(1+µ) W µ,ν Hx,md [n] dν -E sr,µ H [n], (23)
leading to Eq. ( 21). We emphasize that, for µ > 0, Ψ µ [n] differs from the KS determinant. As a result, the expression for the "exchange" energy in Eq. ( 21) contains a correlation contribution, in addition to the short-range exchange energy. Note also that the complementary short-range correlation functional in Eq. ( 20) differs from the conventional one introduced in Eq. ( 13):

E sr,µ c,md [n] = E sr,µ c [n] + Φ KS [n]| Ŵ sr,µ ee |Φ KS [n] -Ψ µ [n]| Ŵ sr,µ ee |Ψ µ [n] . (24)
This also becomes clear when expressing the MD short-range correlation energy,

E sr,µ c,md [n] = 1 µ/(1+µ) W µ,ν c,md [n] dν, ( 25
)
in terms of the corresponding correlation integrand

W µ,ν c,md [n] = W ν Hxc [n] -W µ,ν Hx,md [n] = ∂ Ŵ lr,ν/(1-ν) ee ∂ν Ψ ν/(1-ν) [n] - ∂ Ŵ lr,ν/(1-ν) ee ∂ν Ψ µ [n] , (26)
which differs from the conventional integrand in Eq. ( 17) only in the use of Ψ µ [n] rather than Φ KS [n] in the last (subtracted) term. Note that the correlation energies E c [n] and E sr,µ c,md [n] are related in a simple manner, as seen by inserting Eq. ( 19) into Eq. ( 24) and rearranging,

Ψ µ [n]| T + Ŵee |Ψ µ [n] + E sr,µ c,md [n] = Φ KS [n]| T + Ŵee |Φ KS [n] + E c [n], (27)
where the long-range correlation and its coupling with the short-range interaction is contained in the expectation value on the left-hand side but in the correlation functional on the righthand side. Substitution of the srHxc energy decomposition in Eqs. (20) and ( 21) back into the ground-state energy expression in Eq. ( 9) leads to

E = min n Ψ µ [n]| T + Ŵee + Vne |Ψ µ [n] + E sr,µ c,md [n] . (28)
Since Ψ µ [n Ψ µ ]=Ψ µ , we conclude from Eq. ( 10) that the exact ground-state energy can be re-expressed as

E = Ψ µ | T + Ŵee + Vne |Ψ µ + E sr,µ c,md [n Ψ µ ]. (29)
We emphasize that the expression in Eq. ( 29) is exact when the energy is calculated from the self-consistent wavefunction in Eq. (10). We here introduce an approximation, the range-separated hybrid model with full-range integrals (RSHf), where the energy in Eq. ( 29)

is instead computed from the HF-srDFT wavefunction. A multi-configuration extension is obtained when using the MC-srDFT rather than HF-srDFT wavefunction, defining thus a range-separated multi-configuration hybrid model with full-range integrals (RSMCHf).

In the RSHf and RSMCHf schemes, the wavefunctions are optimized with the conventional short-range exchange-correlation density-functional, while the energy is computed with the alternative separation of exchange and correlation energies. As discussed in Sec. IV, this approach may not be sufficiently accurate when approximate functionals are used, especially when static correlation becomes important. The alternative decomposition of the short-range energy should then be used for the optimization of the wavefunction. However, unlike in srDFT, the minimization over densities in Eq. ( 28) cannot be replaced by a minimization over wavefunctions,

E = min Ψ Ψ| T + Ŵee + Vne |Ψ + E sr,µ c,md [n Ψ ] , (30)
simply because the minimizing wavefunction would be the ground state of a fully-interacting system, leading thus to double counting. On the other hand, invoking the one-to-one correspondence between densities and local potentials, a multi-determinant extension of KS-OEP schemes can be formulated.

D. Multi-determinant range-separated OEP approach

Exact formulation

When using the local potential rather than the density as a basic variable, the exact ground-state energy expression in Eq. ( 28) can be rewritten as 39

E =min v Ψ µ [v]| T + Ŵee + Vne |Ψ µ [v] + E sr,µ c,md [n Ψ µ [v] ] , (31)
where Ψ µ [v] is the ground state of the long-range interacting Hamiltonian with the local potential v:

Ψ µ [v] = arg min Ψ Ψ| T + Ŵ lr,µ ee |Ψ + dr v(r)n Ψ (r) . (32)
Since Ψ µ [v] is the solution to a (linear) eigenvalue equation, the formulation in Eqs. (31) and (32) will be referred to as non-self-consistent. According to Eq. ( 25), the MD short-range correlation energy vanishes as µ → +∞ and the minimizing potential in Eq. ( 31) is simply the nuclear potential. Regular wave-function theory is then recovered. On the other hand, when µ = 0, the energy in Eq. ( 31) reduces to a KS-OEP energy where the exact-exchange (EXX) term 56 is used in conjunction with the standard correlation density functional.

When 0 < µ < +∞, we obtain a rigorous combination of wave-function and KS-OEP density-functional approaches, referred to as srOEP in the following. According to Eqs. ( 11), ( 12) and ( 29), the exact minimizing potential is then given by

v µ (r) = v ne (r) + δE sr,µ Hxc δn(r) [n Ψ µ ], (33)
which can be rewritten as

v µ (r) = v ne (r) + δE sr,µ H δn(r) [n Ψ µ ] + δE sr,µ
x,md δn(r)

[n Ψ µ ] + δE sr,µ c,md δn(r) [n Ψ µ ], (34)
in terms of the srHxc decomposition in Eq. (20). Since the MD short-range correlation energy is here an explicit functional of the density, for which local density approximations have been proposed, 39,55 only the MD short-range exchange part needs to be optimized in Eq. (34), being an implicit functional of the density according to Eq. (21). Note that the srOEP energy in Eq. ( 31) can be rewritten as

E = min v Ψµ [v]| T + Ŵee + Vne | Ψµ [v] + E sr,µ c,md [n Ψµ [v] ] , (35)
where the auxiliary wavefunction Ψµ [v] is obtained for a given local potential v as follows:

Ψµ [v] = arg min Ψ Ψ| T + Ŵ lr,µ ee + Vne |Ψ + E sr,µ H [n Ψ ] + E sr,µ c,md [n Ψ ] + dr v(r) n Ψ (r) . (36)
Indeed, as Ψµ [v] satisfies the self-consistent equation T + Ŵ lr,µ ee + Vne + dr

δE sr,µ H δn(r) [n Ψµ [v] ] n(r) + dr δE sr,µ c,md δn(r) [n Ψµ [v] ] + v(r) n(r) | Ψµ [v] = Ẽµ [v]| Ψµ [v] , (37)
it is clear from Eq. ( 34) that the minimum in Eq. ( 35) is reached for the local potential

ṽµ (r) = δE sr,µ x,md δn(r) [n Ψ µ ] ( 38
)
since Ψµ [ṽ µ ] = Ψ µ . The formulation in Eqs. ( 35) and ( 36) will be referred to as self-consistent.

It is equivalent to the non-self-consistent formulation in Eq. ( 31) as long as there are no restrictions in the form of the optimized potential. This statement holds if an approximate MD short-range correlation density functional is used in conjunction with approximate long-range interacting HF or MCSCF wavefunctions, as considered in the rest of this work.

Since optimized potentials are usually expanded in a finite basis, 56 the non-self-consistent and self-consistent formulations in Eqs. ( 31) and ( 35), respectively, may give different results if the basis set in the first formulation is not sufficiently large to describe both short-range Hartree and MD correlation density-functional potentials accurately. The advantage of the second formulation lies in the fact that the basis set is used only to represent the MD exchange part of the short-range potential, see Eq. ( 38); the remaining Hartree and MD correlation short-range contributions are calculated as functional derivatives, see Eq. ( 37). The drawback with respect to the implementation is related to the computation of the energy gradient needed for optimizing the potential. As discussed in Sec. II E, the gradient requires the calculation of the linear response function for the long-range interacting wavefunction. The self-consistent formulation is thus less trivial to implement, requiring the implementation of second-order functional derivative contributions (kernel) to the linear response equations. 57 All approximate srOEP models introduced in the following are therefore based on the non-self-consistent formulation in Eq. (31).

Approximate formulations

As in srDFT, the exact auxiliary wavefunction in srOEP is multi-determinantal, being the ground state of a long-range interacting system. In the simplest HF-srOEP approach, the minimization in Eq. ( 32) is over single-determinantal wavefunctions Φ,

Φ µ [v] = arg min Φ Φ| T + Ŵ lr,µ ee |Φ + dr v(r)n Φ (r) . (39)
The HF-srOEP energy is then defined as

E srOEP HF = min v Φ µ [v]| T + Ŵee + Vne |Φ µ [v] + E sr,µ c,md [n Φ µ [v] ] . (40)
A multi-configurational extension is obtained by restricting the minimization in Eq. ( 32) to MCSCF wavefunctions that belong to a given active space S M :

Ψ µ M [v] = arg min Ψ∈S M Ψ| T + Ŵ lr,µ ee |Ψ + dr v(r) n Ψ (r) . (41)
This approach leads to the MC-srOEP energy expression

E srOEP MC = min v Ψ µ M [v]| T + Ŵee + Vne |Ψ µ M [v] + E sr,µ c,md [n Ψ µ M [v] ] , (42)
similar to the CAS-DFT energy expression of Refs. 12 and 18, based on the regular Hamiltonian and a complementary correlation functional. However, unlike in CAS-DFT, the correlation functional in the MC-srOEP method is universal in the sense that it does not depend on the active space S M . In addition, the minimization over local potentials (rather than over wavefunctions) ensures that the MCSCF model is applied to a long-range interacting

Hamiltonian. As a result, the active space can be enlarged to the full configuration-interaction (FCI) limit with no risk of double counting correlation effects.

If we now compare MC-srOEP with the MCOEP approach of Weimer et al. 23 , they differ in many respects. First, the MCOEP wavefunction is a linear combination of determinants constructed from KS-OEP-type orbitals. As a result, the optimization of the OEP requires the computation of the linear response of KS determinants related to changes in the potential, for which simple analytical expressions can be derived. In this respect, MC-srOEP is more complicated to implement, as the linear response of an MC long-range interacting wavefunction is required.

Second, an important difference between MC-srOEP and MCOEP methods relates to the active space. In the exact formulation of the MCOEP method, the density constructed from the MCOEP wave function for a fixed active space is equal to the exact ground-state density of the physical system and the exact ground-state energy is recovered. Therefore, the complementary density-functional correlation term, which describes dynamical correlation, depends on the active space. Developing approximate density functionals for this scheme is a difficult task as correlation effects may be double counted.

On the other hand, the "exact" MC-srOEP wavefunction gives the exact energy and density only in the FCI limit in an infinite basis set. Nevertheless, with relatively small µ values and the same active space as in the MCOEP expansion, the "exact" MC-srOEP density and energy should be close to the exact density and energy of the physical system, given that short-range effects are described by the exact MD short-range correlation functional and long-range effects are treated exactly in the given active space. 55 As already pointed out, the advantage of such a scheme is that the complementary MD short-range correlation density-functional does not depend on the active space, making it easier to model. In this work, the MD srLDA functional of Paziani et al. 49 is used.

The last approximation discussed here concerns the srOEP parameterization. Following

Wu and Yang, 56 we introduce an expansion of the potential

v(r) = v ne (r) + v sr,µ ref (r) + t b t g t (r), (43)
where the short-range analogue of the Fermi-Amaldi potential, calculated for a fixed Nelectron density n 0 , is employed as the reference potential:

v sr,µ ref (r) = N -1 N dr ′ n 0 (r ′ )w sr,µ ee (|r -r ′ |). (44)
We use the same basis set {g t } for the expansion of the potential and the molecular orbitals.

This parameterization allows for the use of analytic derivatives in quasi-Newton approaches to perform the optimization of Eqs. ( 40) and (42), and thus determine the potential expansion coefficients {b t }. As a first step, we here present the derivation of the HF-srOEP gradient.

The implementation of the analytical MC-srOEP gradient is in progress and will be presented in a separate paper. The MC-srOEP results presented in Sec. IV B 4 for H 2 were obtained numerically, by finite differences.

E. Analytical HF-srOEP energy gradient

The computation of the HF-srOEP energy in Eq. ( 40) can be performed with quasi-Newton approaches, using the coefficients {b t } in the potential expansion of Eq. ( 43) as variational parameters.

Let v 0 denote the trial potential defined by the initial set of coefficients {b 39) is denoted Φ µ 0 in the following. Variations in the potential coefficients

(0) t }. The associated determinant Φ µ [v 0 ] in Eq. (
v 0 (r) → v 0 (r) + t ǫ t g t (r) (45)
can be interpreted as static perturbations, where the property operators are the Gaussians g t with perturbation strengths ǫ t . We denote by i, j and a, b the occupied and unoccupied realvalued orbitals in Φ µ 0 , respectively. We use a second-quantized exponential parameterization 58 for the determinant Φ in Eq. (39),

|Φ(κ) = e -κ |Φ µ 0 , (46) where κ
= a,i κ ai Êai -Êia , Êai = â † a,α âi,α + â † a,β âi,β . (47)
The HF-srOEP energy gradient can then be expressed in terms of the orbital rotation vector

κ =      . . . κ ai . . .      (48)
as follows:

dE dǫ t 0 = ∂E ∂κ T 0 ∂κ ∂ǫ t 0 , (49)
with, according to Eq. ( 40),

E(κ) = Φ(κ)| T + Ŵee + Vne |Φ(κ) + E sr,µ c,md [n(κ)], (50)
and

n(κ, r) = Φ(κ)|n(r)|Φ(κ) . (51)
Note that the Hellmann-Feynman theorem cannot be applied in this context since the HF-srOEP energy depends implicitly on the potential. By analogy with regular HF theory, [START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF] the energy gradient components can be written as

∂E ∂κ ai 0 = -4 f ai + a|v sr,µ c,md [n Φ µ 0 ]|i , (52)
where f ai is the conventional Fock-matrix element computed with HF-srOEP orbitals, while the MD short-range correlation density-functional potential v sr,µ c,md [n Φ µ 0 ](r) = δE sr,µ c,md /δn(r)[n Φ µ 0 ] is calculated for the HF-srOEP density. As shown in the Appendix, the linear response vector is obtained as follows

E [2]µ ∂κ ∂ǫ t 0 = -g [1] t , (53)
where the gradient property vector is equal to

g [1] t = -4      . . . a|ĝ t |i . . .      . ( 54
)
The long-range analog of the HF Hessian, E [2]µ , is in the canonical HF-srOEP orbital basis equal to 58

E [2]µ ai,bj = 4 δ ab δ ij (ε µ a -ε µ i ) + 4 ab|ij lr,µ -ai|bj lr,µ -ai|jb lr,µ , (55)
where ε µ a and ε µ i are the unoccupied and occupied HF-srOEP orbital energies, respectively. We conclude from Eqs. ( 49) and ( 53) that the HF-srOEP energy gradient can be written as

dE dǫ t 0 = - ∂E ∂κ T 0 E [2]µ -1 g [1] t , (56)
or, equivalently,

dE dǫ t 0 = κ T g [1] t , (57)
where κ fulfills the linear response equation

E [2]µ κ = - ∂E ∂κ 0 . ( 58
)
Note that all components of the HF-srOEP energy gradient are thus computed from one single linear response vector κ. The latter can be obtained straightforwardly from a standard second-order HF wavefunction optimizer 58 when (i) using long-range integrals and substituting the trial srOEP for the nuclear potential in the Hessian and (ii) adding the MD short-range density-functional potential calculated for the HF-srOEP density to the Fock operator in the energy gradient.

We finally mention that, at µ = 0, the long-range integrals are zero and the MD short-range correlation density-functional potential reduces to the conventional correlation densityfunctional potential v c [n](r). As a result, the HF-srOEP determinant becomes the standard KS-OEP determinant Φ KS , the orbital energies reduce to conventional KS-OEP energies ε a and ε i , and the Hessian becomes diagonal:

E [2]0 ai,bj = 4δ ab δ ij (ε a -ε i ) . ( 59
)
We thus obtain from Eqs. ( 52), ( 54) and ( 56) the following analytical expression for the energy gradient

dE dǫ t 0 -→ µ=0 4 a,i a|ĝ t |i ε i -ε a × (f ai + a|v c [n Φ KS ]|i ) . ( 60
)
When the correlation potential is neglected, the KS-EXX energy gradient expression of Yang and Wu 56 is recovered.

F. Summary

In conventional MC-srDFT the KS decomposition of the complementary short-range exchange-correlation density-functional energy is used. Within the local density approximation, 48 the scheme will be referred to as MC-srLDA. The alternative separation of exchange and correlation energies that we investigate in this work relies on the multi-determinantal (MD) long-range interacting wavefunction rather than the KS determinant. As a result, long-range and short-range interactions can be recombined in the energy expression. The latter is thus rewritten as the sum of the expectation value for the regular Hamiltonian and a complementary short-range density-functional correlation energy that is referred to as MD. The long-range MC wavefunction to be inserted into this energy expression cannot be obtained straightfowardly by minimization over the wavefunction parameters otherwise double counting occurs.

Various approximations utilising the alternative MD decomposition are considered in this work. The simplest consists of using the MC-srLDA wavefunction. This approximation is referred to as RSMCHf. We may also consider, for analysis purposes, the single-determinant version of RSMCHf, that we refer to as RSHf and which consists of computing the energy with the HF-srLDA determinant rather than the MC-srLDA wavefunction. A more sophisticated procedure uses short-range OEPs. These can be optimized either at the HF or MC levels, leading to the HF-srOEP and MC-srOEP models. The analytical energy gradient has been derived and implemented for HF-srOEP. The implementation of the analytical MC-srOEP gradient is currently in progress. For analysis purposes, the long-range MC wavefunction can still be computed without reoptimization of the srOEP. This scheme, where a frozen effective potential (FEP) is employed, will be referred to as MC-srFEP. The HF-srOEP potential has been used as the srFEP in the following.

Working equations associated with all these schemes are given in Table I.

III. COMPUTATIONAL DETAILS

The various range-separated DFT schemes listed in Sec. II F have been implemented in a development version of the DALTON2011 program. [START_REF] Dalton | an ab initio electronic structure program[END_REF] The MD srLDA correlation functional of Paziani et al. 49 has been used. MC-srLDA wavefunctions and energies have been computed

with the srLDA exchange-correlation functional of Toulouse et al. 48 For the HF-srOEP and MC-srOEP approaches the minimizations of Eq. ( 40) and Eq. ( 42), respectively, were performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm 56 .

The initial Hessian was taken to be the approximate Hessian expression of Ref. 60. Convergence to a gradient norm below 10 -5 or a maximum absolute change in potential coefficients of less than 10 -12 is typically achieved in less than 20 iterations with this choice. For the HF-srOEP approach the energy gradient required at each iteration is computed according to Eqs. ( 57) and (58).

For the MC-srOEP approach, analytical gradients are not yet implemented (though we have derived their form), however, for analysis purposes in the present work calculations are carried out by determining the required gradient by finite difference. Since this substantially increases the number of energy evaluations to be performed we have considered the fully optimized MC-srOEP approach only for the H 2 molecule. Potential energy curves, equilibrium bond lengths and dissociation energies have been calculated for H 2 , Li 2 , N 2 and H 2 O. All calculations were performed with uncontracted cc-pVTZ basis sets 61,62 for both the orbital and potential expansions. Un-contraction of the basis sets and the use of the same sets for each expansion ensures smooth physically reasonable srOEP potentials are obtained.

The active orbital spaces used in the multi-configuration calculations are 1σ g 1σ u for H 2 and 2σ g 2σ u 1π u 3σ g 1π g 3σ u for N 2 and Li 2 . For H 2 O, the active orbital space is denoted 3.1.2.0, which signifies the number of orbitals in the a 1 .b 1 .b 2 .a 2 symmetries, respectively. The C 2 symmetry axis is along the z axis, and the σ v and σ ′ v mirror planes are σ v (xz) and σ v (yz), respectively.

IV. RESULTS AND DISCUSSION

A. Choice of the µ parameter

In the context of range-separated hybrid functionals, where range separation is used for the exchange energy only, the µ parameter is usually optimized semi-empirically for thermochemistry and other desired properties, leading thus to an average system-independent value in the range 0.4-0.5. [63][64][65][66][67][68][69] Quite recently, Baer et al. 70 proposed to choose µ such that Koopmans' theorem for both neutral and anion is obeyed, as closely as possible.

This procedure, which relies on first principles, enables one to tune the µ parameter for a given system. Let us stress that, in the exact theory, any µ value would provide the same (exact) ground-state energy. The problem of choosing µ occurs in practice because approximate short-range exchange functionals are used. In the context of multi-determinant range-separated DFT, where range separation is used for both exchange and correlation energies, the optimal choice of µ is even more problematic as, in practice, both approximate wavefunctions and density functionals are employed. Ángyán and coworkers 30,35,[71][72][73] use for example µ = 0.5 in their range-separated second-order Møller-Plesset (MP2) or random phase approximation (RPA) calculations on weakly interacting systems. This value has been calibrated in a completely different context, that is the exchange-only range-separated hybrid one, for reproducing atomization energies of small molecules. 67 In the particular case of the homonuclear rare-gas dimers, Goll et al. 31 alternatively proposed to choose for µ the inverse of the van der Waals radius in their range-separated coupled cluster (CC) calculations, so that intra-atomic correlations could be essentially treated in DFT.

In the context of MC-srDFT, Fromager et al. [START_REF] Fromager | [END_REF] investigated the possibility of choosing µ in such a way that static and dynamical correlations could be assigned to the long-range MCSCF and short-range density-functional correlation energies, respectively. As static correlation is usually not a purely long-range effect, even in the simple case of the dissociated H 2 molecule, 37 the authors focused on the dynamical correlation, suggesting that µ should be chosen small enough that the Coulomb hole is essentially treated in DFT. On the other hand it must also be chosen large enough that, in cases where static correlation becomes significant, the wave function can become sufficiently multi-configurational. The authors proposed from these considerations the following prescription: the largest value of µ for which the MC-srDFT wave function is well approximated by a single determinant, in systems where static correlation is not significant, should be considered as optimal. Let us stress that such a prescription does not guarantee that MC-srDFT will perform well when applied to systems with static correlation. It only ensures that the Coulomb hole is described within DFT and that the long-range part of the static correlation is assigned to MCSCF.

Obviously, with this choice, the complementary short-range correlation density-functional is expected to model the short-range part of the static correlation. As discussed further in Sec. IV, this can be problematic when stretching a bond for example. Numerical values for µ were obtained when analyzing long-range correlation effects as µ varies. Two strategies were proposed. The first one is based on the energy and consists of examining the total energy difference ∆E µ c between the MC and HF approximations for a given range-separated scheme. The best µ value is then determined by examining ∆E µ c for systems dominated by dynamical correlation and choosing the µ value at which this quantity falls below a threshold of -1 mE h . The second strategy examines, for the same systems, the natural orbital occupancies within the MC approximation and choosing the µ value at which these deviate from 2 (with a threshold of 10 -4 ). Calculations on a small test set of systems containing light elements all yielded the optimal µ = 0.4 value.

We now investigate whether this value is still optimal when a different separation of exchange and correlation energies is employed. We use for the discussion the H 2 molecule in its equilibrium geometry (R e = 0.741 Å74,75 ) as an example of system that is completely dominated by dynamical correlation. N 2 and Li 2 in their equilibrium geometries (R e =1.097 Å74 and 2.673 Å76,77 , respectively) will then be considered. These systems are interesting as they exhibit, Li 2 in particular, a multi-configurational character already at equilibrium. As discussed in the following, it is of course not as pronounced as in the dissociation limit but it is not negligible. Following Ref. 28 from zero (to within 10 -3 a.u.) for much smaller µ values than at the MC-srLDA level. For H 2 , for example, the RSMCHf and MC-srLDA µ values are 0.25 and 0.5, respectively.

It is tempting to conclude that the prescription of Fromager et al. [START_REF] Fromager | [END_REF] leads to different optimal µ values when considering RSMCHf energies. The situation is more subtle, however.

The ∆E µ c energy difference calculated within RSMCHf contains different correlation effects to the one obtained within the MC-srLDA scheme. In the latter case, ∆E µ c is essentially the purely long-range correlation energy while, within RSMCHf, it also contains the coupling between short-and long-range correlations that arises from the expectation value of the short-range interaction over the long-range correlated MC-srLDA wavefunction (see Eq. ( 29)).

This appears clearly in range-separated second-order density-functional perturbation theory (DFPT2) where a long-range MP2 description is used rather than a MCSCF one. 78 For small µ values, the difference between MC-srLDA and RSMCHf correlation energies can be rationalized by considering the Taylor expansions of the long-and short-range interactions,

w lr,µ ee (r 12 ) = 2 √ π µ - 1 3 µ 3 r 2 12 + O(µ 5 ) , w sr,µ ee (r 12 ) = 1 r 12 + O(µ). (61)
As the deviation of the exact long-range interacting wavefunction from the KS determinant varies as µ 3 (see Appendix B4 of Ref. 46), the long-range correlation energy varies as µ 6

and the coupling between long-and short-range correlations as µ 3 . As a result, the latter is expected to deviate more rapidly from zero as µ increases. This is the reason why a threshold of -1mE h on ∆E µ c will provide a smaller µ value for RSMCHf than for MC-srLDA. As shown in Fig. 1, the same conclusion can be drawn for MC-srFEP and MC-srOEP. This was expected as all three models use the same energy expression. They only differ by the long-range MC wavefunction that is inserted into this expression for the computation of the energy. If we want to follow the prescription of Fromager et al., [START_REF] Fromager | [END_REF] which relies on the analysis of purely long-range correlation effects, one should extract from ∆E µ c the purely long-range correlation energy or increase the threshold. As this analysis is performed on static-correlation free electronic systems, range-separated DFPT2 78 is expected to be a good approximation to RSMCHf, especially for small µ values. In this context, the coupling between long-and short-range correlations can be separated from the purely long-range correlation energy and each term can be computed when varying µ. Results obtained for rare gas atoms are presented in Ref. 78. When µ = 0.4, the coupling term equals -4 and -20 mE h in He and Ne, respectively, while the long-range correlation energy is above -1 mE h in both systems. This suggests that the coupling term is more system-dependent than the purely long-range correlation energy. It was expected as the former is expressed in terms of both long-and short-range integrals. In order to reduce the system-dependency of the µ parameter, one may want to examine purely long-range correlation effects only, leading thus to the optimal µ = 0.4 value for the rare gas atoms 78 .

Returning to H 2 and the RSMCHf model, we can still utilize, as an alternative, the second strategy of Ref. 28 that relies on the analysis of the natural orbital occupancies as µ increases from zero. Since RSMCHf and MC-srLDA wavefunctions are identical by definition, we can simply refer to Ref. 28 and conclude that µ = 0.4 is also optimal in this context. This is illustrated by the MC-srLDA occupancies of H 2 in Fig. 2. Interestingly, similar conclusions can be drawn for H 2 at both MC-srFEP and MC-srOEP levels. This should clearly be investigated on more static-correlation-free systems, once the analytical MC-srOEP energy gradient is implemented. This is left for future work. Note that µ = 0.4 is large enough to assign static correlation in Li 2 and N 2 , or at least a part of it, to the long-range MCSCF as suggested by their natural orbital occupancies in Fig. 2.

In summary, our preliminary calculations suggest that it is relevant to use µ = 0.4 in conjunction with the alternative separation of exchange and correlation energies of Toulouse et al. 39 B. Potential curves of H 2 , N 2 , Li 2 and H 2 O

RSMCHf equilibrium distances and binding energies

In this section, we compare the performance of the RSMCHf and MC-srLDA schemes for the description of the potential energy curves (PECs) of H 2 , N 2 , Li 2 and H 2 O. In the latter case, the symmetric dissociation at the experimental equilibrium angle H-O-H of 104.5 •79 was investigated. PECs as well as equilibrium bond distances and dissociation energies are given in Fig. 3 and Tables II and III, respectively. According to Sec. IV A the µ parameter has been set to 0.4.

For the H 2 molecule, close to the equilibrium geometry, MC-srLDA has a total energy that is too positive, though it recovers more short-range dynamical correlation than standard MCSCF. Although the qualitative shape of the curve is much better than restricted Hartree-Fock (not shown), the energy in the dissociation regime is much too positive. Overall this leads to a dissociation energy which is much too large, as shown in Table II. The RSMCHf approach gives energies close to equilibrium that are slightly too negative, whilst those at dissociation are significantly too negative. The result is a slight underestimation of the dissociation energy. The deviations from the exact curve in the dissociation regime come active natural orbitals as a function of µ for H 2 , N 2 and Li 2 at their experimental equilibrium geometries [74][75][76][77] .

from the complementary MD srLDA correlation functional, as illustrated by the RSMCHf (no src) PEC in Fig. 3, obtained when subtracting the former from the RSMCHf energy. As expected, the RSMCHf (no src) energy, which is equal to the expectation value of the regular Hamiltonian over the MC-srLDA wavefunction, is greater than the pure MCSCF energy for all bond distances. Returning to the equilibrium distance, the too negative MD srLDA correlation energy for H 2 at µ = 0.4 was already observed by Gori-Giorgi and Savin (see the uppermost panel in Fig. 7 of Ref. 55), who obtained an error (about 0.01E h in absolute value) close to that at the RSMCHf level (0.007E h ) when comparing with the "exact" energy of Ref. 75. In spite of these errors, the RSMCHf model clearly improves upon the dissociation energy obtained at the MC-srLDA level, reducing the absolute error from 27 to 7% for H 2 , with similar conclusions for the other systems.

The calculated equilibrium bond distances compare relatively well with experiment at both MC-srLDA and RSMCHf levels (see Table II). While the MC-srLDA model slightly overestimates the bond distance of H 2 by 0.015 Å, the RSMCHf model underestimates it by 0.02 Å. This slight over-binding is induced by the complementary MD srLDA correlation functional, as suggested by the RSMCHf (no src) equilibrium distance of 0.744 Å, which is almost equal to the experimental value. Note that MC-srLDA and HF-srLDA bond distances are almost identical since, for µ = 0.4, the MC-srLDA wave function is well approximated by a single determinant (see Sec. IV A and Fig. 4).

On the other hand, RSHf and RSMCHf equilibrium distances are quite similar but not as close as HF-srLDA and MC-srLDA distances (see Table II). The (slight) difference may be caused by the short/long-range correlation coupling, which contributes significantly to the RSMCHf correlation energy, where it is treated explicitly within the MCSCF, unlike for the MC-srLDA energy, where it is described within DFT (see Sec. IV A). Similar conclusions can be drawn for the other systems when comparing MC-srLDA with RSMCHf equilibrium distances, as shown in Tables II andIII.

The treatment of static correlation within RSMCHf

Regarding the H 2 dissociation limit, the large error in the MC-srLDA energy was interpreted in Ref. 28 as a self-interaction error of the complementary spin-unpolarised srLDA exchange-correlation functional. We emphasize that the complementary short-range correlation energy, whose exact expression is given in Eq. ( 18), is not supposed to be zero at large internuclear distances-instead, it should compensate the short-range Hartree and exchange energy contributions. This was shown by Teale et al., 37 who computed accurately the correlation integrand in Eq. ( 17) for various bond distances. According to Eq. ( 18), the accurate short-range correlation energy is obtained for µ = 0.4 when integrating the correlation integrand in Fig. 6 (f) of Ref. 37 from µ/(1 + µ) = 0.286 to 1. This quantity is obviously not equal to zero. The sum of short-range Hartree, exchange and correlation energies should, on the other hand, vanish in the dissociation limit. The corresponding integrand in Eq. ( 7) is plotted in Fig. 6 (d) of Ref. 37. Pure short-range exchange-correlation density functionals like srLDA are simply unable to compensate the short-range Hartree term. [START_REF] Fromager | [END_REF] The error in the total energy at large distances is significantly reduced when using the RSMCHf model, based on a different decomposition of the short-range exchange-correlation energy. The complementary short-range MD density-functional correlation energy is defined, in RSMCHf, with respect to the long-range interacting wavefunction rather than the KS non-interacting one. Its exact expression is given in Eq. ( 25) in terms of the short-range MD correlation integrand of Eq. ( 26). Since the long-range wavefunction in the dissociation limit of H 2 reduces to the Heitler-London wavefunction for all non-zero µ values, 55 this integrand and, consequently, the short-range MD correlation energy should vanish upon bond stretching. This is an important difference between MC-srDFT and RSMCHf schemes.

While the former is expected to describe much of the static correlation with the short-range correlation functional, the latter treats all static correlation with MCSCF, at least in the dissociation limit.

Comparing the RSMCHf and RSMCHf (no src) PECs of H 2 at large separations, it is clear that the short-range MD correlation energy is not well described within the local density approximation, as expected from the work of Gori-Giorgi and Savin. 55 The error can be interpreted as a self-correlation error of the spin-unpolarized MD srLDA functional since the molecule is correctly dissociated into two neutral hydrogen atoms (see the MC-srLDA natural orbital occupancies in Fig. 4). It is, however, not exclusively due to the functionalindeed, the RSMCHf (no src) PEC deviates slightly from the exact PEC. The deviation comes from the MC-srLDA wavefunction that is used to compute the RSMCHf (no src) energy. The former contains self-interaction errors because of the srLDA exchange-correlation density-functional potential. This problem will be addressed in the following sections.

Finally, we observe for H 2 a slight bump in the intermediate region (R = 2.25 Å) of the RSMCHf PEC. It is even more pronounced for H 2 O and N 2 . No bump appears for Li 2 , possibly because (unlike the other systems) it has a significant multi-configuration character already at equilibrium (see Fig. 4). Interestingly, at R = 2.117 Å, the RSMCHf (no src) energy of H 2 deviates from the exact one by about 0.03 E h , which could be interpreted as the correct value for the short-range MD correlation energy (the RSMCHf energy is relatively close to the exact one). However, from the srOEP calculations of Gori-Giorgi and Savin (the lowest panel in Fig. 7 of Ref. 55), the accurate value of this energy for µ = 0.4 is 0.01-0.015 E h . This difference suggests that the srLDA exchange-correlation density-functional potential, from which the MC-srLDA wavefunction (and hence the RSMCHf energy) is obtained, may not be accurate enough, especially as the wavefunction is strongly multi-configurational in this region, as shown in Fig. 4. Interestingly, errors on the short-range potential and the short-range MD correlation density functional seem to compensate in the RSMCHf energy when R = 2.117 Å but not for the other bond distances. Use of srOEPs rather than srLDA potentials is then a reasonable alternative, as is investigated in the following.

MC-srFEP results

The MC-srFEP model that was introduced in Sec. II F corresponds to a zero-iteration MC-srOEP calculation where the long-range MC wavefunction only is optimized. The initial HF-srOEP potential is simply frozen. Equilibrium bond distances and binding energies obtained at the MC-srFEP level are given in Tables II and III. Full PECs have been plotted in Fig. 3 for H 2 and Li 2 . Convergence problems in the HF-srOEP calculation occurred for 1.6 ≤ R ≤ 3.5 Å in N 2 and for all O-H distances beyond 2.0 Å in H 2 O-that is, when the long-range wavefunction becomes strongly multi-configurational, as shown in Fig. 4. In these cases, we expect the convergence of the srOEP to be manageable only at the long-range MCSCF level, requiring the implementation of analytical gradients, which is currently in progress.

For H 2 and Li 2 , we observe that the MC-srFEP and RSMCHf PECs remain relatively close for all bond distances, the natural orbital occupations being almost identical (see Fig. 4). For H 2 at equilibrium, the MC-srFEP energy is slightly lower than the RSMCHf energy, meaning that the HF-srOEP potential is better than the srLDA one, according to the variational principle in Eq. (42). However, as the bond is stretched and the wavefunction becomes multi-configurational, the MC-srFEP energy becomes higher than the RSMCHf one. This is expected since the HF-srOEP potential is unaffected by long-range correlation. For comparison, the MC-srFEP (no src) PEC has been computed. In this case the frozen potential is calculated at the HF-srOEP (no src) level, that is without the MD srLDA correlation term.

In the dissociation limit, the MC-srFEP (no src) energy is too high which clearly indicates that the HF-srOEP (no src) potential is a poor approximation to the exact short-range potential. Optimization of the srOEP at the long-range MCSCF level should, on the other hand, give essentially the exact solution. 55 It is noteworthy that, in a minimal atomic-orbital basis, the exact energy would actually be obtained at both RSMCHf (no src) and MC-srFEP (no src) levels of theory since the wavefunction is then fixed to 1/ √ 2 |1σ 2 g -|1σ 2 u , where the bonding 1σ g and the anti-bonding 1σ u molecular orbitals are simply linear combinations of the atomic 1s orbitals of each hydrogen atom. 38 In our calculations, as larger basis are used, orbitals can rotate so that different energies can be obtained.

Returning to RSMCHf, the srLDA density-functional potential, even though it may be a crude approximation to the exact potential, includes multi-configuration effects through the density. It is therefore more accurate than the HF-srOEP potential. Note also that, in the intermediate region, the maximum of the bump is at a higher energy with the MC-srFEP model than with the RSMCHf model. In conclusion, MC-srFEP provides better binding energies than RSMCHf but this essentially relies on error compensation.

MC-srOEP results for H 2

We now consider the calculation of the PEC of H 2 at the MC-srOEP level. For comparison, the MC-srOEP (no src) PEC was also calculated. These PECs are shown in Fig. 3. As expected, the MC-srOEP and MC-srFEP PECs are almost identical in the equilibrium region where static correlation is negligible. As in the MC-srLDA scheme, the MC-srOEP wavefunction is essentially a single determinant (see Fig. 4).

Upon bond stretching, the MC-srOEP PEC deviates from the MC-srFEP curve. Interestingly, the separation occurs when the MC-srFEP energy is higher than the RSMCHf one, that is when the HF-srOEP potential becomes less accurate than the srLDA one, as discussed previously. The difference between MC-srOEP and MC-srFEP energies is reflected both in the orbital occupations and in the orbitals (not shown). The MC-srOEP wavefunction has a much more pronounced multi-configurational character than the MC-srFEP wavefunction, as shown in Fig. 4. This explains why the MC-srOEP energy reaches it asymptotic limit at a shorter bond distance (R ≈ 2.2 Å) than does the MC-srFEP energy (R ≈ 3.5 Å). As expected, the exact energy is recovered at the MC-srOEP (no src) level in the dissociation limit. Due to self-correlation errors in the MD srLDA functional, the MC-srOEP energy is too low at the dissociation. This leads to an underestimation of the binding energy, which is actually even more pronounced than for RSMCHf and MC-srFEP.

In the intermediate region, at R = 2.117 Å, the MC-srOEP (no src) energy differs from the exact one by 0.017E h , as expected for the accurate short-range MD correlation energy from the work of Gori-Giorgi and Savin. 55 The MD srLDA correlation energy obtained by us (-0.03E h ) agrees perfectly with their value, see the lowest panel in Fig. 7 of Ref. 55. Note also that the bump observed at the RSMCHf and MC-srFEP levels is significantly reduced at the MC-srOEP level but not completely removed.

In conclusion, our preliminary MC-srOEP calculations on H 2 confirm that the MD srLDA correlation functional of Paziani et al. 49 can be improved upon. The conclusion is basically the same as for MC-srLDA: better short-range functionals should be developed. However, the exact complementary short-range correlation energy to be modelled (namely the MD one) is supposed to vanish in the dissociation limit, unlike the KS short-range correlation energy that is used in conventional MC-srDFT models. It may therefore be easier to develop better short-range MD correlation density functionals. The ab initio calculation of the range-separated adiabatic connection 37 is a valuable tool for such developments.

Analysis of the srOEPs

When applying the OEP method care must be taken to ensure that smooth potentials are obtained. In the present work we have used uncontracted basis sets and chosen the orbital and potential sets identical to help ensure this is the case. In Fig. 5 we present the exchange-correlation contributions to the HF-srOEP potentials. For this figure we have generated the potentials using the smoothing-norm approach of Heaton-Burgess et al. [START_REF] Heaton-Burgess | [END_REF]81 with a smoothing parameter of 10 -5 . This value perturbs the energy by less than 10 -5 a.u. from the unconstrained energies presented in the rest of this work. The result is a a potential that is everywhere smooth. Without the application of this approach the potential very close to the nuclei exhibits a large spike, owing to the fact that this region has essentially no contribution when determining the energy of the system. The potential in other regions of space is essentially unchanged from its unconstrained counterpart.

In the top panel of Fig. 5 the potentials exhibit relatively few features, although as µ increases the positions of the nuclei can be discerned. As would be expected the from the discussion in Sec. II the srOEP potential approaches zero as the value of µ is increased. The other feature of the potentials that is clearly visible is their rate of decay as a function of µ.

For µ = 0.00 a KS-OEP is obtained and exhibits the usual -1/r asymptotic decay. As mu is increased the decay of the potential becomes more rapid, reflecting the form of Eq. (44).

For µ = 0.4 the potential contributions along the bond axis are essentially zero beyond approximately 5 a.u. from each atom. Interestingly, the potential near the atoms is also affected rather strongly by the change in µ, this may again reflect the inclusion of long-range / short-range coupling effects in the partitioning of Toulouse et al. 39 . For the longer bond lengths in the lower two panels of Fig. 5 similar conclusions can be reached, though as the atoms are further separated their positions become clearer in the associated potentials.

For the MC-srOEP approach the exchange-correlation potentials obtained (not shown) are very similar, having only very slightly more negative potentials surrounding the nuclei and slightly more positive potentials further away in the intermediate "shoulder" regions.

Whilst these subtle changes are essential for proper optimization of the MC-srOEP energy, they do not significantly alter the potentials from those of the HF-srOEP approach on the scale shown.

V. CONCLUSIONS

An alternative separation of exchange and correlation energies has been investigated in the context of multi-configuration range-separated DFT. The new decomposition of the shortrange exchange-correlation energy relies on the auxiliary long-range interacting wavefunction rather than the KS determinant. This approach, first proposed by Toulouse et al., 39 has two advantages relative to the traditional KS decomposition. First, the MCSCF part of the energy is now computed with the regular (fully-interacting) Hamiltonian, following CAS-DFT approaches. 12,18 Second, the exact complementary short-range correlation energy vanishes upon dissociation of H 2 , meaning that the static correlation is fully assigned to MCSCF.

The drawback is that, because of double counting, the long-range interacting wave function used to compute the energy cannot be obtained by minimizing the energy expression with respect to the wavefunction parameters. Different approaches that overcome this problem have been investigated. The first simply computes the energy from the long-range interacting MCSCF wavefunction that is optimized with the traditional KS decomposition of the shortrange density-functional energy. A more sophisticated scheme uses short-range OEPs. The resulting combination of OEP techniques with wavefunction theory has been investigated in this work, at the HF and MCSCF levels.

In the HF case, an analytical expression for the energy gradient has been derived and implemented. Calculations have been performed within the short-range local density approximation on H 2 , N 2 , Li 2 and H 2 O. Significant improvements in binding energies are obtained with the new decomposition of the short-range energy, relative to the traditional one. The importance of optimizing the short-range OEP at the MCSCF level when static correlation becomes significant has been demonstrated for H 2 , using a finite-difference gradient. For further assessment of this approach, the analytical gradient is under implementation.

Our preliminary calculations indicate that the local density approximation is not accurate enough for modelling the complementary short-range correlation energy-better functionals may be developed from an accurate calculation of range-separated adiabatic-connection paths.

Work is in progress in this direction.

APPENDIX: HF-SROEP LINEAR RESPONSE EQUATION

The response of the HF-srOEP determinant related to variations ǫ t in the srOEP coefficients is obtained when considering the auxiliary energy

E µ (κ, ǫ t ) = Φ(κ)| T + Ŵ lr,µ ee + V0 |Φ(κ) + t ǫ t Φ(κ)|ĝ t |Φ(κ) = E µ (κ) + t ǫ t Φ(κ)|ĝ t |Φ(κ) , (62)
where V0 = dr v 0 (r) n(r) is the trial srOEP and the Gaussian operator ĝt = dr g t (r) n(r) is analogous to a property operator in response theory 82 . From the variational condition in Eq. ( 39) we obtain

∀ǫ t ∂E µ (κ, ǫ t ) ∂κ κ(ǫt) = 0, (63)
which leads to the linear response equation

d dǫ t ∂E µ (κ, ǫ t ) ∂κ 0 = 0. ( 64
)
Using the Taylor expansion through second order

E µ (κ) = E µ (0) + 1 2 κ T E [2]µ κ + . . . , (65)
where the HF-type Hessian E [2]µ is constructed from the auxiliary long-range interacting Hamiltonian T + Ŵ lr,µ ee + V0 , and rewriting the first-order derivative of the Gaussian expectation value as 58

∂ ∂κ Φ(κ)|ĝ t |Φ(κ) 0 = g [1] t , (66)
where the gradient Gaussian vector expression is given in Eq. ( 54), we finally obtain Eq. ( 53).

20 -Figure 1 .

 201 Figure 1. The quantity ∆E µ c for (from top to bottom) H 2 , N 2 and Li 2 for each method as a function of the parameter µ. See text for further details.

Figure 2 .

 2 Figure 2. Occupation numbers of MC-srLDA (dotted), MC-srFEP (dashed) and MC-srOEP (solid)

Figure 3 .

 3 Figure 3. Potential energy curves (a. u.) of H 2 (upper left panel) and interaction (binding) energies of N 2 (upper right panel), H 2 O (lower left panel) and Li 2 (lower right panel) obtained by means of the new multi-configuration range-separated schemes. The coupling parameter is µ=0.4. The exact interaction energy and PEC curves are taken from Ref. 75 for H 2 and Ref. 74 for N 2 .

Figure 4 .

 4 Figure 4. Occupation numbers of MC-srOEP (solid), MC-srFEP (dotted) and MC-srLDA (dashed) active natural orbitals as a function of bond distance ( Å) (from top to bottom) H 2 , Li 2 and N 2 .The parameter µ is set to 0.4.

Figure 5 .

 5 Figure 5. The exchange-correlation contributions to the srOEPs plotted along the bond axis for the H 2 molecule at the HF-srOEP level. The top panel corresponds to R = 1.3 Å, the middle to R = 2.8 Åand the bottom to R = 4.2 Å. In each panel the potentials are shown for µ = 0.00, 0.10, 0.20, 0.30, 0.40, 0.50, 0.75, 1.00, 2.00, 3.00, 4.00, 5.00 and these may be distinguished by noting that the potential value at z = 0.0 increases with increasing µ.

FIGURE CAPTIONS Figure 1 :

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 :

 2 Figure 2: Occupation numbers of MC-srLDA (dotted), MC-srFEP (dashed) and MC-srOEP (solid) active natural orbitals as a function of µ for H 2 , N 2 and Li 2 at their experimental equilibrium geometries 74-77 .

Figure 3 :

 3 Figure 3: Potential energy curves (a. u.) of H 2 (upper left panel) and interaction (binding) energies of N 2 (upper right panel), H 2 O (lower left panel) and Li 2 (lower right panel) obtained by means of the new multi-configuration range-separated schemes. The range separation parameter is µ=0.4. The exact interaction energy and PEC curves are taken from Ref. 75 for H 2 and Ref. 74 for N 2 .

Figure 4 :

 4 Figure 4: Occupation numbers of MC-srLDA (dashed), MC-srFEP (dotted) and MC-srOEP (solid) active natural orbitals as a function of bond distance ( Å) (from top to bottom) H 2 , Li 2 and N 2 74-77 . The parameter µ is set to 0.4.

Figure 5 :

 5 Figure 5: The exchange-correlation contributions to the srOEPs plotted along the bond axis for the H 2 molecule at the HF-srOEP level. The top panel corresponds to R = 1.3 Å, the middle to R = 2.8 Å and the bottom to R = 4.2 Å. In each panel the potentials are shown for µ = 0.00, 0.10, 0.20, 0.30, 0.40, 0.50, 0.75, 1.00, 2.00, 3.00, 4.00, 5.00 and these may be distinguished by noting that the potential value at z = 0.0 increases with increasing µ.

Table II .

 II Stoyanova et al., Journal of Chemical Physics

		R e	D e
	H 2
	HF-srLDA	0.755 10.39
	RSHf	0.717 11.73
	HF-srOEP	0.715 11.75
	MCSCF	0.755 4.14
	MC-srLDA	0.756 6.05
	RSMCHf	0.724 4.41
	RSMCHf (no src)	0.744 3.91
	MC-srFEP	0.720 4.53
	MC-srFEP (no src) 0.743 4.09
	MC-srOEP	0.722 4.18
	MC-srOEP (no src) 0.743 3.77
	Exp.	0.741 a 4.75 a
	N 2
	HF-srLDA	1.082 30.85
	RSHf	1.051 36.86
	HF-srOEP	1.052 34.55
	MCSCF	1.105 9.19
	MC-srLDA	1.087 16.18
	RSMCHf	1.079 7.96
	RSMCHf (no src)	1.096 6.86
	MC-srFEP	1.077 8.34
	MC-srFEP (no src) 1.095 7.44
	Exp.	1.097 a 9.91 a

a Ref.

74

ACKNOWLEDGMENTS E.F. and A.S. thank ANR (DYQUMA project) as well as Hans Jørgen Aa. Jensen, Andreas Savin and Yann Cornaton for fruitful discussions. T. H. and A. M. T. acknowledge supported by the Norwegian Research Council through the CoE Centre for Theoretical and Computational Chemistry (CTCC) Grant No. 179568/V30 and by the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683. A. M. T. is also grateful for support from the Royal Society University Research Fellowship scheme. A. S. acknowledges support from the CTCC for a visit to Oslo.

 I: Wavefunction, local potential and energy expressions associated with all rangeseparated methods discussed in this work. Single determinantal trial wavefunctions are denoted Φ. Ĥ = T + Ŵee + Vne and Ĥlr,µ = T + Ŵ lr,µ ee + Vne correspond to the physical (fully interacting) and long-range interacting Hamiltonians, respectively. S M denotes the active space used in the long-range MCSCF calculation.

Table III