loos/FarDFT
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
0ed66dcf65
FarDFT/Manuscript/FarDFT.tex

801 lines
48 KiB
TeX
Raw Blame History

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amsmath,amssymb,amsfonts,physics,mhchem}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}
\usepackage{grffile}
\usepackage[
colorlinks=true,
citecolor=blue,
breaklinks=true
]{hyperref}
\urlstyle{same}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\cloclo}[1]{\textcolor{purple}{#1}}
\newcommand{\bruno}[1]{\textcolor{blue}{Bruno: #1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashCM}[1]{\textcolor{purple}{\sout{#1}}}
\newcommand{\cmark}{\color{green}{\text{\ding{51}}}}
\newcommand{\xmark}{\color{red}{\text{\ding{55}}}}
%useful stuff
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\ra}{\rightarrow}
\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}
\newcommand{\cD}{\mathcal{D}}
% operators
\newcommand{\hH}{\Hat{H}}
\newcommand{\hHc}{\Hat{h}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hWee}{\Hat{W}_\text{ee}}
\newcommand{\hGam}[1]{\Hat{\Gamma}^{#1}}
\newcommand{\hgam}[1]{\Hat{\gamma}^{#1}}
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\hVne}{\Hat{V}_\text{ne}}
\newcommand{\vne}{v_\text{ne}}
% functionals, potentials, densities, etc
\newcommand{\F}[2]{F_{#1}^{#2}}
\newcommand{\Ts}[1]{T_\text{s}^{#1}}
\newcommand{\eps}[2]{\varepsilon_{#1}^{#2}}
\newcommand{\Eps}[2]{\mathcal{E}_{#1}^{#2}}
\newcommand{\e}[2]{\epsilon_{#1}^{#2}}
\newcommand{\kin}[2]{t_\text{#1}^{#2}}
\newcommand{\E}[2]{E_{#1}^{#2}}
\newcommand{\bE}[2]{\overline{E}_{#1}^{#2}}
\newcommand{\tE}[2]{\widetilde{E}_{#1}^{#2}}
\newcommand{\be}[2]{\overline{\epsilon}_{#1}^{#2}}
\newcommand{\n}[2]{n_{#1}^{#2}}
\newcommand{\Cx}[1]{C_\text{x}^{#1}}
% energies
\newcommand{\EHF}{E_\text{HF}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\HF}{\text{HF}}
\newcommand{\LDA}{\text{LDA}}
\newcommand{\SD}{\text{S}}
\newcommand{\VWN}{\text{VWN5}}
\newcommand{\eVWN}{\text{eVWN5}}
\newcommand{\SVWN}{\text{SVWN5}}
\newcommand{\LIM}{\text{LIM}}
\newcommand{\MOM}{\text{MOM}}
\newcommand{\Hxc}{\text{Hxc}}
\newcommand{\Ha}{\text{H}}
\newcommand{\ex}{\text{x}}
\newcommand{\co}{\text{c}}
\newcommand{\xc}{\text{xc}}
% matrices
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bw}{\boldsymbol{w}}
\newcommand{\bHc}{\boldsymbol{h}}
\newcommand{\Ex}[2]{\Omega_{#1}^{#2}}
\newcommand{\tEx}[2]{\widetilde{\Omega}_{#1}^{#2}}
% elements
\newcommand{\ew}[1]{w_{#1}}
\newcommand{\eHc}[1]{h_{#1}}
\newcommand{\eJ}[1]{J_{#1}}
\newcommand{\eK}[1]{K_{#1}}
\newcommand{\eF}[2]{F_{#1}^{#2}}
\newcommand{\ON}[2]{f_{#1}^{#2}}
\newcommand{\Det}[2]{\Phi_{#1}^{#2}}
% Numbers
\newcommand{\nEns}{M}
\newcommand{\nEl}{N}
\newcommand{\nOrb}{K}
% Ao and MO basis
\newcommand{\MO}[2]{\phi_{#1}^{#2}}
\newcommand{\cMO}[2]{c_{#1}^{#2}}
\newcommand{\AO}[1]{\chi_{#1}}
\newcommand{\RHH}{R_{\ce{H-H}}}
% units
\newcommand{\IneV}[1]{#1 eV}
\newcommand{\InAU}[1]{#1 a.u.}
\newcommand{\InAA}[1]{#1 \AA}
\newcommand{\kcal}{kcal/mol}
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques, Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LCQ}{Laboratoire de Chimie Quantique, Institut de Chimie, CNRS, Universit\'e de Strasbourg, Strasbourg, France}
\newcommand{\UL}{Instituut-Lorentz, Universiteit Leiden, P.O.~Box 9506, 2300 RA Leiden, The Netherlands}
\newcommand{\VU}{Division of Theoretical Chemistry, Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands}
\begin{document}
\title{Weight Dependence of Local Exchange-Correlation Functionals in Ensemble Density-Functional Theory: Double Excitations in Two-Electron Systems}
\author{Clotilde \surname{Marut}}
\affiliation{\LCPQ}
\author{Bruno \surname{Senjean}}
\affiliation{\UL}
\affiliation{\VU}
\author{Emmanuel \surname{Fromager}}
\affiliation{\LCQ}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\begin{abstract}
Gross--Oliveira--Kohn (GOK) ensemble density-functional theory (GOK-DFT) is a time-\textit{independent} formalism which allows to compute excitation energies via the derivative of the ensemble energy with respect to the weight of each excited state.
Contrary to the time-dependent version of density-functional theory (TD-DFT), double excitations can be easily computed within GOK-DFT.
However, to take full advantage of this formalism, one must have access to a \textit{weight-dependent} exchange-correlation functional in order to model the infamous ensemble derivative contribution to the excitation energies.
In the present article, we discuss the construction of first-rung (\ie, local) weight-dependent exchange-correlation density-functional approximations for two-electron atomic and molecular systems (He and H$_2$) specifically designed for the computation of double excitations within GOK-DFT.
A specific protocol is proposed to obtain accurate energies associated with double excitations.
\end{abstract}
\maketitle
%%%%%%%%%%%%%%%%%%%%
%%% INTRODUCTION %%%
%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Time-dependent density-functional theory (TD-DFT) has been the dominant force in the calculation of excitation energies of molecular systems in the last two decades.\cite{Casida,Ulrich_2012,Loos_2020a}
At a relatively low computational cost (at least compared to the other excited-state \textit{ab initio} methods), TD-DFT can provide accurate transition energies for low-lying excited states of organic molecules (see, for example, Ref.~\onlinecite{Dreuw_2005} and references therein).
Importantly, setting up a TD-DFT calculation for a given system is an almost pain-free process from a user perspective as the only (yet essential) input variable is the choice of the so-called exchange-correlation (xc) functional.
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) formalism 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 which have both the exact same one-electron density.
However, TD-DFT is far from being perfect as, in practice, drastic approximations must be made for the xc functional.
One of its issues 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.
Moreover, because it is so popular, it has been studied in excruciated details by the community, and some researchers have quickly unveiled various theoretical and practical deficiencies of approximate TD-DFT.
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.
One key consequence of this so-called adiabatic approximation (based on the assumption that the density varies slowly with time) is that double excitations are completely absent from the TD-DFT spectra. \cite{Levine_2006,Tozer_2000,Elliott_2011}
Although these double excitations are usually experimentally dark (which means they usually cannot be observed in photo-absorption spectroscopy), these states play, indirectly, a key role in many photochemistry mechanisms. \cite{Boggio-Pasqua_2007} They are, moreover, a real challenge for high-level computational methods. \cite{Loos_2018,Loos_2019,Loos_2020b}
One possible solution to access double excitations within TD-DFT is provided by spin-flip TD-DFT which describes double excitations as single excitations from the lowest triplet state. \cite{Huix-Rotllant_2010,Krylov_2001,Shao_2003,Wang_2004,Wang_2006,Minezawa_2009}
However, spin contamination might be an issue. \cite{Huix-Rotllant_2010}
In order to go beyond the adiabatic approximation, a dressed TD-DFT approach has been proposed by Maitra and coworkers \cite{Maitra_2004,Cave_2004} (see also Refs.~\onlinecite{Mazur_2009,Mazur_2011,Huix-Rotllant_2011,Elliott_2011,Maitra_2012}).
In this approach the xc kernel is made frequency dependent, which allows to treat doubly-excited states. \cite{Romaniello_2009a,Sangalli_2011,Loos_2019}
Maybe surprisingly, another possible way of accessing double excitations is to resort to a time-\textit{independent} formalism. \cite{Yang_2017,Sagredo_2018,Deur_2019}
With a computational cost similar to traditional KS-DFT, DFT for ensembles (eDFT) \cite{Theophilou_1979,Gross_1988a,Gross_1988b,Oliveira_1988} is a viable alternative following such a strategy currently under active development.\cite{Gidopoulos_2002,Franck_2014,Borgoo_2015,Kazaryan_2008,Gould_2013,Gould_2014,Filatov_2015,Filatov_2015b,Filatov_2015c,Gould_2017,Deur_2017,Gould_2018,Gould_2019,Sagredo_2018,Ayers_2018,Deur_2018,Deur_2019,Kraisler_2013,Kraisler_2014,Alam_2016,Alam_2017,Nagy_1998,Nagy_2001,Nagy_2005,Pastorczak_2013,Pastorczak_2014,Pribram-Jones_2014,Yang_2013a,Yang_2014,Yang_2017,Senjean_2015,Senjean_2016,Smith_2016,Senjean_2018}
In the assumption of monotonically decreasing weights, eDFT for excited states has the undeniable advantage to be based on a rigorous variational principle for ground and excited states, the so-called Gross--Oliveria--Kohn (GOK) variational principle. \cite{Gross_1988a}
In short, GOK-DFT (\ie, eDFT for neutral excitations) is the density-based analog of state-averaged wave function methods, and excitation energies can then be easily extracted from the total ensemble energy. \cite{Deur_2019}
Although the formal foundations of GOK-DFT have been set three decades ago, \cite{Gross_1988a,Gross_1988b,Oliveira_1988} its practical developments have been rather slow.
We believe that it is partly due to the lack of accurate approximations for GOK-DFT.
In particular, to the best of our knowledge, an explicitly weight-dependent density-functional approximation for ensemble (eDFA) has never been developed for atoms and molecules.
The present contribution is a small step towards this goal.
When one talks about constructing functionals, the local-density approximation (LDA) is never far away.
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 LDA functionals. \cite{Loos_2014a,Loos_2014b,Loos_2017a}
Electrons restricted to remain on the surface of a $\cD$-sphere (where $\cD$ is the dimensionality of the surface of the sphere) are an example of finite UEGs (FUEGs). \cite{Loos_2011b}
In particular, we combine these FUEGs with the usual infinite UEG (IUEG) to construct a weigh-dependent LDA correlation functional for ensembles, which is specifically designed to compute double excitations within GOK-DFT, and automatically incorporates the infamous 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.
Section \ref{sec:compdet} provides the computational details.
The results of our calculations for two-electron systems are reported and discussed in Sec.~\ref{sec:res}.
Finally, we draw our conclusions in Sec.~\ref{sec:ccl}.
Unless otherwise stated, atomic units are used throughout.
%%%%%%%%%%%%%%%%%%%%
%%% THEORY %%%
%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{sec:theo}
Let us consider a GOK ensemble of $\nEns$ electronic states with individual energies $\E{}{(0)} \le \ldots \le \E{}{(\nEns-1)}$, and (normalised) monotonically decreasing weights $\bw = (\ew{0},\ldots,\ew{M-1})$, \ie, $\sum_{I=0}^{\nEns-1} \ew{I} = 1$, and $\ew{0} \ge \ldots \ge \ew{\nEns-1}$.
The corresponding ensemble energy
\begin{equation}
\E{}{\bw} = \sum_{I=0}^{\nEns-1} \ew{I} \E{}{(I)}
\end{equation}
fulfils the variational principle
as follows\cite{Gross_1988a}
\begin{eqnarray}\label{eq:ens_energy}
\E{}{\bw} = \min_{\hGam{\bw}} \Tr[\hGam{\bw} \hH],
\end{eqnarray}
where $\hH = \hT + \hWee + \hVne$ contains the kinetic, electron-electron and nuclei-electron interaction potential operators, respectively, $\Tr$ denotes the trace and $\hGam{\bw}$ is a trial density matrix of the form
\begin{eqnarray}
\hGam{\bw} = \sum_{I=0}^{\nEns - 1} \ew{I} \dyad*{\overline{\Psi}^{(I)}},
\end{eqnarray}
where $\lbrace \overline{\Psi}^{(I)} \rbrace_{0 \le I \le \nEns-1}$ is a set of $\nEns$ orthonormal trial wave functions.
The lower bound of Eq.~\eqref{eq:ens_energy} is reached when the set of wave functions correspond to the exact eigenstates of $\hH$, \ie, $\lbrace \overline{\Psi}^{(I)} \rbrace_{0 \le I \le \nEns-1} = \lbrace \Psi^{(I)} \rbrace_{0 \le I \le \nEns-1}$.
Multiplet degeneracies can be easily handled by assigning the same weight to the degenerate states.
One of the key feature of the GOK ensemble is that individual excitation energies are extracted from the ensemble energy via differentiation with respect to individual weights:
\begin{equation}
\pdv{\E{}{\bw}}{\ew{I}} = \E{}{(I)} - \E{}{(0)} = \Ex{}{(I)}.
\end{equation}
Turning to GOK-DFT, the extension of the Hohenberg--Kohn theorem to ensembles allows to rewrite the exact variational expression for the ensemble energy as\cite{Gross_1988a}
\begin{equation}
\label{eq:Ew-GOK}
\E{}{\bw} = \min_{\n{}{}} \qty{ \F{}{\bw}[\n{}{}] + \int \vne(\br{}) \n{}{}(\br{}) d\br{} },
\end{equation}
where $\vne(\br{})$ is the external potential and $\F{}{\bw}[\n{}{}]$ is the universal ensemble functional
(the weight-dependent analog of the Hohenberg--Kohn universal functional for ensembles).
In the KS formulation, this functional is decomposed as
\begin{equation}
\F{}{\bw}[\n{}{}]
= \Ts{\bw}[\n{}{}] + \E{\Hxc}{\bw}[\n{}{}]
= \Tr[ \hgam{\bw} \hT ] + \Tr[ \hgam{\bw} \hWee ],
\end{equation}
where
$\Ts{\bw}[\n{}{}]$ is the noninteracting ensemble kinetic energy functional,
\begin{equation}
\hgam{\bw} = \sum_{I=0}^{\nEns-1} \ew{I} \dyad{\Det{I}{\bw}}
\end{equation}
is the density matrix operator, $\lbrace \Det{I}{\bw} \rbrace_{0 \le I \le \nEns-1}$ are single-determinant wave functions (or configuration state functions) built with KS orbitals $\lbrace \MO{p}{\bw}(\br{}) \rbrace_{1 \le p \le \nOrb}$, and
\begin{equation}
\label{eq:exc_def}
\begin{split}
\E{\Hxc}{\bw}[\n{}{}]
& = \E{\Ha}{}[\n{}{}] + \E{\xc}{\bw}[\n{}{}]
\\
& = \frac{1}{2} \iint \frac{\n{}{}(\br{}) \n{}{}(\br{}')}{\abs{\br{}-\br{}'}} d\br{} d\br{}'
+ \int \e{\xc}{\bw}(\n{}{}(\br{})) \n{}{}(\br{}) d\br{}
\end{split}
\end{equation}
is the ensemble Hartree-exchange-correlation (Hxc) functional.
Note that the weight-independent Hartree functional $\E{\Ha}{}[\n{}{}]$ causes the infamous ghost-interaction error \cite{Gidopoulos_2002, Pastorczak_2014, Alam_2016, Alam_2017, Gould_2017} in eDFT, which is supposed to be cancelled by the weight-dependent xc functional $\E{\xc}{\bw}[\n{}{}]$.
From the GOK-DFT ensemble energy expression in Eq.~\eqref{eq:Ew-GOK}, we obtain \cite{Gross_1988b,Deur_2019}
\begin{equation}
\begin{split}
\label{eq:dEdw}
\pdv{\E{}{\bw}}{\ew{I}}
& = \E{}{(I)} - \E{}{(0)}
\\
& = \Eps{I}{\bw} - \Eps{0}{\bw} + \left. \pdv{\E{\xc}{\bw}[\n{}{}]}{\ew{I}} \right|_{\n{}{} = \n{}{\bw}(\br{})},
\end{split}
\end{equation}
where
\begin{align}
\label{eq:nw}
\n{}{\bw}(\br{}) & = \sum_{I=0}^{\nEns-1} \ew{I} \n{\Det{I}{\bw}}{}(\br{}),
\\
\label{eq:nI}
\n{\Det{I}{\bw}}{}(\br{}) & = \sum_{p}^{\nOrb} \ON{p}{(I)} [\MO{p}{\bw}(\br{})]^2
\end{align}
are the ensemble and individual one-electron densities, respectively,
\begin{equation}
\label{eq:KS-energy}
\Eps{I}{\bw} = \sum_{p}^{\nOrb} \ON{p}{(I)} \eps{p}{\bw}
\end{equation}
is the weight-dependent KS energy of state $I$, and $\eps{p}{\bw}$ is the KS orbital energy associated with $\MO{p}{\bw}(\br{})$ [$\ON{p}{(I)}$ being its occupancy for the state $I$].
The latters are determined by solving the ensemble KS equation
\begin{equation}
\label{eq:eKS}
\qty{ \hHc(\br{}) + \fdv{\E{\Hxc}{\bw}[\n{}{\bw}]}{\n{}{}(\br{})}} \MO{p}{\bw}(\br{}) = \eps{p}{\bw} \MO{p}{\bw}(\br{}),
\end{equation}
where $\hHc(\br{}) = -\nabla^2/2 + \vne(\br{})$, and
\begin{equation}
\fdv{\E{\Hxc}{\bw}[\n{}{}]}{\n{}{}(\br{})}
= \fdv{\E{\Ha}{}[\n{}{}]}{\n{}{}(\br{})} + \fdv{\E{\xc}{\bw}[\n{}{}]}{\n{}{}(\br{})}
\end{equation}
is the Hxc potential, with
\begin{subequations}
\begin{align}
\fdv{\E{\Ha}{}[\n{}{}]}{\n{}{}(\br{})}
& = \frac{1}{2} \int \frac{\n{}{}(\br{}')}{\abs{\br{}-\br{}'}} d\br{}',
\\
\fdv{\E{\xc}{\bw}[\n{}{}]}{\n{}{}(\br{})}
& = \left. \pdv{\e{\xc}{\bw{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{}(\br{})} \n{}{}(\br{}) + \e{\xc}{\bw{}}(\n{}{}(\br{})).
\end{align}
\end{subequations}
Equation \eqref{eq:dEdw} is our working equation for computing excitation energies from a practical point of view.
Note that the individual densities $\n{\Det{I}{\bw}}{}(\br{})$ defined in Eq.~\eqref{eq:nI} do not match the \textit{exact} individual-state densities as the non-interacting KS ensemble is expected to reproduce the true interacting ensemble density $\n{}{\bw}(\br{})$ defined in Eq.~\eqref{eq:nw}, and not each individual density.
Nevertheless, these densities can still be extracted in principle exactly from the KS ensemble as shown by Fromager. \cite{Fromager_2020}
In the following, we adopt the usual decomposition, and write down the weight-dependent xc functional as
\begin{equation}
\e{\xc}{\ew{}}(\n{}{}) = \e{\ex}{\ew{}}(\n{}{}) + \e{\co}{\ew{}}(\n{}{}),
\end{equation}
where $\e{\ex}{\ew{}}(\n{}{})$ and $\e{\co}{\ew{}}(\n{}{})$ are the weight-dependent exchange and correlation functionals, respectively.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% COMPUTATIONAL DETAILS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:compdet}
The self-consistent GOK-DFT calculations have been performed in a restricted formalism with the \texttt{QuAcK} software, \cite{QuAcK} which is freely available on \texttt{github}, and where the present weight-dependent functionals have been implemented.
For more details about the self-consistent implementation of GOK-DFT, we refer the interested reader to Ref.~\onlinecite{Loos_2020} where additional technical details can be found.
For all calculations, we use the aug-cc-pVXZ (X = D, T, Q, and 5) Dunning family of atomic basis sets. \cite{Dunning_1989,Kendall_1992,Woon_1994}
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.
Although one should have $0 \le \ew{} \le 1/2$ to ensure the GOK variational principle, 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}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}
\label{sec:res}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hydrogen molecule at equilibrium}
\label{sec:H2}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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 Slater-Dirac (LDA) local exchange functional, \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}
The ensemble energy $\E{}{\ew{}}$ is depicted in Fig.~\ref{fig:Ew_H2} as a function of the weight $0 \le \ew{} \le 1$.
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$.
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{}$.
\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.
\label{fig:Ew_H2}
}
\end{figure}
\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.
\label{fig:Om_H2}
}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Weight-dependent exchange functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Second, in order to remove this spurious curvature of the ensemble energy (which is mostly due to the ghost-interaction error, 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$.
Doing so, we have found that the present weight-dependent exchange functional (denoted as GIC-S in the following as its main purpose is to correct for the ghost-interaction error)
\begin{equation}
\e{\ex}{\ew{},\text{GIC-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 ],
\end{equation}
and
\begin{subequations}
\begin{align}
\alpha & = + 0.575\,178,
&
\beta & = - 0.021\,108,
&
\gamma & = - 0.367\,189,
\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$.
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.
Finally, let us mention that, around $\ew{} = 0$, the behaviour of Eq.~\eqref{eq:Cxw} is linear.
We shall come back to this point later on.
\begin{figure}
\includegraphics[width=\linewidth]{Cxw}
\caption{
$\Cx{\ew{}}/\Cx{\ew{}=0}$ as a function of $\ew{}$ [see Eq.~\eqref{eq:Cxw}] computed with the aug-cc-pVTZ basis set for the \ce{He} atom (blue) and the \ce{H2} molecule at $\RHH = 1.4$ bohr (red) and $\RHH = 3.7$ bohr (green).
\label{fig:Cxw}
}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Weight-independent correlation functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Third, we add up correlation effects via the VWN5 local correlation functional. \cite{Vosko_1980}
For the sake of clarity, the explicit expression of the VWN5 functional is not reported here but it can be found in Ref.~\onlinecite{Vosko_1980}.
The combination of the Slater and VWN5 functionals (SVWN5) yield a highly convex ensemble energy (as shown in Fig.~\ref{fig:Ew_H2}), while the combination of GIC-S and VWN5 (GIC-SVWN5) exhibit a small curvature (the ensemble energy is slightly concave) and improved excitation energies, especially at small weights, where the GIC-SVWN5 excitation energy is almost spot on.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Weight-dependent correlation functional}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fourth, in the spirit of our recent work, \cite{Loos_2020} we have designed a weight-dependent correlation functional.
To build this weight-dependent correlation functional, we consider the singlet ground state and the first singlet doubly-excited state of a two-electron FUEGs which consists of two electrons confined to the surface of a 3-sphere (also known as a glome).\cite{Loos_2009a,Loos_2009c,Loos_2010e}
Notably, these two states have the same (uniform) density $\n{}{} = 2/(2\pi^2 R^3)$, where $R$ is the radius of the 3-sphere onto which the electrons are confined.
Indeed, the orbitals for an electron on a 3-sphere of unit radius are the normalised hyperspherical harmonics $Y_{\ell\mu}$, where $\ell$ is the principal quantum number and $\mu$ is a composite index of the remaining two quantum numbers. \cite{AveryBook, Avery_1993}
As mentioned above, we confine our attention to paramagnetic (or unpolarised) systems, and in particular to the simple two-electron system in which the orbital with $\ell = 0$ is doubly-occupied by one spin-up and one spin-down electron, thus yielding an electron density that is uniform on the 3-sphere.
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
\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}.
\label{eq:eHF_1}
\end{align}
\end{subequations}
Thanks to highly-accurate calculations and the expressions of the HF energies provided by Eqs.~\eqref{eq:eHF_0} and \eqref{eq:eHF_1}, \cite{Loos_2009a,Loos_2009c,Loos_2010e} one can write down, for each state, an accurate analytical expression of the reduced correlation energy \cite{Loos_2013a, Loos_2014a} via the following Pad\'e approximant \cite{Sun_2016,Loos_2020}
\begin{equation}
\label{eq:ec}
\e{\co}{(I)}(\n{}{}) = \frac{a_1^{(I)}}{1 + a_2^{(I)} \n{}{-1/6} + a_3^{(I)} \n{}{-1/3}},
\end{equation}
where $a_2^{(I)}$ and $a_3^{(I)}$ are state-specific fitting parameters, which are provided in Table \ref{tab:OG_func}.
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:
\begin{equation}
\label{eq:ecw}
\e{\co}{\ew{}}(\n{}{}) = (1-\ew{}) \e{\co}{(0)}(\n{}{}) + \ew{} \e{\co}{(1)}(\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.
The data gathered in Table \ref{tab:Ref} are also reported.
}
\label{fig:Ec}
\end{figure}
%%% %%% %%%
%%% TABLE I %%%
\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.
}
\begin{ruledtabular}
\begin{tabular}{lcc}
& \tabc{Ground state} & \tabc{Doubly-excited state} \\
$R$ & \tabc{$I=0$} & \tabc{$I=1$} \\
\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$ \\
\end{tabular}
\end{ruledtabular}
\end{table}
%%% TABLE II %%%
\begin{table}
\caption{
\label{tab:OG_func}
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$} \\
\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$ \\
\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).
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{}{})
\end{equation}
via the following global, state-independent shift:
\begin{equation}
\be{\co}{(I)}(\n{}{}) = \e{\co}{(I)}(\n{}{}) + \e{\co}{\VWN}(\n{}{}) - \e{\co}{(0)}(\n{}{}).
\end{equation}
In the following, we name this weight-dependent correlation functional ``eVWN5'' as it is a natural extension of the VWN5 local correlation functional for ensembles.
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{}{})]
\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{}{})$.
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),
\end{equation}
showing that the weight correction is purely linear in eVWN5 and entirely dependent on the FUEG model.
As shown in Fig.~\ref{fig:Ew_H2}, the GIC-SeVWN5 is slightly less concave than its GIC-SVWN5 counterpart and it also improves (not by much) the excitation energy (see Fig.~\ref{fig:Om_H2}).
For a more qualitative picture, Table \ref{tab:BigTab_H2} reports excitation energies for various methods and basis sets.
In particular, we report the excitation energies obtained with GOK-DFT in the zero-weight limit (\ie, $\ew{} = 0$) and for the equi-ensemble (\ie, $\ew{} = 1/2$).
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}
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 GIC-SeVWN5 functional at zero weight are the most accurate with an improvement of $0.25$ eV as compared to GIC-SVWN5, which is due to the ensemble derivative contribution of the eVWN5 functional.
The GIC-SeVWN5 excitation energies at equi-weights (\ie, $\ew{} = 1/2$) are less satisfactory, but still remains in good agreement with FCI, with again a small improvement as compared to GIC-SVWN5.
It is also important to mention that the GIC-S functional does not alter the MOM excitation energy as the correction vanishes accordingly for $\ew{} = 1$ (\textit{vide supra}).
Note that by construction,
for ensemble energies that are quadratic with respect to the weight (which is almost always the case in this paper),
LIM and MOM can be reduced to a single calculation
at $w = 1/4$ and $w=1/2$, respectively, instead of performing an interpolation between two different calculations.
Finally, although we had to design a system-specific, weight-dependent exchange functional to reach such accuracy, we have not used any high-level reference data (such as FCI) to tune our functional, the only requirement being the linearity of the ensemble energy (obtained with LDA exchange) between $\ew{} = 0$ and $1$.
%%% TABLE III %%%
\begin{table}
\caption{
Excitation energies (in eV) associated with the lowest double excitation of \ce{H2} with $\RHH = 1.4$ bohr for various methods, combinations of xc functionals, and basis sets.
\label{tab:BigTab_H2}
}
\begin{ruledtabular}
\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 \\
\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 & VWN5 & aug-cc-pVDZ & 37.83 & 31.19 & 35.66 & 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 \\
& & aug-cc-pVTZ & 37.61 & 37.04 & 37.28 & 29.34 \\
& & aug-cc-pVQZ & 37.32 & 37.14 & 37.27 & 29.34 \\
\\
S & & aug-cc-pVDZ & 19.44 & 27.35 & 23.54 & 26.60 \\
& & aug-cc-pVTZ & 19.47 & 27.42 & 23.62 & 26.67 \\
& & aug-cc-pVQZ & 19.41 & 27.42 & 23.62 & 26.67 \\
\\
S & VWN5 & aug-cc-pVDZ & 21.04 & 27.76 & 24.40 & 27.10 \\
& & aug-cc-pVTZ & 21.14 & 27.81 & 24.46 & 27.17 \\
& & aug-cc-pVQZ & 21.13 & 27.81 & 24.46 & 27.17 \\
\\
S & eVWN5 & aug-cc-pVDZ & 21.28 & 27.92 & 24.49 & 27.27 \\
& & aug-cc-pVTZ & 21.39 & 27.98 & 24.55 & 27.34 \\
& & aug-cc-pVQZ & 21.38 & 27.97 & 24.55 & 27.34 \\
\\
GIC-S & & aug-cc-pVDZ & 26.83 & 26.51 & 26.53 & 26.60 \\
& & aug-cc-pVTZ & 26.88 & 26.59 & 26.61 & 26.67 \\
& & aug-cc-pVQZ & 26.82 & 26.60 & 26.62 & 26.67 \\
\\
GIC-S & VWN5 & aug-cc-pVDZ & 28.54 & 26.94 & 27.48 & 27.10 \\
& & aug-cc-pVTZ & 28.66 & 27.00 & 27.56 & 27.17 \\
& & aug-cc-pVQZ & 28.64 & 27.00 & 27.56 & 27.17 \\
\\
GIC-S & eVWN5 & aug-cc-pVDZ & 28.78 & 27.10 & 27.56 & 27.27 \\
& & aug-cc-pVTZ & 28.90 & 27.16 & 27.64 & 27.34 \\
& & aug-cc-pVQZ & 28.89 & 27.16 & 27.65 & 27.34 \\
\hline
B & LYP & aug-mcc-pV8Z & & & & 28.42 \\
B3 & LYP & aug-mcc-pV8Z & & & & 27.77 \\
HF & LYP & aug-mcc-pV8Z & & & & 29.18 \\
HF & & aug-mcc-pV8Z & & & & 28.65 \\
\hline
\mc{6}{l}{Accurate\fnm[1]} & 28.75 \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{FCI/aug-mcc-pV8Z calculation from Ref.~\onlinecite{Barca_2018a}.}
\end{table}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hydrogen molecule at stretched geometry}
\label{sec:H2st}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
To investigate the weight dependence of the xc functional in the strong correlation regime, we now consider the \ce{H2} molecule in a stretched geometry ($\RHH = 3.7$ bohr).
For this particular geometry, the doubly-excited state becomes the lowest excited state.
We then follow the same protocol as in Sec.~\ref{sec:H2}, and considering again the aug-cc-pVTZ basis set, we design a GIC-S functional for this system at $\RHH = 3.7$ bohr.
It yields $\alpha = +0.019\,226$, $\beta = -0.017\,996$, and $\gamma = -0.022\,945$ [see Eq.~\eqref{eq:Cxw}].
The weight dependence of $\Cx{\ew{}}$ is illustrated in Fig.~\ref{fig:Cxw} (green curve).
One clearly sees that the correction brought by GIC-S is much more gentle than at $\RHH = 1.4$ bohr, which means that the ensemble energy obtained with the LDA exchange functional is much more linear at $\RHH = 3.7$ bohr.
In other words, the ghost-interaction ``hole'' depicted in Fig.~\ref{fig:Cxw} is thus much more shallow at stretched geometry.
Note that this linearity at $\RHH = 3.7$ bohr was also observed using weight-independent xc functionals in Ref.~\onlinecite{Senjean_2015}.
Table \ref{tab:BigTab_H2st} reports, for the aug-cc-pVTZ basis set (which delivers basis set converged results), the same set of calculations as in Table \ref{tab:BigTab_H2}.
As a reference value, we computed a FCI/aug-cc-pV5Z excitation energy of $8.69$ eV, which compares well with previous studies. \cite{Senjean_2015}
For $\RHH = 3.7$ bohr, it is much harder to get an accurate estimate of the excitation energy, the closest match being reached with HF exchange and eVWN5 correlation at equi-ensemble.
%\bruno{? I don't see it, for me HF is really bad here, especially due to its very hight dependence on the weight ! Maybe you are just referring to MOM ?}.
As expected from the linearity of the ensemble energy, the GIC-S functional coupled or not with a correlation functional yield extremely stable excitation energies as a function of the weight, with only a few tenths of eV difference between the zero- and equi-weights limits.
Nonetheless, the excitation energy is still off by $3$ eV.
The fundamental theoretical reason of such a poor agreement is not clear.
The fact that HF exchange yields better excitation energies hints at the effect of self-interaction error.
%\bruno{I'm a bit surprise that the ensemble correction to the correlation functional does not improve things at all... Is the derivative discontinuity, computed with this functional, almost 0 here ?}
%%% TABLE IV %%%
\begin{table}
\caption{
Excitation energies (in eV) associated with the lowest double excitation of \ce{H2} at $\RHH = 3.7$ bohr obtained with the aug-cc-pVTZ basis set for various methods and combinations of xc functionals.
\label{tab:BigTab_H2st}
}
\begin{ruledtabular}
\begin{tabular}{llcccc}
\mc{2}{c}{xc functional} & \mc{2}{c}{GOK} \\
\cline{1-2} \cline{3-4}
\tabc{x} & \tabc{c} & $\ew{} = 0$ & $\ew{} = 1/2$ & LIM & MOM \\
\hline
HF & & 19.09 & 6.59 & 12.92 & 6.52 \\
HF & VWN5 & 19.40 & 6.54 & 13.02 & 6.49 \\
HF & eVWN5 & 19.59 & 6.72 & 13.11 & \fnm[1] \\
S & & 5.31 & 5.60 & 5.46 & 5.56 \\
S & VWN5 & 5.34 & 5.57 & 5.46 & 5.52 \\
S & eVWN5 & 5.53 & 5.76 & 5.56 & 5.72 \\
GIC-S & & 5.55 & 5.56 & 5.56 & 5.56 \\
GIC-S & VWN5 & 5.58 & 5.53 & 5.57 & 5.52 \\
GIC-S & eVWN5 & 5.77 & 5.72 & 5.66 & 5.72 \\
\hline
B & LYP & & & & 5.28 \\
B3 & LYP & & & & 5.55 \\
HF & LYP & & & & 6.68 \\
\hline
\mc{5}{l}{Accurate\fnm[2]} & 8.69 \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{KS calculation does not converge.}
\fnt[2]{FCI/aug-cc-pV5Z calculation performed with QUANTUM PACKAGE. \cite{QP2}}
\end{table}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Helium atom}
\label{sec:He}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
As a final example, we consider the \ce{He} atom which can be seen as the limiting form of the \ce{H2} molecule for very short bond lengths.
In \ce{He}, the lowest doubly-excited state is an auto-ionising resonance state, extremely high in energy and lies in the continuum. \cite{Madden_1963}
In Ref.~\onlinecite{Burges_1995}, highly-accurate calculations estimate an excitation energy of $2.126$ hartree for this $1s^2 \rightarrow 2s^2$ transition.
Nonetheless, it can be nicely described with a Gaussian basis set containing enough diffuse functions.
Consequently, we consider for this particular example the d-aug-cc-pVQZ basis set which contains two sets of diffuse functions.
The excitation energies associated with this double excitation computed with various methods and combinations of xc functionals are gathered in Table \ref{tab:BigTab_He}.
The parameters of the GIC-S weight-dependent exchange functional (computed with the smaller aug-cc-pVTZ basis) are $\alpha = +1.912\,574$, $\beta = +2.715\,267$, and $\gamma = +2.163\,422$ [see Eq.~\eqref{eq:Cxw}], the curvature of the ensemble energy being more pronounced in \ce{He} than in \ce{H2} (blue curve in Fig.~\ref{fig:Cxw}).
In other words, the ghost-interaction hole is deeper.
The results reported in Table \ref{tab:BigTab_He} evidence this strong weight dependence of the excitation energies for HF or LDA exchange.
The GIC-S exchange functional attenuates significantly this dependence, and when coupled with the eVWN5 weight-dependent correlation functional, the GIC-SeVWN5 excitation energy for $\ew{} = 0$ is only $8$ millihartree (or $0.22$ eV) off the reference value.
%\bruno{But also with GIC-SVWN5, as in the rest of this article, so one could wonder about the usefulness of the eVWN5 functional...}
As in the case of \ce{H2}, the excitation energies obtained at zero-weight are more accurate than at equi-weight, while the opposite conclusion were made in Ref.~\onlinecite{Loos_2020}.
This motivates further the importance of developing weight-dependent functionals that yields linear ensemble energies in order to get rid of the weight-dependency of the excitation energy.
As a final comment, let us stress again that the present protocol does not rely on high-level calculations as the sole requirement for constructing the GIC-S functional is the linearity of the ensemble energy with respect to the weight of the double excitation.
%%% TABLE V %%%
\begin{table}
\caption{
Excitation energies (in hartree) associated with the lowest double excitation of \ce{He} obtained with the d-aug-cc-pVQZ basis set for various methods and combinations of xc functionals.
\label{tab:BigTab_He}
}
\begin{ruledtabular}
\begin{tabular}{llcccc}
\mc{2}{c}{xc functional} & \mc{2}{c}{GOK} \\
\cline{1-2} \cline{3-4}
\tabc{x} & \tabc{c} & $\ew{} = 0$ & $\ew{} = 1/2$ & LIM & MOM \\
\hline
HF & & 1.874 & 2.212 & 2.080 & 2.142 \\
HF & VWN5 & 1.988 & 2.260 & 2.153 & 2.193 \\
HF & eVWN5 & 2.000 & 2.264 & 2.156 & 2.196 \\
S & & 1.062 & 2.056 & 1.547 & 2.030 \\
S & VWN5 & 1.163 & 2.104 & 1.612 & 2.079 \\
S & eVWN5 & 1.174 & 2.108 & 1.615 & 2.083 \\
GIC-S & & 1.996 & 2.044 & 1.988 & 2.030 \\
GIC-S & VWN5 & 2.107 & 2.097 & 2.060 & 2.079 \\
GIC-S & eVWN5 & 2.118 & 2.100 & 2.063 & 2.083 \\
\hline
B & LYP & & & & 2.147 \\
B3 & LYP & & & & 2.150 \\
HF & LYP & & & & 2.171 \\
\hline
\mc{5}{l}{Accurate\fnm[1]} & 2.126 \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{Explicitly-correlated calculations from Ref.~\onlinecite{Burges_1995}.}
\end{table}
%%% TABLE I %%%
%\begin{table}
%\caption{
%Excitation energies (in eV) associated with the lowest double excitation of \ce{HNO} obtained with the aug-cc-pVDZ basis set for various methods and combinations of xc functionals.
%\label{tab:BigTab_H2st}
%}
%\begin{ruledtabular}
%\begin{tabular}{llcccc}
% \mc{2}{c}{xc functional} & \mc{2}{c}{GOK} \\
% \cline{1-2} \cline{3-4}
% \tabc{x} & \tabc{c} & $\ew{} = 0$ & $\ew{} = 1/2$ & LIM & MOM \\
% \hline
% HF & & & & & \\
% HF & VWN5 & & & & \\
% S & & 1.72 & 4.00 & 2.86 & 3.99 \\
% S & VWN5 & & & & \\
% GIC-S & & 3.99 & 3.99 & 3.99 & 3.99 \\
% GIC-S & VWN5 & 4.05 & 4.03 & 4.04 & 4.03 \\
% \hline
% S & PW92 & & & & 4.00\fnm[1] \\
% PBE & PBE & & & & 4.13\fnm[1] \\
% SCAN & SCAN & & & & 4.24\fnm[1] \\
% B97M-V & B97M-V & & & & 4.33\fnm[1] \\
% PBE0 & PBE0 & & & & 4.24\fnm[1] \\
% \hline
% \mc{5}{l}{Theoretical best estimate\fnm[2]} & 4.32 \\
%\end{tabular}
%\end{ruledtabular}
%\fnt[1]{Square gradient minimization (SGM) approach from Ref.~\onlinecite{Hait_2020} obtained with the aug-cc-pVTZ basis set. SGM is theoretically equivalent to MOM.}
%\fnt[2]{Theoretical best estimate from Ref.~\onlinecite{Loos_2019} obtained at the (extrapolated) FCI/aug-cc-pVQZ level.}
%\end{table}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%
%%% CONCLUSION %%%
%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:ccl}
In the light of the results obtained in this study on double excitations computed within the GOK-DFT framework, we believe that the development of more universal weight-dependent exchange and correlation functionals has a bright future, and we hope to be able to report further on this in the near future.
%%%%%%%%%%%%%%%%%%%%%%%%
%%% ACKNOWLEDGEMENTS %%%
%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
PFL thanks Radovan Bast and Anthony Scemama for technical assistance, as well as Julien Toulouse for stimulating discussions on double excitations.
CM thanks the \textit{Universit\'e Paul Sabatier} (Toulouse, France) for a PhD scholarship.
%PFL thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481) for financial support.
This work has also been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''.}
\end{acknowledgements}
%%%%%%%%%%%%%%%%%%%%
%%% BIBLIOGRAPHY %%%
%%%%%%%%%%%%%%%%%%%%
\bibliography{FarDFT}
\end{document}
Reference in New Issue View Git Blame Copy Permalink