put SI as an appendix

Manuscript/srDFT_SC.tex

@@ -71,6 +71,7 @@
 \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]}
@@ -142,6 +143,7 @@
 % 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}}}
@@ -597,8 +599,7 @@ The second approach to approximate the exact on-top pair density consists in usi
 \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. When using $\ntwoextrap(\ntwo,\mu)$ in a functional, we will simply refer it as ``OT'', which stands for "on-top". 
 
-We then define four complementary functionals:
+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}:
@@ -613,18 +614,19 @@ We then define four complementary functionals:
 	\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}:
+%\item[iii)] $\pbeuegns$ which combines a zero spin polarization and the UEG on-top pair density of Eq.~\eqref{eq:def_n2ueg}:
-\begin{equation}
+%\begin{equation}
-	\label{eq:def_pbeueg_iii}
+%	\label{eq:def_pbeueg_iii}
-	\bar{E}^\Bas_{\pbeuegns}  =  \int  \d\br{} \,\denr  \ecmd(\argrpbeuegns),
+%	\bar{E}^\Bas_{\pbeuegns}  =  \int  \d\br{} \,\denr  \ecmd(\argrpbeuegns),
-\end{equation}
+%\end{equation}
-\item[iv)] $\pbeontns$ which combines a zero spin polarization and the on-top pair density of Eq.~\eqref{eq:def_n2extrap}:
+\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 four functionals is tested in the following.
+The performance of each of these functionals is tested in the following. Notice that we did not define a spin-unpolarized PBE-UEG functional which would be significantly inferior to the three other functionals. Indeed, without knowledge of the spin polarization or the on-top pair density, such a functional would be inaccurate for the hydrogen atom for example.
+
 
 
 %DFT: BLACK BOX and not CASSCF 
@@ -729,7 +731,7 @@ More quantitatively, the values of $D_0$ are within chemical accuracy (\ie, an e
 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 $D_0$ 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 used in these functionals: i) the explicit use of the on-top pair density coming from the CASSCF wave function [Eq.~\eqref{eq:def_n2extrap}] is preferable over the use of the UEG on-top pair density [Eq.~\eqref{eq:def_n2ueg}] which is somewhat understandable, and ii) removing the dependency on any kind of spin polarization does not lead to significant loss of accuracy providing that one employs a qualitatively correct on-top pair density. The latter point is crucial as it shows 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. 
 
-Finally, the reader would have noticed that we did not report the performance of the $\pbeuegns$ functional as its performance are significantly inferior than the three other functionals. The main reason behind this comes from the fact that $\pbeuegns$ has no direct or indirect knowledge of the on-top pair density of the system. Therefore, it yields a non-zero correlation energy for the totally dissociated \ce{H10} chain even if the on-top pair density is vanishingly small. This necessary lowers the value of $D_0$. Therefore, from hereon, we simply discard the $\pbeuegns$ functional.
+%Finally, the reader would have noticed that we did not report the performance of the $\pbeuegns$ functional as its performance are significantly inferior than the three other functionals. The main reason behind this comes from the fact that $\pbeuegns$ has no direct or indirect knowledge of the on-top pair density of the system. Therefore, it yields a non-zero correlation energy for the totally dissociated \ce{H10} chain even if the on-top pair density is vanishingly small. This necessary lowers the value of $D_0$. Therefore, from hereon, we simply discard the $\pbeuegns$ functional.
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -821,6 +823,168 @@ Regarding the results of the present approach, the basis-set correction systemat
 
 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 energy missing because of the basis set incompleteness. This behaviour is actually expected as dispersion is 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}
+
+\subsection{Sufficient condition for size consistency}
+
+The basis-set correction is expressed as an integral in real space
+\begin{equation}
+ \label{eq:def_ecmdpbebasis}
+ \begin{aligned}
+ & \efuncdenpbe{\argebasis} = \\ & \int  \text{d}\br{} \,\denr   \ecmd(\argrebasis),
+ \end{aligned}
+\end{equation}
+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{equation}
+ \label{eq:def_ecmdpbebasis}
+ \begin{aligned}
+ & \efuncdenpbeAB{\argebasis} = \\ & \int_{\Omega_\text{A}}  \text{d}\br{} \,\denr   \ecmd(\argrebasis) \\ & + \int_{\Omega_\text{B}}  \text{d}\br{} \,\denr   \ecmd(\argrebasis).
+ \end{aligned}
+\end{equation}
+Therefore, a sufficient condition to obtain size consistency is that all the local quantities $\argrebasis$ are \textit{intensive}, i.e. that they \textit{locally} coincide in the supersystem $\text{A}+\text{B}$ and in each isolated fragment $\text{X}=\text{A}$ or $\text{B}$. Hence, for $\br{} \in \Omega_\text{X}$, we should have
+\begin{subequations}
+\begin{equation}
+n_\text{A+B}(\br{}) = n_\text{X}(\br{}),
+\label{nAB}
+\end{equation}
+\begin{equation}
+\zeta_\text{A+B}(\br{}) = \zeta_\text{X}(\br{}),
+\label{zAB}
+\end{equation}
+\begin{equation}
+s_\text{A+B}(\br{}) = s_\text{X}(\br{}),
+\label{sAB}
+\end{equation}
+\begin{equation}
+n_{2,\text{A+B}}(\br{}) = n_{2,\text{X}}(\br{}),
+\label{n2AB}
+\end{equation}
+\begin{equation}
+\mu_{\text{A+B}}(\br{}) = \mu_{\text{X}}(\br{}),
+\label{muAB}
+\end{equation}
+\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 fulfil in the presence of degeneracies since the system X can be in a different mixed state (i.e. 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, i.e.
+\begin{equation}
+ \ket{\wf{\text{A}+\text{B}}{}} = \ket{\wf{\text{A}}{}} \otimes \ket{\wf{\text{B}}{}}, 
+\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 explicity 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$, i.e. 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 (e.g., coming from atomic p states) can also be avoided by selecting the same member of the ensemble in the supersystem and in the isolated fragement. This applies to the systems treated in this work.
+
+\subsection{Intensivity of the on-top pair density and of the local range-separation parameter}
+
+The on-top pair density can be written in an orthonormal spatial orbital basis $\{\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 partitioned into orbitals localized on the fragment A and orbitals localized on B, i.e. $\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}
+ \label{eq:def_n2}
+ 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}
+ \label{eq:def_n2}
+ n_{2,\text{X}}(\br{})  = \sum_{pqrs \in \Bas_\text{X}} \SO{p}{} \SO{q}{} \Gam{pq}{rs}  \SO{r}{} \SO{s}{}. 
+\end{equation}
+This shows that the on-top pair density is a local intensive quantity.
+
+The local range-separation parameter is defined by 
+\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 parameter 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 gives
+\begin{equation}
+ \label{eq:def_murABsum}
+ \mu_{\text{A}+\text{B}}(\bfr{}) = \mu_{\text{A}}(\bfr{}) + \mu_{\text{B}}(\bfr{}),
+\end{equation}
+with $\mu_{\text{X}}(\bfr{}) = (\sqrt{\pi}/2) f_{\text{X}}(\bfr{},\bfr{})/n_{2,\text{X}}(\br{})$. The local range-separation parameter is thus a local intensive quantity.
+
+In conclusion, 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 considerations}
+
+The computational cost of the present approach is driven by two quantities: the computation of the on-top pair density and the $\murpsibas$ on the real-space grid. Within a blind approach, for each grid point the computational cost is of order $n_{\Bas}^4$ and $n_{\Bas}^6$ for the on-top pair density $n_{2,\wf{\Bas}{}}(\br{})$ and the local range separation parameter $\murpsibas$, respectively. 
+Nevertheless, using CASSCF wave functions to compute these quantities leads to significant simplifications which can substantially reduce the CPU time. 
+
+\subsection{Computation of the on-top pair density for a CASSCF wave function}
+
+Given a generic wave function developed on a basis set with $n_{\Bas}$ basis functions, the evaluation of the on-top pair density is of order $\left(n_{\Bas}\right)^4$. 
+Nevertheless, assuming that the wave function $\Psi^{\Bas}$ is of CASSCF type, a lot of simplifications happen. 
+If the active space is referred as the set of spatial orbitals $\mathcal{A}$ which are labelled by the indices $t,u,v,w$, and the doubly occupied orbitals are the set of spatial orbitals $\mathcal{C}$ labeled by the indices $i,j$, one can write the on-top pair density of a CASSCF wave function as
+\begin{equation}
+ \label{def_n2_good}
+ n_{2,\wf{\Bas}{}}(\br{})  = n_{2,\mathcal{A}}(\br{}) + n_{\mathcal{C}}(\br{}) n_{\mathcal{A}}(\br{}) + \left( n_{\mathcal{C}}(\br{})\right)^2
+\end{equation}
+where 
+\begin{equation}
+ \label{def_n2_act}
+ n_{2,\mathcal{A}}(\br{}) = \sum_{t,u,v,w \, \in \mathcal{A}} 2 \mel*{\wf{}{\Bas}}{ \aic{t_\downarrow}\aic{u_\uparrow}\ai{v_\uparrow}\ai{w_\downarrow}}{\wf{}{\Bas}} \phi_t (\br{})  \phi_u (\br{})  \phi_v (\br{})  \phi_w (\br{})
+\end{equation} 
+is the purely active part of the on-top pair density, 
+\begin{equation}
+ n_{\mathcal{C}}(\br{}) = \sum_{i\, \in \mathcal{C}} \left(\phi_i (\br{}) \right)^2, 
+\end{equation}
+and 
+\begin{equation}
+ n_{\mathcal{A}}(\br{}) = \sum_{t,u\, \in \mathcal{A}} \phi_t (\br{})  \phi_u (\br{}) 
+ \mel*{\wf{}{\Bas}}{ \aic{t_\downarrow}\ai{u_\downarrow}  +   \aic{t_\uparrow}\ai{u_\uparrow}}{\wf{}{\Bas}} 
+\end{equation}
+is the purely active one-body density. 
+Written as in eq. \eqref{def_n2_good}, the leading computational cost is the evaluation of $n_{2,\mathcal{A}}(\br{})$ which, according to eq. \eqref{def_n2_act}, scales as $\left( n_{\mathcal{A}}\right) ^4$ where $n_{\mathcal{A}}$ is the number of active orbitals which is much smaller than the number of basis functions $n_{\Bas}$. Therefore, the final computational scaling of the on-top pair density for a CASSCF wave function over the whole real-space grid is of $\left( n_{\mathcal{A}}\right) ^4 n_G$, where  $n_G$ is the number of grid points. 
+
+\subsection{Computation of $\murpsibas$}
+At a given grid point, the computation of $\murpsibas$ needs the computation of $f_{\wf{}{}}(\bfr{},\bfr{}) $ defined in eq. \eqref{eq:def_f} and the on-top pair density defined in eq. \eqref{eq:def_n2}. In the previous paragraph we gave an explicit form of the on-top pair density in the case of a CASSCF wave function with a computational scaling of $\left( n_{\mathcal{A}}\right)^4$. In the present paragraph we focus on simplifications that one can obtain for the computation of $f_{\wf{}{}}(\bfr{},\bfr{}) $ in the case of a CASSCF wave function. 
+
+One can rewrite $f_{\wf{}{}}(\bfr{},\bfr{}) $ as 
+\begin{equation}
+ \label{eq:f_good}
+ f_{\wf{}{}}(\bfr{},\bfr{})  =  \sum_{r,s \in \Bas} \mathcal{V}_r^s(\bfr{}) \, \mathcal{N}_{r}^s(\bfr{})
+\end{equation}
+where 
+\begin{equation}
+ \mathcal{V}_r^s(\bfr{}) = \sum_{p,q \in \Bas} V_{pq}^{rs} \phi_p(\br{}) \phi_q(\br{})
+\end{equation}
+and 
+\begin{equation}
+ \mathcal{N}_{r}^s(\bfr{}) = \sum_{p,q \in \Bas} \Gam{pq}{rs} \phi_p(\br{}) \phi_q(\br{}) .
+\end{equation}
+\textit{A priori}, for a given grid point, the computational scaling of eq. \eqref{eq:f_good} is of $\left(n_{\Bas}\right)^4$ and the total computational cost over the whole grid scales therefore as $\left(n_{\Bas}\right)^4 n_G$. 
+
+In the case of a CASSCF wave function, it is interesting to notice that $\Gam{pq}{rs}$ vanishes if one index $p,q,r,s$ does not belong 
+to the set of of doubly occupied or active orbitals $\mathcal{C}\cup \mathcal{A}$. Therefore, at a given grid point, the matrix $\mathcal{N}_{r}^s(\bfr{})$ has only at most $\left(n_{\mathcal{A}}+n_{\mathcal{C}}\right)^2$ non-zero elements, which is usually much smaller than $\left(n_{\Bas}\right)^2$. 
+As a consequence, in the case of a CASSCF wave function one can rewrite $f_{\wf{}{}}(\bfr{},\bfr{})$ as 
+\begin{equation}
+ f_{\wf{}{}}(\bfr{},\bfr{}) = \sum_{r,s \in \mathcal{C}\cup\mathcal{A}} \mathcal{V}_r^s(\bfr{}) \, \mathcal{N}_{r}^s(\bfr{}). 
+\end{equation}
+Therefore the final computational cost of $f_{\wf{}{}}(\bfr{},\bfr{})$ is dominated by that of $\mathcal{V}_r^s(\bfr{})$, which scales as $\left(n_{\mathcal{A}}+n_{\mathcal{C}}\right)^2 \left( n_{\Bas} \right)^2 n_G$, which is much weaker than $\left(n_{\Bas}\right)^4 n_G$. 
+
 \bibliography{srDFT_SC}
 
 \end{document}