loos/eDFT_FUEG
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
d4e16e4b0e
eDFT_FUEG/Manuscript/eDFT.tex
Emmanuel Fromager d4e16e4b0e Manu: saving work
2020-02-12 12:04:03 +01:00

1333 lines
59 KiB
TeX
Raw Blame History

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,mhchem,longtable}
\usepackage{mathpazo,libertine}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{RGB}{0, 180, 0}
\usepackage{hyperref}
\hypersetup{
colorlinks=true,
linkcolor=blue,
filecolor=blue,
urlcolor=blue,
citecolor=blue
}
%useful stuff
\newcommand{\cdash}{\multicolumn{1}{c}{---}}
\newcommand{\mc}{\multicolumn}
\newcommand{\mr}{\multirow}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\la}{\lambda}
\newcommand{\si}{\sigma}
% operators
\newcommand{\hHc}{\Hat{h}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hVext}{\Hat{V}_\text{ext}}
\newcommand{\hWee}{\Hat{W}_\text{ee}}
% functionals, potentials, densities, etc
\newcommand{\eps}{\epsilon}
\newcommand{\e}[2]{\eps_\text{#1}^{#2}}
\newcommand{\E}[2]{E_\text{#1}^{#2}}
\newcommand{\bE}[2]{\overline{E}_\text{#1}^{#2}}
\newcommand{\be}[2]{\overline{\eps}_\text{#1}^{#2}}
\newcommand{\bv}[2]{\overline{f}_\text{#1}^{#2}}
\newcommand{\n}[2]{n_{#1}^{#2}}
\newcommand{\DD}[2]{\Delta_\text{#1}^{#2}}
\newcommand{\LZ}[2]{\Xi_\text{#1}^{#2}}
% energies
\newcommand{\EHF}{E_\text{HF}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\Ecat}{E_\text{cat}}
\newcommand{\Eneu}{E_\text{neu}}
\newcommand{\Eani}{E_\text{ani}}
\newcommand{\EPT}{E_\text{PT2}}
\newcommand{\EFCI}{E_\text{FCI}}
% matrices/operator
\newcommand{\br}{\bm{r}}
\newcommand{\bw}{{\bm{w}}}
\newcommand{\bG}{\bm{G}}
\newcommand{\bS}{\bm{S}}
\newcommand{\bGamma}[1]{\bm{\Gamma}^{#1}}
\newcommand{\opGamma}[1]{\hat{\Gamma}^{#1}}
\newcommand{\bHc}{\bm{h}}
\newcommand{\bF}[1]{\bm{F}^{#1}}
\newcommand{\Ex}[1]{\Omega^{#1}}
% elements
\newcommand{\ew}[1]{w_{#1}}
\newcommand{\eG}[1]{G_{#1}}
\newcommand{\eS}[1]{S_{#1}}
\newcommand{\eGamma}[2]{\Gamma_{#1}^{#2}}
\newcommand{\hGamma}[2]{\Hat{\Gamma}_{#1}^{#2}}
\newcommand{\eHc}[1]{h_{#1}}
\newcommand{\eF}[2]{F_{#1}^{#2}}
% Numbers
\newcommand{\Nel}{N}
\newcommand{\Nbas}{K}
% Ao and MO basis
\newcommand{\MO}[2]{\phi_{#1}^{#2}}
\newcommand{\cMO}[2]{c_{#1}^{#2}}
\newcommand{\AO}[1]{\chi_{#1}}
% units
\newcommand{\IneV}[1]{#1~eV}
\newcommand{\InAU}[1]{#1~a.u.}
\newcommand{\InAA}[1]{#1~\AA}
\newcommand{\SI}{\textcolor{blue}{supplementary material}}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LCQ}{Laboratoire de Chimie Quantique, Institut de Chimie, CNRS, Universit\'e de Strasbourg, Strasbourg, France}
%%% added by Manu %%%
\newcommand{\manu}[1]{{\textcolor{blue}{ Manu: #1 }} }
\newcommand{\beq}{\begin{eqnarray}}
\newcommand{\eeq}{\end{eqnarray}}
\newcommand{\bmk}{\bm{\kappa}} % orbital rotation vector
\newcommand{\bmg}{\bm{\Gamma}} % orbital rotation vector
\newcommand{\bxi}{\bm{\xi}}
\newcommand{\bfx}{{\bf{x}}}
\newcommand{\bfr}{{\bf{r}}}
\DeclareMathOperator*{\argmin}{arg\,min}
\newcommand{\blue}[1]{{\textcolor{blue}{#1}}}
%%%%
\begin{document}
\title{Weight-dependent local density-functional approximations for ensembles}
\author{Pierre-Fran\c{c}ois Loos}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Emmanuel Fromager}
\email{fromagere@unistra.fr}
\affiliation{\LCQ}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
We report a first generation of local, weight-dependent correlation density-functional approximations (DFAs) that incorporate information about both ground and excited states in the context of density-functional theory for ensembles (eDFT).
These density-functional approximations for ensembles (eDFAs) are specially designed for the computation of single and double excitations within eDFT, and can be seen as a natural extension of the ubiquitous local-density approximation for ensemble (eLDA).
The resulting eDFAs, based on both finite and infinite uniform electron gas models, automatically incorporate the infamous derivative discontinuity contributions to the excitation energies through their explicit ensemble weight dependence.
Their accuracy is illustrated by computing single and double excitations in one-dimensional many-electron systems in the weak, intermediate and strong correlation regimes.
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Over the last two decades, density-functional theory (DFT) \cite{Hohenberg_1964, Kohn_1965} has become the method of choice for modeling the electronic structure of large molecular systems and materials. \cite{ParrBook}
The main reason is that, within DFT, the quantum contributions to the electronic repulsion energy --- the so-called exchange-correlation (xc) energy --- is rewritten as a functional of the electron density $\n{}{}(\br)$, the latter being a much simpler quantity than the many-electron wave function.
The complexity of the many-body problem is then transferred to the xc functional.
Despite its success, the standard Kohn-Sham (KS) formulation of DFT \cite{Kohn_1965} (KS-DFT) suffers, in practice, from various deficiencies. \cite{Woodcock_2002, Tozer_2003, Tozer_1999, Dreuw_2003, Sobolewski_2003, Dreuw_2004, Tozer_1998, Tozer_2000, Casida_1998, Casida_2000, Tapavicza_2008, Levine_2006}
The description of strongly multiconfigurational ground states (often referred to as ``strong correlation problem'') still remains a challenge.
Another issue, which is partly connected to the previous one, is the description of electronically-excited states.
The standard approach for modeling excited states in DFT is linear response time-dependent DFT (TDDFT). \cite{Casida}
In this case, the electronic spectrum relies on the (unperturbed) ground-state KS picture, which may break down when electron correlation is strong.
Moreover, in exact TDDFT, the xc functional is time dependent.
The simplest and most widespread approximation in state-of-the-art electronic structure programs where TDDFT is implemented consists in neglecting memory effects. \cite{Casida}
In other words, within this so-called adiabatic approximation, the xc functional is assumed to be local in time.
As a result, double electronic excitations are completely absent from the TDDFT spectrum, thus reducing further the applicability of TDDFT. \cite{Maitra_2004,Cave_2004,Mazur_2009,Mazur_2011,Huix-Rotllant_2011,Elliott_2011,Maitra_2012,Sundstrom_2014}
When affordable (i.e., for relatively small molecules), time-independent state-averaged wave function methods \cite{Roos,Andersson_1990,Angeli_2001a,Angeli_2001b,Angeli_2002} can be employed to fix the various issues mentioned above.
The basic idea is to describe a finite ensemble of states (ground and excited) altogether, i.e., with the same set of orbitals.
Interestingly, a similar approach exists in DFT.
Ensemble DFT (eDFT) was proposed at the end of the 80's by Gross, Oliveira and Kohn (GOK), \cite{Gross_1988, Gross_1988a, Oliveira_1988} and is a generalization of Theophilou's variational principle for equi-ensembles. \cite{Theophilou_1979}
In eDFT, the (time-independent) xc functional depends explicitly on the weights assigned to the states that belong to the ensemble of interest.
This weight dependence of the xc functional plays a crucial role in the calculation of excitation energies.
It actually accounts for the infamous derivative discontinuity contribution to energy gaps. \cite{Levy_1995, Perdew_1983}
\alert{Shall we further discuss the derivative discontinuity? Why is it important and where is it coming from?}
Despite its formal beauty and the fact that eDFT can in principle tackle near-degenerate situations and multiple excitations, it has not been given much attention until recently. \cite{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_2018a,Deur_2018b,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,Senjean_2018,Smith_2016}
The main reason is simply the absence of density-functional approximations (DFAs) for ensembles in the literature.
Recent works on this topic are still fundamental and exploratory, as they rely either on simple (but nontrivial) models like the Hubbard dimer \cite{Carrascal_2015,Deur_2017,Deur_2018a,Deur_2018b,Senjean_2015,Senjean_2016,Senjean_2018,Sagredo_2018} or on atoms for which highly accurate or exact-exchange-only calculations have been performed. \cite{Yang_2014,Yang_2017}
In both cases, the key problem, namely the design of weight-dependent DFAs for ensembles (eDFAs), remains open.
A first step towards this goal is presented in this article with the ambition to turn, in the near future, eDFT into a practical computational method for modeling excited states in molecules and extended systems.
\alert{Mention WIDFA?}
In the following, the present methodology is illustrated on \emph{strict} one-dimensional (1D), spin-polarized electronic systems. \cite{Loos_2012, Loos_2013a, Loos_2014a, Loos_2014b}
In other words, the Coulomb interaction used in this work describes particles which are \emph{strictly} restricted to move within a 1D sub-space of three-dimensional space.
Despite their simplicity, 1D models are scrutinized as paradigms for quasi-1D materials \cite{Schulz_1993, Fogler_2005a} such as carbon nanotubes \cite{Bockrath_1999, Ishii_2003, Deshpande_2008} or nanowires. \cite{Meyer_2009, Deshpande_2010}
%Early models of 1D atoms using this interaction have been used to study the effects of external fields upon Rydberg atoms \cite{Burnett_1993, Mayle_2007} and the dynamics of surface-state electrons in liquid helium. \cite{Nieto_2000, Patil_2001}
This description of 1D systems also has interesting connections with the exotic chemistry of ultra-high magnetic fields (such as those in white dwarf stars), where the electronic cloud is dramatically compressed perpendicular to the magnetic field. \cite{Schmelcher_1990, Lange_2012, Schmelcher_2012}
In these extreme conditions, where magnetic effects compete with Coulombic forces, entirely new bonding paradigms emerge. \cite{Schmelcher_1990, Schmelcher_1997, Tellgren_2008, Tellgren_2009, Lange_2012, Schmelcher_2012, Boblest_2014, Stopkowicz_2015}
Atomic units are used throughout.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:eDFT}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{GOK-DFT}
The GOK ensemble energy is defined as follows:
\beq
E^{\bw}=\sum_{K\geq0}w_KE^{(K)},
\eeq
where the $K$th energy level $E^{(K)}$ [$K=0$ refers to the ground state] is the eigenvalue of the electronic Hamiltonian $\hat{H}=\hat{h}
%\sum^N_{i=1}\hat{h}(i)
+\hat{W}_{\rm ee}$, where $\hat{h}\equiv \sum^N_{i=1}\left(-\frac{1}{2}\nabla_{\br_i}^2+v_{\rm ne}(\br_i)\right)$ is the one-electron operator describing kinetic and nuclear attraction energies, and $\hat{W}_{\rm ee}$ is the electron repulsion operator.
The (positive) ensemble weights $w_K$ decrease with increasing index $K$. They are normalized, i.e.
\beq
w_0=1-\sum_{K>0}w_K,
\eeq
so that only the weights $\bw\equiv\left(w_1,w_2,\ldots w_K,\ldots\right)$ assigned to the excited states can vary independently.
For simplicity we will assume in the following that the energies are not degenerate. Note that the theory can be extended to multiplets simply by assigning the same ensemble weight to all degenerate states~\cite{}. In GOK-DFT, the ensemble energy is determined variationally as follows:
\beq\label{eq:var_ener_gokdft}
E^{{\bw}}=
\underset{\opGamma{{\bw}}}{\rm min}\left\{
{\rm
Tr}\left[\opGamma{{\bw}}\hat{h}\right]
+
{E}^{{\bw}}_{\rm
Hx}\left[n_{\opGamma{\bw}}\right]
+
{E}^{{\bw}}_{\rm
c}\left[n_{\opGamma{\bw}}\right]
\right\},
\eeq
where ${\rm
Tr}$ denotes the trace and the trial ensemble density matrix operator reads
\beq
\opGamma{{\bw}}=\sum_{K\geq 0}w_K\ket{\Phi^{(K)}}\bra{\Phi^{(K)}}.
\eeq
The determinants (or configuration state functions) $\Phi^{(K)}$ are all constructed from the same set of (ensemble Kohn--Sham) orbitals that is optimized variationally and the trial ensemble density is simply the weighted sum of the individual densities:
\beq\label{eq:KS_ens_density}
n_{\opGamma{\bw}}(\br)=\sum_{K\geq 0}w_Kn_{\Phi^{(K)}}(\br).
\eeq
As readily seen from Eq.~(\ref{eq:var_ener_gokdft}), both Hartree-exchange and
correlation energies are described with density functionals that are {\it weight-dependent}. We focus here on the (exact) Hx part which is defined as follows:
\beq\label{eq:exact_ens_Hx}
{E}^{{\bw}}_{\rm
Hx}[n]=\sum_{K\geq 0}w_K\bra{\Phi^{(K)}[n]}\hat{W}_{\rm
ee}\ket{\Phi^{(K)}[n]}
\eeq
where the KS wavefunctions fulfill the ensemble density constraint
\beq
\sum_{K\geq 0}w_Kn_{\Phi^{(K)}[n]}(\br)=n(\br).
\eeq
The (approximate) description of the correlation part is discussed in Sec.~\ref{sec:eDFA}.\\
In practice, one is not much interested in ensemble energies but rather in excitation energies or individual energy levels (for geometry optimizations, for example). The latter can be extracted exactly as follows~\cite{}:
\beq\label{eq:indiv_ener_from_ens}
E^{(I)}&=&E^{{\bw}}+\sum_{K>0}\left(\delta_{IK}-w_K\right)\dfrac{\partial
E^{{\bw}}}{\partial w_K},
\eeq
where, according to the {\it variational} ensemble energy expression of Eq.~(\ref{eq:var_ener_gokdft}), the derivative in $w_K$ can be evaluated from the minimizing KS wavefunctions $\Phi^{(K)}=\Phi^{(K),\bw}$ as follows:
\beq\label{eq:deriv_Ew_wk}
&&\dfrac{\partial
E^{{\bw}}}{\partial w_K}=\bra{\Phi^{(K)}}\hat{h}\ket{\Phi^{(K)}}-\bra{\Phi^{(0)}}\hat{h}\ket{\Phi^{(0)}}
\nonumber\\
&&+\Bigg[\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
Hx}\left[n\right]}{\delta
n({\br})}\left(n_{\Phi^{(K)}}(\br)-n_{\Phi^{(0)}}(\br)\right)
%\nonumber\\
%&&
+
%\left.
\dfrac{\partial {E}^{{\bw}}_{\rm
Hx}\left[n\right]}{\partial w_K}
%\right|_{n=n_{\opGamma{\bw}}}
\nonumber\\
&&+\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
c}\left[n\right]}{\delta
n({\br})}\left(n_{\Phi^{(K)}}(\br)-n_{\Phi^{(0)}}(\br)\right)
%\nonumber\\
%&&
+
%\left.
\dfrac{\partial {E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}
%\right|
\Bigg]_{n=n_{\opGamma{\bw}}}
.
\nonumber\\
\eeq
The Hx contribution to Eq.~(\ref{eq:deriv_Ew_wk}) can be rewritten as follows:
\beq\label{eq:_deriv_wk_Hx}
\left.\dfrac{\partial
}{\partial \xi_K}
\left({E}^{{\bxi}}_{\rm
Hx}\left[n^{\bxi,\bxi}\right]
-
{E}^{{\bw}}_{\rm
Hx}\left[n^{\bw,\bxi}\right]
\right)\right|_{\bxi=\bw},
\eeq
where $\bxi\equiv (\xi_1,\xi_2,\ldots,\xi_K,\ldots)$ and
\beq
n^{\bw,\bxi}(\br)=\sum_{K\geq 0}w_Kn_{\Phi^{(K),\bxi}}(\br).
\eeq
Since, for given ensemble weight $\bw$ and $\bxi$ values, the ensemble densities $n^{\bxi,\bxi}$ and $n^{\bw,\bxi}$ are generated from the {\it same} KS potential (which is unique up to a constant), it comes
from the exact expression in Eq.~(\ref{eq:exact_ens_Hx}) that
\beq
{E}^{{\bxi}}_{\rm
Hx}\left[n^{\bxi,\bxi}\right]=\sum_{K\geq 0}\xi_K\bra{\Phi^{(K),\bxi}}\hat{W}_{\rm
ee}\ket{\Phi^{(K),\bxi}}
\eeq
with $\xi_0=1-\sum_{K>0}\xi_K$
and
\beq
{E}^{{\bw}}_{\rm
Hx}\left[n^{\bw,\bxi}\right]=\sum_{K\geq 0}w_K\bra{\Phi^{(K),\bxi}}\hat{W}_{\rm
ee}\ket{\Phi^{(K),\bxi}},
\eeq
thus leading, according to Eqs.~(\ref{eq:deriv_Ew_wk}) and (\ref{eq:_deriv_wk_Hx}), to the simplified expression
\beq\label{eq:deriv_Ew_wk_simplified}
&&\dfrac{\partial
E^{{\bw}}}{\partial w_K}=\bra{\Phi^{(K)}}\hat{H}\ket{\Phi^{(K)}}-\bra{\Phi^{(0)}}\hat{H}\ket{\Phi^{(0)}}
\nonumber\\
&&+\Bigg[
\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
c}\left[n\right]}{\delta
n({\br})}\left(n_{\Phi^{(K)}}(\br)-n_{\Phi^{(0)}}(\br)\right)
%\nonumber\\
%&&
+
%\left.
\dfrac{\partial {E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}
%\right|
\Bigg]_{n=n_{\opGamma{\bw}}}
.
\nonumber\\
\eeq
Since the ensemble energy can be evaluated as follows:
\beq
E^{{\bw}}=\sum_{K\geq 0}w_K\bra{\Phi^{(K)}}\hat{H}\ket{\Phi^{(K)}}+{E}^{{\bw}}_{\rm
c}\left[n_{\opGamma{\bw}}\right],
\eeq
with $\Phi^{(K)}=\Phi^{(K),\bw}$ [note that, when the minimum is reached
in Eq.~(\ref{eq:var_ener_gokdft}), $n_{\opGamma{\bw}}=n^{\bw,\bw}$],
we finally recover from Eqs.~(\ref{eq:KS_ens_density}) and
(\ref{eq:indiv_ener_from_ens}) the {\it exact} expression of Ref.~\cite{} for the $I$th energy level:
\beq\label{eq:exact_ener_level_dets}
E^{(I)}&=&\bra{\Phi^{(I)}}\hat{H}\ket{\Phi^{(I)}}+{E}^{{\bw}}_{\rm
c}\left[n_{\opGamma{\bw}}\right]
\nonumber\\
&&+
\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
c}\left[n_{\opGamma{\bw}}\right]}{\delta
n({\br})}\left(n_{\Phi^{(I)}}(\br)-n_{\opGamma{\bw}}(\br)\right)
\nonumber\\
&&+
\sum_{K>0}\left(\delta_{IK}-w_K\right)\left.\dfrac{\partial {E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}
\right|
_{n=n_{\opGamma{\bw}}}.
\eeq
For implementation purposes, we will use in the rest of this work
(one-electron reduced) density matrices
as basic variables, rather than Slater determinants. If we expand the
(spin-) orbitals [from which the latter are constructed] in an atomic
orbital (AO) basis,
\beq
\varphi_p({\bfx})=\sum_{\mu}c_{{\mu p}}\AO{\mu}(\bfx),
\eeq
then the density matrix elements obtained from the
determinant $\Phi^{(K)}$ can be expressed as follows in the AO basis:
\beq
\Gamma_{\mu\nu}^{(K)}=\sum_{p\in (K)}c_{\mu p}c_{\nu p},
\eeq
where the summation runs over the spin-orbitals that are occupied in
$\Phi^{(K)}$. Note that the density of the $K$th KS state reads
\beq
n_{\bmg^{(K)}}(\br)=\sum_{\sigma=\alpha,\beta}\sum_{\mu\nu}\AO{\mu}({\br,
\sigma})\AO{\nu}(\br,\sigma){\Gamma}^{(K)}_{\mu\nu}.
\eeq
We can then construct the ensemble density matrix
and the ensemble density as follows:
\beq
{\bmg}^{{\bw}}=\sum_{K\geq 0}w_K{\bmg}^{(K)}
\eeq
and
\beq
n_{\bmg^\bw}({\br})=\sum_{\sigma=\alpha,\beta}\sum_{\mu\nu}\AO{\mu}({\br,\sigma})\AO{\nu}(\br,\sigma){\Gamma}^{\bw}_{\mu\nu},
\eeq
respectively. The exact energy level expression in Eq.~(\ref{eq:exact_ener_level_dets}) can be
rewritten as follows:
\beq
E^{(I)}&&={\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]
+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
\bmg^{(I)})
\nonumber\\
&&+{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]
+\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
\nonumber\\
&&+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left. \dfrac{\partial {E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
,
\eeq
where ${\bm
h}\equiv\left\{\langle\AO{\mu}\vert-\frac{1}{2}\nabla_{\br}^2+v_{\rm
ne}(\br)\vert\AO{\nu}\rangle\right\}_{\mu\nu}$ and ${\bm G}\equiv{\bm J}-{\bm K}$ denote
the Coulomb-exchange
integrals.
%%%%%%%%%%%%%%%%%%%%%
\iffalse%%%% Manu's derivation ...
\blue{
\beq
n^{\bw}({\br})&=&\sum_{K\geq 0}\sum_{\sigma=\alpha,\beta}{\tt
w}_Kn^{(K)}({\bfx})
\nonumber\\
&=&
\sum_{K\geq 0}\sum_{\sigma=\alpha,\beta}{\tt
w}_K\sum_{pq}\varphi_p({\bfx})\varphi_q({\bfx})\Gamma_{pq}^{(K)}
\nonumber\\
&=&
\sum_{\sigma=\alpha,\beta}
\sum_{K\geq 0}
{\tt
w}_K\sum_{p\in (K)}\varphi^2_p({\bfx})
\nonumber\\
&=&
\sum_{\sigma=\alpha,\beta}
\sum_{K\geq 0}
{\tt
w}_K
\sum_{\mu\nu}
\sum_{p\in (K)}c_{\mu p}c_{\nu p}\AO{\mu}({\bfx})\AO{\nu}({\bfx})
\nonumber\\
&=&\sum_{\sigma=\alpha,\beta}\sum_{\mu\nu}\AO{\mu}({\bfx})\AO{\nu}({\bfx}){\Gamma}^{\bw}_{\mu\nu}
\eeq
}
\fi%%%%%%%% end
%%%%%%%%%%%%%%%
%\subsection{Hybrid GOK-DFT}
%%%%%%%%%%%%%%%
\subsection{Approximations}
As Hartree and exchange energies cannot be separated in the
one-dimension systems considered in the rest of this work, we will substitute the Hartree--Fock
density-matrix-functional interaction energy,
\beq\label{eq:eHF-dens_mat_func}
W_{\rm
HF}\left[{\bmg}\right]=\frac{1}{2} \Tr(\bmg \, \bG \, \bmg),
\eeq
for the Hx density-functional energy in the variational energy
expression of Eq.~(\ref{eq:var_ener_gokdft}):
\beq
{\bmg}^{\bw}\approx\argmin_{{\bm\gamma}^{\bw}}
\Big\{
{\rm
Tr}\left[{\bm \gamma}^{{\bw}}{\bm h}\right]+W_{\rm
HF}\left[{\bm\gamma}^{\bw}\right]
+
{E}^{{\bw}}_{\rm
c}\left[n_{\bm\gamma^{\bw}}\right]
%+E^{{\bw}}_{\rm c}\left[n_{\hat{\Gamma}^{{\bw}}}\right]
\Big\}.
\nonumber\\
\eeq
Note that this approximation, where the ensemble density matrix is
optimized from a non-local exchange potential [rather than a local one,
as expected from Eq.~(\ref{eq:var_ener_gokdft})] is applicable to real
(three-dimension) systems. As readily seen from
Eq.~(\ref{eq:eHF-dens_mat_func}), {\it ghost-interaction} errors will be
introduced in the ensemble HF interaction energy:
In order to remove ghost interactions from the variational energy
expression used in the first step, we then employ the (in-principle-exact)
expression in Eq.~(\ref{eq:exact_ind_ener_OEP-like}). In this second
step, the response of the individual density matrices to weight
variations (last term on the right-hand side of
Eq.~(\ref{eq:exact_ind_ener_OEP-like})) is neglected. The complete GIC
procedure can be summarized as follows,
and
In order to compute (approximate) energy levels within generalized
GOK-DFT we use a two-step procedure. The first step consists in
optimizing variationally the ensemble density matrix according to
Eq.~(\ref{eq:var_princ_Gamma_ens}) with an approximate Hxc ensemble
functional where (i) the ghost-interaction correction functional $\overline{E}^{{\bw}}_{\rm
Hx}[n]$ in
Eq.~(\ref{eq:exact_GIC}) is
neglected, for simplicity, and (ii) the weight-dependent correlation
energy is described at the local density level of approximation.
At this
level of approximation, the two correlation functionals $\overline{E}^{{\bw}}_{\rm
c}[n]$ and ${E}^{{\bw}}_{\rm
c}[n]$ are actually identical and can be expressed as
\beq\label{eq:eLDA_corr_fun}
{E}^{{\bw}}_{\rm
c}[n]=\int d\br\;n(\br)\epsilon_{c}^{\bw}(n(\br)).
\eeq
More
details about the construction of such a functional will be given in the
following.
\beq
E^{(I)}&&\approx{\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]
+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
\bmg^{(I)})
\nonumber\\
&&+{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]
+\int d\br\,\dfrac{\delta {E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
\nonumber\\
&&+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left. \dfrac{\partial {E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
,
\eeq
thus leading to the final implementable expression [see Eq.~(\ref{eq:eLDA_corr_fun})]
\beq
E^{(I)}&&\approx{\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]
+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
\bmg^{(I)})
\nonumber\\
&&+\int d\br\,
{\epsilon}^{{\bw}}_{\rm
c}(n_{\bmg^{\bw}}(\br))\,n_{\bmg^{(I)}}(\br)
\nonumber\\
&&
+\int d\br\,\left.\dfrac{\partial {\epsilon}^{{\bw}}_{\rm
c}(n)}{\partial n}\right|_{n=n_{\bmg^{\bw}}(\br)}n_{\bmg^{\bw}}(\br)\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
\nonumber\\
&&
+\int d\br\,\sum_{K>0}\left(\delta_{IK}-w_K\right)n_{\bmg^{\bw}}(\br)\left.
\dfrac{\partial {\epsilon}^{{\bw}}_{\rm
c}(n)}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}(\br)}.
\eeq
\blue{$================================$}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory (old)}
\label{sec:eDFT_old}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Kohn--Sham formulation of GOK-DFT}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hybrid GOK-DFT}
\label{sec:geKS}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Since Hartree and exchange energy contributions cannot be separated in
the one-dimensional case, we introduce in the following an alternative
formulation of KS-eDFT where, in complete analogy with the generalized
KS scheme, a HF-like Hartree-exchange energy is employed. This
formulation is in principle exact and applicable to higher dimensions.
When $\bw=0$, the
conventional ground-state universal functional is recovered,
\beq
F^{\bw=0}[n]=F[n]=\underset{\Psi\rightarrow n}{\rm min}
\bra{\Psi}\hat{T}+\hat{W}_{\rm
ee}\ket{\Psi},
\eeq
where the ensemble reduces to a single wavefunction. In the latter case,
the HF-like expression (or a fraction of it, as usually done in
practical calculations) for the Hx energy can be introduced rigorously
into DFT by considering the following decomposition,
\beq\label{eq:generalized_KS-DFT_decomp}
F[n]&=&
\underset{\Phi\rightarrow n}{\rm min}
\bra{\Phi}\hat{T}+\hat{W}_{\rm
ee}\ket{\Phi}+\overline{E}_{\rm c}[n]
\nonumber\\
&=&
\underset{\bmg^\Phi\rightarrow n}{\rm min}
\left\{{\rm
Tr}\left[\bmg^\Phi{\bm t}\right]+W_{\rm
HF}\left[{\bmg}^{\Phi}\right]\right\}+
\overline{E}_{\rm c}[n]
,
\eeq
where ${\bm t}$ is the matrix representation of the one-electron kinetic
energy operator, $\bmg^\Phi$ is the one-electron reduced density
matrix (just referred to as density matrix in the following) obtained
from $\Phi$,
and
\beq
W_{\rm
HF}\left[{\bmg}\right]\equiv\frac{1}{2} \Tr(\bmg \, \bG \, \bmg)
\eeq
is the conventional density-matrix functional HF Hartree-exchange
energy. By analogy with Eq.~(\ref{eq:generalized_KS-DFT_decomp}), we
decompose the ensemble universal functional as follows:
\beq\label{eq:generalized_F_w}
F^{\bw}[n]&=&
\underset{\hat{\Gamma}^{{\bw}}\rightarrow n}{\rm min}\left\{{\rm
Tr}\left[\hat{\Gamma}^{{\bw}}\hat{T}\right]
+W_{\rm
HF}\left[{\bmg}^{\bw}\right]\right\}
+\overline{E}^{{\bw}}_{\rm
Hxc}[n]
\nonumber\\
&=&
\underset{{\bmg}^{{\bw}}\rightarrow n}{\rm min}
\left\{
{\rm Tr}
\left[{\bmg}^{\bw}{\bm t}\right]
+W_{\rm HF}\left[{\bmg}^{\bw}\right]
\right\}+
\overline{E}^{\bw}_{\rm Hxc}[n],
\eeq
where the minimization in Eq.~(\ref{eq:ens_LL_func}) has been restricted
to density matrix operators
\beq
\hat{\Gamma}^{{\bw}}=\sum_{K\geq 0}w_K\vert\Phi^{(K)}\rangle\langle\Phi^{(K)}\vert=\sum_{K\geq 0}w_K\hat{\Gamma}^{(K)}
\eeq
that are constructed from single Slater
determinants $\Phi^{(K)}$. Note that the density matrices
${\bmg}^{(K)}={\bmg}^{\Phi^{(K)}}$ are idempotent and diagonal in the
same spin-orbital basis). On the other hand, the ensemble
density matrix ${\bmg}^{{\bw}}=\sum_{K\geq 0}w_K{\bmg}^{(K)}$, which is a convex combination of the ${\bmg}^{(K)}$
matrices, is {\it not} idempotent, unless ${\bw}=0$.
Using an ensemble is, in this context,
analogous to assigning
fractional occupation numbers (which are determined from the ensemble
weights) to the KS orbitals.\\
Another issue with the use of
ensembles in DFT is the introduction of spurious ghost-interaction errors
(i.e. unphysical interactions between different states) into the
ensemble energy when inserting ${\bmg}^{{\bw}}$ into the HF
density-matrix functional Hx energy $W_{\rm
HF}\left[\bmg\right]$. This type of errors is specific to ensembles
which explains why, in constrast to ground-state DFT [see
Eq.~(\ref{eq:generalized_KS-DFT_decomp})], a complementary ensemble Hx
energy is needed to recover a ghost-interaction-free energy:
\beq\label{eq:exact_GIC}
\overline{E}^{{\bw}}_{\rm
Hx}[n]&=&
{\rm
Tr}\left[\hat{\Gamma}^{{\bw}}[n]\hat{W}_{\rm ee}\right]-W_{\rm
HF}\left[{\bmg}^{\bw}[n]\right]
\nonumber\\
&=&
\sum_{K\geq0}w_KW_{\rm
HF}\left[{\bmg}^{(K)}[n]\right]
-W_{\rm
HF}\left[{\bmg}^{\bw}[n]\right],
\eeq
where ${\bmg}^{\bw}[n]$ is the minimizing ensemble density matrix in
Eq.~(\ref{eq:generalized_F_w}) and, by construction, $\overline{E}^{{\bw}=0}_{\rm
Hx}[n]=0$. Consequently, the ensemble correlation functional can be
expressed as follows [see Eq.~(\ref{eq:generalized_F_w})]:
\beq
\overline{E}^{{\bw}}_{\rm
c}[n]&=&
\overline{E}^{\bw}_{\rm Hxc}[n]-\overline{E}^{{\bw}}_{\rm
Hx}[n]
\nonumber\\
&=&{\rm
Tr}\left[\hat{\gamma}^{{\bw}}[n]\left(\hat{T}+\hat{W}_{\rm
ee}\right)\right]
%\nonumber\\
%&&
-
{\rm
Tr}\left[\hat{\Gamma}^{{\bw}}[n]\left(\hat{T}+\hat{W}_{\rm
ee}\right)\right]
\nonumber\\
&=&
\sum_{K\geq 0}w_K\Bigg(\bra{\Psi^{(K)}[n]}\hat{T}+\hat{W}_{\rm
ee}\ket{\Psi^{(K)}[n]}
\nonumber\\
&&-\bra{\Phi^{(K)}[n]}\hat{T}+\hat{W}_{\rm
ee}\ket{\Phi^{(K)}[n]}\Bigg),
\eeq
where $\hat{\gamma}^{{\bw}}[n]$ and $\hat{\Gamma}^{{\bw}}[n]$ are the minimizing density matrix
operators in Eqs.~(\ref{eq:ens_LL_func}) and
(\ref{eq:generalized_KS-DFT_decomp}), respectively.\\
In eDFT, the ensemble energy $E^{{\bw}}=\sum_{K\geq
0}w_KE^{(K)}$ is obtained variationally as follows:
\beq
E^{{\bw}}=\underset{n}{\rm min}\Big\{
F^{\bw}[n]+\int d\br\,v_{\rm ext}(\br)n(\br)
\Big\}.
\eeq
Combining the latter expression with the decomposition in
Eq.~(\ref{eq:generalized_KS-DFT_decomp}) leads to
\beq
E^{{\bw}}=
\underset{n}{\rm min}\Bigg\{
\underset{{\bmg}^{{\bw}}\rightarrow n}{\rm min}\Big\{
{\rm
Tr}\left[{\bmg}^{{\bw}}{\bm h}\right]+W_{\rm
HF}\left[{\bmg}^{\bw}\right]
+
\overline{E}^{{\bw}}_{\rm
Hxc}\left[n_{\bmg^{\bw}}\right]
%+E^{{\bw}}_{\rm c}\left[n_{\hat{\Gamma}^{{\bw}}}\right]
\Big\}
\Bigg\}
\nonumber\\
\eeq
or, equivalently,
\beq\label{eq:var_princ_Gamma_ens}
E^{{\bw}}=
\underset{{\bmg}^{{\bw}}}{\rm min}\Big\{
{\rm
Tr}\left[{\bmg}^{{\bw}}{\bm h}\right]+W_{\rm
HF}\left[{\bmg}^{\bw}\right]
+
\overline{E}^{{\bw}}_{\rm
Hxc}\left[n_{\bmg^{\bw}}\right]
%+E^{{\bw}}_{\rm c}\left[n_{\hat{\Gamma}^{{\bw}}}\right]
\Big\}
,
\eeq
where $n_{\bmg^{\bw}}$ is the density obtained from the density matrix
${\bmg}^{\bw}$ and ${\bm h}={\bm t}+{\bm v}_{\rm ext}$ is the total one-electron
Hamiltonian matrix representation. When the minimum is reached, the
ensemble energy and its derivatives can be used to extract individual
ground- and excited-state energies as follows:\cite{Deur_2018b}
\beq
E^{(I)}&=&E^{{\bw}}+\sum_{K>0}\left(\delta_{IK}-w_K\right)\dfrac{\partial
E^{{\bw}}}{\partial w_K}.
\eeq
Since, according to the Hellmann--Feynman theorem, the ensemble energy
derivative reads
\beq
\dfrac{\partial E^{{\bw}}}{\partial w_K}&=&{\rm
Tr}\left[\left({\bmg}^{(K)}-{\bmg}^{(0)}\right){\bm h}\right]
\nonumber\\
&&+\Tr\left[\left({\bmg}^{(K)}-{\bmg}^{(0)}\right) \, \bG \, \bmg^{\bw}\right]
\nonumber\\
&&+
\int d\br\,\dfrac{\delta \overline{E}^{{\bw}}_{\rm
Hxc}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\left(n_{\bmg^{(K)}}(\br)-n_{\bmg^{(0)}}(\br)\right)
\nonumber\\
&&
+\left. \dfrac{\partial \overline{E}^{{\bw}}_{\rm
Hxc}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
,
\eeq
we finally obtain from Eqs.~(\ref{eq:var_princ_Gamma_ens}) and (\ref{eq:indiv_ener_from_ens}) the following in-principle-exact expressions for the
energy levels within the ensemble:
\beq
&&E^{(I)}=
{\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]+
\Tr\left[\left(\bmg^{(I)}-\dfrac{1}{2}\bmg^{\bw}\right) \, \bG \, \bmg^{\bw}\right]
\nonumber\\
&&+\overline{E}^{{\bw}}_{\rm
Hxc}\left[n_{\bmg^{\bw}}\right]
+\int d\br\,\dfrac{\delta \overline{E}^{{\bw}}_{\rm
Hxc}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
\nonumber\\
&&
+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left. \dfrac{\partial \overline{E}^{{\bw}}_{\rm
Hxc}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}.
\eeq
%+\Tr(\bmg^{(I)} \, \bG \, \bmg^{\bw})
%-\dfrac{1}{2}\Tr(\bmg^{\bw} \, \bG \, \bmg^{\bw})+...
\alert{
Note that
\beq
\overline{E}^{{\bw}}_{\rm
Hx}\left[n_{\bmg^{\bw}}\right]=
\frac{1}{2} \sum_{L\geq0}w_L \Tr(\bmg^{(L)} \, \bG \, \bmg^{(L)})
-\frac{1}{2}\Tr(\bmg^{\bw} \, \bG \, \bmg^{\bw})
\nonumber\\
\eeq
and
\beq
\left.\dfrac{\partial \overline{E}^{{\bw}}_{\rm
Hx}[n]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}&=&
\frac{1}{2} \Tr(\bmg^{(K)} \, \bG \, \bmg^{(K)})-\frac{1}{2}
\Tr(\bmg^{(0)} \, \bG \, \bmg^{(0)})
\nonumber\\
&&-\Tr\left[\left({\bmg}^{(K)}-{\bmg}^{(0)}\right) \, \bG \, \bmg^{\bw}\right]
+\ldots
\eeq
thus leading to
\beq
&&\overline{E}^{{\bw}}_{\rm
Hx}\left[n_{\bmg^{\bw}}\right]+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left. \dfrac{\partial \overline{E}^{{\bw}}_{\rm
Hx}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
\nonumber\\
&&=
\overline{E}^{{\bw}}_{\rm
Hx}\left[n_{\bmg^{\bw}}\right]+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \, \bmg^{(I)})
-\frac{1}{2} \sum_{L\geq0}w_L \Tr(\bmg^{(L)} \, \bG \, \bmg^{(L)})
\nonumber\\
&&-\Tr\left[\left({\bmg}^{(I)}-{\bmg}^{\bw}\right) \, \bG \, \bmg^{\bw}\right]
+\ldots
\nonumber\\
&&=\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
\bmg^{(I)})-\Tr\left[\left({\bmg}^{(I)}-\dfrac{1}{2}\bmg^{\bw}\right)
\, \bG \, \bmg^{\bw}\right]
+\ldots
\eeq
}
At the eLDA level:
\beq
\overline{E}^{{\bw}}_{\rm
Hxc}\left[n\right]\rightarrow\int d\br\,n(\br)\overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n(\br))
\eeq
\beq
\dfrac{\delta \overline{E}^{{\bw}}_{\rm
Hxc}\left[n\right]}{\delta
n({\br})}\rightarrow \overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n(\br))+n(\br)\left.\dfrac{\partial \overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n)}{\partial n}\right|_{n=n(\br)}
\eeq
\beq
&&E^{(I)}\rightarrow
{\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]+
\Tr\left[\left(\bmg^{(I)}-\dfrac{1}{2}\bmg^{\bw}\right) \, \bG \, \bmg^{\bw}\right]
\nonumber\\
&&
+\int d\br\,
\overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n_{\bmg^{\bw}}(\br))\,n_{\bmg^{(I)}}(\br)
\nonumber\\
&&
+\int d\br\,\left.\dfrac{\partial \overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n)}{\partial n}\right|_{n=n_{\bmg^{\bw}}(\br)}n_{\bmg^{\bw}}(\br)\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
\nonumber\\
&&
+\int d\br\,\sum_{K>0}\left(\delta_{IK}-w_K\right)n_{\bmg^{\bw}}(\br)\left.
\dfrac{\partial \overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n)}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}(\br)}.
\eeq
\alert{
or, equivalently,
\beq
&&E^{(I)}\rightarrow
{\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]+
\Tr\left[\left(\bmg^{(I)}-\dfrac{1}{2}\bmg^{\bw}\right) \, \bG \, \bmg^{\bw}\right]
\nonumber\\
&&
+\int d\br\,
\dfrac{\delta \overline{E}^{{\bw}}_{\rm
Hxc}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\,n_{\bmg^{(I)}}(\br)
\nonumber\\
&&
-\int d\br\,\left.\dfrac{\partial \overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n)}{\partial n}\right|_{n=n_{\bmg^{\bw}}(\br)}\Big(n_{\bmg^{\bw}}(\br)\Big)^2
\nonumber\\
&&
+\int d\br\,\sum_{K>0}\left(\delta_{IK}-w_K\right)n_{\bmg^{\bw}}(\br)\left.
\dfrac{\partial \overline{\epsilon}^{{\bw}}_{\rm
Hxc}(n)}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}(\br)}.
\eeq
}
\subsection{Exact ensemble exchange in hybrid GOK-DFT}
In the exact theory, the minimizing density matrix in
Eq.~(\ref{eq:var_princ_Gamma_ens}) is such that
\beq
{\bmg}^{(K)}[n_{{\bmg}^{{\bw}}}]={\bmg}^{(K)},\hspace{0.2cm}\forall
K\geq0,
\eeq
and therefore
\beq
{\bmg}^{{\bw}}\left[n_{{\bmg}^{{\bw}}}\right]={\bmg}^{{\bw}}.
\eeq
Combining the latter Eqs. with
Eqs. (\ref{eq:exact_GIC}), (\ref{eq:var_princ_Gamma_ens}) leads to
the final ensemble energy expression
\beq\label{eq:exact_Eens_EEXX}
E^{{\bw}}={\rm
Tr}\left[{\bmg}^{{\bw}}{\bm h}\right]+\frac{1}{2} \sum_{L\geq0}w_L
\Tr(\bmg^{(L)} \, \bG \, \bmg^{(L)})
+\overline{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right].
\nonumber\\
\eeq
Note that
\beq
E^{{\bw}}&\neq& \underset{\left\{{\bmg}^{(L)}\right\}_{L\geq 0}}{\rm min}\Big\{
{\rm
Tr}\left[{\bmg}^{{\bw}}{\bm h}\right]+\frac{1}{2} \sum_{L\geq0}w_L
\Tr(\bmg^{(L)} \, \bG \, \bmg^{(L)})
\nonumber\\
&&+\overline{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]
\Bigg\}
\eeq
For $K>0$
\beq\label{eq:XE_EEXX}
&&\dfrac{\partial E^{{\bw}}}{\partial w_K}=
{\rm
Tr}\left[\left({\bmg}^{(K)}-{\bmg}^{(0)}\right){\bm h}\right]
\nonumber\\
&&+\frac{1}{2}\Tr(\bmg^{(K)} \, \bG \, \bmg^{(K)})
-\frac{1}{2}\Tr(\bmg^{(0)} \, \bG \, \bmg^{(0)})
\nonumber\\
&&
+\int d\br\,\dfrac{\delta \overline{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\left(n_{\bmg^{(K)}}(\br)-n_{\bmg^{(0)}}(\br)\right)
+\left. \dfrac{\partial \overline{E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
\nonumber\\
&&+\sum_{L\geq0}w_L{\rm
Tr}\left[\dfrac{\partial\bmg^{(L)}}{\partial w_K}{\bm h}\right]
%\nonumber\\
%&&
+\sum_{L\geq0}w_L
\Tr(\bmg^{(L)} \, \bG \, \dfrac{\partial\bmg^{(L)}}{\partial w_K})
\nonumber\\
&&
+\sum_{L\geq0}w_L\int d\br\,\dfrac{\delta \overline{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}n_{\frac{\partial \bmg^{(L)}}{\partial w_K}}(\br).
\eeq
If we introduce individual Fock matrices
\beq
{\bm F}^{(L)}={\bm h}+\bG \,\bmg^{(L)}+\overline{\bm v}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right],
\eeq
the last three terms can be simply rewritten as
\beq
\sum_{L\geq0}w_L{\rm
Tr}\left[{\bm F}^{(L)}\frac{\partial \bmg^{(L)}}{\partial w_K}\right].
\eeq
According to Eqs.~(\ref{eq:indiv_ener_from_ens}),
(\ref{eq:exact_Eens_EEXX}), and (\ref{eq:XE_EEXX}),
\beq\label{eq:exact_ind_ener_OEP-like}
E^{(I)}&&={\rm
Tr}\left[{\bmg}^{(I)}{\bm h}\right]
+\frac{1}{2} \Tr(\bmg^{(I)} \, \bG \,
\bmg^{(I)})
\nonumber\\
&&+\overline{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]
+\int d\br\,\dfrac{\delta \overline{E}^{{\bw}}_{\rm
c}\left[n_{\bmg^{\bw}}\right]}{\delta
n({\br})}\left(n_{\bmg^{(I)}}(\br)-n_{\bmg^{\bw}}(\br)\right)
\nonumber\\
&&+\sum_{K>0}\left(\delta_{IK}-w_K\right)\left. \dfrac{\partial \overline{E}^{{\bw}}_{\rm
c}\left[n\right]}{\partial w_K}\right|_{n=n_{\bmg^{\bw}}}
\nonumber\\
&&
+\sum_{K>0}\left(\delta_{IK}-w_K\right)\sum_{L\geq0}w_L{\rm
Tr}\left[{\bm F}^{(L)}\frac{\partial \bmg^{(L)}}{\partial w_K}\right]
.
\eeq
\subsection{Ensemble
correlation LDA and ghost interaction correction
}
\alert{Secs. \ref{sec:KS-eDFT}-\ref{sec:E_I} should maybe be moved to an appendix or merged
with the theory section (?)}
%%%%%%%%%%%%%%%%
\section{Implementation}
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{KS-eDFT for excited states}
\label{sec:KS-eDFT}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Here, we explain how to perform a self-consistent KS calculation for ensembles (KS-eDFT) in the context of excited states.
In order to take into account both single and double excitations simultaneously, we consider a three-state ensemble including:
(i) the ground state ($I=0$), (ii) the first singly-excited state ($I=1$), and (iii) the first doubly-excited state ($I=2$).
Generalization to a larger number of states is straightforward and is left for future work.
By definition, the ensemble energy is
\begin{equation}
\E{}{\bw} = (1 - \ew{1} - \ew{2}) \E{}{(0)} + \ew{1} \E{}{(1)} + \ew{2} \E{}{(2)}.
\end{equation}
The $\E{}{(I)}$'s are individual energies, while $\ew{1}$ and $\ew{2}$ are the weights assigned to the single and double excitation, respectively.
To ensure the GOK variational principle, \cite{Gross_1988a} the weights must fulfil the following conditions:
$0 \le \ew{1} \le 1/3$ and $\ew{2} \le \ew{1} \le (1-\ew{2})/2$.
Note that, in order to extract individual energies from a single KS-eDFT calculation [see Subsec.~\ref{sec:E_I}], the weights must remain independent.
By construction, the excitation energies are
\begin{equation}
\label{eq:Ex}
\Ex{(I)} = \pdv{\E{}{(I)}}{\ew{I}} = \E{}{(I)} - \E{}{(0)}.
\end{equation}
In the following, the orbitals $\MO{p}{\bw}(\br)$ are defined as linear combination of basis functions $\AO{\mu}(\br)$, such as
\begin{equation}
\MO{p}{\bw}(\br) = \sum_{\mu=1}^{\Nbas} \cMO{\mu p}{\bw} \, \AO{\mu}(\br).
\end{equation}
Within the self-consistent KS-eDFT process, one is looking for the following weight-dependent density matrix:
\begin{equation}
\label{eq:Gamma}
\bGamma{\bw} = (1-\ew{1}-\ew{2}) \bGamma{(0)} - \ew{1} \bGamma{(1)} - \ew{2} \bGamma{(2)},
\end{equation}
where $\bw = (\ew{1},\ew{2})$ and $\bGamma{(I)}$ is the $I$th-state density matrix with elements
\begin{equation}
\label{eq:eGamma}
\eGamma{\mu\nu}{(I)} = \sum_{i=1}^{\Nel-I} \cMO{\mu i}{\bw} \cMO{\nu i}{\bw} + \sum_{a=\Nel+1}^{\Nel+I} \cMO{\mu a}{\bw} \cMO{\nu a}{\bw}.
\end{equation}
The coefficients $\cMO{\mu p}{\bw}$ used to construct the density matrix $\bGamma{\bw}$ in Eq.~\eqref{eq:Gamma} are obtained by diagonalizing the following Fock matrix
\begin{multline}
\label{eq:F}
\eF{\mu\nu}{\bw}
= \eHc{\mu\nu} + \sum_{\la\si} \eGamma{\la\si}{\bw} \eG{\mu\nu\la\si}
\\
+ \int \left. \fdv{\bE{Hxc}{\bw}[\n{}{}]}{\n{}{}(\br)} \right|_{\n{}{} = \n{}{\bw}(\br)} \AO{\mu}(\br) \AO{\nu}(\br) d\br,
\end{multline}
which itself depends on $\bGamma{\bw}$.
In Eq.~\eqref{eq:F}, $\hHc$ is the core Hamiltonian (including kinetic and electron-nucleus attraction terms), $\eG{\mu\nu\la\si} = (\mu\nu|\la\si) - (\mu\si|\la\nu)$,
\begin{equation}
(\mu\nu|\la\si) = \iint \frac{\AO\mu(\br_1) \AO\nu(\br_1) \AO\la(\br_2) \AO\si(\br_2)}{\abs{\br_1 - \br_2}} d\br_1 d\br_2
\end{equation}
are two-electron repulsion integrals,
$\bE{Hxc}{\bw}[\n{}{}(\br)] = \n{}{}(\br) \be{Hxc}{\bw}[\n{}{}(\br)]$ and $\be{Hxc}{\bw}[\n{}{}(\br)]$ is the weight-dependent Hartree-exchange-correlation functional to be built in the present study.
The one-electron ensemble density is
\begin{equation}
\n{}{\bw}(\br) = \sum_{\mu\nu} \AO{\mu}(\br) \, \eGamma{\mu\nu}{\bw} \, \AO{\nu}(\br),
\end{equation}
with a similar expression for $\n{}{(I)}(\br)$, while the ensemble energy reads
\begin{equation}
\label{eq:Ew}
\E{}{\bw}
= \Tr(\bGamma{\bw} \, \bHc)
+ \frac{1}{2} \Tr(\bGamma{\bw} \, \bG \, \bGamma{\bw})
% \\
% + \int \e{c}{\bw}[\n{}{\bw}(\br)] \n{}{\bw}(\br) d\br.
+ \int \bE{Hxc}{\bw}[\n{}{\bw}(\br)] d\br.
\end{equation}
The self-consistent process described above is carried on until $\max \abs{\bF{\bw} \, \bGamma{\bw} \, \bS - \bS \, \bGamma{\bw} \, \bF{\bw}} < \tau$, where $\tau$ is a user-defined threshold and $\eS{\mu\nu} = \braket{\AO{\mu}}{\AO{\nu}}$ are elements of the overlap matrix $\bS$.
Note that because the second term of the RHS of Eq.~\eqref{eq:Ew} is quadratic in $\bGamma{\bw}$, the weight-dependent energy contains the so-called ghost interaction which makes the ensemble energy non linear. \cite{Gidopoulos_2002, Pastorczak_2014, Alam_2016, Alam_2017, Gould_2017}
Below, we propose a ghost-interaction correction (GIC) in order to minimize this error.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Extracting individual energies}
\label{sec:E_I}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Following Deur and Fromager, \cite{Deur_2018b} it is possible to extract individual energies, $\E{}{(I)}$, from the ensemble energy [see Eq.~\eqref{eq:Ew}] as follows:
\begin{multline}
\E{}{(I)} = \Tr(\bGamma{(I)} \, \bHc) + \frac{1}{2} \Tr(\bGamma{(I)} \, \bG \, \bGamma{(I)})
\\
+ \int \left. \fdv{\bE{Hxc}{\bw}[\n{}{}]}{\n{}{}(\br)} \right|_{\n{}{} = \n{}{\bw}(\br)} \n{}{(I)}(\br) d\br
+ \LZ{Hxc}{} + \DD{Hxc}{(I)}.
\end{multline}
Note that a \emph{single} KS-eDFT calculation is required to extract the three individual energies.
\alert{Mention LIM?}
The (state-independent) Levy-Zahariev shift and the so-called derivative discontinuity are given by
\begin{align}
\LZ{Hxc}{} & = - \int \left. \fdv{\be{Hxc}{\bw}[\n{}{}]}{\n{}{}(\br)} \right|_{\n{}{} = \n{}{\bw}(\br)} \n{}{\bw}(\br)^2 d\br,
\\
\DD{Hxc}{(I)} & = \sum_{J>0} (\delta_{IJ} - \ew{J}) \int \left. \pdv{\be{Hxc}{\bw}[\n{}{}]}{\ew{J}}\right|_{\n{}{} = \n{}{\bw}(\br)} \n{}{\bw}(\br) d\br.
\end{align}
Because the Levy-Zahariev shift is state independent, it does not contribute to excitation energies [see Eq.~\eqref{eq:Ex}].
The only remaining piece of information to define at this stage is the weight-dependent Hartree-exchange-correlation functional $\be{Hxc}{\bw}(\n{}{})$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Density-functional approximations for ensembles}
\label{sec:eDFA}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We decompose the weight-dependent functional as
\begin{equation}
\be{Hxc}{\bw}(\n{}{}) = \be{Hx}{\bw}(\n{}{}) + \be{c}{\bw}(\n{}{}),
\end{equation}
where $\be{Hx}{\bw}(\n{}{})$ is a weight-dependent Hartree-exchange functional designed to correct the ghost interaction \cite{Gidopoulos_2002, Pastorczak_2014, Alam_2016, Alam_2017, Gould_2017} [see Subsec.~\ref{sec:GIC}] and $\be{c}{\bw}(\n{}{})$ is a weight-dependent correlation functional [see Subsec.~\ref{sec:Ec}].
The construction of these two functionals is described below.
Note that, because we consider strict 1D systems, one cannot decompose further the Hartree-exchange contribution as each component diverges independently but their sum is finite. \cite{Astrakharchik_2011, Lee_2011a, Loos_2012, Loos_2013, Loos_2013a}
Most of the standard local and semi-local DFAs rely on the infinite uniform electron gas (UEG) model (also known as jellium). \cite{ParrBook, Loos_2016}
One major drawback of the jellium paradigm, when it comes to develop eDFAs, is that the ground and excited states cannot be easily identified like in a molecule. \cite{Gill_2012, Loos_2012, Loos_2014a, Loos_2014b, Agboola_2015, Loos_2017a}
Moreover, because the infinite UEG model is a metal, it is gapless, which means that both the fundamental and optical gaps are zero.
From this point of view, using finite UEGs \cite{Loos_2011b, Gill_2012} (which have, like an atom, discrete energy levels) to construct eDFAs can be seen as more relevant. \cite{Loos_2014a, Loos_2014b, Loos_2017a}
Here, we propose to construct a weight-dependent eDFA for the calculations of excited states in 1D systems.
As a finite uniform electron gas, we consider the ringium model in which electrons move on a perfect ring (i.e., a circle). \cite{Loos_2012, Loos_2013a, Loos_2014b}
The most appealing feature of ringium (regarding the development of functionals in the context of eDFT) is the fact that both ground- and excited-state densities are uniform.
As a result, the ensemble density will remain constant (and uniform) as the ensemble weights vary.
This is a necessary condition for being able to model derivative discontinuities.
The present weight-dependent eDFA is specifically designed for the calculation of excitation energies within eDFT.
As mentioned previously, we consider a three-state ensemble including the ground state ($I=0$), the first singly-excited state ($I=1$), and the first doubly-excited state ($I=2$) of the (spin-polarized) two-electron ringium system.
All these states have the same (uniform) density $\n{}{} = 2/(2\pi R)$ where $R$ is the radius of the ring where the electrons are confined.
We refer the interested reader to Refs.~\onlinecite{Loos_2012, Loos_2013a, Loos_2014b} for more details about this paradigm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Ghost-interaction correction}
\label{sec:GIC}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The GIC weight-dependent Hartree-exchange functional is defined as
\begin{multline}
\be{Hx}{\bw}(\n{}{\bw}) = (1-\sum_{I>0} \ew{I}) \be{Hx}{}(\n{}{(0)}) + \sum_{I>0} \ew{I} \be{Hx}{}(\n{}{(I)})
\\
- \be{Hx}{(I)}(\n{}{\bw}),
\end{multline}
where
\begin{equation}
\be{Hx}{}(\n{}{}) = \iint \frac{\n{}{}(\br_1) \n{}{}(\br_2) - \n{}{}(\br_1,\br_2)^2}{r_{12}} d\br_1 d\br_2,
\end{equation}
and
\begin{equation}
\n{}{(I)}(\omega) = (\pi R)^{-1} \cos[(I+1) \omega/2]
\end{equation}
is the first-order density matrix with $\omega$ the interelectronic angle.
It yields
\begin{equation}
\be{Hx}{}(\n{}{}) = \n{}{} \qty[ a_1 \ew{1} (\ew{1} - 1) + a_2 \ew{1} \ew{2} + a_3 \ew{2} (\ew{2} - 1)],
\end{equation}
with $a_1 = 2 \ln 2 - 1/3$, $a_2 = 8/3$ and $a_3 = 32/15$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Weight-dependent correlation functional}
\label{sec:Ec}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Based on highly-accurate calculations (see {\SI} for additional details), one can write down, for each state, an accurate analytical expression of the reduced (i.e., per electron) correlation energy \cite{Loos_2013a, Loos_2014a} via the following Pad\'e approximant
\begin{equation}
\label{eq:ec}
\e{c}{(I)}(\n{}{}) = \frac{c_1^{(I)}\,\n{}{}}{\n{}{} + c_2^{(I)} \sqrt{\n{}{}} + c_3^{(I)}},
\end{equation}
where the $c_k^{(I)}$'s are state-specific fitting parameters, which are provided in Table \ref{tab:OG_func}.
The value of $c_1^{(I)}$ is obtained via the exact high-density expansion of the correlation energy. \cite{Loos_2013a, Loos_2014a}
Equation \eqref{eq:ec} provides three state-specific correlation DFAs based on a two-electron system.
Combining these, one can build a three-state weight-dependent correlation eDFA:
\begin{equation}
\label{eq:ecw}
\e{c}{\bw}(\n{}{}) = (1-\ew{1}-\ew{2}) \e{c}{(0)}(\n{}{}) + \ew{1} \e{c}{(1)}(\n{}{}) + \ew{2} \e{c}{(2)}(\n{}{}).
\end{equation}
%%% TABLE 1 %%%
\begin{table*}
\caption{
\label{tab:OG_func}
Parameters of the correlation DFAs defined in Eq.~\eqref{eq:ec}.}
% \begin{ruledtabular}
\begin{tabular}{lcddd}
\hline\hline
State & $I$ & \tabc{$c_1^{(I)}$} & \tabc{$c_2^{(I)}$} & \tabc{$c_3^{(I)}$} \\
\hline
Ground state & $0$ & -0.0137078 & 0.0538982 & 0.0751740 \\
Singly-excited state & $1$ & -0.0238184 & 0.00413142 & 0.0568648 \\
Doubly-excited state & $2$ & -0.00935749 & -0.0261936 & 0.0336645 \\
\hline\hline
\end{tabular}
% \end{ruledtabular}
\end{table*}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{LDA-centered functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to make the two-electron-based eDFA defined in Eq.~\eqref{eq:ecw} more universal and to ``center'' it on the jellium reference (as commonly done in DFT), we propose to \emph{shift} it as follows:
\begin{equation}
\label{eq:becw}
\be{c}{\bw}(\n{}{}) = (1-\ew{1}-\ew{2}) \be{c}{(0)}(\n{}{}) + \ew{1} \be{c}{(1)}(\n{}{}) + \ew{2} \be{c}{(2)}(\n{}{}),
\end{equation}
where
\begin{equation}
\be{c}{(I)}(\n{}{}) = \e{c}{(I)}(\n{}{}) + \e{c}{\text{LDA}}(\n{}{}) - \e{c}{(0)}(\n{}{}).
\end{equation}
The local-density approximation (LDA) correlation functional,
\begin{equation}
\e{c}{\text{LDA}}(\n{}{}) = c_1^\text{LDA} \, F\qty[1,\frac{3}{2},c_3^\text{LDA}, \frac{c_1^\text{LDA}(1-c_3^\text{LDA})}{c_2^\text{LDA}} {\n{}{}}^{-1}],
\end{equation}
specifically designed for 1D systems in Ref.~\onlinecite{Loos_2013} as been used, where $F(a,b,c,x)$ is the Gauss hypergeometric function, \cite{NISTbook} and
\begin{align}
c_1^\text{LDA} & = - \frac{\pi^2}{360},
&
c_2^\text{LDA} & = \frac{3}{4} - \frac{\ln{2\pi}}{2},
&
c_3^\text{LDA} & = 2.408779.
\end{align}
Equation \eqref{eq:becw} can be recast
\begin{equation}
\label{eq:eLDA}
\begin{split}
\be{c}{\bw}(\n{}{})
& = \e{c}{\text{LDA}}(\n{}{})
\\
& + \ew{1} \qty[\e{c}{(1)}(\n{}{})-\e{c}{(0)}(\n{}{})] + \ew{2} \qty[\e{c}{(2)}(\n{}{})-\e{c}{(0)}(\n{}{})],
\end{split}
\end{equation}
which nicely highlights the centrality of the LDA in the present eDFA.
In particular, $\be{c}{(0,0)}(\n{}{}) = \e{c}{\text{LDA}}(\n{}{})$.
Consequently, in the following, we name this correlation functional ``eLDA'' as it is a natural extension of the LDA for ensembles.
This procedure can be theoretically justified by the generalized adiabatic connection formalism for ensembles (GACE) which was originally derived by Franck and Fromager. \cite{Franck_2014}
Within this in-principle-exact formalism, the (weight-dependent) correlation energy of the ensemble is constructed from the (weight-independent) ground-state functional (such as the LDA), yielding Eq.~\eqref{eq:eLDA}.
This is a crucial point as we intend 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).
Finally, we note that, by construction,
\begin{equation}
\left. \pdv{\be{c}{\bw}[\n{}{}]}{\ew{J}}\right|_{\n{}{} = \n{}{\bw}(\br)} = \be{c}{(J)}[\n{}{\bw}(\br)] - \be{c}{(0)}[\n{}{\bw}(\br)].
\end{equation}
\alert{As shown by Gould and Pittalis, comment on density- and and state-driven errors. \cite{Gould_2019}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Having defined the eLDA functional in the previous section [see Eq.~\eqref{eq:eLDA}], we now turn to its validation.
Our testing playground for the validation of the eLDA functional is the ubiquitous ``electrons in a box'' model where $\Nel$ electrons are confined in a 1D box of length $L$, a family of systems that we call $\Nel$-boxium.
In particular, we investigate systems where $L$ ranges from $\pi/8$ to $8\pi$ and $2 \le \Nel \le 7$.
\alert{Comment on the quality of these density: density- and functional-driven errors?}
These inhomogeneous systems have non-trivial electronic structure properties which can be tuned by varying the box length.
For small $L$, the system is weakly correlated, while strong correlation effects dominate in the large-$L$ regime. \cite{Rogers_2017,Rogers_2016}
We use as basis functions the (orthonormal) orbitals of the one-electron system, i.e.
\begin{equation}
\AO{\mu}(x) =
\begin{cases}
\sqrt{2/L} \cos(\mu \pi x/L), & \mu \text{ is odd,}
\\
\sqrt{2/L} \sin(\mu \pi x/L), & \mu \text{ is even,}
\end{cases}
\end{equation}
with $ \mu = 1,\ldots,\Nbas$ and $\Nbas = 30$ for all calculations.
For the self-consistent calculations (such as HF, KS or eKS), the convergence threshold has been set to $\tau = 10^{-7}$.
For KS and eKS calculations, a Gauss-Legendre quadrature is employed to compute numerical integrals.
In order to test the present eLDA functional we have performed various sets of calculations.
To get reference excitation energies for both the single and double excitations, we have performed full configuration interaction (FCI) calculations with the Knowles-Handy FCI program described in Ref.~\onlinecite{Knowles_1989}.
For the single excitations, we have also performed time-dependent HF (TDHF), configuration interaction singles (CIS) and TDLDA calculations. \cite{Dreuw_2005}
For TDLDA, the validity of the Tamm-Dancoff approximation (TDA) has been also tested.
Concerning the eKS calculations, two sets of weight have been tested: the zero-weight limit where $\bw = (0,0)$ and the equi-ensemble (or state-averaged) limit where $\bw = (1/3,1/3)$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In Fig.~\ref{fig:EvsL}, we report the error (in \%) in excitation energies (compared to FCI) for various methods and box sizes in the case of 5-boxium (i.e., $\Nel = 5$).
Similar graphs are obtained for the other $\Nel$ values and they can be found --- alongside the numerical data associated with each method --- in the {\SI}.
In the weakly correlated regime (i.e., small $L$), all methods provide accurate estimates of the excitation energies.
When the box gets larger, they start to deviate.
For the single excitation, TDHF is extremely accurate over the whole range of $L$ values, while CIS is slightly less accurate and starts to overestimate the excitation energy by a few percent at $L=8\pi$.
TDLDA yields larger errors at large $L$ by underestimating the excitation energies.
TDA-TDLDA slightly corrects this trend thanks to error compensation.
Concerning the eLDA functional, our results clearly evidences that the equi-weights [i.e., $\bw = (1/3,1/3)$] excitation energies are much more accurate than the ones obtained in the zero-weight limit [i.e., $\bw = (0,0)$].
This is especially true for the single excitation which is significantly improved by using state-averaged weights.
The effect on the double excitation is less pronounced.
Overall, one clearly sees that, with state-averaged weights, the eLDA functional yields accurate excitation energies for both single and double excitations.
This conclusion is verified for smaller and larger number of electrons (see {\SI}).
\alert{Shall I test the one-electron system for self-interaction?}
%%% FIG 1 %%%
\begin{figure}
\includegraphics[width=\linewidth]{EvsL_5}
\caption{
\label{fig:EvsL}
Error with respect to FCI in single and double excitation energies for 5-boxium for various methods and box length $L$.
Graphs for additional values of $\Nel$ can be found as {\SI}.
}
\end{figure}
%%% %%% %%%
Figure \ref{fig:EvsN} reports the error (in \%) in excitation energies, for the same methods, as a function of $\Nel$ and fixed $L$ (in this case $L=\pi$).
The graphs associated with other $L$ values are reported as {\SI}.
Again, the graph for $L=\pi$ is quite typical and we draw similar conclusions as in the previous paragraph: irrespectively of the number of electrons, the eLDA functional with state-averaged weights is able to accurately model single and double excitations.
As a rule of thumb, we see that eLDA single excitations are of the same quality as the ones obtained in the linear response formalism (such as TDHF or TDLDA), while double excitations only deviates from the FCI values by a few tenth of percent for $L=\pi$, an error of the same order as CIS or TDA-TDLDA.
Even for larger boxes, the discrepancy between FCI and eLDA for double excitations is only a few percent.
%%% FIG 2 %%%
\begin{figure}
\includegraphics[width=\linewidth]{EvsN_1}
\caption{
\label{fig:EvsN}
Error with respect to FCI in single and double excitation energies for $\Nel$-boxium for various methods and number of electrons $\Nel$ at $L=\pi$.
Graphs for additional values of $L$ can be found as {\SI}.
}
\end{figure}
%%% %%% %%%
\alert{Need further discussion on DD and LZ shift. Linearity of energy wrt weights?}
\alert{For small $L$, the single and double excitations are ``pure''. In other words, the excitation is dominated by a single reference Slater determinant.
However, when the box gets larger, there is a strong mixing between different degree of excitations.
In particular, the single and double excitations strongly mix.
This is clearly evidenced if one looks at the weights of the different configurations in the FCI wave function.
In one hand, if one does construct a eDFA with a single state (either single or double), one clearly sees that the results quickly deteriorates when the box gets larger.
On the other hand, building a functional which does mix singles and doubles corrects this by allowing configuration mixing.}
\alert{It might be useful to add eHF results where one switch off the correlation part.
For both zero weight and state-averaged weights?
It would highlight the contribution of the derivative discontinuity.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Concluding remarks}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the present article, we have constructed a weight-dependent three-state DFA in the context of ensemble DFT.
This eDFA delivers accurate excitation energies for both single and double excitations.
Generalization to a larger number of states is straightforward and will be investigated in future work.
Using similar ideas, a three-dimensional version \cite{Loos_2009,Loos_2009c,Loos_2010,Loos_2010d,Loos_2017a} of the present eDFA is currently under development to model excited states in molecules and solids.
Similar to the present excited-state methodology for ensembles, one can easily design a local eDFA for the calculations of the ionization potential, electron affinity, and fundamental gap.
This can be done by constructing DFAs for the one- and three-electron ground state systems, and combining them with the two-electron DFA in complete analogy with Eqs.~\eqref{eq:ec} and \eqref{eq:ecw}.
However, as shown by Senjean and Fromager, \cite{Senjean_2018} one must modify the weights accordingly in order to maintain a constant density.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supplementary material}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
See {\SI} for the additional details about the construction of the functionals, raw data and additional graphs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
E.~F.~thanks the \textit{Agence Nationale de la Recherche} (MCFUNEX project, Grant No.~ANR-14-CE06-0014-01) for funding.
\end{acknowledgements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% REMOVED FROM THE MAIN TEXT by Manu %%%%%%%%%%%%
\iffalse%%%%
Indeed,
\beq
\left[{\bmg}^{{\bw}}\right]^2&=&\sum_{K,L\geq
0}w_Kw_L{\bmg}^{(K)}{\bmg}^{(L)}
\nonumber\\
&=&\sum_{K\geq
0}\left(w_K\right)^2{\bmg}^{(K)}+\sum_{K\neq L\geq
0}w_Kw_L{\bmg}^{(K)}{\bmg}^{(L)}
\nonumber\\
&=&
{\bmg}^{{\bw}}+\sum_{K,L\geq
0}w_K\left(w_L-\delta_{KL}\right){\bmg}^{(K)}{\bmg}^{(L)}
\nonumber\\
&=&{\bmg}^{{\bw}}+w_0{\bmg}^{(0)}\times\sum_{K>0}w_K\left(2{\bmg}^{(K)}-1\right)
\nonumber\\
&&+\sum_{K, L >0
}w_K\left(w_L-\delta_{KL}\right){\bmg}^{(K)}{\bmg}^{(L)}
\nonumber\\
&\neq&{\bmg}^{{\bw}}
.
\eeq
%%%% End -- REMOVED FROM THE MAIN TEXT by Manu %%%%%%%%%%%%
\fi%%%
\bibliography{eDFT}
\end{document}
Reference in New Issue View Git Blame Copy Permalink