making progress on the ensemble energies

Manuscript/FarDFT.tex

@@ -40,7 +40,7 @@
 \newcommand{\hH}{\Hat{H}}
 \newcommand{\hHc}{\Hat{h}}
 \newcommand{\hT}{\Hat{T}}
-\newcommand{\bH}{\bm{H}}
+\newcommand{\bH}{\boldsymbol{H}}
 \newcommand{\hVext}{\Hat{V}_\text{ext}}
 \newcommand{\vext}{v_\text{ext}}
 \newcommand{\hWee}{\Hat{W}_\text{ee}}
@@ -54,6 +54,7 @@
 \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}}
@@ -61,11 +62,6 @@
 % 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}}
 \newcommand{\HF}{\text{HF}}
 \newcommand{\LDA}{\text{LDA}}
 \newcommand{\eLDA}{\text{eLDA}}
@@ -77,14 +73,15 @@
 \newcommand{\xc}{\text{xc}}
 
 % matrices
-\newcommand{\br}{\bm{r}}
+\newcommand{\br}{\boldsymbol{r}}
-\newcommand{\bw}{\bm{w}}
+\newcommand{\bw}{\boldsymbol{w}}
-\newcommand{\bG}{\bm{G}}
+\newcommand{\bG}{\boldsymbol{G}}
-\newcommand{\bS}{\bm{S}}
+\newcommand{\bS}{\boldsymbol{S}}
-\newcommand{\bGamma}[1]{\bm{\Gamma}^{#1}}
+\newcommand{\bGamma}[1]{\boldsymbol{\Gamma}^{#1}}
-\newcommand{\bHc}{\bm{h}}
+\newcommand{\bHc}{\boldsymbol{h}}
-\newcommand{\bF}[1]{\bm{F}^{#1}}
+\newcommand{\bF}[1]{\boldsymbol{F}^{#1}}
 \newcommand{\Ex}[2]{\Omega_{#1}^{#2}}
+\newcommand{\tEx}[2]{\widetilde{\Omega}_{#1}^{#2}}
 
 % elements
 \newcommand{\ew}[1]{w_{#1}}
@@ -231,7 +228,8 @@ and
 	\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{}.
+	& = \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}
 are the noninteracting ensemble kinetic energy functional and ensemble Hartree-exchange-correlation (Hxc) functional, respectively.
@@ -242,14 +240,34 @@ From the GOK-DFT ensemble energy expression in Eq.~\eqref{eq:Ew-GOK}, we obtain
 	\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}},
+	= \Eps{I}{\bw} - \Eps{0}{\bw} + \left. \pdv{\E{\xc}{\bw}[\n{}{}]}{\ew{I}} \right|_{\n{}{} = \n{}{\bw}(\br{})},
 \end{equation}
-where $\Eps{I}{\bw} = \sum_{p}^{\Norb} \ON{p}{(I)} \eps{p}{\bw}$, $\eps{p}{\bw}$ is the $p$th KS orbital energy given by the ensemble KS equation
+where 
 \begin{equation}
+	\n{}{\bw}(\br{}) = \sum_{I=0}^{\Nens-1} \ew{I} \n{}{(I)}(\br{})
+\end{equation}
+is the ensemble density, 
+\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, 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$) given by the ensemble KS equation
+\begin{equation}
+\label{eq:eKS}
 	\qty( \hHc(\br{}) + \fdv{\E{\Hxc}{\bw}[\n{}{}]}{\n{}{}(\br{})}) \MO{p}{\bw}(\br{}) = \eps{p}{\bw} \MO{p}{\bw}(\br{}),
 \end{equation}
-where $\hHc(\br{}) = -\frac{\nabla^2}{2}  + \vext(\br{})$, $\MO{p}{\bw}(\br{})$ is a KS orbital, $\ON{p}{(I)}$ its occupancy for the state $I$, and $\n{}{\bw} = \sum_{I=0}^{\Nens-1} \ew{I} \n{}{(I)}$ is the ensemble density.
+where $\hHc(\br{}) = -\nabla^2/2  + \vext(\br{})$, and
-Equation \eqref{eq:dEdw} is our working equation for computing excitation energies.
+\begin{equation}
+\begin{split}
+	\fdv{\E{\Hxc}{\bw}[\n{}{}]}{\n{}{}(\br{})} 
+	& = \fdv{\E{\Ha}{\bw}[\n{}{}]}{\n{}{}(\br{})} + \fdv{\E{\xc}{\bw}[\n{}{}]}{\n{}{}(\br{})}
+	\\
+	& = \frac{1}{2} \int \frac{\n{}{}(\br{}')}{\abs{\br{}-\br{}'}} d\br{}' 
+	+ \left. \pdv{\e{\xc}{\bw{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{}(\br{})} \n{}{}(\br{}) + \e{\xc}{\bw{}}[\n{}{}(\br{})] 
+\end{split}
+\end{equation}
+is the Hxc potential.
+Equation \eqref{eq:dEdw} is our working equation for computing excitation energies from a practical point of view.
 
 %%%%%%%%%%%%%%%%%%
 %%% FUNCTIONAL %%%
@@ -324,12 +342,9 @@ and we then obtain
 We can now combine these two exchange functionals to create a weight-dependent exchange functional for a two-state ensemble
 \begin{equation}
 \label{eq:exw}
-\begin{split}
 	\e{\ex}{\ew{}}(\n{}{}) 
-	& = (1-\ew{}) \e{\ex}{(0)}(\n{}{}) + \ew{} \e{\ex}{(1)}(\n{}{})
+	= (1-\ew{}) \e{\ex}{(0)}(\n{}{}) + \ew{} \e{\ex}{(1)}(\n{}{})
-	\\
+	= \Cx{\ew{}} \n{}{1/3}
-	& = \Cx{\ew{}} \n{}{1/3}
-\end{split}
 \end{equation}
 with
 \begin{equation}
@@ -487,6 +502,9 @@ In the case of a homogeneous system (or equivalently within the LDA), substituti
 \label{sec:res}
 Here, we consider as testing ground the minimal-basis \ce{H2} molecule.
 We select STO-3G as minimal basis, and study the behaviour of the total energy of \ce{H2} as a function of the internuclear distance $\RHH$ (in bohr).
+This minimal-basis example is quite pedagogical as the molecular orbitals are fixed by symmetry.
+Therefore, there is no density-driven error and the only error that we are going to see is the functional-driven error (and this is what we want to study).
+
 The bonding and antibonding orbitals of the \ce{H2} molecule are given by
 \begin{subequations}
 \begin{align}
@@ -519,14 +537,14 @@ with
 \end{subequations}
 Note that, in the HF case, there is no self-interaction error as $\eJ{pp} = \eK{pp}$.
 We also define the HF excitation energy as $\Ex{\HF}{(1)} = \E{\HF}{(1)} - \E{\HF}{(0)}$.
-The HF orbital energies are
+%The HF orbital energies are
-\begin{subequations}
+%\begin{subequations}
-\begin{align}
+%\begin{align}
-	\eps{1}{\HF} & = \eHc{1} + 2\eJ{11} - \eK{11},
+%	\eps{1}{\HF} & = \eHc{1} + 2\eJ{11} - \eK{11},
-	\\
+%	\\
-	\eps{2}{\HF} & = \eHc{2} + 2\eJ{12} - \eK{12}.
+%	\eps{2}{\HF} & = \eHc{2} + 2\eJ{12} - \eK{12}.
-\end{align}
+%\end{align}
-\end{subequations}
+%\end{subequations}
 
 As reference results, we consider CID (configuration interaction with doubles) computed in the same (minimal) basis set.
 The CID energies of the ground state and doubly-excited states are provided by the eigenvalues of the following CID matrix:
@@ -568,14 +586,17 @@ with
 	\n{}{(1)}(\br{}) & = 2 \MO{2}{2}(\br{}),
 \end{align}
 Note that, contrary to the HF case, self-interaction is present in LDA.
-The KS orbital energies are given by
+%The KS orbital energies are given by
-\begin{subequations}
+%\begin{subequations}
-\begin{align}
+%\begin{align}
-	\eps{1}{\LDA} & = \eHc{1} + 2\eJ{11} + \ldots,
+%	\eps{1}{\LDA} 
-	\\
+%	& = \eHc{1} + 2\eJ{11} 
-	\eps{2}{\LDA} & = \eHc{2} + 2\eJ{12} + \ldots.
+%	+ \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}[\n{}{}]}{\n{}{}} \right|_{\n{}{} = \n{}{(0)}(\br{})} \n{}{(0)}(\br{}) d\br{},
-\end{align}
+%	\\
-\end{subequations}
+%	\eps{2}{\LDA} & = \eHc{2} + 2\eJ{12} 
+%	+ \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}[\n{}{}]}{\n{}{}} \right|_{\n{}{} = \n{}{(0)}(\br{})} \n{}{(1)}(\br{}) d\br{}.
+%\end{align}
+%\end{subequations}
 
 At the eLDA, we have
 \begin{subequations}
@@ -599,13 +620,13 @@ Interestingly here, there is a strong connection between the LDA and eLDA excita
 \end{split}
 \end{equation}
 The KS orbital energies are given by
-\begin{subequations}
+%\begin{subequations}
-\begin{align}
+%\begin{align}
-	\eps{1}{\eLDA} & = \eHc{1} + 2\eJ{11} + \ldots,
+%	\eps{1}{\eLDA} & = \eHc{1} + 2\eJ{11} + \ldots,
-	\\
+%	\\
-	\eps{2}{\eLDA} & = \eHc{2} + 2\eJ{12} + \ldots.
+%	\eps{2}{\eLDA} & = \eHc{2} + 2\eJ{12} + \ldots.
-\end{align}
+%\end{align}
-\end{subequations}
+%\end{subequations}
 
 
 These equations can be combined to define three ensemble energies
@@ -629,13 +650,13 @@ Similar energies than the ones given in Eqs.~\eqref{eq:EwHF}, \eqref{eq:EwLDA} a
 	\n{}{\ew{}} = (1-\ew{}) \n{}{(0)} + \ew{} \n{}{(1)}.
 \end{equation}
 (This is what one would do in practice, \ie, by performing a KS ensemble calculation.)
-We will label these energies as $\bE{}{\ew{}}$ to avoid confusion.
+We will label these energies as $\tE{}{\ew{}}$ to avoid confusion.
 \begin{widetext}
 For HF, we have
 \begin{equation}
 \label{eq:bEwHF}
 \begin{split}
-	\bE{\HF}{\ew{}} 
+	\tE{\HF}{\ew{}} 
 	& = \titou{\int \hHc(\br{}) \n{}{\ew{}}(\br{}) d\br{}}
 	+ \frac{1}{2} \iint \frac{\n{}{\ew{}}(\br{})\n{}{\ew{}}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}'
 	\\
@@ -648,7 +669,7 @@ In the case of the LDA, it reads
 \begin{equation}
 \label{eq:bEwLDA}
 \begin{split}
-	\bE{\LDA}{\ew{}} 
+	\tE{\LDA}{\ew{}} 
 	& = \titou{\int \hHc(\br{}) \n{}{\ew{}}(\br{}) d\br{}}
 	+ \iint \frac{\n{}{\ew{}}(\br{})\n{}{\ew{}}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}' 
 	+ \int \e{\xc}{\LDA}[\n{}{\ew{}}(\br{})] \n{}{\ew{}}(\br{}) d\br{}
@@ -665,7 +686,7 @@ For eLDA, the ensemble energy can be decomposed as
 \begin{equation}
 \label{eq:bEweLDA}
 \begin{split}
-	\bE{\eLDA}{\ew{}} 
+	\tE{\eLDA}{\ew{}} 
 	& = \titou{\int \hHc(\br{}) \n{}{\ew{}}(\br{}) d\br{}}
 	+ \iint \frac{\n{}{\ew{}}(\br{})\n{}{\ew{}}(\br{}')}{\abs{\br{} - \br{}'}} d\br{} d\br{}' 
 	+ \int \be{\xc}{\ew{}}[\n{}{\ew{}}(\br{})] \n{}{\ew{}}(\br{}) d\br{}
@@ -694,13 +715,64 @@ This would be, for example, the case with the exact xc functional.
 
 Extracting excitation energies from Eqs.~\eqref{eq:bEwHF}, \eqref{eq:bEwLDA} and \eqref{eq:bEweLDA} is more tricky.
 To do so, we will employ Eq.~\eqref{eq:dEdw}.
-The derivative discontinuity, modelled by the last term of the RHS of Eq.~\eqref{eq:dEdw} and only non-zero in the case of an explicitly weight-dependent functional, is straightforward to compute in our case [see Eq.~\eqref{eq:dexcdw}].
+The two first terms are
-
 \begin{align}
-	\Eps{0}{\ew{}} & = 2 \eHc{1} + \ldots,
+	\Eps{0}{\ew{}} & = 2(1-\ew{}) \eps{1}{\ew{}},
-	\\
+	&
-	\Eps{1}{\ew{}} & = 2 \eHc{2} + \ldots.
+	\Eps{1}{\ew{}} & = 2 \ew{} \eps{2}{\ew{}},
 \end{align}
+where the HF, LDA and eLDA weight-dependent orbital energies are
+\begin{subequations}
+\begin{align}
+	\eps{1}{\ew{},\HF} 
+	& = \eHc{1} + (1-\ew{})(2\eJ{11} - \eK{11}) + \ew{}(2\eJ{12} - \eK{12}),
+	\\
+	\eps{2}{\ew{},\HF} 
+	& = \eHc{2} +  (1-\ew{})(2\eJ{12} - \eK{12}) + \ew{}(2\eJ{22} - \eK{22}), 
+\end{align}
+\end{subequations}
+
+\begin{subequations}
+\begin{align}
+\begin{split}
+	\eps{1}{\ew{},\LDA} 
+	& = \eHc{1} + 2(1-\ew{}) \eJ{11} + 2\ew{} \eJ{12}
+	\\
+	& + \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{(0)}(\br{}) d\br{},
+\end{split}
+	\\
+\begin{split}
+	\eps{2}{\ew{},\LDA} 
+	& = \eHc{2} +  2(1-\ew{}) \eJ{12} + 2 \ew{} \eJ{22}
+	\\
+	& + \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{(1)}(\br{}) d\br{},
+\end{split}
+\end{align}
+\end{subequations}
+
+\begin{subequations}
+\begin{align}
+\begin{split}
+	\eps{1}{\ew{},\eLDA} 
+	& = \eHc{1} + (1-\ew{})(2\eJ{11} - \eK{11}) + \ew{}(2\eJ{12} - \eK{12}) 
+	\\
+	& + \frac{1}{2} \int \left. \fdv{\bE{\xc}{\ew{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{(0)}(\br{}) d\br{},
+\end{split}
+	\\
+\begin{split}
+	\eps{2}{\ew{},\eLDA} & = \eHc{2} +  (1-\ew{})(2\eJ{12} - \eK{12}) + \ew{}(2\eJ{22} - \eK{22}) 
+	\\
+	& + \frac{1}{2} \int \left. \fdv{\bE{\xc}{\ew{}}(\n{}{})}{\n{}{}} \right|_{\n{}{} = \n{}{\ew{}}(\br{})} \n{}{(1)}(\br{}) d\br{},
+\end{split}
+\end{align}
+\end{subequations}
+respectively.
+
+The derivative discontinuity is modelled by the last term of the RHS of Eq.~\eqref{eq:dEdw}.
+Note that this contribution is only non-zero in the case of an explicitly weight-dependent functional [see Eq.~\eqref{eq:dexcdw}].
+
+
+
 
 %%%%%%%%%%%%%%%%%%
 %%% CONCLUSION %%%