Done with 1st draft

Manuscript/FarDFT.bib

@@ -1,13 +1,37 @@
 %% This BibTeX bibliography file was created using BibDesk.
 %% http://bibdesk.sourceforge.net/
 
-%% Created for Pierre-Francois Loos at 2020-04-10 15:29:54 +0200 
+%% Created for Pierre-Francois Loos at 2020-04-10 22:33:34 +0200 
 
 
 %% Saved with string encoding Unicode (UTF-8) 
 
 
 
+@article{Loos_2020a,
+	Author = {P. F. Loos and A. Scemama and D. Jacquemin},
+	Date-Added = {2020-04-10 22:11:02 +0200},
+	Date-Modified = {2020-04-10 22:11:55 +0200},
+	Doi = {10.1021/acs.jpclett.9b03652},
+	Journal = {J. Phys. Chem. Lett.},
+	Pages = {974},
+	Title = {The Quest for Highly-Accurate Excitation Energies: a Computational Perspective},
+	Volume = {11},
+	Year = {2020},
+	Bdsk-Url-1 = {https://doi.org/10.1021/acs.jpclett.9b03652}}
+
+@article{Loos_2020b,
+	Author = {P. F. Loos and F. Lipparini and M. Boggio-Pasqua and A. Scemama and D. Jacquemin},
+	Date-Added = {2020-04-10 22:09:54 +0200},
+	Date-Modified = {2020-04-10 22:15:37 +0200},
+	Doi = {10.1021/acs.jctc.9b01216},
+	Journal = {J. Chem. Theory Comput.},
+	Pages = {1711},
+	Title = {A mountaineering strategy to excited states: highly-accurate energies and benchmarks for medium size molecules,},
+	Volume = {16},
+	Year = {2020},
+	Bdsk-Url-1 = {https://doi.org/10.1021/acs.jctc.9b01216}}
+
 @article{Hait_2020,
 	Author = {D. Hait and M. Head-Gordon},
 	Date-Added = {2020-04-10 15:18:47 +0200},
@@ -17,17 +41,20 @@
 	Pages = {1699--1710},
 	Title = {Excited state orbital optimization via minimizing the square of the gradient: General approach and application to singly and doubly excited states via density functional theory},
 	Volume = {16},
-	Year = {2020}}
+	Year = {2020},
+	Bdsk-Url-1 = {https://doi.org/10.1021/acs.jctc.9b01127}}
 
 @article{Madden_1963,
 	Author = {R. P. Madden and K. Codling},
 	Date-Added = {2020-04-10 15:13:19 +0200},
-	Date-Modified = {2020-04-10 15:14:37 +0200},
+	Date-Modified = {2020-04-10 22:33:34 +0200},
 	Doi = {10.1103/PhysRevLett.10.516},
 	Journal = {Phys. Rev. Lett.},
+	Pages = {516},
 	Title = {New Autoionizing Atomic Energy Levels in He, Ne, and Ar},
 	Volume = {10},
-	Year = {1963}}
+	Year = {1963},
+	Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevLett.10.516}}
 
 @article{Becke_1988a,
 	Author = {A. D. Becke},
@@ -5123,12 +5150,14 @@
 @article{Loos_2018,
 	Author = {P. F. Loos and A. Scemama and A. Blondel and Y. Garniron and M. Caffarel and D. Jacquemin},
 	Date-Added = {2018-10-24 22:38:52 +0200},
-	Date-Modified = {2018-10-24 22:38:52 +0200},
+	Date-Modified = {2020-04-10 22:19:17 +0200},
+	Doi = {10.1021/acs.jctc.8b00406},
 	Journal = {J. Chem. Theory Comput.},
 	Pages = {4360},
 	Title = {A Mountaineering Strategy to Excited States: Highly-Accurate Reference Energies and Benchmarks},
 	Volume = {14},
-	Year = {2018}}
+	Year = {2018},
+	Bdsk-Url-1 = {https://doi.org/10.1021/acs.jctc.8b00406}}
 
 @article{Lorentzon_1995,
 	Author = {Lorentzon, Johan and Malmqvist, Per-Ake and Fiilscher, Markus},

Manuscript/FarDFT.tex

@@ -134,10 +134,11 @@
 	\affiliation{\LCPQ}
 
 \begin{abstract}
-Gross-Oliveira-Kohn (GOK) ensemble density-functional theory (GOK-DFT) is a time-independent formalism which allows to compute excitation energies via the derivative of the ensemble energy with respect to the weight of each excited state.
+Gross--Oliveira--Kohn (GOK) ensemble density-functional theory (GOK-DFT) is a time-\textit{independent} formalism which allows to compute excitation energies via the derivative of the ensemble energy with respect to the weight of each excited state.
 Contrary to the time-dependent version of density-functional theory (TD-DFT), double excitations can be easily computed within GOK-DFT.
-However, to take full advantage of this formalism, one must have access to a \textit{weight-dependent} exchange-correlation functional in order to model the infamous derivative discontinuity contribution to the excitation energies.
+However, to take full advantage of this formalism, one must have access to a \textit{weight-dependent} exchange-correlation functional in order to model the infamous ensemble derivative contribution to the excitation energies.
 In the present article, we discuss the construction of first-rung (\ie, local) weight-dependent exchange-correlation density-functional approximations for two-electron atomic and molecular systems (He and H$_2$) specifically designed for the computation of double excitations within GOK-DFT.
+A specific protocol is proposed to obtain accurate energies associated with double excitations.
 \end{abstract}
 
 \maketitle
@@ -146,7 +147,7 @@ In the present article, we discuss the construction of first-rung (\ie, local) w
 %%% INTRODUCTION %%%
 %%%%%%%%%%%%%%%%%%%%
 \section{Introduction}
-Time-dependent density-functional theory (TD-DFT) has been the dominant force in the calculation of excitation energies of molecular systems in the last two decades.\cite{Casida,Ulrich_2012}
+Time-dependent density-functional theory (TD-DFT) has been the dominant force in the calculation of excitation energies of molecular systems in the last two decades.\cite{Casida,Ulrich_2012,Loos_2020a}
 At a relatively low computational cost (at least compared to the other excited-state \textit{ab initio} methods), TD-DFT can provide accurate transition energies for low-lying excited states of organic molecules (see, for example, Ref.~\onlinecite{Dreuw_2005} and references therein).
 Importantly, setting up a TD-DFT calculation for a given system is an almost pain-free process from a user perspective as the only (yet essential) input variable is the choice of the so-called exchange-correlation (xc) functional.
 
@@ -161,7 +162,7 @@ For example, TD-DFT has problems with charge-transfer \cite{Tozer_1999,Dreuw_200
 The development of range-separated hybrids provides an effective solution to this problem. \cite{Tawada_2004,Yanai_2004} 
 From a practical point of view, the TD-DFT xc kernel is usually considered as static instead of being frequency dependent.
 One key consequence of this so-called adiabatic approximation (based on the assumption that the density varies slowly with time) is that double excitations are completely absent from the TD-DFT spectra. \cite{Levine_2006,Tozer_2000,Elliott_2011}
-Although these double excitations are usually experimentally dark (which means they usually cannot be observed in photo-absorption spectroscopy), these states play, indirectly, a key role in many photochemistry mechanisms. \cite{Boggio-Pasqua_2007,Loos_2019}
+Although these double excitations are usually experimentally dark (which means they usually cannot be observed in photo-absorption spectroscopy), these states play, indirectly, a key role in many photochemistry mechanisms. \cite{Boggio-Pasqua_2007} They are, moreover, a real challenge for high-level computational methods. \cite{Loos_2018,Loos_2019,Loos_2020b}
 
 One possible solution to access double excitations within TD-DFT is provided by spin-flip TD-DFT which describes double excitations as single excitations from the lowest triplet state. \cite{Huix-Rotllant_2010,Krylov_2001,Shao_2003,Wang_2004,Wang_2006,Minezawa_2009}
 However, spin contamination might be an issue. \cite{Huix-Rotllant_2010}
@@ -170,8 +171,8 @@ In this approach the xc kernel is made frequency dependent, which allows to trea
 
 Maybe surprisingly, another possible way of accessing double excitations is to resort to a time-\textit{independent} formalism. \cite{Yang_2017,Sagredo_2018,Deur_2019}
 With a computational cost similar to traditional KS-DFT, DFT for ensembles (eDFT) \cite{Theophilou_1979,Gross_1988a,Gross_1988b,Oliveira_1988} is a viable alternative following such a strategy currently under active development.\cite{Gidopoulos_2002,Franck_2014,Borgoo_2015,Kazaryan_2008,Gould_2013,Gould_2014,Filatov_2015,Filatov_2015b,Filatov_2015c,Gould_2017,Deur_2017,Gould_2018,Gould_2019,Sagredo_2018,Ayers_2018,Deur_2018,Deur_2019,Kraisler_2013,Kraisler_2014,Alam_2016,Alam_2017,Nagy_1998,Nagy_2001,Nagy_2005,Pastorczak_2013,Pastorczak_2014,Pribram-Jones_2014,Yang_2013a,Yang_2014,Yang_2017,Senjean_2015,Senjean_2016,Smith_2016,Senjean_2018}
-In the assumption of monotonically decreasing weights, eDFT for excited states has the undeniable advantage to be based on a rigorous variational principle for ground and excited states,  the so-called Gross-Oliveria-Kohn (GOK) variational principle.  \cite{Gross_1988a} 
+In the assumption of monotonically decreasing weights, eDFT for excited states has the undeniable advantage to be based on a rigorous variational principle for ground and excited states,  the so-called Gross--Oliveria--Kohn (GOK) variational principle.  \cite{Gross_1988a} 
-In short, GOK-DFT (\ie, eDFT for excited states) is the density-based analog of state-averaged wave function methods, and excitation energies can then be easily extracted from the total ensemble energy. \cite{Deur_2019}
+In short, GOK-DFT (\ie, eDFT for neutral excitations) is the density-based analog of state-averaged wave function methods, and excitation energies can then be easily extracted from the total ensemble energy. \cite{Deur_2019}
 Although the formal foundations of GOK-DFT have been set three decades ago, \cite{Gross_1988a,Gross_1988b,Oliveira_1988} its practical developments have been rather slow.
 We believe that it is partly due to the lack of accurate approximations for GOK-DFT.
 In particular, to the best of our knowledge, an explicitly weight-dependent density-functional approximation for ensemble (eDFA) has never been developed for atoms and molecules.
@@ -182,7 +183,7 @@ The LDA, as we know it, is based on the uniform electron gas (UEG) also known as
 Although the Hohenberg--Kohn theorems \cite{Hohenberg_1964} are here to provide firm theoretical grounds to DFT, modern KS-DFT rests largely on the presumed similarity between this hypothetical UEG and the electronic behaviour in a real system. \cite{Kohn_1965}
 However, Loos and Gill have recently shown that there exists other UEGs which contain finite numbers of electrons (more like in a molecule), \cite{Loos_2011b,Gill_2012} and that they can be exploited to construct LDA functionals. \cite{Loos_2014a,Loos_2014b,Loos_2017a}
 Electrons restricted to remain on the surface of a $\cD$-sphere (where $\cD$ is the dimensionality of the surface of the sphere) are an example of finite UEGs (FUEGs). \cite{Loos_2011b}
-In particular, we combine these FUEGs with the usual infinite UEG (IUEG) to construct a weigh-dependent LDA correlation functional for ensembles, which is specifically designed to compute double excitations within GOK-DFT, and automatically incorporates the infamous derivative discontinuity contribution to the excitation energies through its explicit ensemble weight dependence. \cite{Levy_1995, Perdew_1983}
+In particular, we combine these FUEGs with the usual infinite UEG (IUEG) to construct a weigh-dependent LDA correlation functional for ensembles, which is specifically designed to compute double excitations within GOK-DFT, and automatically incorporates the infamous ensemble derivative contribution to the excitation energies through its explicit ensemble weight dependence. \cite{Levy_1995, Perdew_1983}
 
 The paper is organised as follows.
 In Sec.~\ref{sec:theo}, the theory behind GOK-DFT is presented.
@@ -225,7 +226,7 @@ Turning to GOK-DFT, the extension of the Hohenberg--Kohn theorem to ensembles al
 	\E{}{\bw} = \min_{\n{}{}} \qty{ \F{}{\bw}[\n{}{}] + \int \vne(\br{}) \n{}{}(\br{}) d\br{} },
 \end{equation}
 where $\vne(\br{})$ is the external potential and $\F{}{\bw}[\n{}{}]$ is the universal ensemble functional
-(the weight-dependent analog of the Hohenberg-Kohn universal functional for ensembles).
+(the weight-dependent analog of the Hohenberg--Kohn universal functional for ensembles).
 In the KS formulation, this functional is decomposed as
 \begin{equation}
 	\F{}{\bw}[\n{}{}] 
@@ -412,19 +413,7 @@ We shall come back to this point later on.
 
 Third, we add up correlation effects via the VWN5 local correlation functional. \cite{Vosko_1980}
 For the sake of clarity, the explicit expression of the VWN5 functional is not reported here but it can be found in Ref.~\onlinecite{Vosko_1980}.
-The combination of the Slater and VWN5 functionals (SVWN5) yield a highly convex ensemble energy (as shown in Fig.~\ref{fig:Ew_H2}), while the combination of the GIC-S and VWN5 functionals (GIC-SVWN5) exhibit a small curvature (the ensemble energy is slightly concave) and improved excitation energies, especially at small weights, where the SVWN5 excitation energy is almost spot on. 
+The combination of the Slater and VWN5 functionals (SVWN5) yield a highly convex ensemble energy (as shown in Fig.~\ref{fig:Ew_H2}), while the combination of GIC-S and VWN5 (GIC-SVWN5) exhibit a small curvature (the ensemble energy is slightly concave) and improved excitation energies, especially at small weights, where the GIC-SVWN5 excitation energy is almost spot on. 
-
-%%%%%%%%%%%%%%%%%%
-%%% FUNCTIONAL %%%
-%%%%%%%%%%%%%%%%%%
-%\section{Functional}
-%\label{sec:func}
-%The present work deals with the explicit construction of the (reduced) LDA xc functional $\e{\xc}{\bw}[\n{}{}]$ defined in Eq.~\eqref{eq:exc_def}.
-%The generalisation to a larger number of states (in particular the inclusion of the first singly-excited state) is trivial and left for future work.
-
-%The construction of these two functionals is described below.
-%Extension to spin-polarised systems will be reported in future work.
-
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsubsection{Weight-dependent correlation functional}
@@ -448,49 +437,7 @@ The reduced (\ie, per electron) Hartree-Fock (HF) energies for these two states
 	\label{eq:eHF_1}
 \end{align}
 \end{subequations}
-%These two energies can be conveniently decomposed as 
+
-%\begin{equation}
-%	\e{\HF}{(I)}(\n{}{}) = \kin{s}{(I)}(\n{}{}) + \e{\Ha}{(I)}(\n{}{}) + \e{\ex}{(I)}(\n{}{}),
-%\end{equation}
-%with
-%\begin{subequations}
-%\begin{align}
-%	\kin{s}{(0)}(\n{}{}) & = 0,
-%	&
-%	\kin{s}{(1)}(\n{}{}) & = \frac{3\pi^{2}}{2} \qty(\frac{\n{}{}}{\pi})^{2/3},
-%	\\
-%	\e{\Ha}{(0)}(\n{}{}) & = \frac{8}{3} \qty(\frac{\n{}{}}{\pi})^{1/3},
-%	&
-%	\e{\Ha}{(1)}(\n{}{}) & = \frac{352}{105} \qty(\frac{\n{}{}}{\pi})^{1/3},
-%	\\
-%	\e{\ex}{(0)}(\n{}{}) & = - \frac{4}{3} \qty(\frac{\n{}{}}{\pi})^{1/3},
-%	&
-%	\e{\ex}{(1)}(\n{}{}) & = - \frac{176}{105} \qty(\frac{\n{}{}}{\pi})^{1/3}.	
-%\end{align}
-%\end{subequations}
-%
-%In analogy with the conventional Dirac exchange functional, \cite{Dirac_1930} we write down the exchange functional of each individual state as
-%\begin{equation}
-%	\e{\ex}{(I)}(\n{}{}) = \Cx{(I)} \n{}{1/3},
-%\end{equation}
-%and we then obtain
-%\begin{align}
-%	 \Cx{(0)} & = - \frac{4}{3} \qty( \frac{1}{\pi} )^{1/3},
-%	 &
-%	 \Cx{(1)} & = - \frac{176}{105} \qty( \frac{1}{\pi} )^{1/3}.
-%\end{align}
-%We can now combine these two exchange functionals to create a weight-dependent exchange functional for a two-state ensemble
-%\begin{equation}
-%\label{eq:exw}
-%	\e{\ex}{\ew{}}(\n{}{}) 
-%	= (1-\ew{}) \e{\ex}{(0)}(\n{}{}) + \ew{} \e{\ex}{(1)}(\n{}{})
-%	= \Cx{\ew{}} \n{}{1/3}
-%\end{equation}
-%with
-%\begin{equation}
-%	\Cx{\ew{}} = (1-\ew{}) \Cx{(0)} + \ew{} \Cx{(1)}.
-%\end{equation}
-%Conveniently, the weight dependence of the exchange functional can be transferred to the $\Cx{}$ coefficient; an expected feature from a theoretical point of view, yet a nice property from a more practical aspect.
 Thanks to highly-accurate calculations and the expressions of the HF energies provided by Eqs.~\eqref{eq:eHF_0} and \eqref{eq:eHF_1}, \cite{Loos_2009a,Loos_2009c,Loos_2010e} one can write down, for each state, an accurate analytical expression of the reduced correlation energy \cite{Loos_2013a, Loos_2014a} via the following Pad\'e approximant \cite{Sun_2016,Loos_2020}
 \begin{equation}
 \label{eq:ec}
@@ -624,7 +571,7 @@ They can then be obtained via GOK-DFT ensemble calculations by performing a line
 The results gathered in Table \ref{tab:BigTab_H2} show that the GOK-DFT excitation energies obtained with the GIC-SeVWN5 functional at zero weight are the most accurate with an improvement of $0.25$ eV as compared to GIC-SVWN5, which is due to the ensemble derivative contribution of the eVWN5 functional.
 The GIC-SeVWN5 excitation energies at equi-weights (\ie, $\ew{} = 1/2$) are less satisfactory, but still remains in good agreement with FCI, with again a small improvement as compared to GIC-SVWN5.
 It is also important to mention that the GIC-S functional does not alter the MOM excitation energy as the correction vanishes accordingly for $\ew{} = 1$ (\textit{vide supra}).
-Finally, note that, although we had to design a system-specific, weight-dependent exchange functional to reach such accuracy, we have not used any high-level reference data (such as FCI) to tune our functional, the only requirement being the linearity of the ensemble energy (obtained with the Slater exchange functional) between $\ew{} = 0$ and $1$.
+Finally, note that, although we had to design a system-specific, weight-dependent exchange functional to reach such accuracy, we have not used any high-level reference data (such as FCI) to tune our functional, the only requirement being the linearity of the ensemble energy (obtained with LDA exchange) between $\ew{} = 0$ and $1$.
 
 %%% TABLE I %%%
 \begin{table}
@@ -696,6 +643,7 @@ For this particular geometry, the doubly-excited state becomes the lowest excite
 We then follow the same protocol as in Sec.~\ref{sec:H2}, and considering again the aug-cc-pVTZ basis set, we design a GIC-S functional for this system at $\RHH = 3.7$ bohr.
 It yields $\alpha = +0.019\,226$, $\beta = -0.017\,996$, and $\gamma = -0.022\,945$ [see Eq.~\eqref{eq:Cxw}].
 The weight dependence of $\Cx{\ew{}}$ is illustrated in Fig.~\ref{fig:Cxw} (green curve). 
+
 One clearly sees that the correction brought by GIC-S is much more gentle than at $\RHH = 1.4$ bohr, which means that the ensemble energy obtained with the LDA exchange functional is much more linear at $\RHH = 3.7$ bohr.
 In other words, the ghost-interaction ``hole'' depicted in Fig.~\ref{fig:Cxw} is thus much more shallow at stretched geometry.
 Table \ref{tab:BigTab_H2st} reports, for the aug-cc-pVTZ basis set (which delivers basis set converged results), the same set of calculations as in Table \ref{tab:BigTab_H2}.
@@ -704,7 +652,7 @@ For $\RHH = 3.7$ bohr, it is much harder to get an accurate estimate of the exci
 The GIC-S functional coupled or not with a correlation functional yield extremely stable excitation energies as a function of the weight, with only a few tenths of eV difference between the zero- and equi-weights limits.
 Nonetheless, the excitation energy is still off by 3 eV.
 The fundamental theoretical reason of such a poor agreement is not clear.
-The fact that HF exchange yields better excitation energy hints at the effect of self-interaction error.
+The fact that HF exchange yields better excitation energies hints at the effect of self-interaction error.
 
 %%% TABLE I %%%
 \begin{table}
@@ -747,9 +695,9 @@ Excitation energies (in eV) associated with the lowest double excitation of \ce{
 
 As a final example, we consider the \ce{He} atom which can be seen as the limiting form of the \ce{H2} molecule for very short bond lengths.
 In \ce{He}, the lowest doubly-excited state is an auto-ionising resonance state, extremely high in energy and lies in the continuum. \cite{Madden_1963}
-In Ref.~\onlinecite{Burges_1995}, highly-accurate calculations estimate an excitation energy of $2.126$ hartree.
+In Ref.~\onlinecite{Burges_1995}, highly-accurate calculations estimate an excitation energy of $2.126$ hartree for this $1s^2 \rightarrow 2s^2$ transition.
 Nonetheless, it can be nicely described with a Gaussian basis set containing enough diffuse functions.
-Consequently,  we considered for this particular example the d-aug-cc-pVQZ basis set which contains two sets of diffuse functions.
+Consequently,  we consider for this particular example the d-aug-cc-pVQZ basis set which contains two sets of diffuse functions.
 The excitation energies associated with this double excitation computed with various methods and combinations of xc functionals are gathered in Table \ref{tab:BigTab_He}.
 
 The parameters of the GIC-S weight-dependent exchange functional (computed with the smaller aug-cc-pVTZ basis) are $\alpha = +1.912\,574$, $\beta = +2.715\,267$, and $\gamma = +2.163\,422$ [see Eq.~\eqref{eq:Cxw}], the curvature of the ensemble energy being more pronounced in \ce{He} than in \ce{H2} (blue curve in Fig.~\ref{fig:Cxw}).
@@ -758,7 +706,7 @@ The results reported in Table \ref{tab:BigTab_He} evidence this strong weight de
 
 The GIC-S exchange functional attenuates significantly this dependence, and when coupled with the eVWN5 weight-dependent correlation functional, the GIC-SeVWN5 excitation energy for $\ew{} = 0$ is only $8$ millihartree (or $0.22$ eV) off the reference value.
 As in the case of \ce{H2}, the excitation energies obtained at zero-weight are more accurate than at equi-weight.
-As a final comment, let us stress that the present protocol does not rely on high-level calculations as the sole requirement for constructing the GIC-S functional is the linearity of the ensemble energy with respect to the weight of the double excitation.
+As a final comment, let us stress again that the present protocol does not rely on high-level calculations as the sole requirement for constructing the GIC-S functional is the linearity of the ensemble energy with respect to the weight of the double excitation.
 
 
 %%% TABLE I %%%