revising manuscript

This commit is contained in:
Pierre-Francois Loos 2020-04-27 23:04:16 +02:00
parent 96baa3f9c5
commit 2bd087c89f
6 changed files with 4124 additions and 6736 deletions

10423
FarDFT.nb

File diff suppressed because it is too large Load Diff

Binary file not shown.

Binary file not shown.

View File

@ -84,6 +84,8 @@
\newcommand{\xc}{\text{xc}}
% matrices
\newcommand{\bO}{\boldsymbol{0}}
\newcommand{\bthird}{\boldsymbol{1/3}}
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bw}{\boldsymbol{w}}
\newcommand{\bHc}{\boldsymbol{h}}
@ -92,6 +94,7 @@
% elements
\newcommand{\ew}[1]{w_{#1}}
\newcommand{\eW}{\xi}
\newcommand{\eHc}[1]{h_{#1}}
\newcommand{\eJ}[1]{J_{#1}}
\newcommand{\eK}[1]{K_{#1}}
@ -106,6 +109,8 @@
% Ao and MO basis
\newcommand{\MO}[2]{\phi_{#1}^{#2}}
\newcommand{\HOMO}[1]{\MO{\text{HOMO}}{#1}}
\newcommand{\LUMO}[1]{\MO{\text{LUMO}}{#1}}
\newcommand{\cMO}[2]{c_{#1}^{#2}}
\newcommand{\AO}[1]{\chi_{#1}}
@ -143,15 +148,15 @@
\affiliation{\LCPQ}
\begin{abstract}
\titou{Gross--Oliveira--Kohn (GOK) ensemble density-functional theory (GOK-DFT)
Gross--Oliveira--Kohn (GOK) ensemble density-functional theory (GOK-DFT)
is a time-\textit{independent} extension of density-functional theory (DFT) which
allows to compute excited-state
energies via the derivatives of the ensemble energy with
respect to the ensemble weights.}
respect to the ensemble weights.
Contrary to the time-dependent version of DFT (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 ensemble derivative contribution to the excitation energies.
In the present article, we discuss the construction of first-rung (\textit{i.e.}, 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.
\titou{In the spirit of optimally-tuned range-separated hybrid functionals,} a specific protocol is proposed to obtain accurate energies associated with double excitations.
In the spirit of optimally-tuned range-separated hybrid functionals, a two-step system-dependent procedure is proposed to obtain accurate energies associated with double excitations.
\end{abstract}
\maketitle
@ -162,26 +167,26 @@ In the present article, we discuss the construction of first-rung (\textit{i.e.}
\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_1995,Ulrich_2012,Loos_2020a}
At a moderate 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).
\titou{Importantly, within the widely-used adiabatic approximation, setting up a TD-DFT calculation for a given system is an
Importantly, within the widely-used adiabatic approximation, 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
ground-state exchange-correlation (xc) functional.}
ground-state exchange-correlation (xc) functional.
Similar to density-functional theory (DFT), \cite{Hohenberg_1964,Kohn_1965,ParrBook} TD-DFT is an in-principle exact theory which formal foundations relie on the Runge-Gross theorem. \cite{Runge_1984}
The Kohn-Sham (KS) formulation of TD-DFT transfers the
complexity of the many-body problem to the xc functional thanks to a
judicious mapping between a time-dependent non-interacting reference
system and its interacting analog \titou{which have both
exactly the same one-electron density.}
system and its interacting analog which have both
exactly the same one-electron density.
\titou{However, TD-DFT is far from being perfect as, in practice, drastic approximations must be made.
However, TD-DFT is far from being perfect as, in practice, drastic approximations must be made.
First, within the linear-response approximation, the electronic spectrum relies on the (unperturbed) pure-ground-state KS picture, \cite{Runge_1984, Casida_1995, Casida_2012} which may not be adequate in certain situations (such as strong correlation).
Second, the time dependence of the functional is usually treated at the local approximation level within the standard adiabatic approximation.
In other words, memory effects are absent from the xc functional which is assumed to be local in time
(the xc energy is in fact an xc action, not an energy functional). \cite{Vignale_2008}
Third and more importantly in the present context, a major issue of TD-DFT actually originates directly from the choice of the xc functional, and more specifically, the possible (not to say likely) substantial variations in the quality of the excitation energies for two different choices of xc functionals.}
Third and more importantly in the present context, a major issue of TD-DFT actually originates directly from the choice of the xc functional, and more specifically, the possible (not to say likely) substantial variations in the quality of the excitation energies for two different choices of xc functionals.
\titou{Because its popularity, approximate TD-DFT has been studied in excruciated details by the community, and some researchers have quickly unveiled various theoretical and practical deficiencies.}
Because its popularity, approximate TD-DFT has been studied in excruciated details by the community, and some researchers have quickly unveiled various theoretical and practical deficiencies.
For example, TD-DFT has problems with charge-transfer \cite{Tozer_1999,Dreuw_2003,Sobolewski_2003,Dreuw_2004,Maitra_2017} and Rydberg \cite{Tozer_1998,Tozer_2000,Casida_1998,Casida_2000,Tozer_2003} excited states (the excitation energies are usually drastically underestimated) due to the wrong asymptotic behaviour of the semi-local xc functional.
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.
@ -199,17 +204,17 @@ In the assumption of monotonically decreasing weights, eDFT for excited states h
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.
\titou{In particular, to the best of our knowledge, albeit several attempts have been made, \cite{Nagy_1996,Paragi_2001} an explicitly
In particular, to the best of our knowledge, although several attempts have been made, \cite{Nagy_1996,Paragi_2001} an explicitly
weight-dependent density-functional approximation for ensembles (eDFA)
has never been developed for atoms and molecules from first principles.
The present contribution paves the way towards this goal.}
The present contribution paves the way towards this goal.
When one talks about constructing functionals, the local-density
approximation (LDA) is never far away.
%\manu{too ``oral'' style I think}. Let's be fun Manu!
The LDA, as we know it, is based on the uniform electron gas (UEG) also known as jellium, an hypothetical infinite substance where an infinite number of electrons ``bathe'' in a (uniform) positively-charged jelly. \cite{Loos_2016}
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 \titou{ground-state functionals as shown in Refs.~\onlinecite{Loos_2014a,Loos_2014b,Loos_2017a}, where the authors proposed generalised LDA exchange and correlation functionals.}
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 ground-state functionals as shown in Refs.~\onlinecite{Loos_2014a,Loos_2014b,Loos_2017a}, where the authors proposed generalised LDA exchange and correlation functionals.
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}
%\manu{It goes much too fast here. One should make a clear distinction
@ -218,9 +223,9 @@ Electrons restricted to remain on the surface of a $\cD$-sphere (where $\cD$ is
%to ringium. In the present work we extend the approach to glomium. As we
%did in our previous work we should motivate the use of FUEGs for
%developing weight-dependent functionals.}
\titou{Very recently, \cite{Loos_2020} two of the present authors have taken advantages of these FUEGs to construct a local, weight-dependent correlation functional specifically designed for one-dimensional many-electron systems.
Very recently, \cite{Loos_2020} two of the present authors have taken advantages of these FUEGs to construct a local, weight-dependent correlation functional specifically designed for one-dimensional many-electron systems.
Unlike any standard functional, this first-rung functional incorporates derivative discontinuities thanks to its natural weight dependence, and has shown to deliver accurate excitation energies for both single and double excitations.
Extending this methodology to more realistic (atomic and molecular) systems, we combine here 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}}
Extending this methodology to more realistic (atomic and molecular) systems, we combine here 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.
@ -459,91 +464,92 @@ is the Hxc potential, with
%%%%% end stuff removed by Manu %%%%%%
\fi%%%%
\section{Some thougths illustrated with the Hubbard dimer model}
The definition of an ensemble density functional relies on the concavity
of the ensemble energy with respect to the external potential. In the
case of the Hubbard dimer, the singlet triensemble non-interacting
energy (which contains both singly- and doubly-excited states) reads
\beq
\begin{split}
\mathcal{E}_{\rm KS}^{\bw}\left(\Delta
v\right)=&(1-\ew{1}-\ew{2})\mathcal{E}_0\left(\Delta
v\right)+\ew{1}\mathcal{E}_1\left(\Delta
v\right)
\\
&+\ew{2}\mathcal{E}_2\left(\Delta
v\right),
\end{split}
\eeq
where $\mathcal{E}_0\left(\Delta
v\right)=2\varepsilon_0\left(\Delta
v\right)$, $\mathcal{E}_1\left(\Delta
v\right)=0$, $\mathcal{E}_2\left(\Delta
v\right)=-2\varepsilon_0\left(\Delta
v\right)$, and
\beq
\varepsilon_0\left(\Delta
v\right)=-\sqrt{t^2+\dfrac{\Delta v^2}{4}},
\eeq
thus leading to
\beq
\mathcal{E}_{\rm KS}^{\bw}\left(\Delta
v\right)=-2\left(1-\ew{1}-2\ew{2}\right)\sqrt{t^2+\dfrac{\Delta
v^2}{4}}.
\eeq
If we ignore the single excitation ($\ew{1}=0$) and denote
$\ew{}=\ew{2}$, the ensemble energy becomes
\beq
\mathcal{E}_{\rm KS}^{\ew{}}\left(\Delta
v\right)=-2(1-2\ew{})\sqrt{t^2+\dfrac{\Delta
v^2}{4}}.
\eeq
As readily seen, it is concave only if $\ew{}\leq 1/2$. Outside the
usual range of weight values, it is convex, thus preventing any density
to be ensemble non-interacting $v$-representable. This statement is
based on the Legendre--Fenchel transform expression of the
non-interacting ensemble kinetic energy functional:
\beq
T^{\ew{}}_{\rm s}(n)=\sup_{\Delta
v}\left\{\mathcal{E}_{\rm KS}^{\ew{}}\left(\Delta
v\right)+\Delta
v\times(n-1)\right\}.
\eeq
In this simple example, ignoring the single excitation is fine. However,
considering $1/2\leq \ew{}\leq 1$ is meaningless. Of course, if we
employ approximate ground-state-based density-functional potentials and
manage to converge the KS wavefunctions, one may obtain something
interesting. But I have no idea how meaningful such a solution is.\\
In the interacting case, the bi-ensemble (with the double excitation
only) energy reads
\beq
%\section{Some thougths illustrated with the Hubbard dimer model}
%
%The definition of an ensemble density functional relies on the concavity
%of the ensemble energy with respect to the external potential. In the
%case of the Hubbard dimer, the singlet triensemble non-interacting
%energy (which contains both singly- and doubly-excited states) reads
%\beq
%\begin{split}
E^{\ew{}}\left(\Delta
v\right)&=&(1-\ew{})E_0\left(\Delta
v\right)+\ew{}E_2\left(\Delta
v\right)
\nonumber
\\
&=&(1-\ew{})E_0\left(\Delta
v\right)+\ew{}\Big(2U-E_0\left(\Delta
v\right)-E_1\left(\Delta
v\right)\Big)
\nonumber
\\
&=&(1-2\ew{})E_0\left(\Delta
v\right)-\ew{}E_1\left(\Delta
v\right)+2U\ew{}.
%\mathcal{E}_{\rm KS}^{\bw}\left(\Delta
%v\right)=&(1-\ew{1}-\ew{2})\mathcal{E}_0\left(\Delta
%v\right)+\ew{1}\mathcal{E}_1\left(\Delta
%v\right)
%\\
%&+\ew{2}\mathcal{E}_2\left(\Delta
%v\right),
%\end{split}
\eeq
In the vicinity of the symmetric regime ($\Delta
v=0$), the excited-state energy is $E_1\left(\Delta
v\right)\approx U$. In this case, the ensemble energy is concave if
$\ew{}\leq 1/2$. One should check if $(1-2\ew{})E_0\left(\Delta
v\right)-\ew{}E_1\left(\Delta
v\right)$ remains concave away from this regime (I see no reason why it
should be).
%\eeq
%where $\mathcal{E}_0\left(\Delta
%v\right)=2\varepsilon_0\left(\Delta
%v\right)$, $\mathcal{E}_1\left(\Delta
%v\right)=0$, $\mathcal{E}_2\left(\Delta
%v\right)=-2\varepsilon_0\left(\Delta
%v\right)$, and
%\beq
%\varepsilon_0\left(\Delta
%v\right)=-\sqrt{t^2+\dfrac{\Delta v^2}{4}},
%\eeq
%thus leading to
%\beq
%\mathcal{E}_{\rm KS}^{\bw}\left(\Delta
%v\right)=-2\left(1-\ew{1}-2\ew{2}\right)\sqrt{t^2+\dfrac{\Delta
%v^2}{4}}.
%\eeq
%If we ignore the single excitation ($\ew{1}=0$) and denote
%$\ew{}=\ew{2}$, the ensemble energy becomes
%\beq
%\mathcal{E}_{\rm KS}^{\ew{}}\left(\Delta
%v\right)=-2(1-2\ew{})\sqrt{t^2+\dfrac{\Delta
%v^2}{4}}.
%\eeq
%As readily seen, it is concave only if $\ew{}\leq 1/2$. Outside the
%usual range of weight values, it is convex, thus preventing any density
%to be ensemble non-interacting $v$-representable. This statement is
%based on the Legendre--Fenchel transform expression of the
%non-interacting ensemble kinetic energy functional:
%\beq
%T^{\ew{}}_{\rm s}(n)=\sup_{\Delta
%v}\left\{\mathcal{E}_{\rm KS}^{\ew{}}\left(\Delta
%v\right)+\Delta
%v\times(n-1)\right\}.
%\eeq
%In this simple example, ignoring the single excitation is fine. However,
%considering $1/2\leq \ew{}\leq 1$ is meaningless. Of course, if we
%employ approximate ground-state-based density-functional potentials and
%manage to converge the KS wavefunctions, one may obtain something
%interesting. But I have no idea how meaningful such a solution is.\\
%
%In the interacting case, the bi-ensemble (with the double excitation
%only) energy reads
%\beq
%%\begin{split}
%E^{\ew{}}\left(\Delta
%v\right)&=&(1-\ew{})E_0\left(\Delta
%v\right)+\ew{}E_2\left(\Delta
%v\right)
%\nonumber
%\\
%&=&(1-\ew{})E_0\left(\Delta
%v\right)+\ew{}\Big(2U-E_0\left(\Delta
%v\right)-E_1\left(\Delta
%v\right)\Big)
%\nonumber
%\\
%&=&(1-2\ew{})E_0\left(\Delta
%v\right)-\ew{}E_1\left(\Delta
%v\right)+2U\ew{}.
%%\end{split}
%\eeq
%In the vicinity of the symmetric regime ($\Delta
%v=0$), the excited-state energy is $E_1\left(\Delta
%v\right)\approx U$. In this case, the ensemble energy is concave if
%$\ew{}\leq 1/2$. One should check if $(1-2\ew{})E_0\left(\Delta
%v\right)-\ew{}E_1\left(\Delta
%v\right)$ remains concave away from this regime (I see no reason why it
%should be).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% COMPUTATIONAL DETAILS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -556,13 +562,37 @@ For all calculations, we use the aug-cc-pVXZ (X = D, T, Q, and 5) Dunning family
Numerical quadratures are performed with the \texttt{numgrid} library \cite{numgrid} using 194 angular points (Lebedev grid) and a radial precision of $10^{-6}$. \cite{Becke_1988b,Lindh_2001}
This study deals only with spin-unpolarised systems, \ie, $\n{\uparrow}{} = \n{\downarrow}{} = \n{}{}/2$ (where $\n{\uparrow}{}$ and $\n{\downarrow}{}$ are the spin-up and spin-down electron densities).
Moreover, we restrict our study to the case of a two-state ensemble (\ie, $\nEns = 2$) where both the ground state ($I=0$ with weight $1 - \ew{}$) and the first doubly-excited state ($I=1$ with weight $\ew{}$) are considered.
\titou{To ensure the GOK variational principle, one should then have $0 \le \ew{} \le 1/2$.
However, we will sometimes ``violate'' this variational constraint.
Indeed, the limit $\ew{} = 1$ is of particular interest as it corresponds to a genuine saddle point of the KS equations, and match perfectly the results obtained with the maximum overlap method (MOM) developed by Gilbert, Gill and coworkers. \cite{Gilbert_2008,Barca_2018a,Barca_2018b}
Although the range $1/2 < \ew{} \leq 1$ stands a little bit beyond the theory discussed previously, we look at these solutions for analysis purposes mainly.
These solutions of the density matrix operator functional in Eq.~\eqref{eq:min_KS_DM} correspond to stationary points rather than minimising ones.
Applying GOK-DFT in this range of weights would simply consists in switching the ground and excited states if true minimisations of the ensemble energy were performed.}
Moreover, we restrict our study to the case of a three-state ensemble (\ie, $\nEns = 3$) where the ground state ($I=0$ with weight $1 - \ew{1} - \ew{2}$), the first singly-excited state ($I=1$ with weight $\ew{1}$), as well as the first doubly-excited state ($I=2$ with weight $\ew{2}$) are considered.
To ensure the GOK variational principle, one should then have $0 \le \ew{2} \le 1/3$ and $\ew{2} \le \ew{1} \le (1 - \ew{2})/2$.
Taking a two-electron system as an example, the individual one-electron densities reads
\begin{subequations}
\begin{align}
\n{}{(0)} & = 2 \HOMO{2},
\\
\n{}{(1)} & = \HOMO{2} + \LUMO{2},
\\
\n{}{(2)} & = 2 \LUMO{2},
\end{align}
\end{subequations}
and they can be combined to produce the ensemble density
\begin{equation}
\label{eq:nw1w2}
\n{}{(\ew{1},\ew{2})}
= (1 - \ew{1} - \ew{2}) \n{}{(0)}
+ \ew{1} \n{}{(1)} + \ew{2} \n{}{(2)}.
\end{equation}
Equation \eqref{eq:nw1w2} can be conveniently recast as a single-weight quantity
\begin{equation}
\n{}{\eW} = (1 - \eW) \n{}{(0)} + \eW \n{}{(2)},
\end{equation}
with $\eW = \ew{1}/2 + \ew{2}$ and $0 \le \eW \le 1/2$.
Unless otherwise stated, we set the same weight to the two excited states (\ie, $\ew{1} = \ew{2}$), and we consider the zero-weight limit $\eW = 0$ (\ie, $\ew{1} = \ew{2} = 0$), and the equiweight ensemble $\eW = 1/2$ (\ie, $\ew{1} = \ew{2} = 1/3$).
Nonetheless, we will sometimes ``violate'' the GOK variational principle in order to build our weight-dependent functionals.
Indeed, the limit $\ew{2} = 1$ (which corresponds to a pure excited state) is of particular interest as it is a genuine saddle point of the KS equations, and match perfectly the results obtained with the maximum overlap method (MOM) developed by Gilbert, Gill and coworkers. \cite{Gilbert_2008,Barca_2018a,Barca_2018b}
%Although the range $1/2 < \ew{} \leq 1$ stands a little bit beyond the theory discussed previously, we look at these solutions for analysis purposes mainly.
%These solutions of the density matrix operator functional in Eq.~\eqref{eq:min_KS_DM} correspond to stationary points rather than minimising ones.
%Applying GOK-DFT in this range of weights would simply consists in switching the ground and excited states if true minimisations of the ensemble energy were performed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -579,38 +609,34 @@ Applying GOK-DFT in this range of weights would simply consists in switching the
\subsubsection{Weight-independent exchange functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
First, we compute the ensemble energy of the \ce{H2} molecule (at equilibrium bond length, \ie, $\RHH = 1.4$ bohr) using the aug-cc-pVTZ basis set and the weight-independent \titou{LDA Slater exchange functional (\ie, no correlation functional is employed)}, \cite{Dirac_1930, Slater_1951} which is explicitly given by
First, we compute the ensemble energy of the \ce{H2} molecule (at equilibrium bond length, \ie, $\RHH = 1.4$ bohr) using the aug-cc-pVTZ basis set and the weight-independent LDA Slater exchange functional (\ie, no correlation functional is employed), \cite{Dirac_1930, Slater_1951} which is explicitly given by
\begin{align}
\e{\ex}{\text{S}}(\n{}{}) & = \Cx{} \n{}{1/3},
&
\Cx{} & = -\frac{3}{4} \qty(\frac{3}{\pi})^{1/3}.
\end{align}
In the case of \ce{H2}, the ensemble is composed by the $1\sigma_g^2$ ground state and the lowest doubly-excited state of configuration $1\sigma_u^2$, which has an autoionising resonance nature. \cite{Bottcher_1974}
\manu{At equilibrium, I expect the singly-excited configuration
$1\sigma_g2\sigma_g$ to be lower in energy. From the point of view of
GOK-DFT I do not see how we can reach the doubly-excited state while
ignoring the singly-excited one. One can always argue that we explore
stationary points (and not minima) but an obvious and important question that any
referee working on GOK-DFT would ask is: How would your results
be changed if you were incorporating the single excitation in your
ensemble? In one way or another
we have to look at this, even within the simplest weight-independent
approximation.}
The ensemble energy $\E{}{\ew{}}$ is depicted in Fig.~\ref{fig:Ew_H2} as a function of the weight $0 \le \ew{} \le 1$.
In the case of \ce{H2}, the ensemble is composed by the $1\sigma_g^2$ ground state, the lowest singly-excited state $1\sigma_g 1\sigma_u$ and the lowest doubly-excited state of configuration $1\sigma_u^2$ (which has an autoionising resonance nature \cite{Bottcher_1974}).
%\manu{At equilibrium, I expect the singly-excited configuration
%$1\sigma_g2\sigma_g$ to be lower in energy. From the point of view of
%GOK-DFT I do not see how we can reach the doubly-excited state while
%ignoring the singly-excited one. One can always argue that we explore
%stationary points (and not minima) but an obvious and important question that any
%referee working on GOK-DFT would ask is: How would your results
%be changed if you were incorporating the single excitation in your
%ensemble? In one way or another
%we have to look at this, even within the simplest weight-independent
%approximation.}
The ensemble energy $\E{}{\eW{}}$ is depicted in Fig.~\ref{fig:Ew_H2} as a function of the composite weight $0 \le \eW{} \le 1/2$.
Because this exchange functional does not depend on the ensemble weight, there is no contribution from the ensemble derivative term [last term in Eq.~\eqref{eq:dEdw}].
As anticipated, $\E{}{\ew{}}$ is far from being linear, which means that the excitation energy obtained via the derivative of the ensemble energy varies significantly with the weight of the double excitation (see Fig.~\ref{fig:Om_H2}).
Taking as a reference the full configuration interaction (FCI) value of $28.75$ eV obtained with the aug-mcc-pV8Z basis set, \cite{Barca_2018a} one can see that the excitation energy varies by more than $8$ eV from $\ew{} = 0$ to $1/2$.
As anticipated, $\E{}{\eW{}}$ is far from being linear, which means that the excitation energy obtained via the derivative of the ensemble energy varies significantly with the weight of the double excitation (see Fig.~\ref{fig:Om_H2}).
Taking as a reference the full configuration interaction (FCI) value of $28.75$ eV obtained with the aug-mcc-pV8Z basis set, \cite{Barca_2018a} one can see that the excitation energy varies by more than $8$ eV from $\eW{} = 0$ to $1/2$.
Note that the exact xc ensemble functional would yield a perfectly linear energy and, hence, the same value of the excitation energy independently of $\ew{}$.
\titou{The eDFT purist is going to surprised to see that we left out the singly-excited states $\sigma_g \sigma_u$ from the ensemble as this state is lower in energy than the doubly-excited state of configuration $1\sigma_u^2$.
As we wish to stick with a restricted formalism, the single excitation is naturally left out of the ensemble.
However, as a sanity check, we have tried to introduce the single excitations as well.}
\begin{figure}
\includegraphics[width=\linewidth]{Ew_H2}
\caption{
\ce{H2} at equilibrium bond length: deviation from linearity of the ensemble energy $\E{}{\ew{}}$ (in hartree) as a function of the weight of the double excitation $\ew{}$ for various functionals and the aug-cc-pVTZ basis set.
\titou{See main text for the definition of the various functionals.}
\ce{H2} at equilibrium bond length: deviation from linearity of the ensemble energy $\E{}{\eW}$ (in hartree) as a function of the weight of the double excitation $\eW = \ew{1}/2 + \ew{2}$ and $\ew{1} = \ew{2}$ for various functionals and the aug-cc-pVTZ basis set.
See main text for the definition of the various functional's acronyms.
\label{fig:Ew_H2}
}
\end{figure}
@ -618,7 +644,7 @@ However, as a sanity check, we have tried to introduce the single excitations as
\begin{figure}
\includegraphics[width=\linewidth]{Om_H2}
\caption{
\ce{H2} at equilibrium bond length: error (with respect to FCI) in excitation energy (in eV) as a function of the weight of the double excitation $\ew{}$ for various functionals and the aug-cc-pVTZ basis set.
\ce{H2} at equilibrium bond length: error (with respect to FCI) in excitation energy (in eV) \titou{for the double excitation} as a function of the composite weight $\eW$ for various functionals and the aug-cc-pVTZ basis set.
\label{fig:Om_H2}
}
\end{figure}
@ -630,41 +656,43 @@ However, as a sanity check, we have tried to introduce the single excitations as
Second, in order to remove this spurious curvature of the ensemble
energy (which is mostly due to the ghost-interaction error, \cite{Loos_2020} but not only),
one can easily reverse-engineer (for this particular system, geometry, and basis set) a local exchange functional to make $\E{}{\ew{}}$ as linear as possible for $0 \le \ew{} \le 1$.
\manu{Something that seems important to me: you may require linearity in
the range $0\leq \ew{}\leq 1/2$. The excitation energy you would obtain
is simply the one of LIM, right? I suspect that by considering the
endpoint $\ew{}=1$ you change the excitation energy substantially. How
different are the results? At first sight, it seems like MOM gives you
the excitation energy that drives the
parameterization of the functional. Regarding the excitation energy, the
parameterized functional does not bring any additional information,
right? Maybe I miss something. Of course it gives ideas about how to
construct functionals. Maybe we need to elaborate more on this. For
example, its combination with correlation functionals (as done in the
following) is very interesting. It should be introduced as a kind of
two-step procedure.}
Doing so, we have found that the present weight-dependent exchange functional \titou{(denoted as CC-S for ``curvature-corrected'' Slater functional)}
one can easily reverse-engineer (for this particular system, geometry, and basis set) a local exchange functional to make $\E{}{\eW}$ as linear as possible for $0 \le \eW \le 1/2$.
%\manu{Something that seems important to me: you may require linearity in
%the range $0\leq \ew{}\leq 1/2$. The excitation energy you would obtain
%is simply the one of LIM, right? I suspect that by considering the
%endpoint $\ew{}=1$ you change the excitation energy substantially. How
%different are the results? At first sight, it seems like MOM gives you
%the excitation energy that drives the
%parameterization of the functional. Regarding the excitation energy, the
%parameterized functional does not bring any additional information,
%right? Maybe I miss something. Of course it gives ideas about how to
%construct functionals. Maybe we need to elaborate more on this. For
%example, its combination with correlation functionals (as done in the
%following) is very interesting. It should be introduced as a kind of
%two-step procedure.}
Doing so, we have found that the present weight-dependent exchange functional (denoted as CC-S for ``curvature-corrected'' Slater functional)
\begin{equation}
\e{\ex}{\ew{},\text{CC-S}}(\n{}{}) = \Cx{\ew{}} \n{}{1/3},
\e{\ex}{\eW,\text{CC-S}}(\n{}{}) = \Cx{\eW} \n{}{1/3},
\end{equation}
with
\begin{equation}
\label{eq:Cxw}
\frac{\Cx{\ew{}}}{\Cx{}} = 1 - \ew{} (1 - \ew{})\qty[ \alpha + \beta (\ew{} - 1/2) + \gamma (w - 1/2)^2 ],
\frac{\Cx{\eW}}{\Cx{}} = 1 - \eW (1 - \eW)\qty[ \alpha + \beta (\eW - 1/2) + \gamma (\eW - 1/2)^2 ],
\end{equation}
and
\begin{subequations}
\begin{align}
\alpha & = + 0.575\,178,
\alpha & = + 0.573\,919,
&
\beta & = - 0.021\,108,
\beta & = - 0.000\,347,
&
\gamma & = - 0.367\,189,
\gamma & = - 0.267\,634,
\end{align}
\end{subequations}
makes the ensemble almost perfectly linear (see Fig.~\ref{fig:Ew_H2}), and the excitation energy much more stable (with respect to $\ew{}$) and closer to the FCI reference (see Fig.~\ref{fig:Om_H2}).
The parameters $\alpha$, $\beta$, and $\gamma$ entering Eq.~\eqref{eq:Cxw} have been obtained via a least-square fit of the non-linear component of the ensemble energy computed between $\ew{} = 0$ and $\ew{} = 1$ by steps of $0.025$.
makes the ensemble almost perfectly linear (see Fig.~\ref{fig:Ew_H2}), and the excitation energy much more stable (with respect to $\eW$) and closer to the FCI reference (see Fig.~\ref{fig:Om_H2}).
The parameters $\alpha$, $\beta$, and $\gamma$ entering Eq.~\eqref{eq:Cxw} have been obtained via a least-square fit of the non-linear component of the ensemble energy computed between $\eW = 0$ and $\eW = 1/2$ by steps of $0.025$.
Note that this procedure makes, by construction, the ensemble energy also linear with respect to $\ew{1}$ and $\ew{2}$.
As readily seen from Eq.~\eqref{eq:Cxw} and graphically illustrated in Fig.~\ref{fig:Cxw} (red curve), the weight-dependent correction does not affect the two ghost-interaction-free limits at $\ew{} = 0$ and $\ew{} = 1$, as $\Cx{\ew{}}$ reduces to $\Cx{}$ in these two limits.
Maybe surprisingly, one would have noticed that we are not only using data from $0 \le \ew{} \le 1/2$, but we also consider ensemble energies for $1/2 < \ew{} \le 1$, which is deterred by the GOK variational principle. \cite{Gross_1988a}
However, it is important to ensure that the weight-dependent functional does not alter the $\ew{} = 1$ limit, which is a genuine saddle point of the KS equations, as mentioned above.
@ -706,14 +734,17 @@ As mentioned above, we confine our attention to paramagnetic (or unpolarised) sy
Note that the present paradigm is equivalent to the conventional IUEG model in the thermodynamic limit. \cite{Loos_2011b}
We refer the interested reader to Refs.~\onlinecite{Loos_2011b,Loos_2017a} for more details about this paradigm.
The reduced (\ie, per electron) Hartree-Fock (HF) energies for these two states are
The reduced (\ie, per electron) Hartree-Fock (HF) energies for these three states are
\begin{subequations}
\begin{align}
\e{\HF}{(0)}(\n{}{}) & = \frac{4}{3} \qty(\frac{\n{}{}}{\pi})^{1/3},
\label{eq:eHF_0}
\\
\e{\HF}{(1)}(\n{}{}) & = \frac{3\pi^{2}}{2} \qty(\frac{\n{}{}}{\pi})^{2/3} + \frac{176}{105} \qty(\frac{\n{}{}}{\pi})^{1/3}.
\e{\HF}{(1)}(\n{}{}) & = \frac{3\pi^{2}}{4} \qty(\frac{\n{}{}}{\pi})^{2/3} + \frac{16}{10} \qty(\frac{\n{}{}}{\pi})^{1/3}.
\label{eq:eHF_1}
\\
\e{\HF}{(2)}(\n{}{}) & = \frac{3\pi^{2}}{2} \qty(\frac{\n{}{}}{\pi})^{2/3} + \frac{176}{105} \qty(\frac{\n{}{}}{\pi})^{1/3}.
\label{eq:eHF_2}
\end{align}
\end{subequations}
@ -726,17 +757,17 @@ where $a_2^{(I)}$ and $a_3^{(I)}$ are state-specific fitting parameters, which a
The value of $a_1^{(I)}$ is obtained via the exact high-density expansion of the correlation energy. \cite{Loos_2011b}
Equation \eqref{eq:ec} is depicted in Fig.~\ref{fig:Ec} for each state alongside the data gathered in Table \ref{tab:Ref}.
Combining these, we build a two-state weight-dependent correlation functional:
Combining these, we build a three-state weight-dependent correlation functional:
\begin{equation}
\label{eq:ecw}
\e{\co}{\ew{}}(\n{}{}) = (1-\ew{}) \e{\co}{(0)}(\n{}{}) + \ew{} \e{\co}{(1)}(\n{}{}).
\e{\co}{\bw}(\n{}{}) = (1-\ew{1}-\ew{2}) \e{\co}{(0)}(\n{}{}) + \ew{1} \e{\co}{(1)}(\n{}{}) + \ew{2} \e{\co}{(2)}(\n{}{})
\end{equation}
%%% FIG 1 %%%
\begin{figure}
\includegraphics[width=0.8\linewidth]{fig1}
\caption{
Reduced (i.e., per electron) correlation energy $\e{\co}{(I)}$ [see Eq.~\eqref{eq:ec}] as a function of $R = 1/(\pi^2 \n{}{})^{1/3}$ for the ground state ($I=0$), and the first doubly-excited state ($I=1$) of the (spin-unpolarised) two-electron FUEG.
Reduced (i.e., per electron) correlation energy $\e{\co}{(I)}$ [see Eq.~\eqref{eq:ec}] as a function of $R = 1/(\pi^2 \n{}{})^{1/3}$ for the ground state ($I=0$), the first singly-excited state ($I=1$), and the first doubly-excited state ($I=2$) of the two-electron FUEG.
The data gathered in Table \ref{tab:Ref} are also reported.
}
\label{fig:Ec}
@ -747,29 +778,30 @@ Combining these, we build a two-state weight-dependent correlation functional:
\begin{table}
\caption{
\label{tab:Ref}
$-\e{\co}{(I)}$ as a function of the radius of the glome $R = 1/(\pi^2 \n{}{})^{1/3}$ for the ground state ($I=0$), and the first doubly-excited state ($I=1$) of the (spin-unpolarised) two-electron FUEG.
$-\e{\co}{(I)}$ as a function of the radius of the glome $R = 1/(\pi^2 \n{}{})^{1/3}$ for the ground state ($I=0$), the first singly-excited state ($I=1$), and the first doubly-excited state ($I=2$) of the two-electron FUEG.
}
\begin{ruledtabular}
\begin{tabular}{lcc}
& \tabc{Ground state} & \tabc{Doubly-excited state} \\
$R$ & \tabc{$I=0$} & \tabc{$I=1$} \\
\begin{tabular}{lccc}
& \tabc{Ground state} & \tabc{Single excitation} & \tabc{Double excitation} \\
$R$ & \tabc{$I=0$} & \tabc{$I=1$} & \tabc{$I=2$} \\
\hline
$0$ & $0.023\,818$ & $0.014\,463$ \\
$0.1$ & $0.023\,392$ & $0.014\,497$ \\
$0.2$ & $0.022\,979$ & $0.014\,523$ \\
$0.5$ & $0.021\,817$ & $0.014\,561$ \\
$1$ & $0.020\,109$ & $0.014\,512$ \\
$2$ & $0.017\,371$ & $0.014\,142$ \\
$5$ & $0.012\,359$ & $0.012\,334$ \\
$10$ & $0.008\,436$ & $0.009\,716$ \\
$20$ & $0.005\,257$ & $0.006\,744$ \\
$50$ & $0.002\,546$ & $0.003\,584$ \\
$100$ & $0.001\,399$ & $0.002\,059$ \\
$150$ & $0.000\,972$ & $0.001\,458$ \\
$0$ & $0.023\,818$ & 0.028\,281 & $0.014\,463$ \\
$0.1$ & $0.023\,392$ & 0.027\,886 & $0.014\,497$ \\
$0.2$ & $0.022\,979$ & 0.027\,499 & $0.014\,523$ \\
$0.5$ & $0.021\,817$ & 0.026\,394 & $0.014\,561$ \\
$1$ & $0.020\,109$ & 0.024\,718 & $0.014\,512$ \\
$2$ & $0.017\,371$ & 0.021\,901 & $0.014\,142$ \\
$5$ & $0.012\,359$ & 0.016\,295 & $0.012\,334$ \\
$10$ & $0.008\,436$ & 0.011\,494 & $0.009\,716$ \\
$20$ & $0.005\,257$ & 0.007\,349 & $0.006\,744$ \\
$50$ & $0.002\,546$ & 0.003\,643 & $0.003\,584$ \\
$100$ & $0.001\,399$ & 0.002\,025 & $0.002\,059$ \\
$150$ & $0.000\,972$ & 0.001\,414 & $0.001\,458$ \\
\end{tabular}
\end{ruledtabular}
\end{table}
%%% TABLE II %%%
\begin{table}
\caption{
@ -777,26 +809,26 @@ Combining these, we build a two-state weight-dependent correlation functional:
Parameters of the correlation functionals for each individual state defined in Eq.~\eqref{eq:ec}.
The values of $a_1$ are obtained to reproduce the exact high density correlation energy of each individual state, while $a_2$ and $a_3$ are fitted on the numerical values reported in Table \ref{tab:Ref}.}
\begin{ruledtabular}
\begin{tabular}{lll}
& \tabc{Ground state} & \tabc{Doubly-excited state} \\
& \tabc{$I=0$} & \tabc{$I=1$} \\
\begin{tabular}{llll}
& \tabc{Ground state} & \tabc{Single excitation} & \tabc{Double excitation} \\
& \tabc{$I=0$} & \tabc{$I=1$} & \tabc{$I=2$} \\
\hline
$a_1$ & $-0.023\,818\,4$ & $-0.014\,463\,3$ \\
$a_2$ & $+0.005\,409\,94$ & $-0.050\,601\,9$ \\
$a_3$ & $+0.083\,076\,6$ & $+0.033\,141\,7$ \\
$a_1$ & $-0.023\,818\,4$ & $-0.028\,281\,4$ & $-0.014\,463\,3$ \\
$a_2$ & $+0.005\,409\,94$ & $+0.002\,739\,25$ & $-0.050\,602\,0$ \\
$a_3$ & $+0.083\,076\,6$ & $+0.066\,491\,4$ & $+0.033\,141\,7$ \\
\end{tabular}
\end{ruledtabular}
\end{table}
%%% %%% %%% %%%
Because our intent is to incorporate into standard functionals (which are ``universal'' in the sense that they do not depend on the number of electrons) information about excited states that will be extracted from finite systems (whose properties may depend on the number of electrons), we employ a simple embedding scheme where the two-electron FUEG (the impurity) is embedded in the IUEG (the bath).
Because our intent is to incorporate into standard functionals (which are ``universal'' in the sense that they do not depend on the number of electrons) information about excited states that will be extracted from finite systems (whose properties may depend on the number of electrons), we employ a simple ``embedding'' scheme where the two-electron FUEG (the impurity) is embedded in the IUEG (the bath).
As explained further in Ref.~\onlinecite{Loos_2020}, this embedding procedure can be theoretically justified by the generalised adiabatic connection formalism for ensembles originally derived by Franck and Fromager. \cite{Franck_2014}
The weight-dependence of the correlation functional is then carried exclusively by the impurity [\ie, the functional defined in \eqref{eq:ecw}], while the remaining effects are produced by the bath (\ie, the usual LDA correlation functional).
Consistently with such a strategy, Eq.~\eqref{eq:ecw} is ``centred'' on its corresponding weight-independent VWN5 LDA reference
\begin{equation}
\label{eq:becw}
\e{\co}{\ew{},\eVWN}(\n{}{}) = (1-\ew{}) \be{\co}{(0)}(\n{}{}) + \ew{} \be{\co}{(1)}(\n{}{})
\e{\co}{\bw{},\eVWN}(\n{}{}) = (1-\ew{1}-\ew{2}) \be{\co}{(0)}(\n{}{}) + \ew{1} \be{\co}{(1)}(\n{}{}) + \ew{2} \be{\co}{(2)}(\n{}{})
\end{equation}
via the following global, state-independent shift:
\begin{equation}
@ -806,15 +838,21 @@ In the following, we name this weight-dependent correlation functional ``eVWN5''
Also, Eq.~\eqref{eq:becw} can be recast as
\begin{equation}
\label{eq:eLDA}
\e{\co}{\ew{},\eVWN}(\n{}{}) = \e{\co}{\VWN}(\n{}{}) + \ew{} \qty[\e{\co}{(1)}(\n{}{}) - \e{\co}{(0)}(\n{}{})]
\begin{split}
\e{\co}{\bw{},\eVWN}(\n{}{})
& = \e{\co}{\VWN}(\n{}{})
\\
& + \ew{1} \qty[\e{\co}{(1)}(\n{}{}) - \e{\co}{(0)}(\n{}{})]
+ \ew{2} \qty[\e{\co}{(2)}(\n{}{}) - \e{\co}{(0)}(\n{}{})]
\end{split}
\end{equation}
which nicely highlights the centrality of the LDA in the present weight-dependent density-functional approximation for ensembles.
In particular, $\e{\co}{\ew{}=0,\eVWN}(\n{}{}) = \e{\co}{\VWN}(\n{}{})$.
In particular, $\e{\co}{(0,0),\eVWN}(\n{}{}) = \e{\co}{\VWN}(\n{}{})$.
We note also that, by construction, we have
\begin{equation}
\label{eq:dexcdw}
\pdv{\e{\co}{\ew{},\eVWN}(\n{}{})}{\ew{}}
= \e{\co}{(1)}(n) - \e{\co}{(0)}(n),
\pdv{\e{\co}{\bw{},\eVWN}(\n{}{})}{\ew{I}}
= \e{\co}{(I)}(n) - \e{\co}{(0)}(n),
\end{equation}
showing that the weight correction is purely linear in eVWN5 and entirely dependent on the FUEG model.
@ -828,15 +866,24 @@ in the zero-weight limit (\ie, $\ew{} = 0$) and for the equi-ensemble
For comparison purposes, we also report the linear interpolation method (LIM) excitation energies, \cite{Senjean_2015,Senjean_2016}
a pragmatic way of getting weight-independent
excitation energies defined as
\begin{equation}
\Ex{\LIM}{} = 2 (\E{}{\ew{}=1/2} - \E{}{\ew{}=0}),
\end{equation}
as well as the MOM excitation energies. \cite{Gilbert_2008,Barca_2018a,Barca_2018b}
\begin{subequations}
\begin{align}
\Ex{\LIM}{(1)} & = 2 \qty[\E{}{\bw{}=(1/2,0)} - \E{}{\bw{}=(0,0)}],
\\
\Ex{\LIM}{(2)} & = 3 \qty[\E{}{\bw{}=(1/3,1/3)} - \E{}{\bw{}=(1/2,0)}] + \frac{1}{2} \Ex{\LIM}{(1)},
\end{align}
\end{subequations}
which require three independent calculations, as well as the MOM excitation energies \cite{Gilbert_2008,Barca_2018a,Barca_2018b}
\begin{subequations}
\begin{align}
\Ex{\MOM}{(1)} & = \E{}{\bw{}=(1,0)} - \E{}{\bw{}=(0,0)},
\\
\Ex{\MOM}{(2)} & = \E{}{\bw{}=(0,1)} - \E{}{\bw{}=(0,0)}.
\end{align}
\end{subequations}
which also require three separate calculations.
We point out that the MOM excitation energy is obtained by the difference in energy between the doubly-excited state at $\ew{} = 1$ and the ground state at $\ew{} = 0$.
They can then be obtained via GOK-DFT ensemble calculations by performing a linear interpolation between $\ew{} = 0$ and $\ew{} = 1$, \ie,
\begin{equation}
\Ex{\MOM}{} = \E{}{\ew{}=1} - \E{}{\ew{}=0}.
\end{equation}
The results gathered in Table \ref{tab:BigTab_H2} show that the GOK-DFT excitation energies obtained with the CC-SeVWN5 functional at zero weight are the most accurate with an improvement of $0.25$ eV as compared to CC-SVWN5, which is due to the ensemble derivative contribution of the eVWN5 functional.
The CC-SeVWN5 excitation energies at equi-weights (\ie, $\ew{} = 1/2$) are less satisfactory, but still remain in good agreement with FCI, with again a small improvement as compared to CC-SVWN5.
@ -859,17 +906,17 @@ Excitation energies (in eV) associated with the lowest double excitation of \ce{
\begin{tabular}{llccccc}
\mc{2}{c}{xc functional} & & \mc{2}{c}{GOK} \\
\cline{1-2} \cline{4-5}
\tabc{x} & \tabc{c} & Basis & $\ew{} = 0$ & $\ew{} = 1/2$ & LIM & MOM \\
\tabc{x} & \tabc{c} & Basis & $\bw{} = \bO$ & $\bw{} = \bthird$ & LIM & MOM \\
\hline
HF & & aug-cc-pVDZ & 38.52 & 30.86 & 34.55 & 28.65 \\
& & aug-cc-pVTZ & 38.58 & 35.82 & 35.68 & 28.65 \\
& & aug-cc-pVQZ & 39.12 & 35.94 & 35.64 & 28.65 \\
HF & & aug-cc-pVDZ & 38.52 & 30.86 & 34.00 & 28.65 \\
& & aug-cc-pVTZ & 38.58 & 35.82 & 35.80 & 28.65 \\
& & aug-cc-pVQZ & 39.12 & 35.94 & 35.82 & 28.65 \\
\\
HF & VWN5 & aug-cc-pVDZ & 37.83 & 31.19 & 35.66 & 29.17 \\
HF & VWN5 & aug-cc-pVDZ & 37.83 & 31.19 & 34.91 & 29.17 \\
& & aug-cc-pVTZ & 37.35 & 36.98 & 37.23 & 29.17 \\
& & aug-cc-pVQZ & 37.07 & 37.07 & 37.21 & 29.17 \\
\\
HF & eVWN5 & aug-cc-pVDZ & 38.09 & 31.34 & 35.74 & 29.34 \\
HF & eVWN5 & aug-cc-pVDZ & 38.09 & 31.35 & 35.00 & 29.34 \\
& & aug-cc-pVTZ & 37.61 & 37.04 & 37.28 & 29.34 \\
& & aug-cc-pVQZ & 37.32 & 37.14 & 37.27 & 29.34 \\
\\

Binary file not shown.

Binary file not shown.