\documentclass[aip,jcp,reprint,noshowkeys]{revtex4-1}
%\documentclass[aip,jcp,noshowkeys]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,mhchem,longtable}

\usepackage{mathpazo,libertine}
\usepackage[normalem]{ulem}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{RGB}{0, 180, 0}
\newcommand{\beurk}[1]{\textcolor{darkgreen}{#1}}
\newcommand{\trash}[1]{\textcolor{red}{\sout{#1}}}
\usepackage{xspace}

\usepackage{hyperref}
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    filecolor=blue,      
    urlcolor=blue,
    citecolor=blue
}
\newcommand{\cdash}{\multicolumn{1}{c}{---}}
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\mr}{\multirow}
\newcommand{\SI}{\textcolor{blue}{supporting information}}

% second quantized operators
\newcommand{\psix}[1]{\hat{\Psi}\left({\bf X}_{#1}\right)}
\newcommand{\psixc}[1]{\hat{\Psi}^{\dagger}\left({\bf X}_{#1}\right)}
\newcommand{\ai}[1]{\hat{a}_{#1}}
\newcommand{\aic}[1]{\hat{a}^{\dagger}_{#1}}
\newcommand{\vijkl}[0]{V_{ij}^{kl}}
\newcommand{\phix}[2]{\phi_{#1}(\bfr{#2})}
\newcommand{\phixprim}[2]{\phi_{#1}(\bfr{#2}')}


%operators 
\newcommand{\elemm}[3]{{\ensuremath{\bra{#1}{#2}\ket{#3}}\xspace}}
\newcommand{\ovrlp}[2]{{\ensuremath{\langle #1|#2\rangle}\xspace}}

%\newcommand{\ket}[1]{{\ensuremath{|#1\rangle}\xspace}}
%\newcommand{\bra}[1]{{\ensuremath{\langle #1|}\xspace}}

% 

	
% energies
\newcommand{\Ec}{E_\text{c}}
\newcommand{\EPT}{E_\text{PT2}}
\newcommand{\EsCI}{E_\text{sCI}}
\newcommand{\EDMC}{E_\text{DMC}}
\newcommand{\EexFCI}{E_\text{exFCI}}
\newcommand{\EexFCIbasis}{E_\text{exFCI}^{\Bas}}
\newcommand{\EexFCIinfty}{E_\text{exFCI}^{\infty}}
\newcommand{\EexDMC}{E_\text{exDMC}}
\newcommand{\Ead}{\Delta E_\text{ad}}
\newcommand{\efci}[0]{E_{\text{FCI}}^{\Bas}}
\newcommand{\emodel}[0]{E_{\model}^{\Bas}}
\newcommand{\emodelcomplete}[0]{E_{\model}^{\infty}}
\newcommand{\efcicomplete}[0]{E_{\text{FCI}}^{\infty}}
\newcommand{\ecccomplete}[0]{E_{\text{CCSD(T)}}^{\infty}}
\newcommand{\ecc}[0]{E_{\text{CCSD(T)}}^{\Bas}}
\newcommand{\efuncbasisFCI}[0]{\bar{E}^\Bas[\denFCI]}
\newcommand{\efuncbasisfci}[0]{\bar{E}^\Bas[\denfci]}
\newcommand{\efuncbasis}[0]{\bar{E}^\Bas[\den]}
\newcommand{\efuncden}[1]{\bar{E}^\Bas[#1]}
\newcommand{\efuncdenpbe}[1]{\bar{E}_{\text{X}}^\Bas[#1]}
\newcommand{\ecompmodel}[0]{\bar{E}^\Bas[\denmodel]}
\newcommand{\ecmubis}[0]{\bar{E}_{\text{c,md}}^{\text{sr}}[\denr;\,\mu]}
\newcommand{\ecmubisldapbe}[0]{\bar{E}_{\text{c,md}\,\text{PBE}}^{\text{sr}}[\denr;\,\mu]}
\newcommand{\ecmuapprox}[0]{\bar{E}_{\text{c,md-}\mathcal{X}}^{\text{sr}}[\den;\,\mu]}
\newcommand{\ecmuapproxmur}[0]{\bar{E}_{\text{c,md-}\mathcal{X}}^{\text{sr}}[\den;\,\mur]}
\newcommand{\ecmuapproxmurfci}[0]{\bar{E}_{\text{c,md-}\mathcal{X}}^{\text{sr}}[\denfci;\,\mur]}
\newcommand{\ecmuapproxmurmodel}[0]{\bar{E}_{\text{c,md-}\mathcal{X}}^{\text{sr}}[\denmodel;\,\mur]}
\newcommand{\ecompmodellda}[0]{\bar{E}_{\text{LDA}}^{\Bas,\wf{}{\Bas}}[\denmodel]}
\newcommand{\ecompmodelldaval}[0]{\bar{E}_{\text{LDA, val}}^{\Bas,\wf{}{\Bas}}[\den]}
\newcommand{\ecompmodelpbe}[0]{\bar{E}_{\text{PBE}}^{\Bas,\wf{}{\Bas}}[\den]}
\newcommand{\ecompmodelpbeval}[0]{\bar{E}_{\text{PBE, val}}^{\Bas,\wf{}{\Bas}}[\den]}
\newcommand{\emulda}[0]{\bar{\varepsilon}^{\text{sr},\text{unif}}_{\text{c,md}}\left(\denr;\mu({\bf r};\wf{}{\Bas})\right)}
\newcommand{\emuldamodel}[0]{\bar{\varepsilon}^{\text{sr},\text{unif}}_{\text{c,md}}\left(\denmodelr;\mu({\bf r};\wf{}{\Bas})\right)}
\newcommand{\emuldaval}[0]{\bar{\varepsilon}^{\text{sr},\text{unif}}_{\text{c,md}}\left(\denval ({\bf r});\murval;\wf{}{\Bas})\right)}
\newcommand{\ecmd}[0]{\varepsilon^{\text{c,md}}_{\text{PBE}}}
\newcommand{\psibasis}[0]{\Psi^{\basis}}

%pbeuegxiHF
\newcommand{\pbeuegxihf}{\text{PBE-UEG-}\zeta\text{-HF}^\Bas}
\newcommand{\argpbeuegxihf}[0]{\den,\zeta,s,\ntwo_{\text{UEG}},\mu_{\text{HF}}^{\basis}}
\newcommand{\argrpbeuegxihf}[0]{\den(\br{}),\zeta(\br{}),s(\br{}),\ntwo_{\text{UEG}}(\br{}),\mu_{\text{HF}}^{\basis}(\br{})}
%pbeuegxiCAS
\newcommand{\pbeuegxicas}{\text{PBE-UEG-}\zeta\text{-CAS}^\Bas}
\newcommand{\argpbeuegxicas}[0]{\den,\zeta,s,\ntwo_{\text{UEG}},\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeuegxicas}[0]{\den(\br{}),\zeta(\br{}),s(\br{}),\ntwo_{\text{UEG}}(\br{}),\mu_{\text{CAS}}^{\basis}(\br{})}
%pbeuegXiCAS
\newcommand{\pbeuegXicas}{\text{PBE-UEG-}\tilde{\zeta}\text{-CAS}^\Bas}
\newcommand{\argpbeuegXicas}[0]{\den,\tilde{\zeta},s,\ntwo_{\text{UEG}},\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeuegXicas}[0]{\den(\br{}),\tilde{\zeta}(\br{}),s(\br{}),\ntwo_{\text{UEG}}(\br{}),\mu_{\text{CAS}}^{\basis}(\br{})}
%pbeontxiCAS
\newcommand{\pbeontxicas}{\text{PBE-ONT-}\zeta\text{-CAS}^\Bas}
\newcommand{\argpbeontxicas}[0]{\den,\zeta,s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeontxicas}[0]{\den(\br{}),\zeta(\br{}),s(\br{}),\ntwoextrapcas(\br{}),\mu_{\text{CAS}}^{\basis}(\br{})}
%pbeontXiCAS
\newcommand{\pbeontXicas}{\text{PBE-ONT-}\tilde{\zeta}\text{-CAS}^\Bas}
\newcommand{\argpbeontXicas}[0]{\den,\tilde{\zeta},s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeontXicas}[0]{\den(\br{}),\tilde{\zeta}(\br{}),s(\br{}),\ntwoextrapcas(\br{}),\mu_{\text{CAS}}^{\basis}(\br{})}
%pbeont0xiCAS
\newcommand{\pbeontnscas}{\text{PBE-ONT-}0\tilde{\zeta}\text{-CAS}^\Bas}
\newcommand{\argpbeontnscas}[0]{\den,0,s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeontnscas}[0]{\den(\br{}),0,s(\br{}),\ntwoextrapcas(\br{}),\mu_{\text{CAS}}^{\basis}(\br{})}

%%%%%% arguments 

\newcommand{\argepbe}[0]{\den,\zeta,s}
\newcommand{\argebasis}[0]{\den,\zeta,s,\ntwo,\mu_{\Psi^{\basis}}}
\newcommand{\argecmd}[0]{\den,\zeta,s,\ntwo,\mu}
\newcommand{\argepbeueg}[0]{\den,\zeta,s,\ntwo_{\text{UEG}},\mu_{\Psi^{\basis}}}
\newcommand{\argepbeontxicas}[0]{\den,\zeta,s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argepbeuegXihf}[0]{\den,\tilde{\zeta},s,\ntwo_{\text{UEG}},\mu_{\Psi^{\basis}}}
\newcommand{\argrebasis}[0]{\denr,\zeta(\br{}),s,\ntwo(\br{}),\mu_{\Psi^{\basis}}(\br{})}
\newcommand{\argrebasisab}[0]{\denr,\zeta(\br{}),s,\ntwo(\br{}),\mu_{\Psi^{\basis}}(\br{})}


% numbers
\newcommand{\rnum}[0]{{\rm I\!R}}
\newcommand{\bfr}[1]{{\bf r}_{#1}}
\newcommand{\dr}[1]{\text{d}\bfr{#1}}
\newcommand{\rr}[2]{\bfr{#1}, \bfr{#2}}
\newcommand{\rrrr}[4]{\bfr{#1}, \bfr{#2},\bfr{#3},\bfr{#4}  }



% effective interaction
\newcommand{\twodm}[4]{\elemm{\Psi}{\psixc{#4}\psixc{#3} \psix{#2}\psix{#1}}{\Psi}}
\newcommand{\murpsi}[0]{\mu({\bf r};\wf{}{\Bas})}
\newcommand{\ntwo}[0]{n^{(2)}}
\newcommand{\ntwohf}[0]{n^{(2),\text{HF}}}
\newcommand{\ntwophi}[0]{n^{(2)}_{\phi}}
\newcommand{\ntwoextrap}[0]{\tilde{n}^{(2)}_{\psibasis}}
\newcommand{\ntwoextrapcas}[0]{\tilde{n}^{(2)\,\basis}_{\text{CAS}}}
\newcommand{\mur}[0]{\mu({\bf r})}
\newcommand{\murr}[1]{\mu({\bf r}_{#1})}
\newcommand{\murval}[0]{\mu_{\text{val}}({\bf r})}
\newcommand{\murpsival}[0]{\mu_{\text{val}}({\bf r};\wf{}{\Bas})}
\newcommand{\murrval}[1]{\mu_{\text{val}}({\bf r}_{#1})}
\newcommand{\weeopmu}[0]{\hat{W}_{\text{ee}}^{\text{lr},\mu}}


\newcommand{\wbasis}[0]{W_{\wf{}{\Bas}}(\bfr{1},\bfr{2})}
\newcommand{\wbasiscoal}[0]{W_{\wf{}{\Bas}}(\bfr{},\bfr{})}
\newcommand{\wbasisval}[0]{W_{\wf{}{\Bas}}^{\text{val}}(\bfr{1},\bfr{2})}
\newcommand{\fbasis}[0]{f_{\wf{}{\Bas}}(\bfr{1},\bfr{2})}
\newcommand{\fbasisval}[0]{f_{\wf{}{\Bas}}^{\text{val}}(\bfr{1},\bfr{2})}
 \newcommand{\ontop}[2]{ n^{(2)}_{#1}({\bf #2}_1)}
 \newcommand{\twodmrpsi}[0]{ \ntwo_{\wf{}{\Bas}}(\rrrr{1}{2}{2}{1})}
 \newcommand{\twodmrdiagpsi}[0]{ \ntwo_{\wf{}{\Bas}}(\rr{1}{2})}
 \newcommand{\twodmrdiagpsival}[0]{ \ntwo_{\wf{}{\Bas},\,\text{val}}(\rr{1}{2})}
 \newcommand{\gammamnpq}[1]{\Gamma_{mn}^{pq}[#1]}
 \newcommand{\gammamnkl}[0]{\Gamma_{mn}^{kl}}
 \newcommand{\gammaklmn}[1]{\Gamma_{kl}^{mn}[#1]}
%\newcommand{\wbasiscoal}[1]{W_{\wf{}{\Bas}}({\bf r}_{#1})}
\newcommand{\ontoppsi}[1]{ n^{(2)}_{\wf{}{\Bas}}(\bfr{#1},\barr{#1},\barr{#1},\bfr{#1})}
\newcommand{\wbasiscoalval}[1]{W_{\wf{}{\Bas}}^{\text{val}}({\bf r}_{#1})}
\newcommand{\ontoppsival}[1]{ n^{(2)}_{\wf{}{\Bas}}^{\text{val}}(\bfr{#1},\barr{#1},\barr{#1},\bfr{#1})}



\newcommand{\ex}[4]{$^{#1}#2_{#3}^{#4}$}
\newcommand{\ra}{\rightarrow}
\newcommand{\De}{D_\text{e}}

% MODEL 
\newcommand{\model}[0]{\mathcal{Y}}

% densities 
\newcommand{\denmodel}[0]{\den_{\model}^\Bas}
\newcommand{\denmodelr}[0]{\den_{\model}^\Bas ({\bf r})}
\newcommand{\denfci}[0]{\den_{\psifci}}
\newcommand{\denFCI}[0]{\den^{\Bas}_{\text{FCI}}}
\newcommand{\denhf}[0]{\den_{\text{HF}}^\Bas}
\newcommand{\denrfci}[0]{\denr_{\psifci}}
\newcommand{\dencipsir}[0]{{n}_{\text{CIPSI}}^\Bas({\bf r})}
\newcommand{\dencipsi}[0]{{n}_{\text{CIPSI}}^\Bas}
\newcommand{\den}[0]{{n}}
\newcommand{\denval}[0]{{n}^{\text{val}}}
\newcommand{\denr}[0]{{n}({\bf r})}
\newcommand{\onedmval}[0]{\rho_{ij,\sigma}^{\text{val}}}

% wave functions
\newcommand{\psifci}[0]{\Psi^{\Bas}_{\text{FCI}}}
\newcommand{\psimu}[0]{\Psi^{\mu}}
% operators
\newcommand{\weeopbasis}[0]{\hat{W}_{\text{ee}}^\Bas}
\newcommand{\kinop}[0]{\hat{T}}

\newcommand{\weeopbasisval}[0]{\hat{W}_{\text{ee}}^{\Basval}}
\newcommand{\weeop}[0]{\hat{W}_{\text{ee}}}


% units
\newcommand{\IneV}[1]{#1 eV}
\newcommand{\InAU}[1]{#1 a.u.}
\newcommand{\InAA}[1]{#1 \AA}


% methods
\newcommand{\UEG}{\text{UEG}}
\newcommand{\LDA}{\text{LDA}}
\newcommand{\PBE}{\text{PBE}}
\newcommand{\FCI}{\text{FCI}}
\newcommand{\CCSDT}{\text{CCSD(T)}}
\newcommand{\lr}{\text{lr}}
\newcommand{\sr}{\text{sr}}

\newcommand{\Nel}{N}
\newcommand{\V}[2]{V_{#1}^{#2}}


\newcommand{\n}[2]{n_{#1}^{#2}}
\newcommand{\E}[2]{E_{#1}^{#2}}
\newcommand{\bE}[2]{\Bar{E}_{#1}^{#2}}
\newcommand{\bEc}[1]{\Bar{E}_\text{c}^{#1}}
\newcommand{\e}[2]{\varepsilon_{#1}^{#2}}
\newcommand{\be}[2]{\Bar{\varepsilon}_{#1}^{#2}}
\newcommand{\bec}[1]{\Bar{e}^{#1}}
\newcommand{\wf}[2]{\Psi_{#1}^{#2}}
\newcommand{\W}[2]{W_{#1}^{#2}}
\newcommand{\w}[2]{w_{#1}^{#2}}
\newcommand{\hn}[2]{\Hat{n}_{#1}^{#2}}
\newcommand{\rsmu}[2]{\mu_{#1}^{#2}}
\newcommand{\SO}[2]{\phi_{#1}(\br{#2})}

\newcommand{\modX}{\text{X}}
\newcommand{\modY}{\text{Y}}

% basis sets 
\newcommand{\setdenbasis}{\mathcal{N}_{\Bas}}
\newcommand{\Bas}{\mathcal{B}}
\newcommand{\basis}{\mathcal{B}}
\newcommand{\Basval}{\mathcal{B}_\text{val}}
\newcommand{\Val}{\mathcal{V}}
\newcommand{\Cor}{\mathcal{C}}

% operators
\newcommand{\hT}{\Hat{T}}
\newcommand{\hWee}[1]{\Hat{W}_\text{ee}^{#1}}
\newcommand{\f}[2]{f_{#1}^{#2}}
\newcommand{\Gam}[2]{\Gamma_{#1}^{#2}}

% coordinates
\newcommand{\br}[1]{{\mathbf{r}_{#1}}}
\newcommand{\bx}[1]{\mathbf{x}_{#1}}
\newcommand{\dbr}[1]{d\br{#1}}
\newcommand{\PBEspin}{PBEspin}
\newcommand{\PBEueg}{PBE-UEG-{$\tilde{\zeta}$}}

\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LCT}{Laboratoire de Chimie Th\'eorique, Universit\'e Pierre et Marie Curie, Sorbonne Universit\'e, CNRS, Paris, France}

\begin{document}	

\title{Mixing density functional theory and wave function theory for strong correlation: the best of both worlds}

\begin{abstract}
bla bla bla youpi tralala 
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%
The main goal of quantum chemistry is to propose reliable theoretical tools to describe the rich area of chemistry. 
The accurate computation of the electronic structure of molecular systems plays a central role in the development of methods in quantum chemistry, 
but despite intense developments, no definitive solution to that problem have been found. 
The theoretical challenge to be overcome falls back in the category of the quantum many-body problem due the intrinsic quantum nature 
of the electrons and the coulomb repulsion between them, inducing the so-called electronic correlation problem. 
Tackling this problem translate to solving the Schroedinger equation for a $N$~-~electron system, and two roads have emerged to approximate the solution to this formidably complex mathematical problem: the wave function theory (WFT) and density functional theory (DFT). 
Although both WFT and DFT spring from the same problem, their formalisms are very different as the former deals with the complex 
$N$~-~body wave function whereas the latter handles the much simpler one~-~body density. 
The computational cost of DFT is very appealing as in its Kohn-Sham (KS) formulation it can be recast in a mean-field procedure. 
Therefore, although constant efforts are performed to reduce the computational cost of WFT, DFT remains still the workhorse of quantum chemistry. 

From the theoretician point of view, the complexity of description of a given chemical system can be roughly 
categorized by the strength of the electronic correlation appearing in its electronic structure.  
Weakly correlated systems, such as closed-shell organic molecules near their equilibrium geometry, are typically dominated by the avoidance effects when electron are near the electron coalescence point, which are often called short-range correlation effects, 
or far from each other, typically dispersion forces. The theoretical description of weakly correlated systems is one of the more concrete achievement 
of quantum chemistry, and the main remaining issue for these systems is to push the limit in terms of the size of the chemical systems that can be treated.  
The case of the so-called strongly correlated systems, which are ubiquitous in chemistry, is much more problematic as they exhibits 
a much more exotic electronic structure. 
Transition metals containing systems, low-spin open shell systems, covalent bond breaking or excited states 
have all in common that they cannot be even qualitatively described by a single electronic configuration. 
It is now clear that the usual approximations in KS-DFT fails in giving an accurate description of these situations and WFT has become 
the standard for the treatment of strongly correlated systems. 
From the theoretical point of view, the complexity of the strong correlation problem is, at least, two-fold: 
i) the presence of near degeneracies and/or strong interactions among a primary set of electronic configurations 
(the size of which can potentially scale exponentially in some cases) determines the qualitative description of the wave function, 
ii) the quantitative description of the systems must take into account weak correlation effects which requires to take into account many 
other electronic configurations with typically much smaller weights in the wave function.  
Fulfilling these two objectives is a rather complicated task, specially if one adds the requirement of size-extensivity and additivity of the computed energy in the case of non interacting fragments, which is a very desirable property for any approximated method. 

%To tackle this complicated problem, many methods have been proposed and an exhaustive review of the zoology of methods for strong correlation goes beyond the scope and purpose of this article. 

To tackle this problem, many WFT methods have emerged which can be categorized in two branches: the single-reference (SR) 
and multi-reference (MR) methods. 
The SR methods rely on a single electronic configuration as a zeroth-order wave function, typically Hartree-Fock (HF). 
Then the electron correlation is introduced by increasing the rank of multiple hole-particle excitations, 
preferably treated in a coupled-cluster (CC) fashion for the sake of compactness of the wave function and extensivity of the computed energies. 
The advantage of these approaches rely on the rather straightforward way to improve the level of accuracy, 
which consists in increasing the rank of the excitation operators used to generate the CC wave function. 
Despite its appealing elegant simplicity, the computational cost of the CC methods increase drastically with the rank of the excitation 
operators, even if promising alternative approaches have been proposed using stochastic techniques\cite{Thom-PRL-10,ScoTho-JCP-17,SpeNeuVigFraTho-JCP-18,DeuEmiShePie-PRL-17,DeuEmiMagShePie-JCP-18,DeuEmiYumShePie-JCP-19} or symmetry-broken approaches\cite{QiuHenZhaScu-JCP-17,QiuHenZhaScu-JCP-18,GomHenScu-JCP-19}. 
In the MR approaches, the zeroth order wave function consists in a linear combination of Slater determinants which are supposed to concentrate most of strong interactions and near degeneracies inherent in the structure of the Hamiltonian for a strongly correlated system. The usual approach to build such a zeroth-order wave function is to perform a complete active space self consistent field (CASSCF) whose variational property prevent any divergence, and which can provide extensive energies. Of course, the choice of the active space is rather a subtle art and the CASSCF results might strongly depend on the level of chemical/physical knowledge of the user.  
On top of this zeroth-order wave function, weak correlation is introduced by the addition of other configurations through either configuration interaction\cite{WerKno-JCP-88,KnoWer-CPL-88} (MRCI) or perturbation theory (MRPT) and even coupled cluster (MRCC), which have their strengths and weaknesses, 
The advantage of MRCI approaches rely essentially in their simple linear parametrisation for the wave function together with the variational property of their energies, whose inherent drawback is the lack of size extensivity of their energies unless reaching the FCI limit. On the other hand, MRPT and MRCC can provide extensive energies but to the price of rather complicated formalisms, and these approaches might be subject to divergences and/or convergence problems due to the non linearity of the parametrisation for MRCC or a too poor choice of the zeroth-order Hamiltonian. 
A natural alternative is to combine MRCI and MRPT, which falls in the category of selected CI (SCI) which goes back to the late 60's and who has received a revival of interest and applications during the last decade \cite{BenErn-PhysRev-1969,WhiHac-JCP-1969,HurMalRan-1973,EvaDauMal-ChemPhys-83,Cim-JCP-1985,Cim-JCC-1987,IllRubRic-JCP-88,PovRubIll-TCA-92,BunCarRam-JCP-06,AbrSheDav-CPL-05,MusEngels-JCC-06,BytRue-CP-09,GinSceCaf-CJC-13,CafGinScemRam-JCTC-14,GinSceCaf-JCP-15,CafAplGinScem-arxiv-16,CafAplGinSce-JCP-16,SchEva-JCP-16,LiuHofJCTC-16,HolUmrSha-JCP-17,ShaHolJeaAlaUmr-JCTC-17,HolUmrSha-JCP-17,SchEva-JCTC-17,PerCle-JCP-17,OhtJun-JCP-17,Zim-JCP-17,LiOttHolShaUmr-JCP-2018,ChiHolOttUmrShaZim-JPCA-18,SceBenJacCafLoo-JCP-18,LooSceBloGarCafJac-JCTC-18,GarSceGinCaffLoo-JCP-18,SceGarCafLoo-JCTC-18,GarGinMalSce-JCP-16,LooBogSceCafJac-JCTC-19}, and among which the CI perturbatively selected iteratively (CIPSI) can be considered as a pioneer. The main idea of the CIPSI and other related SCI algorithms is to iteratively select the most important Slater determinants thanks to perturbation theory in order to build a MRCI zeroth-order wave function which automatically concentrate the strongly interacting part of the wave function. On top of this MRCI zeroth-order wave function, a rather simple MRPT approach is used to recover the missing weak correlation and the process is iterated until reaching a given stopping criterion. It is important to notice that in the SCI algorithms, neither the SCI or the MRPT are size extensive \text{per se}, but the extensivity property is almost recovered by approaching the FCI limit. 
When the SCI are affordable, their clear advantage are that they provide near FCI wave functions and energies, whatever the level of knowledge of the user on the specific physical/chemical problem considered. The drawback of SCI is certainly their \textit{intrinsic} exponential scaling due to their linear parametrisation. Nevertheless, such an exponential scaling is lowered by the smart selection of the zeroth-order wave function together with the MRPT calculation. 

Besides the difficulties of accurately describing the electronic structure within a given basis set, a crucial component of the limitations of applicability of WFT concerns the slow convergence of the energies and properties with respect to the quality of the basis set. As initially shown by the seminal work of Hylleraas\cite{Hyl-ZP-29} and further developed by Kutzelnigg \textit{et. al.}\cite{Kut-TCA-85,KutKlo-JCP-91, NogKut-JCP-94}, the main convergence problem originates from the divergence of the coulomb interaction at the electron coalescence point, which induces a discontinuity in the first-derivative of the wave function (the so-called electron-electron cusp). Describing such a discontinuity with an incomplete basis set is impossible and as a consequence, the convergence of the computed energies and properties can be strongly affected. To attenuate this problem, extrapolation techniques has been developed, either based on the Hylleraas's expansion of the coulomb operator\cite{HalHelJorKloKocOlsWil-CPL-98}, or more recently based on perturbative arguments\cite{IrmHulGru-arxiv-19}. A more rigorous approach to tackle the basis set convergence problem has been proposed by the so-called R12 and F12 methods\cite{Ten-TCA-12,TenNog-WIREs-12,HatKloKohTew-CR-12, KonBisVal-CR-12, GruHirOhnTen-JCP-17, MaWer-WIREs-18} which introduce a function explicitly depending on the interelectronic coordinates ensuring the correct cusp condition in the wave function, and the resulting correlation energies converge much faster than the usual WFT. For instance, using the explicitly correlated version of coupled cluster with single, double and perturbative triple substitution (CCSD(T)) in a triple-$\zeta$ quality basis set is equivalent to a quintuple-$\zeta$ quality of the usual CCSD(T) method\cite{TewKloNeiHat-PCCP-07}, although inherent computational overhead are introduced by the auxiliary basis sets needed to resolve the rather complex three- and four-electron integrals involved in the F12 theory. 

An alternative point of view is to leave the short-range correlation effects to DFT and to use WFT to deal only with the long-range and/or strong-correlation effects. A rigorous approach to do so is the range-separated DFT (RSDFT) formalism (see Ref.~\onlinecite{TouColSav-PRA-04} and references therein) which rely on a splitting of the coulomb interaction in terms of the interelectronic distance thanks to a range-separation parameter $\mu$. The advantage of such approach is at least two-folds: i) the DFT part deals only with the short-range part of the coulomb interaction, and therefore the usual semi-local approximations to the unknown exchange-correlation functional are more suited to that correlation regime, ii) the WFT part deals with a smooth non divergent interaction, which removes the cusp condition and therefore the basis set convergence is much faster\cite{FraMusLupTou-JCP-15}. 
Therefore, a number of approximate RS-DFT schemes have been developed within single-reference \cite{AngGerSavTou-PRA-05, GolWerSto-PCCP-05, TouGerJanSavAng-PRL-09,JanHenScu-JCP-09, TouZhuSavJanAng-JCP-11, MusReiAngTou-JCP-15} or multi-reference \cite{LeiStoWerSav-CPL-97, FroTouJen-JCP-07, FroCimJen-PRA-10, HedKneKieJenRei-JCP-15, HedTouJen-JCP-18, FerGinTou-JCP-18} WFT approaches. Nevertheless, there are still some open issues in RSDFT, such as the dependence of the quality of the results on the value of the range separation $\mu$ which can be seen as an empirical parameter, together with the remaining self-interaction errors.  

Following this path, a very recent solution to the basis set convergence problem has been proposed by some of the preset authors\cite{GinPraFerAssSavTou-JCP-18} where they proposed to use RSDFT to take into account only the correlation effects outside a given basis set. The key idea in such a work is to realize that as a wave function developed in an incomplete basis set is cusp-less, it could also come from a Hamiltonian with a non divergent electron-electron interaction. Therefore, the authors proposed a mapping with RSDFT through the introduction of an effective non-divergent interaction representing the usual coulomb interaction projected in an incomplete basis set. First applications to weakly correlated molecular systems have been successfully carried recently\cite{LooPraSceTouGin-JCPL-19} together with the first attempt to generalize this approach to excited states\cite{exicted}. 
The goal of the present work is to push the development of this new theory toward the description of strongly correlated systems. 
The paper is organized as follows: in section \ref{sec:theory} we recall the mathematical framework of the basis set correction and we propose a practical extension for strongly correlated systems. Two key aspect are discussed: the extensivity of the correlation energies together with the $S_z$ independence of the results. 
Then in section \ref{sec:results} we discuss the potential energy surfaces (PES) of N$_2$, F$_2$ and H$_{10}$ up to full dissociation as a prototype of strongly correlated problems. Finally, we conclude in section \ref{sec:conclusion}

%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theory}
%%%%%%%%%%%%%%%%%%%%%%%%
The theoretical framework of the basis set correction has been derived in details in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18}, so we recall briefly the main equations involved for the present study. 
First in section \ref{sec:basic} we recall the basic mathematical framework of the present theory by introducing the density functional complementary to a basis set $\Bas$. Then in section \ref{sec:wee} we introduce an effective non divergent interaction in a basis set $\Bas$, which leads us to the definition of an effective range separation parameter varying in space in section \ref{sec:mur}. Thanks to the range separation parameter, we make a mapping with a specific class of RSDFT functionals and propose practical approximations for the unknown density functional complementary to a basis set $\Bas$, for which new approximations for the strong correlation regime are given in section \ref{sec:functional}. 
\subsection{Basic formal equations}
\label{sec:basic}
The exact ground state energy $E_0$ of a $N-$electron system can be obtained by an elegant mathematical framework connecting WFT and DFT, that is the Levy-Lieb constrained search formalism which reads  
\begin{equation}
 \label{eq:levy}
 E_0 = \min_{\denr} \bigg\{ F[\denr] + (v_{\text{ne}} (\br{}) |\denr) \bigg\},
\end{equation}
where $(v_{ne}(\br{})|\denr)$ is the nuclei-electron interaction for a given density $\denr$ and $F[\denr]$ is the so-called Levy-Liev universal density functional 
\begin{equation}
 \label{eq:levy_func}
 F[\denr] = \min_{\Psi \rightarrow \denr} \elemm{\Psi}{\kinop +\weeop }{\Psi}.
\end{equation}
The minimizing density $n_0$ of equation \eqref{eq:levy} is the exact ground state density. 
Nevertheless, in practical calculations the minimization is performed over the set $\setdenbasis$ which are the densities representable in a basis set $\Bas$, we assume from thereon that the densities used in the equations belong to $\setdenbasis$. 

In the present context it is important to notice that in order to recover the \textit{exact} ground state energy, the wave functions $\Psi$ involved in the definition of eq. \eqref{eq:levy_func} must be developed in a complete basis set. 
An important step proposed originally by some of the present authors in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18} 
was to propose to split the minimization in the definition of $F[\denr]$  using $\wf{}{\Bas}$  which are wave functions developed in $\basis$ 
\begin{equation}
 \label{eq:def_levy_bas}
 F[\denr] = \min_{\wf{}{\Bas} \rightarrow \denr} \elemm{\wf{}{\Bas}}{\kinop +\weeop}{\wf{}{\Bas}} + \efuncden{\denr}, 
\end{equation}
which leads to the following definition of $\efuncden{\denr}$ which is the the density functional complementary to the basis set $\Bas$    
\begin{equation}
 \begin{aligned}
 \efuncden{\denr} =& \min_{\Psi \rightarrow \denr} \elemm{\Psi}{\kinop +\weeop }{\Psi} \\ 
                  &- \min_{\Psi^{\Bas} \rightarrow \denr} \elemm{\wf{}{\Bas}}{\kinop +\weeop}{\wf{}{\Bas}}. 
 \end{aligned}
\end{equation}
Therefore thanks to eq. \eqref{eq:def_levy_bas} one can properly connect the DFT formalism with the basis set error in WFT calculations. In other terms, the existence of $\efuncden{\denr}$ means that the correlation effects not taken into account in $\basis$ can be formulated as a density functional. 

Assuming that the density $\denFCI$ associated to the ground state FCI wave function $\psifci$ is a good approximation of the exact density, one obtains the following approximation for the exact ground state density (see equations 12-15 of Ref. \onlinecite{GinPraFerAssSavTou-JCP-18})
\begin{equation}
 \label{eq:e0approx}
 E_0 = \efci + \efuncbasisFCI
\end{equation}
where $\efci$ is the ground state FCI energy within $\Bas$. As it was originally shown in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18} and further emphasized in Ref. \onlinecite{LooPraSceTouGin-JCPL-19,GinSceTouLoo-JCP-19}, the main role of $\efuncbasisFCI$ is to correct for the basis set incompleteness errors, a large part of which originates from the lack of cusp in any wave function developed in an incomplete basis set. 
The whole purpose of this paper is to determine approximations for $\efuncbasisFCI$ which are suited for treating strong correlation regimes. The two requirement for such conditions are that i) it can be defined for multi-reference wave functions, ii) it must provide size extensive energies, iii) it is invariant of the $S_z$ component of a given spin multiplicity. 

\subsection{Definition of an effective interaction within $\Bas$}
\label{sec:wee}
As it was originally shown by Kato\cite{kato}, the cusp in the exact wave function originates from the divergence of the coulomb interaction at the coalescence point. Therefore, a cusp less wave function  $\wf{}{\Bas}$ could also be obtained from a Hamiltonian with a non divergent electron-electron interaction. In other words, the incompleteness of a finite basis set can be understood as the removal of the divergence of the usual coulomb interaction at the electron coalescence point. 

As it was originally derived in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18} (see section D and annexes), one can obtain an effective non divergent interaction, here referred as $\wbasis$, which reproduces the expectation value of the coulomb operator over a given wave function $\wf{}{\Bas}$. As we are interested in the behaviour at the coalescence point, we focus on the opposite spin part of the electron-electron interaction. 

More specifically, we define the effective interaction associated to a given wave function $\wf{}{\Bas}$ as
\begin{equation}
 \label{eq:wbasis}
 \wbasis =
    \begin{cases}
      \fbasis /\twodmrdiagpsi,    & \text{if $\twodmrdiagpsi \ne 0$,}
\\
       \infty,                                                                                          & \text{otherwise,}
    \end{cases}
\end{equation}
where $\twodmrdiagpsi$ is the opposite spin two-body density associated to $\wf{}{\Bas}$ 
\begin{equation}
 \twodmrdiagpsi = \sum_{pqrs} \SO{p}{1} \SO{q}{2} \Gam{pq}{rs} \SO{r}{1} \SO{s}{2}, 
\end{equation}
$\Gam{pq}{rs} = 2 \mel*{\wf{}{\Bas}}{ \aic{r_\downarrow}\aic{s_\uparrow}\ai{q_\uparrow}\ai{p_\downarrow}}{\wf{}{\Bas}}$ its associated two-body tensor, $\SO{p}{}$ are the spatial orthonormal orbitals,  
\begin{equation}
        \label{eq:fbasis}
        \fbasis
        = \sum_{pqrstu \in \Bas}  \SO{p}{1} \SO{q}{2} \V{pq}{rs} \Gam{rs}{tu} \SO{t}{1} \SO{u}{2},
\end{equation}
and $\V{pq}{rs}=\langle pq | rs \rangle$ are the usual two-electron Coulomb integrals.
With such a definition, one can show that $\wbasis$ satisfies 
\begin{equation}
 \int \int \dr{1} \dr{2} \wbasis \twodmrdiagpsi = \int \int \dr{1} \dr{2} \frac{\twodmrdiagpsi}{|\br{1}-\br{2}|}. 
\end{equation}
As it was shown in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18}, the effective interaction $\wbasis$ is necessary finite at coalescence for an incomplete basis set, and tends to the regular coulomb interaction in the limit of a complete basis set for any choice of wave function $\psibasis$, that is 
\begin{equation}
 \label{eq:cbs_wbasis}
 \lim_{\Bas \rightarrow \text{CBS}} \wbasis = \frac{1}{|\br{1}-\br{2}|}\quad \forall\,\psibasis. 
\end{equation}
The condition of equation \eqref{eq:cbs_wbasis} is fundamental as it guarantees the good behaviour of all the theory in the limit of a complete basis set. 

\subsection{Definition of a range-separation parameter varying in real space}
\label{sec:mur}
As the effective interaction within a basis set $\wbasis$ is non divergent, one can fit such a function with a long-range interaction defined in the framework of RSDFT which depends on the range-separation parameter $\mu$
\begin{equation}
 \label{eq:weelr}
 w_{ee}^{\lr}(\mu;r_{12}) = \frac{\text{erf}\big(\mu \,r_{12} \big)}{r_{12}}. 
\end{equation}
As originally proposed in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18}, we introduce a range-separation parameter $\murpsi$ varying in real space
\begin{equation}
 \label{eq:def_mur}
 \murpsi = \frac{\sqrt{\pi}}{2} \wbasiscoal 
\end{equation}
such that 
\begin{equation}
 w_{ee}^{\lr}(\murpsi;0) = \wbasiscoal \quad \forall \, \br{}. 
\end{equation}
Because of the very definition of $\wbasis$, one has the following properties at the CBS limit (see \eqref{eq:cbs_wbasis}) 
\begin{equation}
 \label{eq:cbs_mu}
  \lim_{\Bas \rightarrow \text{CBS}} \murpsi = \infty\quad \forall \,\psibasis, 
\end{equation}
which is fundamental to guarantee the good behaviour of the theory at the CBS limit. 

\subsection{Generic form and properties of the approximations for $\efuncden{\denr}$ }
\subsubsection{Generic form of the approximated functionals}
As originally proposed and motivated in Ref. \onlinecite{GinPraFerAssSavTou-JCP-18}, we approximate the complementary basis set functional $\efuncden{\denr}$ by using the so-called multi-determinant correlation functional (ECMD) introduced by Toulouse and co-workers\cite{TouGorSav-TCA-05}. 
Following the recent work of some of the present authors\cite{LooPraSceTouGin-JCPL-19}, we propose to use a PBE-like functional which uses the total density $\denr$, the spin polarisation $\zeta(\br{}) = n_{\alpha}(\br{}) - n_{\beta}(\br{})$, reduced density gradient $s(\br{}) = \nabla \denr/\denr^{4/3}$ and the on-top pair density $\ntwo(\br{})$. In the present work, all the density-related quantities are computed with the same wave function $\psibasis$ used to define $\murpsi$. 
Therefore, a given approximation X of $\efuncden{\denr}$ have the following generic form
\begin{equation}
 \begin{aligned}
 \label{eq:def_ecmdpbebasis}
 \efuncdenpbe{\argebasis} =  &\int  d\br{} \,\denr  \\ & \ecmd(\argrebasis)
 \end{aligned}
\end{equation}
where $\ecmd(\argecmd)$ is the ECMD correlation energy density defined as
\begin{equation}
 \label{eq:def_ecmdpbe}
 \ecmd(\argecmd) = \frac{\varepsilon_{\text{c,PBE}}(\argepbe)}{1+ \mu^3 \beta(\argepbe)}
\end{equation}
with 
\begin{equation}
 \label{eq:def_beta}
 \beta(\argebasis) = \frac{3}{2\sqrt{\pi}(1 - \sqrt{2})}\frac{\varepsilon_{\text{c,PBE}}(\argepbe)}{\ntwo/\den}, 
\end{equation}
and where $\varepsilon_{\text{c,PBE}}(\argepbe)$ is the usual PBE correlation energy density\cite{PerBurErn-PRL-96}. Before introducing the different flavour of approximated functionals that we will use here (see  \ref{sec:def_func}), we would like to give some motivations based on physical requirements for the such a choice of functional form. 

The actual functional form of $\ecmd(\argecmd)$ have been originally proposed by some of the present authors in the context of RSDFT~\cite{FerGinTou-JCP-18} in order to fulfill the two following limits
\begin{equation}
   \lim_{\mu \rightarrow 0} \ecmd(\argecmd) = \varepsilon_{\text{c,PBE}}(\argepbe),   
\end{equation}
which can be qualified as the weak correlation regime, and the large $\mu$ limit
\begin{equation}
 \label{eq:lim_mularge}
    \ecmd(\argecmd) = \frac{1}{\mu^3} \ntwo + o(\frac{1}{\mu^5}),
\end{equation}
which, as it was previously shown\cite{TouColSav-PRA-04, GoriSav-PRA-06,PazMorGorBac-PRB-06} by various authors, is the exact expression for the ECMD in the limit of large $\mu$, provided that $\ntwo$ is the \textit{exact} on-top pair density of the system. 
In the context of RSDFT, some of the present authors have illustrated in Ref.~\onlinecite{FerGinTou-JCP-18} that the on-top pair density involved in eq. \eqref{eq:def_ecmdpbe} plays a crucial role when reaching the strong correlation regime. The importance of the on-top pair density in the strong correlation regime have been also acknowledged by Pernal and co-workers\cite{GritMeePer-PRA-18} and Gagliardi and co-workers\cite{CarTruGag-JPCA-17}.  
Also, $\ecmd(\argecmd) $ vanishes when $\ntwo$ vanishes 
\begin{equation}
 \label{eq:lim_n2}
  \lim_{\ntwo \rightarrow 0} \ecmd(\argecmd) = 0
\end{equation}
which is exact for systems with a vanishing on-top pair density, such as the totally dissociated H$_2$ which is the archetype of strongly correlated systems.  
Of course, as all RSDFT functionals the function $\ecmd(\argecmd)$ vanishes when $\mu \rightarrow \infty$ 
\begin{equation}
 \label{eq:lim_muinf}
  \lim_{\mu \rightarrow \infty} \ecmd(\argecmd) = 0. 
\end{equation}

\subsubsection{Properties of approximated functionals}
Within the definition of \eqref{eq:def_mur} and \eqref{eq:def_ecmdpbebasis}, any approximated complementary basis set functionals $\efuncdenpbe{\argecmd}$ satisfies two important properties. 
Because of the properties \eqref{eq:cbs_mu} and \eqref{eq:lim_muinf}, $\efuncdenpbe{\argecmd}$ vanishes when reaching the complete basis set limit, whatever the wave function $\psibasis$ used to define the range separation parameter $\mu_{\Psi^{\basis}}$: 
\begin{equation}
 \label{eq:lim_ebasis}
 \lim_{\basis \rightarrow \text{CBS}} \efuncdenpbe{\argecmd} = 0\quad \forall\, \psibasis, 
\end{equation}
which guarantees an unaltered limit when reaching the CBS limit. 
Also, the $\efuncdenpbe{\argecmd}$ vanishes for systems with vanishing on-top pair density, which guarantees the good limit in the case of stretched H$_2$ and for one-electron system. 
Such a property is guaranteed independently by i) the definition of the effective interaction  $\wbasis$ (see equation \eqref{eq:wbasis}) together with the condition \eqref{eq:lim_muinf}, ii) the fact that the $\ecmd(\argecmd)$ vanishes when the on-top pair density vanishes (see equation \eqref{eq:lim_n2}). 

\subsection{Requirements for the approximated functionals in the strong correlation regime}
\subsubsection{Requirements: separability of the energies and $S_z$ invariance}
An important requirement for any electronic structure method is the extensivity of the energy, \textit{i. e.} the additivity of the energies in the case of non interacting fragments, which is particularly important to avoid any ambiguity in computing interaction energies. 
When two subsystems $A$ and $B$ dissociate in closed shell systems, as in the case of weak interactions for instance, a simple HF wave function leads to extensive energies. 
When the two subsystems dissociate in open shell systems, such as in covalent bond breaking, it is well known that the HF approach fail and an alternative is to use a CASSCF wave function which, provided that the active space has been properly chosen, leads to additives energies.  
Another important requirement is the independence of the energy with respect to the $S_z$ component of a given spin state.  
Such a property is also important in the context of covalent bond breaking where the ground state of the super system $A+B$ is in general of low spin while the ground states of the fragments $A$ and $B$ are in high spin which can have multiple $S_z$ components.

\subsubsection{Condition for the functional $\efuncdenpbe{\argebasis}$ to obtain $S_z$ invariance}
A sufficient condition to achieve $S_z$ invariance is to eliminate all dependency to $S_z$, which in the case of $\ecmd(\argecmd)$ is the spin polarisation $\zeta(\br{})$ involved in the correlation energy density $\varepsilon_{\text{c,PBE}}(\argepbe)$ (see equation \eqref{eq:def_ecmdpbe}). 
As originally shown by Perdew and co-workers\cite{PerSavBur-PRA-95}, the dependence on the spin polarisation in the KS-DFT framework can be removed by the rewriting the spin polarisation of a single Slater determinant with only the on-top pair density and the total density. In other terms, the spin density dependence usually introduced in the correlation functionals of KS-DFT tries to mimic the effect of the on-top pair density. 
Based on this reasoning, a similar approach has been used in the context of multi configurational DFT in order to remove the $S_z$ dependency. 
In practice, these approaches introduce the effective spin polarisation 
\begin{equation}
 \label{eq:def_effspin}
 \tilde{\zeta}(n,\ntwo_{\psibasis}) =  
%   \begin{cases}
                       \sqrt{ n^2 - 4 \ntwo_{\psibasis} }
%                              0              & \text{otherwise.}
%   \end{cases}
\end{equation}
which uses the on-top pair density $\ntwo_{\psibasis}$ of a given wave function $\psibasis$. 

The advantages of this approach are at least two folds: i) the effective spin  polarisation $\tilde{\zeta}$ is $S_z$ invariant, ii) it introduces an indirect dependency on the on-top pair density of the wave function $\psibasis$ which usually improves the treatment of strong correlation. 
Nevertheless, the use of $\tilde{\zeta}$ presents several disadvantages as it can become complex when $n^2 - 4 \ntwo_{\psibasis}<0$ and also 
the formula of equation \eqref{eq:def_effspin} is exact only when the density $n$ and on-top pair density $\ntwo_{\psibasis}$ are obtained from a single determinant\cite{PerSavBur-PRA-95}, but it is applied to multi configurational wave functions.  

An alternative to eliminate the $S_z$ dependency would be to simply set $\zeta(\br{})=0$, but this would lower the accuracy of the usual correlation functional, such as the PBE correlation functional used here $\varepsilon_{\text{c,PBE}}(\argepbe)$. Nevertheless, as the spin polarisation usually tries to mimic the on-top pair density and the function $\ecmd(\argecmd)$ explicitly depends on the on-top pair density (see equations \eqref{eq:def_ecmdpbe} and \eqref{eq:def_beta}), we propose here to use the $\ecmd$ functional with \textit{a zero spin polarisation}. This ensures a $S_z$ invariance and, as will be numerically shown, very weakly affect the accuracy of the functional. 

\subsubsection{Conditions on $\psibasis$ for the extensivity}
In the case of the present basis set correction, as $\efuncdenpbe{\argebasis}$ is an integral over $\mathbb{R}^3$ of local quantities, in the case of non overlapping fragments $A\ldots B$ it can be written as the sum of two local contributions: one coming from the integration over the region of the sub-system $A$ and the other one from the region of the sub-system $B$. 
Therefore, a sufficient condition for the extensivity is that these quantities coincide in the isolated systems and in the subsystem of the super system $A\ldots B$. 
As $\efuncdenpbe{\argebasis}$ depends only on quantities which are properties of the wave function $\psibasis$, a sufficient condition for the extensivity of these quantities is that the function factorise in the limit of non-interacting fragments, that is $\Psi_{A\ldots B}^{\basis} = \Psi_A^{\basis} \Psi_B^{\basis}$. 
In the case where the two subsystems $A$ and $B$ dissociate in closed shell systems, a simple HF wave function ensures this property, but when one or several covalent bonds are broken, the use of a properly chosen CASSCF wave function is sufficient to recover this property, as will be numerically illustrated in section \ref{sec:separability}. 
The condition for the active space involved in the CASSCF wave function is that it has to lead to extensive energies in the limit of dissociated fragments. 


\subsection{Different types of approximations for the functional}
\subsubsection{Definition of the protocol to design functionals}
\label{sec:def_func}
As the present work proposes to investigate how different physical quantities impact the description of correlation, we propose here a general protocol and a corresponding nomenclature in order to make things as clear as possible. 
%
Here we propose to investigate the dependency of the functionals $\efuncdenpbe{\argebasis}$ on: i) the wave function $\psibasis$ used to determine the $\murpsi$ and the various density related quantities, ii) the flavour of on-top pair density used, iii) the type of spin polarisation used. 
Therefore, we propose to use the following notations: PBE-"on-top"-"spin polarisation"-$\psibasis$. 

Regarding the spin polarisation that enters into $\varepsilon_{\text{c,PBE}}(\argepbe)$, we will use three different types of formula: i) the usual spin polarisation $\zeta = n_{\alpha} - n_{\beta}$ which \textit{is not} $S_z$ invariant, ii) $\tilde{\zeta}$ defined in equation \eqref{eq:def_effspin} which \textit{is} $S_z$ invariant, and iii) a \textit{zero} spin polarization which of course \textit{is} $S_z$ invariant. 

For the wave function $\psibasis$, we will use either i) a simple RHF/ROHF wave function, ii) a minimal CASSCF leading to additive energies in the case of dissociated covalent bonds. 

Regarding the approximation to the \textit{exact} on-top pair density, we use two different approximations. The first one is based on the uniform electron gas (UEG) and reads 
\begin{equation}
 \label{eq:def_n2ueg}
 \ntwo_{\text{UEG}}(n,\zeta,\br{}) = n(\br{})^2\big(1-\zeta(\br{})\big)g_0\big(n(\br{})\big)
\end{equation}
where the pair-distribution function $g_0(n)$ is taken from equation (46) of Ref. \onlinecite{GorSav-PRA-06}. The approximation of equation \eqref{eq:def_n2ueg} depends on the density and some spin polarisation. Notice that, when using a CASSCF wave function and $\tilde{\zeta}$ as spin polarization, the $\ntwo_{\text{UEG}}$ will depend indirectly on the on-top pair density as $\tilde{\zeta}$ depends on the on-top pair density. 

Another approach consists in taking advantage of the on-top pair density of the wave function $\psibasis$. Following the work of some of the previous authors\cite{FerGinTou-JCP-18,GinSceTouLoo-JCP-19} we introduce the extrapolated on-top pair density $\ntwoextrap$ as
\begin{equation}
 \ntwoextrap(\ntwo_{\psibasis},\mu,\br{}) = \ntwo_{\wf{}{\Bas}}(\br{}) \bigg( 1 + \frac{2}{\sqrt{\pi}\murpsi} \bigg)^{-1}
\end{equation}
which directly follows from the large-$\mu$ extrapolation of the exact on-top pair density proposed by Gori-Giorgi and Savin\cite{GorSav-PRA-06}. 
When using $\ntwoextrap(\ntwo,\mu,\br{})$ in a functional, we will refer simply refer it as "ont". 

\subsubsection{Definition of a hierarchy of functionals}
Within the convention proposed in the section \ref{sec:def_func}, the PBE-UEG-$\zeta$-HF is the functional which was introduced in Ref. \onlinecite{LooPraSceTouGin-JCPL-19} and which reads
\begin{equation}
 \label{eq:def_pbeueg}
 \begin{aligned}
 \pbeuegxihf &\equiv  \int  d\br{} \,\denr  \\ & \ecmd\big(\argrpbeuegxihf\big)
 \end{aligned}
\end{equation}
Therefore, such a functional uses a HF wave function to define; i) the $\murpsi$, ii) the total density, reduced density gradients, usual spin polarisation $\zeta$ and uses the UEG-like on-top pair density with the usual spin polarisation $\zeta$. 
Of course, because of the use of an HF wave function as $\psibasis$, the density related quantities are extensive only in the case of dissociation in closed shell system. Also, one can notice that changing the spin polarisation from $\zeta$ to $\tilde{\zeta}$ does not change the results as by definition, $\tilde{\zeta} = \zeta$ for a single Slater determinant.  

By changing the definition of $\psibasis=\text{HF}$ to $\psibasis=\text{CASSCF}$ on obtains the PBE-UEG-$\zeta$-CAS which reads 
\begin{equation}
 \label{eq:def_pbeueg}
 \begin{aligned}
 \pbeuegxicas &\equiv  \int  d\br{} \,\denr  \\ & \ecmd\big(\argrpbeuegxicas\big)
 \end{aligned}
\end{equation}
where the density, reduced density gradients, usual spin polarisation and UEG on-top pair density are computed from a CASSCF wave function. Therefore, the $\murpsi$, density, reduced density gradient are extensive in the case of dissociated covalent bonding. Nevertheless, the use of the regular spin polarisation $\zeta$ leads to non $S_z$ invariance. 

One can change the spin polarisation to the effective spin polarisation $\tilde{\zeta}$ to obtain the PBE-UEG-$\tilde{\zeta}$-CAS which is $S_z$ invariant, and therefore this functional will reads to 
\begin{equation}
 \label{eq:def_pbeueg}
 \begin{aligned}
 \pbeuegXicas =  &\int  d\br{} \,\denr  \\ & \ecmd(\argrpbeuegXicas). 
 \end{aligned}
\end{equation}
One can also change the flavour of the on-top pair density by taking advantage of the on-top pair density $\ntwo_{\wf{}{\Bas}}(\br{})$ computed with $\psibasis$. 
Therefore, one can define the PBE-ONT-$\zeta$-CAS as 
\begin{equation}
 \label{eq:def_pbeueg}
 \begin{aligned}
 \pbeontXicas =  &\int  d\br{} \,\denr  \\ & \ecmd(\argrpbeontXicas). 
 \end{aligned}
\end{equation}
Such a functional can be further improved by using the $S_z$ invariant effective spin polarisation $\tilde{\zeta}$ to give the PBE-ONT-$\tilde{\zeta}$-CAS. 



%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}
\subsection{Numerical tests of extensivity}
The first numerical results investigated are the numerical tests of extensivity of the various functionals. 
As mentioned before, when considering a super system $A+B$ dissociating into non interacting fragments $A\ldots B$, there are two different situations regarding the extensivity of the energy: when the subsystems $A$ and $B$ dissociate in closed or open shell systems. 
Therefore, we shall consider two systems $A$ and $B$ and compare the sum of the energies obtained with the super system $A\ldots B$ in the limit of non interactive fragments. The error to additivity for a given method $Y$ is therefore defined as $E_Y(A) + E_Y(B) - E_Y(A\ldots B)$. 

\subsubsection{Dissociation to closed shell ground states}
We begin our study by giving numerical evidence for the extensivity of the present basis set correction for systems dissociating in closed shell systems. 
In these cases, the use a HF wave function is sufficient to guarantee the extensivity of the basis set correction, and therefore we use the simple $\pbeuegxihf$ functional. The system under study is $A=\text{F}_2$ at experimental equilibrium geometry (F-F=1.411 angstroms) and $B=\text{Ne}$. 
We report in table \ref{tab:extensiv_closed} the error to additivity for the HF energy and for $\pbeuegxihf$ using the aug-cc-pvdz basis set and using a He core to define the $\mu_{\text{HF}}^{\basis}(\br{})$ and the frozen core densities. 
The numbers in table \ref{tab:extensiv_closed} clearly show that when HF energies are additive, the $\pbeuegxihf$ is also additive. 
Also, the error to additivity using the usual spin polarisation $\zeta$ and the extrapolated on-top pair density are much lowered compared to that using UEG on-top pair density, highlighting the important role played by the on-top pair density of the CASSCF wave function. 

\begin{table*}
\caption{Total energies (in Hartree) for HF and $E$ in aug-cc-pvdz for the Ne atom, F$_2$ (with F-F=1.411 angstroms) and the super non interacting system Ne--F$_2$.   }
\begin{tabular}{lcc}
%\hline
   System            &  HF                     &   $\pbeuegxihf$          \\
\hline                                                                
Ne                   & -128.4963497306184      &   -0.1039022285466806    \\
F$_2$                & -198.698792752661       &   -0.1596345827582842    \\
Ne $\ldots$ F$_2$    & -201.554497420371       &   -0.2635368113049532    \\
 \hline                                                                
Error to additivity  & 3.4   $\times 10^{-13}$ & 1.1 $\times 10^{-14}$               \\
\end{tabular}
\label{tab:extensiv_closed}
\end{table*}

\subsubsection{Dissociation to open shell ground states}
The system studied to investigate the extensivity in the case of dissociation to open shell systems is the completely dissociated N$_2$ molecule which imply the breaking of three covalent bonds. 
As the HF wave function does not lead to extensive energy, it is clear that it cannot be used as $\psibasis$ and therefore for N$_2$ we use a minimal valence CASSCF(6,6) involving the three bonding orbitals ($\sigma$, $\pi_x$, $\pi_y$) and corresponding anti-bonding orbitals and a ROHF wave function for the N atom. 
The numerical results for the extensivity of the various flavours of functionals are given in table \ref{tab:extensiv_open}. From these numbers, one can clearly notice that only the functionals using the effective spin polarisation $\tilde{\zeta}$ are size extensive, whatever the type of on-top pair density used. 

\begin{table*}
\caption{Total energies (in Hartree) for N$_2$ in the aug-cc-pvdz basis set.  }
\begin{tabular}{lccccc}
%\hline
   System            &  ROHF/CASSCF(6,6)       &   $\pbeuegxicas$           &  $\pbeuegXicas$        &   $\pbeontxicas$     &  $\pbeontXicas$        \\
\hline                                                                                                                            
N                    & -128.496349730618       &   -0.0230740500348705      & -0.0230740500348705    &  -0.0247392466968251 &  -0.0247392466968251     \\
N$\ldots$N           & -198.698792752661       &   -0.0691133629633014      & -0.0461481000697329    &  -0.0509457188492165 &  -0.0494784933936403     \\
 \hline                                                                                                                                               
Error to additivity  & 1.0   $\times 10^{-13}$ &    0.02296                 & 8.0 $\times 10^{-15}$  &   0.0015             &   9.9 $\times 10^{-15}$  \\
\end{tabular}
\label{tab:extensiv_open}
\end{table*}

\label{sec:results}
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
        \includegraphics[width=\linewidth]{data/N2/DFT_avdzE_relat.pdf}
        \includegraphics[width=\linewidth]{data/N2/DFT_avdzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{data/N2/DFT_avdzE_error.pdf}
        \caption{
        N$_2$, aug-cc-pvdz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:N2_avdz}}
\end{figure}

\begin{figure}
        \includegraphics[width=\linewidth]{data/N2/DFT_avtzE_relat.pdf}
        \includegraphics[width=\linewidth]{data/N2/DFT_avtzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{data/N2/DFT_avtzE_error.pdf}\\
%       \includegraphics[width=\linewidth]{fig2c}
        \caption{
        N$_2$, aug-cc-pvtz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:N2_avtz}}
\end{figure}

\begin{table*}
\caption{Dissociation energy ($D_0$) computed at different levels in various basis sets.   }
\begin{ruledtabular}
\begin{tabular}{lcccc}

%\hline
   System/basis             &  FCI                    &   FCI+$\pbeuegXicas$   &   FCI+$\pbeontXicas$   &     FCI+$\pbeontnscas$    \\
\hline                                                                                                                 
N$_2$, aug-cc-pvdz          &   321.4$/ $42.8            &   355.6$/$8.6            &   355.0$/$9.2            &     354.0$/$10.2            \\
N$_2$, aug-cc-pvtz          &   347.8$/$16.4            &   361.0$/$3.2            &   362.7$/$1.5            &     362.4$/$1.8             \\
 \hline                                                                
 & \multicolumn{4}{c}{Estimated exact} \\
 & \multicolumn{4}{c}{364.2 } \\
\hline                                                                                                                 
F$_2$, aug-cc-pvdz          &   49.2$/$11.5             &    54.1$/$6.6            &   54.3$/$6.4             &     54.1$/$6.7              \\
F$_2$, aug-cc-pvtz          &   58.6$/$2.1              &    60.6$/$0.1            &   60.9$/$-0.2            &     60.9$/$-0.2             \\
 \hline                                                                
 & \multicolumn{4}{c}{Estimated exact} \\
 & \multicolumn{4}{c}{60.7 } \\
\hline                                                                                                                 
H$_{10}$, cc-pvdz           &  622.1$/$43.3             &   642.6$/$22.8           &   649.2$/$16.2           &     649.5$/$15.9            \\
H$_{10}$, cc-pvtz           &  655.2$/$10.2             &   661.9$/$3.5            &   666.0$/$-0.6           &     666.0$/$-0.6            \\
H$_{10}$, cc-pvqz           &  661.2$/$4.2              &   664.1$/$1.3            &   666.4$/$-1.0           &     666.5$/$-1.1            \\
 \hline                                                                
 & \multicolumn{4}{c}{Estimated exact} \\
 & \multicolumn{4}{c}{665.4 } \\
\end{tabular}
\end{ruledtabular}

\label{tab:extensiv_closed}
\end{table*}

\begin{figure}
        \includegraphics[width=\linewidth]{data/F2/DFT_avdzE_relat.pdf}
        \includegraphics[width=\linewidth]{data/F2/DFT_avdzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{fig2c}
        \caption{
        F$_2$, aug-cc-pvdz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:F2_avdz}}
\end{figure}

\begin{figure}
        \includegraphics[width=\linewidth]{data/F2/DFT_avtzE_relat.pdf}
        \includegraphics[width=\linewidth]{data/F2/DFT_avtzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{fig2c}
        \caption{
        F$_2$, aug-cc-pvtz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:F2_avtz}}
\end{figure}


\begin{figure}
%        \includegraphics[width=\linewidth]{data/H10/DFT_avdzE_relat.pdf}
        \includegraphics[width=\linewidth]{data/H10/DFT_vdzE_relat.pdf}\\
        \includegraphics[width=\linewidth]{data/H10/DFT_vdzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{data/H10/DFT_vdzE_error.pdf}\\
%       \includegraphics[width=\linewidth]{fig2c}
        \caption{
        H$_{10}$, cc-pvdz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:H10_vdz}}
\end{figure}


\begin{figure}
        \includegraphics[width=\linewidth]{data/H10/DFT_vtzE_relat.pdf}\\
        \includegraphics[width=\linewidth]{data/H10/DFT_vtzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{fig2c}
        \caption{
        H$_{10}$, cc-pvtz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:H10_vtz}}
\end{figure}


\begin{figure}
        \includegraphics[width=\linewidth]{data/H10/DFT_vqzE_relat.pdf}\\
        \includegraphics[width=\linewidth]{data/H10/DFT_vqzE_relat_zoom.pdf}
%       \includegraphics[width=\linewidth]{fig2c}
        \caption{
        H$_{10}$, cc-pvqz: Comparison between the near FCI and corrected near FCI energies and the estimated exact one. 
        \label{fig:H10_vqz}}
\end{figure}


\section{Conclusion}
\label{sec:conclusion}



\bibliography{srDFT_SC}

\end{document}
