loos/srDFT_SC
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge branch 'master' of https://git.irsamc.ups-tlse.fr/loos/srDFT_SC

Browse Source
This commit is contained in:
Emmanuel Giner 2020-03-25 15:42:05 +01:00
commit dd83164d7c
5 changed files with 90 additions and 1825 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
Revised_Manuscript/Makefile Executable file
View File

@ -0,0 +1,51 @@
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
#
# Make file simple pour LaTeX
#
# Compilation simple du document
# Compilation complete du document
#
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# fichier principal
TEX=srDFT_SC
# compilateur latex
CC=pdflatex
# defaut = aide :
default: simple
help:
@echo
@echo " Compilation du rapport, regles :"
@echo
@echo " simple : compilation une seule fois"
@echo " all : compile et met la biblio et les references a jour (long)"
@echo " purge : effacement des fichiers autres que .tex et .pdf (.aux .log .out .toc ..."
@echo
# compilation simple
simple: $(TEX).tex
$(CC) $(TEX)
# compilation complete, on passe plusieurs fois pour les references croisee plus bibtex pour la
# biblio
all: $(TEX).tex
$(CC) $(TEX)
$(CC) $(TEX)
@bibtex $(TEX)
$(CC) $(TEX)
$(CC) $(TEX)
@echo
@echo
@echo "fin ;)"
@echo
# efface les fichiers .log .aux .toc .bbl .blg
.PHONY: purge
purge:
@rm -vf *.aux $(TEX).bbl $(TEX).blg $(TEX).log $(TEX).out $(TEX).toc $(TEX).mtc* $(TEX).lof $(TEX).lot $(TEX).maf
@echo "done"

7
Revised_Manuscript/srDFT_SC.tex
View File

@ -747,6 +747,8 @@ For diatomics with the aug-cc-pVDZ and aug-cc-pVTZ basis sets,~\cite{KenDunHar-J
For the three diatomics, we performed an additional exFCI calculation with the aug-cc-pVQZ basis set at the equilibrium geometry to obtain reliable estimates of the FCI/CBS dissociation energy. For the three diatomics, we performed an additional exFCI calculation with the aug-cc-pVQZ basis set at the equilibrium geometry to obtain reliable estimates of the FCI/CBS dissociation energy.
In the case of the \ce{H10} chain, the approximation to the FCI energies together with the estimated exact potential energy curves are obtained from the data of Ref.~\onlinecite{h10_prx} where the authors performed MRCI+Q calculations with a minimal valence active space as reference (see below for the description of the active space). In the case of the \ce{H10} chain, the approximation to the FCI energies together with the estimated exact potential energy curves are obtained from the data of Ref.~\onlinecite{h10_prx} where the authors performed MRCI+Q calculations with a minimal valence active space as reference (see below for the description of the active space).
\alert{We note that, even though we use near-FCI energies in this work, the DFT-based basis-set correction could also be applied to any approximation to FCI such as multireference perturbation theory, similarly to what was done for weakly correlated systems for which the basis-set correction was applied to CCSD(T) calculations~\cite{LooPraSceTouGin-JCPL-19}.}
Regarding the complementary functional, we first perform full-valence CASSCF calculations with the GAMESS-US software~\cite{gamess} to obtain the wave function $\psibasis$. Then, all density-related quantities involved in the functional [density $n(\br{})$, effective spin polarization $\tilde{\zeta}(\br{})$, reduced density gradient $s(\br{})$, and on-top pair density $n_2(\br{})$] together with the local range-separation function $\mu(\br{})$ are calculated with this full-valence CASSCF wave function. The CASSCF calculations are performed with the following active spaces: (10e,10o) for \ce{H10}, (10e,8o) for \ce{N2}, (12e,8o) for \ce{O2}, and (14e,8o) for \ce{F2}. We note that, instead of using CASSCF wave functions for $\psibasis$, one could of course use the same selected-CI wave functions used for calculating the energy but the calculations of $n_2(\br{})$ and $\mu(\br{})$ would then be more costly. Regarding the complementary functional, we first perform full-valence CASSCF calculations with the GAMESS-US software~\cite{gamess} to obtain the wave function $\psibasis$. Then, all density-related quantities involved in the functional [density $n(\br{})$, effective spin polarization $\tilde{\zeta}(\br{})$, reduced density gradient $s(\br{})$, and on-top pair density $n_2(\br{})$] together with the local range-separation function $\mu(\br{})$ are calculated with this full-valence CASSCF wave function. The CASSCF calculations are performed with the following active spaces: (10e,10o) for \ce{H10}, (10e,8o) for \ce{N2}, (12e,8o) for \ce{O2}, and (14e,8o) for \ce{F2}. We note that, instead of using CASSCF wave functions for $\psibasis$, one could of course use the same selected-CI wave functions used for calculating the energy but the calculations of $n_2(\br{})$ and $\mu(\br{})$ would then be more costly.
Also, as the frozen-core approximation is used in all our selected-CI calculations, we use the corresponding valence-only complementary functionals (see Subsec.~\ref{sec:FC}). Therefore, all density-related quantities exclude any contribution from the 1s core orbitals, and the range-separation function follows the definition given in Eq.~\eqref{eq:def_mur_val}. Also, as the frozen-core approximation is used in all our selected-CI calculations, we use the corresponding valence-only complementary functionals (see Subsec.~\ref{sec:FC}). Therefore, all density-related quantities exclude any contribution from the 1s core orbitals, and the range-separation function follows the definition given in Eq.~\eqref{eq:def_mur_val}.
@ -816,7 +818,6 @@ The estimated exact energies are based on a fit of experimental data and obtaine
\caption{ \caption{
Potential energy curves of the \ce{F2} molecule calculated with exFCI and basis-set corrected exFCI using the aug-cc-pVDZ (top) and aug-cc-pVTZ (bottom) basis sets. Potential energy curves of the \ce{F2} molecule calculated with exFCI and basis-set corrected exFCI using the aug-cc-pVDZ (top) and aug-cc-pVTZ (bottom) basis sets.
The estimated exact energies are based on a fit of experimental data and obtained from Ref.~\onlinecite{LieCle-JCP-74a}. The estimated exact energies are based on a fit of experimental data and obtained from Ref.~\onlinecite{LieCle-JCP-74a}.
The estimated exact energies are based on a fit of the non-relativistic valence-only CEEIS data extracted from Ref.~\onlinecite{BytNagGorRue-JCP-07}.
\label{fig:F2}} \label{fig:F2}}
\end{figure*} \end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -827,7 +828,7 @@ The \ce{N2}, \ce{O2} and \ce{F2} molecules are complementary to the \ce{H10} sys
We report in Figs.~\ref{fig:N2}, \ref{fig:O2}, and \ref{fig:F2} the potential energy curves of \ce{N2}, \ce{O2}, and \ce{F2} computed at various approximation levels using the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The atomization energies for each level of theory with different basis sets are reported in Table \ref{tab:d0}. We report in Figs.~\ref{fig:N2}, \ref{fig:O2}, and \ref{fig:F2} the potential energy curves of \ce{N2}, \ce{O2}, and \ce{F2} computed at various approximation levels using the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The atomization energies for each level of theory with different basis sets are reported in Table \ref{tab:d0}.
Just as in \ce{H10}, the accuracy of the atomization energies is globally improved by adding the basis-set correction and it is remarkable that $\pbeontXi$ and $\pbeontns$ provide again very similar results. The latter observation confirms that the dependence on the on-top pair density allows one to remove the dependence of any kind of spin polarization for a quite wide range of covalent bonds and also for an open-shell system like \ce{O2}. More quantitatively, an error below 1.0 mHa compared to the estimated exact valence-only atomization energy is found for \ce{N2}, \ce{O2}, and \ce{F2} with the aug-cc-pVTZ basis set using the $\pbeontns$ functional, whereas such a feat is far from being reached within the same basis set at the near-FCI level. In the case of \ce{F2} it is clear that the addition of diffuse functions in the double- and triple-$\zeta$ basis sets strongly improves the accuracy of the results, which could have be anticipated due to the strong breathing-orbital effect induced by the ionic valence-bond forms in this molecule. \cite{HibHumByrLen-JCP-94} Just as in \ce{H10}, the accuracy of the atomization energies is globally improved by adding the basis-set correction and it is remarkable that $\pbeontXi$ and $\pbeontns$ provide again very similar results. The latter observation confirms that the dependence on the on-top pair density allows one to remove the dependence of any kind of spin polarization for a quite wide range of covalent bonds and also for an open-shell system like \ce{O2}. More quantitatively, an error below 1.0 mHa compared to the estimated exact valence-only atomization energy is found for \ce{N2}, \ce{O2}, and \ce{F2} with the aug-cc-pVTZ basis set using the $\pbeontns$ functional, whereas such a feat is far from being reached within the same basis set at the near-FCI level. In the case of \ce{F2} it is clear that the addition of diffuse functions in the double- and triple-$\zeta$ basis sets strongly improves the accuracy of the results, which could have been anticipated due to the strong breathing-orbital effect induced by the ionic valence-bond forms in this molecule. \cite{HibHumByrLen-JCP-94}
It should be also noticed that when reaching the aug-cc-pVQZ basis set for \ce{N2}, the accuracy of the atomization energy slightly deteriorates for the $\pbeontXi$ and $\pbeontns$ functionals, but it remains nevertheless more accurate than the estimated FCI atomization energy and very close to chemical accuracy. It should be also noticed that when reaching the aug-cc-pVQZ basis set for \ce{N2}, the accuracy of the atomization energy slightly deteriorates for the $\pbeontXi$ and $\pbeontns$ functionals, but it remains nevertheless more accurate than the estimated FCI atomization energy and very close to chemical accuracy.
\manu{ \manu{
@ -923,7 +924,7 @@ where the left-hand-side quantities are for the supersystem and the right-hand-s
\ket*{\wf{\text{A}+\text{B}}{}} = \ket*{\wf{\text{A}}{}} \otimes \ket*{\wf{\text{B}}{}}, \ket*{\wf{\text{A}+\text{B}}{}} = \ket*{\wf{\text{A}}{}} \otimes \ket*{\wf{\text{B}}{}},
\label{PsiAB} \label{PsiAB}
\end{equation} \end{equation}
where $\otimes$ is the antisymmetric tensor product. In this case, it is easy to shown that Eqs.~(\ref{nAB})-(\ref{sAB}) are valid, as well known, and it remains to show that Eqs.~(\ref{n2AB}) and~(\ref{muAB}) are also valid. Before showing this, we note that even though we do not explicitly consider the case of degeneracies, the lack of size consistency which could arise from spin-multiplet degeneracies can be avoided by the same strategy used for imposing the energy independence on $S_z$, \ie, by using the effective spin polarization $\tilde{\zeta}(n(\br{}),n_{2}(\br{}))$ or a zero spin polarization $\zeta(\br{}) = 0$. Moreover, the lack of size consistency which could arise from spatial degeneracies (\eg, coming from atomic $p$ states) can also be avoided by selecting the same member of the ensemble in the supersystem and in the isolated fragment. This applies to the systems treated in this work. where $\otimes$ is the antisymmetric tensor product. In this case, it is easy to shown that Eqs.~(\ref{nAB})-(\ref{sAB}) are valid, as well known, and it remains to show that Eqs.~(\ref{n2AB}) and~(\ref{muAB}) are also valid. Before showing this, we note that even though we do not explicitly consider the case of degeneracies, the lack of size consistency which could arise from spin-multiplet degeneracies can be avoided by the same strategy used for imposing the energy independence on $S_z$, \ie, by using the effective spin polarization $\tilde{\zeta}(n(\br{}),n_{2}(\br{}))$ or a zero spin polarization $\zeta(\br{}) = 0$. Moreover, \alert{for the systems treated in this work}, the lack of size consistency which could arise from spatial degeneracies (coming from atomic $p$ states) can also be avoided by selecting the \alert{same state} in the supersystem and in the isolated fragment. \alert{For example, for the F$_2$ molecule, the CASSCF wave function dissociates into the atomic configuration $\text{p}_\text{x}^2 \text{p}_\text{y}^2 \text{p}_\text{z}^1$ for each fragment, and we thus choose the same configuration for the calculation of the isolated atom. The same argument applies to the N$_2$ and O$_2$ molecules. For other systems, it may not be always possible to do so.}
\subsection{Intensivity of the on-top pair density and the local range-separation function} \subsection{Intensivity of the on-top pair density and the local range-separation function}

968
Revised_Manuscript/srDFT_SC_jt.tex
View File

@ -1,968 +0,0 @@
\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amsmath,amssymb,amsfonts,physics,mhchem,xspace}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}
\usepackage[
colorlinks=true,
citecolor=blue,
breaklinks=true
]{hyperref}
\urlstyle{same}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{HTML}{009900}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\jt}[1]{\textcolor{purple}{#1}}
\newcommand{\manu}[1]{\textcolor{darkgreen}{#1}}
\newcommand{\toto}[1]{\textcolor{brown}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashJT}[1]{\textcolor{purple}{\sout{#1}}}
\newcommand{\trashMG}[1]{\textcolor{darkgreen}{\sout{#1}}}
\newcommand{\trashAS}[1]{\textcolor{brown}{\sout{#1}}}
\newcommand{\MG}[1]{\manu{(\underline{\bf MG}: #1)}}
\newcommand{\JT}[1]{\juju{(\underline{\bf JT}: #1)}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\AS}[1]{\toto{(\underline{\bf AS}: #1)}}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\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}')}
\newcommand{\CBS}{\text{CBS}}
% 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{}}^\Bas[#1]}
\newcommand{\efuncdenpbeAB}[1]{\bar{E}_{\text{A}+\text{B}}^\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]{\bar{\varepsilon}_{\text{c,md}}^{\text{sr},\text{PBE}}}
\newcommand{\psibasis}[0]{\Psi^{\basis}}
\newcommand{\BasFC}{\mathcal{V}}
%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{\pbeuegxi}{\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{\pbeuegXi}{\text{PBE-UEG}}
\newcommand{\argpbeuegXi}[0]{\den,\tilde{\zeta},s,\ntwo^{\text{UEG}},\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeuegXi}[0]{\den(\br{}),\tilde{\zeta}(\br{}),s(\br{}),\ntwo^{\text{UEG}}(\br{}),\mu_{\text{}}(\br{})}
%pbeontxiCAS
%\newcommand{\pbeontxi}{\text{SPBE-OT}}
%\newcommand{\argpbeontxi}[0]{\den,\zeta,s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
%\newcommand{\argrpbeontxi}[0]{\den(\br{}),\zeta(\br{}),s(\br{}),\ntwoextrapcas(\br{}),\mu_{\text{CAS}}^{\basis}(\br{})}
%pbeontXiCAS
\newcommand{\pbeontXi}{\text{PBE-OT}}
\newcommand{\argpbeontXi}[0]{\den,\tilde{\zeta},s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeontXi}[0]{\den(\br{}),\tilde{\zeta}(\br{}),s(\br{}),\ntwoextrapcas(\br{}),\mu_{\text{}}^{}(\br{})}
%pbeont0xiCAS
\newcommand{\pbeuegns}{\text{SU-PBE-UEG}}
\newcommand{\argpbeuegns}[0]{\den,0,s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeuegns}[0]{\den(\br{}),0,s(\br{}),\ntwo^{\text{UEG}}(\br{}),\mu_{\text{}}^{}(\br{})}
%pbeont0xiCAS
\newcommand{\pbeontns}{\text{SU-PBE-OT}}
\newcommand{\argpbeontns}[0]{\den,0,s,\ntwoextrapcas,\mu_{\text{CAS}}^{\basis}}
\newcommand{\argrpbeontns}[0]{\den(\br{}),0,s(\br{}),\ntwoextrapcas(\br{}),\mu_{\text{}}^{}(\br{})}
%%%%%% arguments
\newcommand{\argepbe}[0]{\den,\zeta,s}
\newcommand{\argebasis}[0]{\den,\zeta,\ntwo,\mu}
\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(\br{}),\ntwo(\br{}),\mu(\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]{\mel{\Psi}{\psixc{#4}\psixc{#3} \psix{#2}\psix{#1}}{\Psi}}
\newcommand{\murpsi}[0]{\mu_{\wf{}{\Bas}}({\bf r})}
\newcommand{\murpsibas}[0]{\mu_{\wf{}{\Bas}}({\bf r})}
\newcommand{\ntwo}[0]{n_{2}}
\newcommand{\ntwohf}[0]{n_2^{\text{HF}}}
\newcommand{\ntwophi}[0]{n_2^{{\phi}}}
\newcommand{\ntwoextrap}[0]{\mathring{n}_{2}^{\text{}}}
\newcommand{\ntwoextrapcas}[0]{\mathring{n}_2^{\text{}}}
\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_{\wf{}{\Bas}}^{\text{val}}({\bf r})}
\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]{ n_{2,{\wf{}{\Bas}}}(\rr{1}{2})}
\newcommand{\twodmrdiagpsival}[0]{n_{2,\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}{W_{\wf{}{\Bas}}^{\text{val}}({\bf r},{\bf r})}
\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{\denbas}[0]{\den^\Bas}
\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
\renewcommand{\d}{\text{d}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hWee}[1]{\Hat{W}_\text{ee}^{#1}}
\newcommand{\f}[2]{f_{#1}^{#2}}
\newcommand{\Gam}[2]{\Gamma_{#1}^{#2}}
\newcommand{\isEquivTo}[1]{\underset{#1}{\sim}}
% 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{\LCT}{Laboratoire de Chimie Th\'eorique (UMR 7616), Sorbonne Universit\'e, CNRS, Paris, France}
\newcommand{\ISCD}{Institut des Sciences du Calcul et des Donn\'ees, Sorbonne Universit\'e, Paris, France}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\IUF}{Institut Universitaire de France, Paris, France}
\begin{document}
\title{A density-based basis-set correction for strongly correlated molecular systems}
\author{Emmanuel Giner}
\email{emmanuel.giner@lct.jussieu.fr}
\affiliation{\LCT}
\author{Barth\'el\'emy Pradines}
\affiliation{\LCT}
\affiliation{\ISCD}
\author{Anthony Scemama}
\affiliation{\LCPQ}
\author{Pierre-Fran\c{c}ois Loos}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Julien Toulouse}
\email{toulouse@lct.jussieu.fr}
\affiliation{\LCT}
\affiliation{\IUF}
\begin{abstract}
We extend to strongly correlated molecular systems the recently introduced basis-set incompleteness correction based on density-functional theory (DFT) [E. Giner \textit{et al.}, \href{https://doi.org/10.1063/1.5052714}{J. Chem. Phys. \textbf{149}, 194301 (2018)}]. This basis-set correction relies on a mapping between wave-function calculations in a finite basis set and range-separated DFT (RSDFT) through the definition of an effective non-divergent interaction corresponding to the electron-electron Coulomb interaction projected in the finite basis set. This enables the use of RSDFT-type complementary density functionals to recover the dominant part of the short-range correlation effects missing in this finite basis set. To study both weak and strong correlation regimes we consider the potential energy curves of the \ce{H10}, \ce{N2}, \ce{O2}, and \ce{F2} molecules up to the dissociation limit, and we explore various approximations of complementary functionals fulfilling two key properties: spin-multiplet degeneracy (\ie, independence of the energy with respect to the spin projection $S_z$) and size consistency. Specifically, we investigate the dependence of the functional on different types of on-top pair densities and spin polarizations. The key result of this study is that the explicit dependence on the on-top pair density allows one to completely remove the dependence on any form of spin polarization without any significant loss of accuracy. Quantitatively, we show that the basis-set correction reaches chemical accuracy on atomization energies with triple-$\zeta$ quality basis sets for most of the systems studied here. Also, the present basis-set incompleteness correction provides smooth potential energy curves along the whole range of internuclear distances.
\end{abstract}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%
The general goal of quantum chemistry is to provide reliable theoretical tools to explore the rich area of chemistry. More specifically, developments in quantum chemistry primarily aim at accurately computing the electronic structure of molecular systems. Despite intense developments, no definitive solution to this problem has been found. The theoretical challenge to tackle belongs to the quantum many-body problem, due to the intrinsic quantum nature of the electrons and the Coulomb repulsion between them. This so-called electronic correlation problem corresponds to finding a solution to the Schr\"odinger equation for a $N$-electron system, and two main roads have emerged to approximate this solution: wave-function theory (WFT) \cite{Pop-RMP-99} and density-functional theory (DFT). \cite{Koh-RMP-99} Although both WFT and DFT spring from the same Schr\"odinger equation, they rely on very different formalisms, as the former deals with the complicated $N$-electron wave function whereas the latter focuses on the much simpler one-electron density. In its Kohn-Sham (KS) formulation, \cite{KohSha-PR-65} the computational cost of DFT is very appealing since it is a simple mean-field procedure. Therefore, although continued efforts have been made to reduce the computational cost of WFT, DFT still remains the workhorse of quantum computational chemistry.
The difficulty of obtaining a reliable theoretical description of a given chemical system can be roughly categorized by the strength of the electronic correlation. The so-called weakly correlated systems, such as closed-shell organic molecules near their equilibrium geometry, are typically dominated by correlation effects which do not affect the qualitative mean-field picture of the system. These weak-correlation effects can be either short range (near the electron-electron coalescence points) \cite{HatKloKohTew-CR-12} or long range (London dispersion interactions). \cite{AngDobJanGou-BOOK-20} The theoretical description of weakly correlated systems is one of the most concrete achievement of quantum chemistry, and the main remaining challenge for these systems is to push the limit of the chemical system size that can be treated. The case of the so-called strongly correlated systems, which are ubiquitous in chemistry, is more problematic as they exhibit a much more complex electronic structure. For example, transition metal complexes, low-spin open-shell systems, covalent bond breaking situations have all in common that they cannot be even qualitatively described by a single electronic configuration. It is now clear that the usual semilocal density-functional approximations of KS DFT fail to accurately describe these situations \cite{GorSeiSav-PCCP-08,GagTruLiCarHoyBa-ACR-17} and WFT is king for the treatment of strongly correlated systems.
In practice, WFT uses a finite one-electron basis set. The exact solution of the Schr\"odinger equation within this basis set is then provided by full configuration interaction (FCI) which consists in a linear-algebra eigenvalue problem with a dimension scaling exponentially with the system size. Due to this exponential growth of the FCI computational cost, introducing approximations is necessary, with at least two difficulties for strongly correlated systems: i) the qualitative description of the wave function is determined by a primary set of electronic configurations (whose size can scale exponentially in many cases) among which near degeneracies and/or strong interactions appear in the Hamiltonian matrix; ii) the quantitative description of the system requires also to account for weak correlation effects which involve many other electronic configurations with typically much smaller weights in the wave function. Addressing simultaneously these two issues is a rather complicated task for a given approximate WFT method, especially if one adds the requirement of satisfying formal properties, such as spin-multiplet degeneracy (\ie, independence of the energy with respect to the spin projection $S_z$) and size consistency.
Beside the difficulties of accurately describing the molecular electronic structure within a given basis set, a crucial limitation of WFT methods is the slow convergence of the energy (and related properties) with respect to the size of the one-electron basis set. As initially shown by the seminal work of Hylleraas \cite{Hyl-ZP-29} and further developed by Kutzelnigg and coworkers, \cite{Kut-TCA-85,KutKlo-JCP-91, NogKut-JCP-94} the main convergence problem originates from the divergence of the electron-electron Coulomb interaction at the coalescence point, which induces a discontinuity in the first derivative of the exact wave function (the so-called electron-electron cusp). Describing such a discontinuity with an incomplete one-electron basis set is impossible and, as a consequence, the convergence of the computed energies and properties are strongly affected. To alleviate this problem, extrapolation techniques have been developed, either based on a partial-wave expansion analysis, \cite{HelKloKocNog-JCP-97,HalHelJorKloKocOlsWil-CPL-98} or more recently based on perturbative arguments. \cite{IrmHulGru-PRL-19,IrmGru-JCP-2019} A more rigorous approach to tackle the basis-set convergence problem is provided by the so-called explicitly correlated F12 (or R12) methods \cite{Ten-TCA-12,TenNog-WIREs-12,HatKloKohTew-CR-12, KonBisVal-CR-12, GruHirOhnTen-JCP-17, MaWer-WIREs-18} which introduce a geminal function depending explicitly on the interelectronic distance. This ensures a correct representation of the Coulomb correlation hole around the electron-electron coalescence point, and leads to a much faster convergence of the energy than usual WFT methods. For instance, using the explicitly correlated version of coupled cluster with singles, doubles, and perturbative triples [CCSD(T)] in a triple-$\zeta$ basis set is equivalent to using a quintuple-$\zeta$ basis set with the usual CCSD(T) method, \cite{TewKloNeiHat-PCCP-07} although a computational overhead is introduced by the auxiliary basis set needed to compute the three-electron integrals involved in F12 theory. \cite{BarLoo-JCP-17} In addition to the computational cost, a possible drawback of F12 theory is its rather complex formalism which requires non-trivial developments for adapting it to a new method. For strongly correlated systems, several multi-reference methods have been extended to explicit correlation (see, for example, Refs.~\onlinecite{Ten-CPL-07,ShiWer-JCP-10,TorKniWer-JCP-11,DemStanMatTenPitNog-PCCP-12,GuoSivValNee-JCP-17}), including approaches based on the so-called universal F12 theory which are potentially applicable to any electronic-structure computational methods. \cite{TorVal-JCP-09,KonVal-JCP-11,HauMaoMukKlo-CPL-12,BooCleAlaTew-JCP-12}
An alternative way to improve the convergence towards the complete-basis-set (CBS) limit is to treat the short-range correlation effects within DFT and to use WFT methods to deal only with the long-range and/or strong correlation effects. A rigorous approach achieving this mixing of DFT and WFT is range-separated DFT (RSDFT) (see Ref.~\onlinecite{TouColSav-PRA-04} and references therein) which relies on a decomposition of the electron-electron Coulomb interaction in terms of the interelectronic distance thanks to a range-separation parameter $\mu$. The advantage of this approach is at least two-fold: i) the DFT part deals primarily with the short-range part of the Coulomb interaction, and consequently the usual semilocal density-functional approximations are more accurate than for standard KS DFT; ii) the WFT part deals only with a smooth non-divergent interaction, and consequently the wave function has no electron-electron cusp \cite{GorSav-PRA-06} and the basis-set convergence is much faster. \cite{FraMusLupTou-JCP-15} A number of approximate RSDFT schemes have been developed involving single-reference \cite{AngGerSavTou-PRA-05, GolWerSto-PCCP-05, TouGerJanSavAng-PRL-09,JanHenScu-JCP-09, TouZhuSavJanAng-JCP-11, MusReiAngTou-JCP-15,KalTou-JCP-18,KalMusTou-JCP-19} and multi-reference \cite{LeiStoWerSav-CPL-97, FroTouJen-JCP-07, FroCimJen-PRA-10, HedKneKieJenRei-JCP-15, HedTouJen-JCP-18, FerGinTou-JCP-18} WFT methods. Nevertheless, there are still some open issues in RSDFT, such as remaining fractional-charge and fractional-spin errors in the short-range density functionals \cite{MusTou-MP-17} or the dependence of the quality of the results on the value of the range-separation parameter $\mu$.
Building on the development of RSDFT, a possible solution to the basis-set convergence problem has been recently proposed by some of the present authors~\cite{GinPraFerAssSavTou-JCP-18} in which RSDFT functionals are used to recover only the correlation effects outside a given basis set. The key point here is to realize that a wave function developed in an incomplete basis set is cuspless and could also originate from a Hamiltonian with a non-divergent long-range electron-electron interaction. Therefore, a mapping with RSDFT can be performed through the introduction of an effective non-divergent interaction representing the usual electron-electron Coulomb interaction projected in an incomplete basis set. First applications to weakly correlated molecular systems have been successfully carried out, \cite{LooPraSceTouGin-JCPL-19} together with extensions of this approach to the calculations of excitation energies \cite{GinSceTouLoo-JCP-19} and ionization potentials. \cite{LooPraSceGinTou-JCTC-20} The goal of the present work is to further develop this approach for the description of strongly correlated systems.
The paper is organized as follows. In Sec.~\ref{sec:theory}, we recall the mathematical framework of the basis-set correction and we present its extension for strongly correlated systems. In particular, our focus is primarily set on imposing two key formal properties: spin-multiplet degeneracy and size consistency.
Then, in Sec.~\ref{sec:results}, we apply the method to the calculation of the potential energy curves of the \ce{H10}, \ce{N2}, \ce{O2}, and \ce{F2} molecules up to the dissociation limit. Finally, we conclude in Sec.~\ref{sec:conclusion}.
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theory}
%%%%%%%%%%%%%%%%%%%%%%%%
As the theory behind the present basis-set correction has been exposed in details in Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}, we only briefly recall the main equations and concepts needed for this study in Secs.~\ref{sec:basic}, \ref{sec:wee}, and \ref{sec:mur}. More specifically, in Sec.~\ref{sec:basic}, we recall the basic mathematical framework of the present theory by introducing the complementary functional to a basis set. Section \ref{sec:wee} introduces the effective non-divergent interaction in the basis set, which leads us to the definition of the effective \textit{local} range-separation function in Sec.~\ref{sec:mur}. Then, Sec.~\ref{sec:functional} exposes the new approximate RSDFT-based complementary correlation functionals. The generic form of such functionals is exposed in Sec.~\ref{sec:functional_form}, their properties in the context of the basis-set correction are discussed in Sec.~\ref{sec:functional_prop}, and the specific requirements for strong correlation are discussed in Sec.~\ref{sec:requirements}. Finally, the actual functionals used in this work are introduced in Sec.~\ref{sec:def_func}.
\subsection{Basic theory}
\label{sec:basic}
The exact ground-state energy $E_0$ of a $N$-electron system can, in principle, be obtained in DFT by a minimization over $N$-representable one-electron densities $\denr$
\begin{equation}
\label{eq:levy}
E_0 = \min_{\den} \bigg\{ F[\den] + \int \d \br{} v_{\text{ne}} (\br{}) \denr \bigg\},
\end{equation}
where $v_\text{ne}(\br{})$ is the nuclei-electron potential, and $F[\den]$ is the universal Levy-Lieb density functional written with the constrained search formalism as~\cite{Lev-PNAS-79,Lie-IJQC-83}
\begin{equation}
\label{eq:levy_func}
F[\den] = \min_{\Psi \to \den} \mel{\Psi}{\kinop +\weeop}{\Psi},
\end{equation}
where $\kinop$ and $\weeop$ are the kinetic and electron-electron Coulomb operators, and the notation $\Psi \to \den$ means that the wave function $\Psi$ yields the density $\den$.
The minimizing density $n_0$ in Eq.~\eqref{eq:levy} is the exact ground-state density. Nevertheless, in practical calculations, the accessible densities are necessarily restricted to the set of densities ``representable in a basis set $\Bas$'', \ie, densities coming from wave functions expandable in the $N$-electron Hilbert space generated by the one-electron basis set $\Bas$. In the following, we always consider only such representable-in-$\Bas$ densities. With this restriction, Eq.~\eqref{eq:levy} then gives an upper bound $E_0^\Bas$ of the exact ground-state energy. Since the density has a faster convergence with the size of the basis set than the wave function, this restriction is a rather weak one and we can consider that $E_0^\Bas$ is an acceptable approximation to the exact ground-state energy, \ie, $E_0^\Bas \approx E_0$.
In the present context, it is important to notice that the wave functions $\Psi$ defined in Eq.~\eqref{eq:levy_func} are not restricted to a finite basis set, \ie, they should be expanded in a complete basis set. In Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}, it was then proposed to decompose $F[\den]$ as, for a representable-in-$\Bas$ density $\den$,
\begin{equation}
\label{eq:def_levy_bas}
F[\den] = \min_{\wf{}{\Bas} \to \den} \mel*{\wf{}{\Bas}}{\kinop +\weeop}{\wf{}{\Bas}} + \efuncden{\den},
\end{equation}
where $\wf{}{\Bas}$ are wave functions expandable in the $N$-electron Hilbert space generated by $\basis$, and
\begin{equation}
\begin{aligned}
\efuncden{\den} = \min_{\Psi \to \den} \mel*{\Psi}{\kinop +\weeop }{\Psi}
- \min_{\Psi^{\Bas} \to \den} \mel*{\wf{}{\Bas}}{\kinop +\weeop}{\wf{}{\Bas}}
\end{aligned}
\end{equation}
is the complementary density functional to the basis set $\Bas$.
Introducing the decomposition in Eq.~\eqref{eq:def_levy_bas} back into Eq.~\eqref{eq:levy} yields
\begin{multline}
\label{eq:E0basminPsiB}
E_0^\Bas = \min_{\Psi^{\Bas}} \bigg\{ \mel*{\wf{}{\Bas}}{\kinop +\weeop}{\wf{}{\Bas}} + \efuncden{\den_{{\Psi^{\Bas}}}}
\\
+ \int \d \br{} v_{\text{ne}} (\br{}) \den_{\Psi^{\Bas}}(\br{}) \bigg\},
\end{multline}
where the minimization is only over wave functions $\wf{}{\Bas}$ restricted to the basis set $\basis$ and $\den_{{\Psi^{\Bas}}}(\br{})$ refers to the density generated from $\wf{}{\Bas}$. Therefore, thanks to Eq.~\eqref{eq:E0basminPsiB}, one can properly combine a WFT calculation in a finite basis set with a density functional (hereafter referred to as complementary functional) accounting for the correlation effects that are not included in the basis set.
As a simple non-self-consistent version of this approach, we can approximate the minimizing wave function $\Psi_0^{\Bas}$ in Eq.~\eqref{eq:E0basminPsiB} by the ground-state FCI wave function $\psifci$ within $\Bas$, and we then obtain the following approximation for the exact ground-state energy [see Eqs.~(12)--(15) of Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}]
\begin{equation}
\label{eq:e0approx}
E_0 \approx E_0^\Bas \approx \efci + \efuncbasisFCI,
\end{equation}
where $\efci$ and $n_\text{FCI}^\Bas$ are the ground-state FCI energy and density, respectively. As it was originally shown in Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18} and further emphasized in Refs.~\onlinecite{LooPraSceTouGin-JCPL-19,GinSceTouLoo-JCP-19}, the main role of $\efuncbasisFCI$ is to correct for the basis-set incompleteness error, a large part of which originating from the lack of electron-electron cusp in the wave function expanded in an incomplete basis set. The whole purpose of this work is to determine approximations for $\efuncbasisFCI$ which are suitable for strongly correlated molecular systems. Two key requirements for this purpose are i) spin-multiplet degeneracy, and ii) size consistency.
\subsection{Effective interaction in a finite basis}
\label{sec:wee}
As originally shown by Kato, \cite{Kat-CPAM-57} the electron-electron cusp of the exact wave function originates from the divergence of the Coulomb interaction at the coalescence point. Therefore, a cuspless wave function $\wf{}{\Bas}$ could also be obtained from a Hamiltonian with a non-divergent electron-electron interaction.
In other words, the impact of the basis set incompleteness can be understood as the removal of the divergence of the usual electron-electron Coulomb interaction.
As originally derived in Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18} (see Sec.~II D~and Appendices), one can obtain an effective non-divergent electron-electron interaction, here referred to as $\wbasis$, which reproduces the expectation value of the electron-electron Coulomb interaction operator over a given wave function $\wf{}{\Bas}$. As we are interested in the behavior at the coalescence point, we focus on the opposite-spin part of the electron-electron interaction. More specifically, the effective electron-electron interaction associated to a given wave function $\wf{}{\Bas}$ is defined as
\begin{equation}
\label{eq:wbasis}
\wbasis =
\begin{cases}
\fbasis /\twodmrdiagpsi, & \text{if $\twodmrdiagpsi \ne 0$,}
\\
\infty, & \text{otherwise,}
\end{cases}
\end{equation}
where
\begin{equation}
\twodmrdiagpsi = \sum_{pqrs \in \Bas} \SO{p}{1} \SO{q}{2} \Gam{pq}{rs} \SO{r}{1} \SO{s}{2},
\end{equation}
is the opposite-spin pair density associated with $\wf{}{\Bas}$, and $\Gam{pq}{rs} = 2 \mel*{\wf{}{\Bas}}{ \aic{r_\downarrow}\aic{s_\uparrow}\ai{q_\uparrow}\ai{p_\downarrow}}{\wf{}{\Bas}}$ its associated tensor in a basis of spatial orthonormal orbitals $\{\SO{p}{}\}$,
\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}= \braket{pq}{rs}$ are the usual two-electron Coulomb integrals.
With such a definition, one can show that $\wbasis$ satisfies
\begin{multline}
\frac{1}{2}\iint \dr{1} \dr{2} \wbasis \twodmrdiagpsi =
\\
\frac{1}{2} \iint \dr{1} \dr{2} \frac{\twodmrdiagpsi}{\abs{\br{1}-\br{2}}}.
\end{multline}
As shown in Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}, the effective interaction $\wbasis$ is necessarily finite at coalescence for an incomplete basis set, and tends to the usual Coulomb interaction in the CBS limit for any choice of wave function $\psibasis$, \ie,
\begin{equation}
\label{eq:cbs_wbasis}
\lim_{\Bas \to \text{CBS}} \wbasis = \frac{1}{\abs{\br{1}-\br{2}}},\quad \forall\,\psibasis.
\end{equation}
The condition in Eq.~\eqref{eq:cbs_wbasis} is fundamental as it guarantees the correct behavior of the theory in the CBS limit.
\subsection{Local range-separation function}
\label{sec:mur}
\subsubsection{General definition}
The effective interaction within a finite basis, $\wbasis$, is bounded and resembles the long-range interaction used in RSDFT
\begin{equation}
\label{eq:weelr}
w_\text{ee}^{\lr}(\mu;r_{12}) = \frac{\text{erf}\big(\mu \,r_{12} \big)}{r_{12}},
\end{equation}
where $\mu$ is the range-separation parameter. As originally proposed in Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}, we make the correspondence between these two interactions by using the local range-separation function
\begin{equation}
\label{eq:def_mur}
\murpsi = \frac{\sqrt{\pi}}{2} \wbasiscoal,
\end{equation}
such that the two interactions coincide at the electron-electron coalescence point for each $\br{}$
\begin{equation}
w_\text{ee}^{\lr}(\murpsi;0) = \wbasiscoal, \quad \forall \, \br{}.
\end{equation}
Because of the very definition of $\wbasis$, one has the following property in the CBS limit [see Eq.~\eqref{eq:cbs_wbasis}]
\begin{equation}
\label{eq:cbs_mu}
\lim_{\Bas \to \text{CBS}} \murpsi = \infty, \quad \forall \,\psibasis,
\end{equation}
which is again fundamental to guarantee the correct behavior of the theory in the CBS limit.
\subsubsection{Frozen-core approximation}
\label{sec:FC}
As all WFT calculations in this work are performed within the frozen-core approximation, we use a ``valence-only'' (or no-core) version of the various quantities needed for the complementary functional introduced in Ref.~\onlinecite{LooPraSceTouGin-JCPL-19}. We partition the basis set as $\Bas = \Cor \bigcup \BasFC$, where $\Cor$ and $\BasFC$ are the sets of core and ``valence'' (\ie, non-core) orbitals, respectively, and define the valence-only local range-separation function as
\begin{equation}
\label{eq:def_mur_val}
\murpsival = \frac{\sqrt{\pi}}{2} \wbasiscoalval{},
\end{equation}
where
\begin{equation}
\label{eq:wbasis_val}
\wbasisval =
\begin{cases}
\fbasisval /\twodmrdiagpsival, & \text{if $\twodmrdiagpsival \ne 0$,}
\\
\infty, & \text{otherwise,}
\end{cases}
\end{equation}
is the valence-only effective interaction and
\begin{gather}
\label{eq:fbasis_val}
\fbasisval
= \sum_{pq\in \Bas} \sum_{rstu \in \BasFC} \SO{p}{1} \SO{q}{2} \V{pq}{rs} \Gam{rs}{tu} \SO{t}{1} \SO{u}{2},
\\
\label{eq:twordm_val}
\twodmrdiagpsival
= \sum_{pqrs \in \BasFC} \SO{p}{1} \SO{q}{2} \Gam{pq}{rs} \SO{r}{1} \SO{s}{2}.
\end{gather}
One would note the restrictions of the sums to the set $\BasFC$ in Eqs.~\eqref{eq:fbasis_val} and \eqref{eq:twordm_val}.
It is also noteworthy that, with the present definition, $\wbasisval$ still tends to the usual Coulomb interaction as $\Bas \to \CBS$.
\subsection{General form of the complementary functional}
\label{sec:functional}
\subsubsection{Generic approximate form}
\label{sec:functional_form}
As originally proposed and motivated in Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}, we approximate the complementary functional $\efuncden{\den}$ by using the so-called correlation energy functional with multideterminant reference (ECMD) introduced by Toulouse \textit{et al.}.\cite{TouGorSav-TCA-05,Tou-THESIS-05} Following the recent work in Ref.~\onlinecite{LooPraSceTouGin-JCPL-19}, we propose to consider a Perdew-Burke-Ernzerhof (PBE)-like functional which uses the one-electron density $\denr$, the spin polarization $\zeta(\br{})=[n_\uparrow(\br{})-n_\downarrow(\br{})]/\denr$ (where $n_\uparrow(\br{})$ and $n_\downarrow(\br{})$ are the spin-up and spin-down densities), the reduced density gradient $s(\br{}) = \nabla \denr/\denr^{4/3}$, and the on-top pair density $\ntwo(\br{})\equiv \ntwo(\br{},\br{})$. In the present work, all these quantities are computed with the same wave function $\psibasis$ used to define $\mur \equiv\murpsi$.
Therefore, $\efuncden{\den}$ has the following generic form
\begin{multline}
\label{eq:def_ecmdpbebasis}
\efuncdenpbe{\argebasis} =
\\
\int \d\br{} \,\denr \ecmd(\argrebasis),
\end{multline}
where
\begin{equation}
\label{eq:def_ecmdpbe}
\ecmd(\argecmd) = \frac{\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)}{1+ \beta(\argepbe,\ntwo) \; \mu^3}
\end{equation}
is the correlation energy per particle, with
\begin{equation}
\label{eq:def_beta}
\beta(\argepbe,\ntwo) = \frac{3}{2\sqrt{\pi}(1 - \sqrt{2})}\frac{\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)}{\ntwo/\den},
\end{equation}
where $\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)$ is the usual PBE correlation energy per particle. \cite{PerBurErn-PRL-96} Before introducing the different flavors of approximate functionals that we will use here (see Sec.~\ref{sec:def_func}), we would like to give some motivations for this choice of functional form.
The form of $\ecmd(\argecmd)$ in Eq.~\eqref{eq:def_ecmdpbe} has been originally proposed in Ref.~\onlinecite{FerGinTou-JCP-18} in the context of RSDFT. In the $\mu\to 0$ limit, it reduces to the usual PBE correlation functional, \ie,
\begin{equation}
\lim_{\mu \to 0} \ecmd(\argecmd) = \varepsilon_{\text{c}}^{\text{PBE}}(\argepbe),
\end{equation}
which is relevant in the weak-correlation (or high-density) limit. In the large-$\mu$ limit, it behaves as
\begin{equation}
\label{eq:lim_mularge}
\ecmd(\argecmd) \isEquivTo{\mu\to\infty} \frac{2\sqrt{\pi}(1 - \sqrt{2})}{3 \mu^3} \frac{\ntwo}{n},
\end{equation}
which is the exact large-$\mu$ behavior of the exact ECMD correlation energy. \cite{PazMorGorBac-PRB-06,FerGinTou-JCP-18} Of course, for a specific system, the large-$\mu$ behavior will be exact only if one uses for $n_2$ the \textit{exact} on-top pair density of this system. This large-$\mu$ limit in Eq.~\eqref{eq:lim_mularge} is relevant in the strong-correlation (or low-density) limit. 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 indeed 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 recently acknowledged by Gagliardi and coworkers \cite{CarTruGag-JPCA-17} and Pernal and coworkers.\cite{GritMeePer-PRA-18}
Note also that $\ecmd(\argecmd)$ vanishes when $\ntwo$ vanishes, \ie,
\begin{equation}
\label{eq:lim_n2}
\lim_{\ntwo \to 0} \ecmd(\argecmd) = 0,
\end{equation}
which is expected for systems with a vanishing on-top pair density.
Finally, the function $\ecmd(\argecmd)$ vanishes when $\mu \to \infty$ like all RSDFT short-range functionals, \ie,
\begin{equation}
\label{eq:lim_muinf}
\lim_{\mu \to \infty} \ecmd(\argecmd) = 0.
\end{equation}
\subsubsection{Two limits where the complementary functional vanishes}
\label{sec:functional_prop}
Within the definitions of Eqs.~\eqref{eq:def_mur} and \eqref{eq:def_ecmdpbebasis}, any approximate complementary functional $\efuncdenpbe{\argebasis}$ satisfies two important properties.
First, thanks to the properties in Eqs.~\eqref{eq:cbs_mu} and~\eqref{eq:lim_muinf}, $\efuncdenpbe{\argebasis}$ vanishes in the CBS limit, independently of the type of wave function $\psibasis$ used to define the local range-separation function $\mu(\br{})$ in a given basis set $\Bas$,
\begin{equation}
\label{eq:lim_ebasis}
\lim_{\basis \to \text{CBS}} \efuncdenpbe{\argebasis} = 0, \quad \forall\, \psibasis.
\end{equation}
Second, $\efuncdenpbe{\argebasis}$ correctly vanishes for systems with uniformly vanishing on-top pair density, such as one-electron systems and for the stretched H$_2$ molecule,
\begin{equation}
\lim_{n_2 \to 0} \efuncdenpbe{\argebasis} = 0.
\end{equation}
This property is doubly guaranteed by i) the choice of setting $\wbasis = \infty$ for a vanishing pair density [see Eq.~\eqref{eq:wbasis}], which leads to $\mu(\br{}) \to \infty$ and thus a vanishing $\ecmd(\argecmd)$ [see Eq.~\eqref{eq:lim_muinf}], and ii) the fact that $\ecmd(\argecmd)$ vanishes anyway when the on-top pair density vanishes [see Eq.~\eqref{eq:lim_n2}].
\subsection{Requirements on the complementary functional for strong correlation}
\label{sec:requirements}
An important requirement for any electronic-structure method is size consistency, \ie, the additivity of the energies of non-interacting fragments, which is mandatory to avoid any ambiguity in computing interaction energies. When two subsystems \ce{A} and \ce{B} dissociate in closed-shell systems, as in the case of weak intermolecular interactions for instance, spin-restricted Hartree-Fock (RHF) is size-consistent. When the two subsystems dissociate in open-shell systems, such as in covalent bond breaking, it is well known that the RHF approach fails and an alternative is to use a complete-active-space self-consistent-field (CASSCF) wave function which, provided that the active space has been properly chosen, leads to additive energies.
Another important requirement is spin-multiplet degeneracy, \ie, the independence of the energy with respect to the $S_z$ component of a given spin state, which is also a property of any exact wave function. Such a property is also important in the context of covalent bond breaking where the ground state of the supersystem $\ce{A + B}$ is generally of lower spin than the corresponding ground states of the fragments (\ce{A} and \ce{B}) which can have multiple $S_z$ components.
\subsubsection{Spin-multiplet degeneracy}
A sufficient condition to achieve spin-multiplet degeneracy is to eliminate all dependencies on $S_z$. In the case of the function $\ecmd(\argecmd)$, this means removing the dependence on the spin polarization $\zeta(\br{})$ originating from the PBE correlation functional $\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)$ [see Eq.~\eqref{eq:def_ecmdpbe}].
To do so, it has been proposed to replace the dependence on the spin polarization by the dependence on the on-top pair density. Most often, it is done by introducing an effective spin polarization~\cite{MosSan-PRA-91,BecSavSto-TCA-95,Sav-INC-96a,Sav-INC-96,MieStoSav-MP-97,TakYamYam-CPL-02,TakYamYam-IJQC-04,GraCre-MP-05,TsuScuSav-JCP-10,LimCarLuoMaOlsTruGag-JCTC-14,GarBulHenScu-JCP-15,GarBulHenScu-PCCP-15,CarTruGag-JCTC-15,GagTruLiCarHoyBa-ACR-17} (see, also, Refs.~\onlinecite{PerSavBur-PRA-95,StaDav-CPL-01})
\begin{equation}
\label{eq:def_effspin}
\tilde{\zeta}(n,n_{2}) = \sqrt{ 1 - 2 \; n_{2}/n^2 }
\end{equation}
expressed as a function of the density $n$ and the on-top pair density $n_2$ calculated from a given wave function. The advantage of this approach is that this effective spin polarization $\tilde{\zeta}$ is independent from $S_z$ since the on-top pair density is $S_z$-independent. Nevertheless, the use of $\tilde{\zeta}$ in Eq.~\eqref{eq:def_effspin} presents some disadvantages since this expression was derived for a single-determinant wave function. Hence, it does not appear justified to use it for a multideterminant wave function. More particularly, it may happen, in the multideterminant case, that $1 - 2 \; n_{2}/n^2 < 0 $ which results in a complex-valued effective spin polarization. \cite{BecSavSto-TCA-95}
Therefore, following other authors, \cite{MieStoSav-MP-97,LimCarLuoMaOlsTruGag-JCTC-14,GarBulHenScu-JCP-15} we use the following definition
\begin{equation}
\label{eq:def_effspin-0}
\tilde{\zeta}(n,n_{2}) =
\begin{cases}
\sqrt{ 1 - 2 \; n_{2}/n^2 }, & \text{if } n^2 \ge 2 n_{2},
\\
0, & \text{otherwise.}
\end{cases}
\end{equation}
An alternative way to eliminate the $S_z$ dependence is to simply set $\zeta=0$, \ie, to resort to the spin-unpolarized functional. This lowers the accuracy for open-shell systems at $\mu=0$, \ie, for the usual PBE correlation functional $\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)$. Nevertheless, we argue that, for sufficiently large $\mu$, it is a viable option. Indeed, the purpose of introducing the spin polarization in semilocal density-functional approximations is to mimic the exact on-top pair density, \cite{PerSavBur-PRA-95} but our functional $\ecmd(\argecmd)$ already explicitly depends on the on-top pair density [see Eqs.~\eqref{eq:def_ecmdpbe} and \eqref{eq:def_beta}]. The dependencies on $\zeta$ and $n_2$ can thus be expected to be largely redundant. Consequently, we propose here to test the use of $\ecmd$ with \textit{a zero spin polarization}. This ensures its $S_z$ independence and, as will be numerically demonstrated, very weakly affects the complementary functional accuracy.
\subsubsection{Size consistency}
Since $\efuncdenpbe{\argebasis}$ is computed via a single integral over $\mathbb{R}^3$ [see Eq.~\eqref{eq:def_ecmdpbebasis}] which involves only local quantities [$n(\br{})$, $\zeta(\br{})$, $s(\br{})$, $n_2(\br{})$, and $\mu(\br{})$], in the case of non-overlapping fragments $\text{A}+\text{B}$, it can be written as the sum of two local contributions: one coming from the integration over the region of subsystem \ce{A} and the other one from the region of subsystem \ce{B}. Therefore, a sufficient condition for size consistency is that these quantities locally coincide in the isolated fragments and in the supersystem $\text{A}+\text{B}$. Since these local quantities are calculated from the wave function $\psibasis$, a sufficient condition is that the wave function is multiplicatively separable in the limit of non-interacting fragments, \ie, $\ket*{\Psi_{\text{A}+\text{B}}^{\basis}} = \ket*{\Psi_{\ce{A}}^{\basis}} \otimes \ket*{\Psi_{\ce{B}}^{\basis}}$. We refer the interested reader to Appendix~\ref{app:sizeconsistency} for a detailed proof and discussion of the latter statement.
In the case where the two subsystems \ce{A} and \ce{B} dissociate in closed-shell systems, a simple RHF wave function ensures this property, but when one or several covalent bonds are broken, a properly chosen CASSCF wave function can be used to recover this property. The underlying active space must however be chosen in such a way that it leads to size-consistent energies in the limit of dissociated fragments.
\subsection{Actual approximations used for the complementary functional}
\label{sec:def_func}
As the present work focuses on the strong-correlation regime, we propose here to investigate only approximate functionals which are $S_z$ independent and size-consistent in the case of covalent bond breaking. Therefore, the wave functions $\psibasis$ used throughout this paper are CASSCF wave functions in order to ensure size consistency of all local quantities. The difference between the different flavors of functionals are only due to the types of spin polarization and on-top pair density used.
Regarding the spin polarization that enters into the function $\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)$, two different types of $S_z$-independent formulations are considered: i) the \textit{effective} spin polarization $\tilde{\zeta}$ defined in Eq.~\eqref{eq:def_effspin-0} and calculated from the CASSCF wave function, and ii) a \textit{zero} spin polarization. In the latter case, the functional is referred as to ``SU'' which stands for ``spin unpolarized''.
Regarding the on-top pair density entering in Eq.~\eqref{eq:def_beta}, 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) \approx n^2\big(1-\zeta^2\big)g_0(n),
\end{equation}
where the pair-distribution function $g_0(n)$ is taken from Eq.~(46) of Ref.~\onlinecite{GorSav-PRA-06}. As the spin polarization appears in Eq.~\eqref{eq:def_n2ueg}, we use the effective spin polarization $\tilde{\zeta}$ of Eq.~\eqref{eq:def_effspin-0} in order to ensure $S_z$ independence. Thus, $\ntwo^{\text{UEG}}$ will depend indirectly on the on-top pair density of the CASSCF wave function through $\tilde{\zeta}$. When using $\ntwo^{\text{UEG}}(\br{}) \equiv \ntwo^{\text{UEG}}(n(\br{}),\tilde{\zeta}(\br{}))$ in a functional, we will refer to it as ``UEG''.
The second approach to approximate the exact on-top pair density consists in using directly the on-top pair density of the CASSCF wave function. Following the work of some of the present authors, \cite{FerGinTou-JCP-18,GinSceTouLoo-JCP-19} we introduce the extrapolated on-top pair density
\begin{equation}
\label{eq:def_n2extrap}
\ntwoextrap(\ntwo,\mu) = \qty( 1 + \frac{2}{\sqrt{\pi}\mu} )^{-1} \; \ntwo,
\end{equation}
which directly follows from the large-$\mu$ extrapolation of the exact on-top pair density derived by Gori-Giorgi and Savin\cite{GorSav-PRA-06} in the context of RSDFT. Thus, the extrapolated on-top pair density $\ntwoextrap$ depends on the local range-separation function $\mu$. When using $\ntwoextrap(\br{}) \equiv \ntwoextrap(\ntwo(\br{}),\mu(\br{}))$ in a functional, we will simply refer it as ``OT'', which stands for "on-top".
We then define three complementary functionals:
\begin{itemize}
\item[i)] $\pbeuegXi$ which combines the effective spin polarization of Eq.~\eqref{eq:def_effspin-0} and the UEG on-top pair density defined in Eq.~\eqref{eq:def_n2ueg}:
\begin{multline}
\label{eq:def_pbeueg_i}
\bar{E}^\Bas_{\pbeuegXi}
\\
= \int \d\br{} \,\denr \ecmd(\argrpbeuegXi),
\end{multline}
\item[ii)] $\pbeontXi$ which combines the effective spin polarization of Eq.~\eqref{eq:def_effspin-0} and the on-top pair density of Eq.~\eqref{eq:def_n2extrap}:
\begin{equation}
\label{eq:def_pbeueg_ii}
\bar{E}^\Bas_{\pbeontXi} = \int \d\br{} \,\denr \ecmd(\argrpbeontXi),
\end{equation}
%\item[iii)] $\pbeuegns$ which combines a zero spin polarization and the UEG on-top pair density of Eq.~\eqref{eq:def_n2ueg}:
%\begin{equation}
% \label{eq:def_pbeueg_iii}
% \bar{E}^\Bas_{\pbeuegns} = \int \d\br{} \,\denr \ecmd(\argrpbeuegns),
%\end{equation}
\item[iii)] $\pbeontns$ which combines a zero spin polarization and the on-top pair density of Eq.~\eqref{eq:def_n2extrap}:
\begin{equation}
\label{eq:def_pbeueg_iv}
\bar{E}^\Bas_{\pbeontns} = \int \d\br{} \,\denr \ecmd(\argrpbeontns).
\end{equation}
\end{itemize}
The performance of each of these functionals is tested in the following. Note that we did not define a spin-unpolarized version of the PBE-UEG functional because it would have been significantly inferior (in terms of performance) compared to the three other functionals. Indeed, because to the lack of knowledge on the spin polarization or on the accurate on-top pair density, such a functional would be inaccurate in particular for open-shell systems. This assumption has been numerically confirmed by calculations.
%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}
\label{sec:results}
\begin{figure*}
\includegraphics[width=0.45\linewidth]{data/H10/DFT_vdzE_relat.eps}
\includegraphics[width=0.45\linewidth]{data/H10/DFT_vdzE_relat_zoom.eps}
\includegraphics[width=0.45\linewidth]{data/H10/DFT_vtzE_relat.eps}
\includegraphics[width=0.45\linewidth]{data/H10/DFT_vtzE_relat_zoom.eps}
\includegraphics[width=0.45\linewidth]{data/H10/DFT_vqzE_relat.eps}
\includegraphics[width=0.45\linewidth]{data/H10/DFT_vqzE_relat_zoom.eps}
\caption{
Potential energy curves of the H$_{10}$ chain with equally-spaced atoms calculated with MRCI+Q and basis-set corrected MRCI+Q using the cc-pVDZ (top, labelled vdz), cc-pVTZ (middle, labelled as ``vtz''), and cc-pVQZ (bottom, labelled as ``vqz'') basis sets.
The MRCI+Q energies and the estimated exact energies have been extracted from Ref.~\onlinecite{h10_prx}.
% \alert{The reported energy is divided by the number of atoms in the chain.}
\label{fig:H10}}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table*}
\caption{Atomization energies (in mHa) and associated errors (in square brackets) with respect to the estimated exact values computed at different levels of theory with various basis sets.}
\begin{ruledtabular}
\begin{tabular}{lrdddd}
System & \tabc{Basis set} & \tabc{MRCI+Q\fnm[1]} & \tabc{(MRCI+Q)+$\pbeuegXi$} & \tabc{(MRCI+Q)+$\pbeontXi$} & \tabc{(MRCI+Q)+$\pbeontns$} \\
\hline
\ce{H10} & cc-pVDZ & 622.1 [43.3] & 642.6 [22.8] & 649.2 [16.2] & 649.5 [15.9] \\
& cc-pVTZ & 655.2 [10.2] & 661.9 [3.5] & 666.0 [-0.6] & 666.0 [-0.6] \\
& cc-pVQZ & 661.2 [4.2] & 664.1 [1.3] & 666.4 [-1.0] & 666.5 [-1.1] \\[0.1cm]
%\cline{2-6}
& &\multicolumn{4}{c}{Estimated exact:\fnm[1] 665.4} \\[0.2cm]
\hline
% & & \tabc{exFCI} & \tabc{exFCI+$\pbeuegXi$} & \tabc{exFCI+$\pbeontXi$} & \tabc{exFCI+$\pbeontns$}\\
%\hline
%\ce{C2} & aug-cc-pVDZ & 204.6 [29.5] & 218.0 [16.1] & 217.4 [16.7] & 217.0 [17.1] \\
% & aug-cc-pVTZ & 223.4 [10.9] & 228.1 [6.0] & 228.6 [5.5] & 226.5 [5.6] \\[0.1cm]
%\cline{2-6}
% & & \multicolumn{4}{c}{Estimated exact:\fnm[2] 234.1} \\[0.2cm]
%\hline
& & \tabc{exFCI} & \tabc{exFCI+$\pbeuegXi$} & \tabc{exFCI+$\pbeontXi$} & \tabc{exFCI+$\pbeontns$}\\
\hline
\ce{N2} & aug-cc-pVDZ & 321.9 [40.8] & 356.2 [6.5] & 355.5 [7.2] & 354.6 [8.1] \\
& aug-cc-pVTZ & 348.5 [14.2] & 361.5 [1.2] & 363.5 [-0.5] & 363.2 [-0.3] \\
& aug-cc-pVQZ & 356.6 [6.1 ] & 362.8 [-0.1] & 364.2 [-1.5] & 364.3 [-1.6] \\[0.1cm]
& & \multicolumn{4}{c}{Estimated exact:\fnm[2] 362.7} \\[0.2cm]
\hline
& & \tabc{exFCI} & \tabc{exFCI+$\pbeuegXi$} & \tabc{exFCI+$\pbeontXi$} & \tabc{exFCI+$\pbeontns$}\\
\hline
\ce{O2} & aug-cc-pVDZ & 171.4 [20.5] & 187.6 [4.3] & 187.6 [4.3] & 187.1 [4.8] \\
& aug-cc-pVTZ & 184.5 [7.4] & 190.3 [1.6] & 191.2 [0.7] & 191.0 [0.9] \\
& aug-cc-pVQZ & 188.3 [3.6] & 190.3 [1.6] & 191.0 [0.9] & 190.9 [1.0] \\[0.1cm]
& & \multicolumn{4}{c}{Estimated exact:\fnm[2] 191.9} \\[0.2cm]
\hline
& & \tabc{exFCI} & \tabc{exFCI+$\pbeuegXi$} & \tabc{exFCI+$\pbeontXi$} & \tabc{exFCI+$\pbeontns$}\\
\hline
\ce{F2} & aug-cc-pVDZ & 49.6 [12.6] & 54.8 [7.4] & 54.9 [7.3] & 54.8 [7.4] \\
& aug-cc-pVTZ & 59.3 [2.9] & 61.2 [1.0] & 61.5 [0.7] & 61.5 [0.7] \\
& aug-cc-pVQZ & 60.1 [2.1] & 61.0 [1.2] & 61.3 [0.9] & 61.3 [0.9] \\[0.1cm]
\cline{2-6}\\[-0.3cm]
& & \tabc{CEEIS\fnm[3]} & \tabc{CEEIS\fnm[3]+$\pbeuegXi$} & \tabc{CEEIS\fnm[3]+$\pbeontXi$} & \tabc{CEEIS\fnm[3]+$\pbeontns$}\\
\cline{2-6}\\[-0.3cm]
& cc-pVDZ & 43.7 [18.5] & 51.0 [11.2] & 51.0 [11.2] & 50.7 [11.5] \\
& cc-pVTZ & 56.3 [5.9] & 59.2 [3.0] & 59.6 [2.6] & 59.5 [2.7] \\
& cc-pVQZ & 59.9 [2.3] & 61.3 [0.9] & 61.6 [0.6] & 61.6 [0.6] \\[0.1cm]
& & \multicolumn{4}{c}{Estimated exact:\fnm[2] 62.2} \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{From Ref.~\onlinecite{h10_prx}.}
\fnt[2]{From the CEEIS valence-only non-relativistic calculations of Ref.~\onlinecite{BytLaiRuedenJCP05}.}
\fnt[3]{From the CEEIS valence-only non-relativistic calculations of Ref.~\onlinecite{BytNagGorRue-JCP-07}.}
\label{tab:d0}
\end{table*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Computational details}
We present potential energy curves of small molecules up to the dissociation limit
to investigate the performance of the basis-set correction in regimes of both weak and strong correlation.
The considered systems are the \ce{H10} linear chain with equally-spaced atoms, and the \ce{N2}, \ce{O2}, and \ce{F2} diatomics.
The computation of the ground-state energy in Eq.~\eqref{eq:e0approx} in a given basis set requires approximations to the FCI energy $\efci$ and to the basis-set correction $\efuncbasisFCI$.
For diatomics with the aug-cc-pVDZ and aug-cc-pVTZ basis sets,~\cite{KenDunHar-JCP-92} energies are obtained using frozen-core selected-CI calculations (using the CIPSI algorithm) followed by the extrapolation scheme proposed by Holmes \textit{et al.} (see Refs.~\onlinecite{HolUmrSha-JCP-17, SceGarCafLoo-JCTC-18, LooSceBloGarCafJac-JCTC-18, SceBenJacCafLoo-JCP-18, LooBogSceCafJac-JCTC-19, QP2} for more detail). All these calculations are performed with the latest version of \textsc{Quantum Package}, \cite{QP2} and will be labeled as exFCI in the following. In the case of \ce{F2}, we also use the correlation energy extrapolated by intrinsic scaling (CEEIS) \cite{BytNagGorRue-JCP-07} as an estimate of the FCI correlation energy with the cc-pVXZ (X $=$ D, T, and Q) basis sets.~\cite{Dun-JCP-89} The estimated exact potential energy curves are obtained from experimental data \cite{LieCle-JCP-74a} for the \ce{N2} and \ce{O2} molecules, and from CEEIS calculations in the case of \ce{F2}. For all geometries and basis sets, the error with respect to the exact FCI energies are estimated to be of the order of $0.5$~mHa.
For the three diatomics, we performed an additional exFCI calculation with the aug-cc-pVQZ basis set at the equilibrium geometry to obtain reliable estimates of the FCI/CBS dissociation energy.
In the case of the \ce{H10} chain, the approximation to the FCI energies together with the estimated exact potential energy curves are obtained from the data of Ref.~\onlinecite{h10_prx} where the authors performed MRCI+Q calculations with a minimal valence active space as reference (see below for the description of the active space).
Regarding the complementary functional, we first perform full-valence CASSCF calculations with the GAMESS-US software~\cite{gamess} to obtain the wave function $\psibasis$. Then, all density-related quantities involved in the functional [density $n(\br{})$, effective spin polarization $\tilde{\zeta}(\br{})$, reduced density gradient $s(\br{})$, and on-top pair density $n_2(\br{})$] together with the local range-separation function $\mu(\br{})$ are calculated with this full-valence CASSCF wave function. The CASSCF calculations are performed with the following active spaces: (10e,10o) for \ce{H10}, (10e,8o) for \ce{N2}, (12e,8o) for \ce{O2}, and (14e,8o) for \ce{F2}. We note that, instead of using CASSCF wave functions for $\psibasis$, one could of course use the same selected-CI wave functions used for calculating the energy but the calculations of $n_2(\br{})$ and $\mu(\br{})$ would then be more costly. Another strategy would be to use for $\psibasis$ size-consistent truncated versions of the selected-CI wave functions but we did not explore this possibility in this work.
Also, as the frozen-core approximation is used in all our selected-CI calculations, we use the corresponding valence-only complementary functionals (see Subsec.~\ref{sec:FC}). Therefore, all density-related quantities exclude any contribution from the 1s core orbitals, and the range-separation function follows the definition given in Eq.~\eqref{eq:def_mur_val}.
It should be stressed that the computational cost of the basis-set correction (see Appendix~\ref{app:computational}) is negligible compared to the cost the selected-CI calculations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
\includegraphics[width=0.45\linewidth]{data/N2/DFT_avdzE_relat.eps}
\includegraphics[width=0.45\linewidth]{data/N2/DFT_avdzE_relat_zoom.eps}
\includegraphics[width=0.45\linewidth]{data/N2/DFT_avtzE_relat.eps}
\includegraphics[width=0.45\linewidth]{data/N2/DFT_avtzE_relat_zoom.eps}
\caption{
Potential energy curves of the \ce{N2} molecule calculated with exFCI and basis-set corrected exFCI using the aug-cc-pVDZ (top, labelled as ``avdz'') and aug-cc-pVTZ (bottom, labelled as ``avtz'') basis sets. The estimated exact energies are based on a fit of experimental data and obtained from Ref.~\onlinecite{LieCle-JCP-74a}.
\label{fig:N2}}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{H$_{10}$ chain}
The \ce{H10} chain with equally-spaced atoms is a prototype of strongly correlated systems as it consists in the simultaneous breaking of 10 interacting covalent $\sigma$ bonds.
As it is a relatively small system, benchmark calculations at near-CBS values are available (see Ref.~\onlinecite{h10_prx} for a detailed study of this system).
We report in Fig.~\ref{fig:H10} the potential energy curves computed using the cc-pVXZ (X $=$ D, T, and Q) basis sets for different levels of approximation, and the corresponding atomization energies are reported in Table \ref{tab:d0}.
As a general trend, the addition of the basis-set correction globally improves
the quality of the potential energy curves, independently of the approximation level of $\efuncbasis$. Also, no erratic behavior is found when stretching the bonds, which shows that the present procedure (\ie, the determination of the range-separation function and the definition of the functionals) is robust when reaching the strong-correlation regime.
In other words, smooth potential energy curves are obtained with the present basis-set correction.
More quantitatively, the values of the atomization energies are within chemical accuracy (\ie, an error below $1.4$ mHa) with the cc-pVTZ basis set when using the $\pbeontXi$ and $\pbeontns$ functionals, whereas such an accuracy is not yet reached at the standard MRCI+Q/cc-pVQZ level of theory.
Analyzing more carefully the performance of the different types of approximate functionals, the results show that $\pbeontXi$ and $\pbeontns$ are very similar (the maximal difference on the atomization energy being 0.3 mHa), and that they give slightly more accurate results than $\pbeuegXi$. These findings provide two important clues on the role of the different physical ingredients included in these functionals: i) the explicit use of the on-top pair density originating from the CASSCF wave function [see Eq.~\eqref{eq:def_n2extrap}] is preferable over the use of the UEG on-top pair density [see Eq.~\eqref{eq:def_n2ueg}] which is somewhat understandable, and ii) removing the dependence on any kind of spin polarization does not lead to a significant loss of accuracy providing that one employs a qualitatively correct on-top pair density. The latter point is crucial as it confirms that the spin polarization in density-functional approximations essentially plays the same role as the on-top pair density. This could have significant implications for the construction of more robust families of density-functional approximations within DFT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure*}
% \includegraphics[width=0.45\linewidth]{data/C2/DFT_avdzE_relat.pdf}
% \includegraphics[width=0.45\linewidth]{data/C2/DFT_avdzE_relat_zoom.pdf}
% \includegraphics[width=0.45\linewidth]{data/C2/DFT_avtzE_relat.pdf}
% \includegraphics[width=0.45\linewidth]{data/C2/DFT_avtzE_relat_zoom.pdf}
% \caption{
% Potential energy curves of the \ce{C2} molecule calculated with exFCI and basis-set corrected exFCI using the aug-cc-pVDZ (top) and aug-cc-pVTZ (bottom) basis sets. The estimated exact energies are based on fit of experimental data and obtained from Ref.~\onlinecite{LieCle-JCP-74a}.
% \label{fig:C2}}
%\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
\includegraphics[width=0.45\linewidth]{data/O2/DFT_avdzE_relat.eps}
\includegraphics[width=0.45\linewidth]{data/O2/DFT_avdzE_relat_zoom.eps}
\includegraphics[width=0.45\linewidth]{data/O2/DFT_avtzE_relat.pdf}
\includegraphics[width=0.45\linewidth]{data/O2/DFT_avtzE_relat_zoom.pdf}
\caption{
Potential energy curves of the \ce{O2} molecule calculated with exFCI and basis-set corrected exFCI using the aug-cc-pVDZ (top, labelled as ``avdz'') and aug-cc-pVTZ (bottom, labelled as ``avtz'') basis sets. The estimated exact energies are based on a fit of experimental data and obtained from Ref.~\onlinecite{LieCle-JCP-74a}.
\label{fig:O2}}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure*}
% \includegraphics[width=0.45\linewidth]{data/F2/DFT_vdzE_relat.pdf}
% \includegraphics[width=0.45\linewidth]{data/F2/DFT_vdzE_relat_zoom.pdf}
% \includegraphics[width=0.45\linewidth]{data/F2/DFT_vtzE_relat.pdf}
% \includegraphics[width=0.45\linewidth]{data/F2/DFT_vtzE_relat_zoom.pdf}
% \includegraphics[width=0.45\linewidth]{data/F2/DFT_vqzE_relat.pdf}
% \includegraphics[width=0.45\linewidth]{data/F2/DFT_vqzE_relat_zoom.pdf}
% \caption{
% Potential energy curves of the \ce{F2} molecule calculated with CEEIS$^1$ and basis-set corrected CEEIS$^1$ using the cc-pVDZ (top), cc-pVTZ (middle) and cc-pVQZ (bottom) basis sets.
% The estimated exact energies are based on a fit of the valence-only extrapolated CEEIS data extracted from Ref.~\onlinecite{BytNagGorRue-JCP-07}.
% \label{fig:F2}}
%\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
\includegraphics[width=0.45\linewidth]{data/F2/DFT_avdzE_relat.pdf}
\includegraphics[width=0.45\linewidth]{data/F2/DFT_avdzE_relat_zoom.pdf}
\includegraphics[width=0.45\linewidth]{data/F2/DFT_avtzE_relat.pdf}
\includegraphics[width=0.45\linewidth]{data/F2/DFT_avtzE_relat_zoom.pdf}
\caption{
Potential energy curves of the \ce{F2} molecule calculated with exFCI and basis-set corrected exFCI using the aug-cc-pVDZ (top, labelled as ``avdz''), aug-cc-pVTZ (bottom, labelled as ``avtz'') basis sets.
The estimated exact energies are based on a fit of the non-relativistic valence-only CEEIS data extracted from Ref.~\onlinecite{BytNagGorRue-JCP-07}.
\label{fig:F2}}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Dissociation of diatomics}
The \ce{N2}, \ce{O2} and \ce{F2} molecules are complementary to the \ce{H10} system for the present study. The level of strong correlation in these diatomics also increases while stretching the bonds, similarly to the case of \ce{H10}, but in addition these molecules exhibit more important and versatile types of weak correlations due to the larger number of electrons. Indeed, the short-range correlation effects are known to play a strong differential effect on the computation of the atomization energy at equilibrium, while the shape of the curve far from the equilibrium geometry is governed by dispersion interactions which are medium to long-range weak-correlation effects. \cite{AngDobJanGou-BOOK-20} The dispersion interactions in \ce{H10} play a minor role on the potential energy curve due to the much smaller number of near-neighbor electron pairs compared to \ce{N2}, \ce{O2} or \ce{F2}. Also, \ce{O2} has a triplet ground state and is therefore a good candidate for checking the spin-polarization dependence of the various functionals proposed here.
We report in Figs.~\ref{fig:N2}, \ref{fig:O2}, and \ref{fig:F2} the potential energy curves of \ce{N2}, \ce{O2}, and \ce{F2} computed at various approximation levels using the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The atomization energies for each level of theory with different basis sets are reported in Table \ref{tab:d0}.
Just as in \ce{H10}, the accuracy of the atomization energies is globally improved by adding the basis-set correction and it is remarkable that $\pbeontXi$ and $\pbeontns$ provide again very similar results. The latter observation confirms that the dependence on the on-top pair density allows one to remove the dependence of any kind of spin polarization for a quite wide range of covalent bonds and also for an open-shell system like \ce{O2}. More quantitatively, an error below 1.0 mHa compared to the estimated exact valence-only atomization energy is found for \ce{N2}, \ce{O2}, and \ce{F2} with the aug-cc-pVTZ basis set using the $\pbeontns$ functional, whereas such a feat is far from being reached within the same basis set at the near-FCI level. In the case of \ce{F2} it is clear that the addition of diffuse functions in the double- and triple-$\zeta$ basis sets strongly improves the accuracy of the results, which could have be anticipated due to the strong breathing-orbital effect induced by the ionic valence-bond forms in this molecule. \cite{HibHumByrLen-JCP-94}
It should be also noticed that when reaching the aug-cc-pVQZ basis set for \ce{N2}, the accuracy of the atomization energy slightly deteriorates for the $\pbeontXi$ and $\pbeontns$ functionals, but it remains nevertheless more accurate than the estimated FCI atomization energy and very close to chemical accuracy.
Regarding now the performance of the basis-set correction along the whole potential energy curve, it is interesting to notice that it fails to provide a noticeable improvement far from the equilibrium geometry. Acknowledging that the weak-correlation effects in these regions are dominated by dispersion interactions which are long-range effects, the failure of the present approximations for the complementary functional can be understood easily. Indeed, the whole scheme designed here is based on the physics of correlation near the electron-electron coalescence point: the local range-separation function $\mu(\br{})$ is based on the value of the effective electron-electron interaction at coalescence and the ECMD functionals are suited for short-range correlation effects. Therefore, the failure of the present basis-set correction to describe dispersion interactions is theoretically expected.
We hope to report further on this in the near future.
\section{Conclusion}
\label{sec:conclusion}
In the present paper we have extended the recently proposed DFT-based basis-set correction to strongly correlated systems. We have applied the method to the \ce{H10}, \ce{N2}, \ce{O2}, and \ce{F2} molecules up to the dissociation limit at near-FCI level in increasingly large basis sets, and investigated how the basis-set correction affects the convergence toward the CBS limit of the potential energy curves of these molecular systems.
The density-based basis-set correction relies on three aspects: i) the definition of an effective non-divergent electron-electron interaction obtained from the expectation value over a wave function $\psibasis$ of the Coulomb electron-electron interaction projected into an incomplete basis set $\basis$; ii) the fit of this effective interaction with the long-range interaction used in RSDFT; and iii) the use of a short-range, complementary functional borrowed from RSDFT. In the present paper, we investigated i) and iii) in the context of strong correlation and focused on potential energy curves and atomization energies. More precisely, we proposed a new scheme to design functionals fulfilling spin-multiplet degeneracy and size consistency. To fulfill such requirements we proposed to use CASSCF wave functions leading to size-consistent energies, and we developed functionals using only $S_z$-independent density-like quantities.
The development of new $S_z$-independent and size-consistent functionals has lead us to investigate the role of two related quantities: the spin polarization and the on-top pair density. One important result of the present study is that by using functionals \textit{explicitly} depending on the on-top pair density, one can eschew its spin-polarization dependence without loss of accuracy. This avoids the commonly used effective spin polarization originally proposed in Ref.~\onlinecite{BecSavSto-TCA-95} which has the disadvantage of possibly becoming complex-valued in the multideterminant case. From a more fundamental aspect, this confirms that, in a DFT framework, the spin polarization mimics the role of the on-top pair density.
Consequently, we believe that one could potentially develop new families of density-functional approximations where the spin polarization is abandoned and replaced by the on-top pair density.
Regarding the results of the present approach, the basis-set correction systematically improves the near-FCI calculations in a given basis set. More quantitatively, it is shown that with only triple-$\zeta$ quality basis sets chemically accurate atomization energies are obtained for all systems whereas the uncorrected near-FCI results are far from this accuracy within the same basis set.
Also, it is shown that the basis-set correction gives substantial differential contribution to potential energy curves close to the equilibrium geometries, but at long internuclear distances it cannot recover the dispersion interaction energy missing because of the basis-set incompleteness. This behavior is actually expected as dispersion interactions are of long-range nature and the present approach is designed to recover only short-range correlation effects.
\appendix
\section{Size consistency of the basis-set correction}
\label{app:sizeconsistency}
\subsection{Sufficient condition for size consistency}
The basis-set correction is expressed as an integral in real space
\begin{multline}
\label{eq:def_ecmdpbebasisAnnex}
\efuncdenpbe{\argebasis} = \\ \int \text{d}\br{} \,\denr \ecmd(\argrebasis),
\end{multline}
where all the local quantities $\argrebasis$ are obtained from the same wave function $\Psi$. In the limit of two non-overlapping and non-interacting dissociated fragments $\text{A}+\text{B}$, this integral can be rewritten as the sum of the integral over the region $\Omega_\text{A}$ and the integral over the region $\Omega_\text{B}$
\begin{multline}
\label{eq:def_ecmdpbebasisAB}
\efuncdenpbeAB{\argebasis} =
\\
\int_{\Omega_\text{A}} \text{d}\br{} \,\denr \ecmd(\argrebasis)
\\
+ \int_{\Omega_\text{B}} \text{d}\br{} \,\denr \ecmd(\argrebasis).
\end{multline}
Therefore, a sufficient condition to obtain size consistency is that all the local quantities $\argrebasis$ are \textit{intensive}, \ie, they \textit{locally} coincide in the supersystem $\text{A}+\text{B}$ and in each isolated fragment $\text{X}=\text{A}$ or $\text{B}$. Hence, we must have, for $\br{} \in \Omega_\text{X}$,
\begin{subequations}
\begin{gather}
n_\text{A+B}(\br{}) = n_\text{X}(\br{}),
\label{nAB}
\\
\zeta_\text{A+B}(\br{}) = \zeta_\text{X}(\br{}),
\label{zAB}
\\
s_\text{A+B}(\br{}) = s_\text{X}(\br{}),
\label{sAB}
\\
n_{2,\text{A+B}}(\br{}) = n_{2,\text{X}}(\br{}),
\label{n2AB}
\\
\mu_{\text{A+B}}(\br{}) = \mu_{\text{X}}(\br{}),
\label{muAB}
\end{gather}
\end{subequations}
where the left-hand-side quantities are for the supersystem and the right-hand-side quantities for an isolated fragment. Such conditions can be difficult to fulfill in the presence of degeneracies since the system X can be in a different mixed state (\ie, ensemble) in the supersystem $\text{A}+\text{B}$ and in the isolated fragment. \cite{Sav-CP-09} Here, we will consider the simple case where the wave function of the supersystem is multiplicatively separable, \ie,
\begin{equation}
\ket*{\wf{\text{A}+\text{B}}{}} = \ket*{\wf{\text{A}}{}} \otimes \ket*{\wf{\text{B}}{}},
\label{PsiAB}
\end{equation}
where $\otimes$ is the antisymmetric tensor product. In this case, it is easy to shown that Eqs.~(\ref{nAB})-(\ref{sAB}) are valid, as well known, and it remains to show that Eqs.~(\ref{n2AB}) and~(\ref{muAB}) are also valid. Before showing this, we note that even though we do not explicitly consider the case of degeneracies, the lack of size consistency which could arise from spin-multiplet degeneracies can be avoided by the same strategy used for imposing the energy independence on $S_z$, \ie, by using the effective spin polarization $\tilde{\zeta}(n(\br{}),n_{2}(\br{}))$ or a zero spin polarization $\zeta(\br{}) = 0$. Moreover, the lack of size consistency which could arise from spatial degeneracies (\eg, coming from atomic $p$ states) can also be avoided by selecting the same member of the ensemble in the supersystem and in the isolated fragment. This applies to the systems treated in this work.
\subsection{Intensivity of the on-top pair density and the local range-separation function}
The on-top pair density can be written in an orthonormal spatial orbital basis set $\{\SO{p}{}\}$ as
\begin{equation}
\label{eq:def_n2}
n_{2{}}(\br{}) = \sum_{pqrs \in \Bas} \SO{p}{} \SO{q}{} \Gam{pq}{rs} \SO{r}{} \SO{s}{},
\end{equation}
with $\Gam{pq}{rs} = 2 \mel*{\wf{}{}}{ \aic{r_\downarrow}\aic{s_\uparrow}\ai{q_\uparrow}\ai{p_\downarrow}}{\wf{}{}}$. As the summations run over all orbitals in the basis set $\Bas$, $n_{2{}}(\br{})$ is invariant to orbital rotations and can thus be expressed in terms of localized orbitals. For two non-overlapping fragments $\text{A}+\text{B}$, the basis set can then be partitioned into orbitals localized on the fragment A and orbitals localized on B, \ie, $\Bas=\Bas_\text{A}\cup \Bas_\text{B}$. Therefore, we see that the on-top pair density of the supersystem $\text{A}+\text{B}$ is additively separable
\begin{equation}
n_{2,\text{A}+\text{B}}(\br{}) = n_{2,\text{A}}(\br{}) + n_{2,\text{B}}(\br{}),
\end{equation}
where $n_{2,\text{X}}(\br{})$ is the on-top pair density of the fragment X
\begin{equation}
n_{2,\text{X}}(\br{}) = \sum_{pqrs \in \Bas_\text{X}} \SO{p}{} \SO{q}{} \Gam{pq}{rs} \SO{r}{} \SO{s}{},
\end{equation}
in which the elements $\Gam{pq}{rs}$ with orbital indices restricted to the fragment X are $\Gam{pq}{rs} = 2 \mel*{\wf{\text{A}+\text{B}}{}}{ \aic{r_\downarrow}\aic{s_\uparrow}\ai{q_\uparrow}\ai{p_\downarrow}}{\wf{\text{A}+\text{B}}{}} = 2 \mel*{\wf{\text{X}}{}}{ \aic{r_\downarrow}\aic{s_\uparrow}\ai{q_\uparrow}\ai{p_\downarrow}}{\wf{\text{X}}{}}$, owing to the multiplicative structure of the wave function [see Eq.~\eqref{PsiAB}]. This shows that the on-top pair density is a local intensive quantity.
The local range-separation function is defined as, for $n_{2}(\br{}) \not=0$,
\begin{equation}
\label{eq:def_murAnnex}
\mur = \frac{\sqrt{\pi}}{2} \frac{f(\bfr{},\bfr{})}{n_{2}(\br{})},
\end{equation}
where
\begin{equation}
\label{eq:def_f}
f(\bfr{},\bfr{}) = \sum_{pqrstu\in \Bas} \SO{p}{ } \SO{q}{ } \V{pq}{rs} \Gam{rs}{tu} \SO{t}{ } \SO{u}{ }.
\end{equation}
Again, $f(\bfr{},\bfr{})$ is invariant to orbital rotations and can be expressed in terms of orbitals localized on the fragments A and B. In the limit of infinitely separated fragments, the Coulomb interaction vanishes between A and B and therefore any two-electron integral $\V{pq}{rs}$ involving orbitals on both A and B vanishes. We thus see that the quantity $f(\bfr{},\bfr{})$ of the supersystem $\text{A}+\text{B}$ is additively separable
\begin{equation}
\label{eq:def_fa+b}
f_{\text{A}+\text{B}}(\bfr{},\bfr{}) = f_{\text{A}}(\bfr{},\bfr{}) + f_{\text{B}}(\bfr{},\bfr{}),
\end{equation}
with
\begin{equation}
\label{eq:def_fX}
f_\text{X}(\bfr{},\bfr{}) = \sum_{pqrstu\in \Bas_\text{X}} \SO{p}{ } \SO{q}{ } \V{pq}{rs} \Gam{rs}{tu} \SO{t}{ } \SO{u}{ }.
\end{equation}
So, $f(\bfr{},\bfr{})$ is a local intensive quantity.
As a consequence, the local range-separation function of the supersystem $\text{A}+\text{B}$ is
\begin{equation}
\label{eq:def_murAB}
\mu_{\text{A}+\text{B}}(\bfr{}) = \frac{\sqrt{\pi}}{2} \frac{f_{\text{A}}(\bfr{},\bfr{}) + f_{\text{B}}(\bfr{},\bfr{})}{n_{2,\text{A}}(\br{}) + n_{2,\text{B}}(\br{})},
\end{equation}
which implies
\begin{equation}
\label{eq:def_murABsum}
\mu_{\text{A}+\text{B}}(\bfr{}) = \mu_{\text{X}}(\bfr{}) \;\; \text{if} \;\; \bfr{} \in \Omega_\text{X},
\end{equation}
where $\mu_{\text{X}}(\bfr{}) = (\sqrt{\pi}/2) f_{\text{X}}(\bfr{},\bfr{})/n_{2,\text{X}}(\br{})$.
The local range-separation function is thus a local intensive quantity.
We can therefore conclude that, if the wave function of the supersystem $\text{A}+\text{B}$ is multiplicative separable, all local quantities used in the basis-set correction functional are intensive and therefore the basis-set correction is size consistent.
\section{Computational cost of the basis-set correction for a CASSCF wave function}
\label{app:computational}
The computational cost of the basis-set correction is determined by the calculation of the on-top pair density $n_{2}(\br{})$ and the local range-separation function $\mur$ on the real-space grid. For a general multideterminant wave function, the computational cost is of order $O(N_\text{grid}N_{\Bas}^4)$ where $N_\text{grid}$ is the number of grid points and $N_{\Bas}$ is the number of basis functions.\cite{LooPraSceTouGin-JCPL-19} For a CASSCF wave function, a significant reduction of the scaling of the computational cost can be achieved.
\subsection{Computation of the on-top pair density}
For a CASSCF wave function $\Psi$, the occupied orbitals can be partitioned into a set of active orbitals $\mathcal{A}$ and a set of inactive (doubly occupied) orbitals $\mathcal{I}$. The CASSCF on-top pair density can then be written as
\begin{equation}
\label{def_n2_good}
n_{2}(\br{}) = n_{2,\mathcal{A}}(\br{}) + n_{\mathcal{A}}(\br{}) n_{\mathcal{I}}(\br{}) + \frac{n_{\mathcal{I}}(\br{})^2}{2},
\end{equation}
where
\begin{subequations}
\begin{align}
\label{def_n2_act}
n_{2,\mathcal{A}}(\br{}) & = \sum_{pqrs \in \mathcal{A}} \SO{p}{} \SO{q}{} \Gam{pq}{rs} \SO{r}{} \SO{s}{},
\\
n_{\mathcal{A}}(\br{}) & = \sum_{pq \in\mathcal{A}} \phi_p (\br{}) \phi_q (\br{})
\mel*{\wf{}{}}{ \aic{p_\uparrow}\ai{q_\uparrow} + \aic{p_\downarrow}\ai{q_\downarrow} }{\wf{}{}},
\\
n_{\mathcal{I}}(\br{}) & = 2 \sum_{p \in \mathcal{I}} \phi_p (\br{})^2
\end{align}
\end{subequations}
are the purely active part of the on-top pair density, the active part of the density, and the inactive part of the density, respectively.
The leading computational cost is the evaluation of $n_{2,\mathcal{A}}(\br{})$ on the grid which, according to Eq.~\eqref{def_n2_act}, scales as $O(N_\text{grid} N_\mathcal{A}^4)$ where $N_{\mathcal{A}}$ is the number of active orbitals which is much smaller than the number of basis functions $N_{\Bas}$.
\subsection{Computation of the local range-separation function}
In addition to the on-top pair density, the computation of $\mur$ needs the computation of $f(\bfr{},\bfr{})$ [see Eq.~\eqref{eq:def_f}] at each grid point. It can be factorized as
\begin{equation}
\label{eq:f_good}
f(\bfr{},\bfr{}) = \sum_{rs \in \Bas} V^{rs}(\bfr{}) \, \Gamma_{rs}(\bfr{}),
\end{equation}
where
\begin{subequations}
\begin{align}
V^{rs}(\bfr{}) & = \sum_{pq \in \Bas} V_{pq}^{rs} \phi_p(\br{}) \phi_q(\br{}),
\\
\Gamma_{rs}(\bfr{}) & = \sum_{pq \in \Bas} \Gam{rs}{pq} \phi_p(\br{}) \phi_q(\br{}) .
\end{align}
\end{subequations}
For a general multideterminant wave function, the computational cost of $f(\bfr{},\bfr{})$ thus scales as $O(N_\text{grid}N_{\Bas}^4)$.
In the case of a CASSCF wave function, $\Gam{rs}{pq}$ vanishes if one index $p,q,r,s$ does not belong to the set of inactive or active occupied orbitals $\mathcal{I}\cup \mathcal{A}$. Therefore, at a given grid point, the number of non-zero elements $\Gamma_{rs}(\bfr{})$ is only at most $(N_{\mathcal{I}}+N_{\mathcal{A}})^2$, which is usually much smaller than $N_{\Bas}^2$. As a consequence, one can also restrict the sum in the calculation of
\begin{equation}
f(\bfr{},\bfr{}) = \sum_{rs \in \mathcal{I}\cup\mathcal{A}} V^{rs}(\bfr{}) \, \Gamma_{rs}(\bfr{}).
\end{equation}
The overall computational cost is dominated by that of $V^{rs}(\bfr{})$, which scales as $O(N_\text{grid}(N_{\mathcal{I}}+N_{\mathcal{A}})^2 N_{\Bas}^2)$, which is much smaller than $O(N_\text{grid}N_{\Bas}^4)$.
\bibliography{srDFT_SC}
\end{document}

753
Revised_Manuscript/srDFT_SC_old.tex
View File

@ -1,753 +0,0 @@
\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}

136
response/response.tex
View File

@ -29,132 +29,66 @@
\sffamily \sffamily
\section*{\sffamily Curing basis-set convergence of wave-function theory using density-functional theory: a systematically improvable approach} \section*{\sffamily A basis-set error correction based on density-functional theory for strongly correlated molecular systems}
We would like to thank the reviewers for their carefull reading of our manuscript, and we reply in the present document to their remarks and criticisms. We would like to thank the reviewers for their carefull reading of our manuscript, and we reply in the present document to their remarks and criticisms.
\subsection*{Comments of reviewer 1 and reply} \subsection*{Comments of reviewer 1 and reply}
\subsubsection*{Comments of reviewer 1} \subsubsection*{Comments of reviewer 1}
The authors propose a method to account for short‐range electron correlation using DFT techniques. Reviewer 1 Evaluations:\\
The idea makes sense. I have long been frustrated (and intrigued) by the fact that DFT handles short‐ Recommendation: Optional revision\\
range correlation effortlessly while wavefunction theory (WFT) requires a large basis set and a lot of New Potential Energy Surface: No\\
computational work. The main problem of all combined WFT/DFT methods is double counting. The Overall Rating (required): Top 5-25\% - significant, novel, and impactful contribution of broad interest \\
authors partially eliminate it by replacing the Coulomb operator with its long‐range component, tailored
to the basis set (which can describe correlation only on the scale of the distance between the closest The manuscript presents an extension of the density-functional-based basis-set
nodes and larger.) correction scheme, by some of the current authors, to strongly correlated systems. This
I have only a few comments. basis-set correction scheme is a beautiful and powerful combination of density-functional
\begin{itemize} and wave-function theoretical concepts and practices.
\item The method proposed is simpler than using explicit interelectronic (F12) coordinates. I wonder how All the basis-set corrections are applied to high-quality approximations of FCI
much the latter cost computationally. The theory is certainly complicated and there are still a few loose results that are computationally costly. A natural question arises, what would be the
ends in R12 or F12 theory but does it slow down the calculations significantly? The authors could add a effect of these basis-set corrections if the computational method were of lower quality.
comment on this. F12 has the advantage that it eliminates double counting rigorously but I am not sure The reader would benefit, if there were at least a discussion of this issue in the
how expensive it is. The method described in the paper is virtually free on the cost scale of a large manuscript.
correlated calculation but adds some of the disadvantages of DFT, mainly noise from numerical A few very minor remarks:\\
quadrature. - The reference numbering is not always consecutive\\
\item Page 1, right column, bottom. The authors repeat the misleading claim that “DFT abandons the - “have be” $\to$ “have been” (in the 4 th sentence of the 3 rd paragraph of IIIC)\\
complex ... wavefunction for the simple one‐body density.” This is valid for orbital‐free DFT which is not - In the caption of Fig. 4, the last two sentences are seemingly contradictory.
yet working. It should be toned down.
\item P. 3, above equ. (4) change “give” to “gives”.
\item P. 4, right column, line 12 from below. The phrase “projected in a basis” is perhaps better if replaced
by “projected to a basis”.
\item I suspect that the CIPSI energy for N + at the AQZ level in Table II (‐54.020414) has a problem. The
energy difference between QMC and CIPSI is generally in the 5 th digit, as pointed out in the paper. For
this figure, the deviation is in the 3 rd digit. Please check it again.
\item I suggest that the authors switch to the IUPAC (and IUPAP) abbreviation for the atomic unit of energy,
E h (and mE h ).
\item In its current form, the method is only applicable to atoms, as the local range separation parameter
depends on the distance from the nucleus. The authors should indicate how it could be generalized to
molecules.
\end{itemize}
\subsection*{Reply to reviewer 1} \subsection*{Reply to reviewer 1}
\begin{itemize} The basis-set correction can indeed be applied to any approximation to FCI such as multireference perturbation theory, similarly to what was done for weakly correlated systems in P.-F. Loos, B. Pradines, A. Scemama, J. Toulouse, and E. Giner, J. Phys. Chem. Lett. 10, 2931 (2019) in which the basis-set correction was applied to CCSD(T) calculations. As for weakly correlated systems, we expect a similar speed-up of the basis convergence and we will report on this in a forthcoming paper. We have added the following sentence in the Computational Details section:\\
\item[] "The main problem of all combined WFT/DFT methods is double counting. The
authors partially eliminate it by replacing the Coulomb operator with its long‐range component, tailored
to the basis set (which can describe correlation only on the scale of the distance between the closest
nodes and larger.)" \\
In the case of standard range-separated DFT (RS-DFT), there is by construction no double counting of the electronic correlation effects as it relies on a clean splitting of the interaction, which is not the case for most of the combined WFT/DFT methods. The practical problems in RS-DFT are the approximated functionals or the ansatz used for the wave function part. Nevertheless, in the present work, as the effective interaction in the basis set is not perfectly represented by the fit that we proposed here (see equation (30) page 7), some double counting can still be present.
\item "The method proposed is simpler than using explicit interelectronic (F12) coordinates. I wonder how
much the latter cost computationally. The theory is certainly complicated and there are still a few loose
ends in R12 or F12 theory but does it slow down the calculations significantly? The authors could add a
comment on this. F12 has the advantage that it eliminates double counting rigorously but I am not sure
how expensive it is. The method described in the paper is virtually free on the cost scale of a large
correlated calculation but adds some of the disadvantages of DFT, mainly noise from numerical
quadrature.
" \\
The authors are aware of several F12 theories but are not expert in that field, so we preferred not to comment too much on the computational cost of F12, specially on the latest developments. Regarding the avoidance of double counting, the reviewer could refer to the very recent study on the practical effect of the truncation of the one-particle basis set in approximated CCSD-F12 approaches: {\it We note that, even though we use near-FCI energies in this work, the DFT-based basis-set correction could also be applied to any approximation to FCI such as multireference perturbation theory, similarly to what was done for weakly correlated systems for which the basis-set correction was applied to CCSD(T) calculations [51].}
M. K. Kesharwani, N Sylvetsky, A. K\"ohn, D. P. Tew, and J. M. L. Martin, Journal of Chemical Physics 149, 154109 (2018).
Nonetheless, even if our approach is conceptually simpler than F12 theory, as it is at its early stage of development, it is hard to compare the computational cost with F12 theories considering the amount of efforts involved in that field.
\item "Page 1, right column, bottom. The authors repeat the misleading claim that “DFT abandons the
complex ... wavefunction for the simple one‐body density.” This is valid for orbital‐free DFT which is not
yet working. It should be toned down. "\\
From a purely formal point of view, we agree that only orbital-free DFT abandons the wave function, as the Kohn-Sham Slater determinant is still needed for the computation of the kinetic part in KS-DFT.
Nonetheless, from a WFT perspective (which is the overall point of view of the paper), KS-DFT is drastically simpler than WFT as it does not rely on the accurate description of the complex and non local two-body density matrix. This is why we stated that "DFT abandons the complex ... wave function for the simple one-body density", and also the sentence in page 1 does not refer to a specific variant of DFT as it refers to the formal implication of the Hohenberg-Kohn theorem.
\item "P. 3, above equ. (4) change “give” to “gives”." \\
The typo was corrected.
\item "P. 4, right column, line 12 from below. The phrase “projected in a basis” is perhaps better if replaced
by “projected to a basis”. \\
We thank the author for the stylistic proposal. After inquiry, we found that the most adequate sentence is "projected on" (which is equivalent to "projected onto") which was changed in the text.
\item "I suspect that the CIPSI energy for N + at the AQZ level in Table II (‐54.020414) has a problem. The
energy difference between QMC and CIPSI is generally in the 5 th digit, as pointed out in the paper. For
this figure, the deviation is in the 3 rd digit. Please check it again." \\
We also found strange the discrepancy between our CIPSI energy and the \textit{i}-FCIQMC result in the case of N$^+$ in the aug-cc-pVQZ basis set. Nonetheless, we double checked our result and pushed our calculations quite far (since $|E_{PT2}| < 0.1 mH$) and obtained essentially the same energy and a variational energy lower than the one in \textit{i}-FCIQMC. As a consequence, considering that the initiator approximation of FCIQMC necessary introduces a bias and that we obtained a lower variational energy than the \textit{i}-FCIQMC, we believe that the problem is on the side of \textit{i}-FCIQMC.
\item "I suggest that the authors switch to the IUPAC (and IUPAP) abbreviation for the atomic unit of energy,
E h (and mE h )."\\
We believe that the present abbreviations for the atomic unit of the energies reported here do not introduce any inconvenience for the clarity of the exposition, therefore we keep the present abbreviations.
\item "In its current form, the method is only applicable to atoms, as the local range separation parameter
depends on the distance from the nucleus. The authors should indicate how it could be generalized to
molecules." \\
I think that the reviewer 1 misunderstood some aspects of the procedure proposed here, as there are no restrictions in its mathematical formulation for the special cases of atomic systems.
Indeed, the definition of the effective electron-electron interaction projected on a basis (see equations (21) to (27) page 5, or appendix A for a more detailed derivation) relies on an expectation value of the Coulomb interaction over a very general wave function which can represent an atom or a molecule. The derived effective interaction (see equation 27) is a general function from ${\rm I\!R}^{6}\rightarrow {\rm I\!R}$, which can be computed in any physical system.
We chose to study atomic systems as these systems are basically free from static correlation effects, exhibit large basis set errors for the calculations of energy differences and are well understood from a physical point of view (which enable us for instance the detailed analysis of the last section). Incoming works will assess the validity of this approach on molecular systems.
\end{itemize}
We have also corrected the three minor problems seen by the referee.
\subsection*{Comments of reviewer 2 and reply} \subsection*{Comments of reviewer 2 and reply}
\subsubsection*{Comments of reviewer 2} \subsubsection*{Comments of reviewer 2}
Reviewer 2 Evaluations:\\ Reviewer 2 Evaluations:\\
Recommendation: Optional revision\\ Recommendation: Optional revision\\
New Potential Energy Surface: No\\ New Potential Energy Surface: No\\
Overall Rating (required): Top 5 \% - highly significant, novel, and impactful contribution of broad interest\\
Reviewer 2 (Comments to the Author):\\ This is a very high-quality manuscript, in which a DFT-based correction for finite basis-set errors in wavefunction theory is extended to the case of strongly correlated molecular systems, by looking at molecular dissociation curves. The beauty of the theory behind this DFT basis-set correction is that it can be applied to any wave function method. The authors use their very smart idea (published in previous works) of a suitable mapping of the interaction projected in the finite basis to a long-range only interaction (erf(mu*r)/r), determining the value of the local range separated parameter mu by matching the interaction at r=0.
Here the authors further extend the DFT part of their theory by using explictilyt the on-top (OT) pair density dependence of the short-range functionals of multideterminant range separated DFT to correct the basis set error. They also show that the OT dependence can fully eliminate the spin-polarization dependence.
The idea is very neat and well explained, and the results are excellent. I recommend publication with the minor revisions below.
In the present manuscript by Giner et al. the authors present a novel method to correct for basis set incompleteness errors in Full CI (FCI) calculations. Their method is based on the local mapping of an effective (basis set dependent) electron-electron interaction to a range separation (RS) parameter as used in RS density functionals. By combining the proposed local mapping and the derived formalism it is possible to introduce a perturbative basis set correction for FCI calculations. The obtained results for spherical systems demonstrate convincingly that their novel method achieves a rapid basis set convergence. 1) The authors should comment more the results of Table 1 for atomization energies. Although their theory is guaranteed to converge to the exact result for the complete basis set limit, the results of table 1 show that for atomization energies this does not always occur in a monotinic way. There are cases (H$_{10}$ and N$_2$ with the OT functionals, O$_2$ with all the functionals and F$_2$ with augmented basis with all the functionals) in which the results are worse for a quadruple-zeta basis than for triple zeta. Could this be due to the fact that when augmenting the basis the OT pair density from the CASSCF wavefunction does not necessarily improve but might oscillate around the Coulomb cusp creating artefacts in the extrapolated OT? After all the basis sets are optimized on the energy and do not really care about the cusp. Or is the effect due to a different accuracy of the correction for the molecule and the atoms? It would be nice if the authors could say something on that.
The manuscript is very well written and the work is of great interest to the electronic structure theory community. I recommend publication of the manuscript.
I only have a few optional suggestions and questions: 2) Appendix A 1: eqs. A3, as discussed by the authors may break down for the case of degenearcy of the isolated fragments. The authors comment that the spin degeneracy can be avoided by considering the functionals without spin polarization and that the spatial degeneracy can be also cured by considering for the fragment the same element of the ensemble as for the suprasystem.
\begin{itemize} However I still wonder if this is ok. For the spin I agree with the authors, but for the spatial degeneracy I still have some doubts. Could the authors comment at least qualitatively on what would happen for C$_2$ and B$_2$ for example? Is the DFT-based basis set correction going to give the same result for the spherically-averaged C or B atom than for the atoms with the occupied p orbital oritented alsong the bond axis? Or does the DFT-based correction show on this the same problems of standard DFT?
\item The present work does not draw any comparison to F12-like approaches in the results section. I suggest for this or future studies to compare results obtained using the present approach to results obtained using FCI + [2]R12 (as proposed by Kong and Valeev [J. Chem. Phys. 135, 214105 (2011)]), or FCI + CT (as proposed in Yanai and Shiozaki [J. Chem. Phys. 136, 084107 (2012)]). It seems that the present approach is computationally simpler and more efficient than F12-like approaches (please correct me if I'm wrong). Therefore it would be interesting to compare these approaches to each other in terms of accuracy.
\item The authors have investigated atoms only. Can the present approach be directly transferred to non-spherical systems, or do the authors expect difficulties? Presenting results for binding curves would be desirable for this or future studies.
\item Typo:\\
p.5 : "which means that its does not..." -> "which means that it does not..."
\end{itemize}
\subsection*{Reply to reviewer 2} \subsection*{Reply to reviewer 2}
\begin{itemize} \begin{itemize}
\item "The present work does not draw any comparison to F12-like approaches in the results section. I suggest for this or future studies to compare results obtained using the present approach to results obtained using FCI + [2]R12 (as proposed by Kong and Valeev [J. Chem. Phys. 135, 214105 (2011)]), or FCI + CT (as proposed in Yanai and Shiozaki [J. Chem. Phys. 136, 084107 (2012)]). It seems that the present approach is computationally simpler and more efficient than F12-like approaches (please correct me if I'm wrong). Therefore it would be interesting to compare these approaches to each other in terms of accuracy." \\ \item[1]
We completely agree with the reviewer 2 that our approach is directly comparable to those proposed by Kong and Valeev, as both approaches can be seen as a post-WFT calculation which aims to account for the incompleteness of the basis set used in quantum chemistry. More generally, we indeed plan to perform a comparative study with F12 approaches which will help us to better understand the strengths and limitations of our new approach. \\
Regarding the method proposed by Yanai and Shiozaki, it is, to our understanding, not straightforwardly comparable to the present work as it deals with a two-body Hamiltonian which is modified by the presence of Slater-geminals. Therefore, it introduces a coupling between the correlation effects in the basis set (at whatever level of treatment selected) and those introduced by the explicit correlation factor. It is then more of a "perturb then diagonalize" approach, whether our approach if closer to a "diagonalize then perturb", as the one of Kong and Valeev. Current developments will try to make the procedure self consistent by carrying the minimization over all possible densities instead of doing the approximation of the FCI density as it is done here (see equation (14)).
\item "The authors have investigated atoms only. Can the present approach be directly transferred to non-spherical systems, or do the authors expect difficulties? Presenting results for binding curves would be desirable for this or future studies." \\ \item[2] Our DFT-based basis set correction does not generally preserve spatial degeneracy for arbitrary states or ensembles. So, in that regard, it has the same problem than standard DFT. Specifically, for the example of C$_2$ and B$_2$, it will not give the same result for the spherically-averaged C or B atom as for the atoms with the occupied p orbital oritented along the bond axis. However, what we write in the Appendix is that for the systems treated in this work the CASSCF wave function dissociates into fragments in simple identified pure states and we can thus choose to perform the calculation on the isolated atom with the same pure state in order to preserve size-consistency. Of course, this may not be always possible for other more complicated systems. We have clarified this point in the Appendix. The new paragraph is now:\\
As mentioned to reviewer 1, we investigated atomic systems only in order to focus our attention on the slow convergence of the dynamical correlation with the basis set (no static correlation involved). Also, these systems are quite well understood and the spherical symmetry has allowed us to perform a simple physical analysis of the behavior of our approach (the last section of the paper). {\it Moreover, for the systems treated in this work, the lack of size consistency which could arise from spatial degeneracies (coming from atomic $p$ states) can also be avoided by selecting the same state in the supersystem and in the isolated fragment. For example, for the F$_2$ molecule, the CASSCF wave function dissociates into the atomic configuration $\text{p}_\text{x}^2 \text{p}_\text{y}^2 \text{p}_\text{z}^1$ for each fragment, and we thus choose the same configuration for the calculation of the isolated atom. The same argument applies to the N$_2$ and O$_2$ molecules. For other systems, it may not be always possible to do so.}
Of course, the approach can be transferred to non-spherical systems since the effective interaction is derived from the expectation value of the Coulomb interaction on a general wave function which can represent molecular systems, and since the functionals used in RS-DFT are perfectly suited for the treatment of molecules. We have already performed calculations on homo- and hetero-nuclear molecules (atomization energies of the G2 sets of molecules) and the presentation and analysis of the results will be the subject of a forthcoming paper.
\item Typo:\\
p.5 : "which means that its does not..." -> "which means that it does not..." \\
The typo was fixed.
\end{itemize} \end{itemize}
\end{document} \end{document}