From 2e25c74b8d629033751ea1ad72e81bec1fb30784 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 27 Nov 2019 17:02:22 +0100 Subject: [PATCH] making progress on the ensemble energies --- Manuscript/FarDFT.tex | 184 +++++++++++++++++++++++++++++------------- 1 file changed, 128 insertions(+), 56 deletions(-) diff --git a/Manuscript/FarDFT.tex b/Manuscript/FarDFT.tex index eef09af..d4df3dd 100644 --- a/Manuscript/FarDFT.tex +++ b/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{\bw}{\bm{w}} -\newcommand{\bG}{\bm{G}} -\newcommand{\bS}{\bm{S}} -\newcommand{\bGamma}[1]{\bm{\Gamma}^{#1}} -\newcommand{\bHc}{\bm{h}} -\newcommand{\bF}[1]{\bm{F}^{#1}} +\newcommand{\br}{\boldsymbol{r}} +\newcommand{\bw}{\boldsymbol{w}} +\newcommand{\bG}{\boldsymbol{G}} +\newcommand{\bS}{\boldsymbol{S}} +\newcommand{\bGamma}[1]{\boldsymbol{\Gamma}^{#1}} +\newcommand{\bHc}{\boldsymbol{h}} +\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,15 +240,35 @@ 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. -Equation \eqref{eq:dEdw} is our working equation for computing excitation energies. - +where $\hHc(\br{}) = -\nabla^2/2 + \vext(\br{})$, and +\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{}{}) - \\ - & = \Cx{\ew{}} \n{}{1/3} -\end{split} + = (1-\ew{}) \e{\ex}{(0)}(\n{}{}) + \ew{} \e{\ex}{(1)}(\n{}{}) + = \Cx{\ew{}} \n{}{1/3} \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 -\begin{subequations} -\begin{align} - \eps{1}{\HF} & = \eHc{1} + 2\eJ{11} - \eK{11}, - \\ - \eps{2}{\HF} & = \eHc{2} + 2\eJ{12} - \eK{12}. -\end{align} -\end{subequations} +%The HF orbital energies are +%\begin{subequations} +%\begin{align} +% \eps{1}{\HF} & = \eHc{1} + 2\eJ{11} - \eK{11}, +% \\ +% \eps{2}{\HF} & = \eHc{2} + 2\eJ{12} - \eK{12}. +%\end{align} +%\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 -\begin{subequations} -\begin{align} - \eps{1}{\LDA} & = \eHc{1} + 2\eJ{11} + \ldots, - \\ - \eps{2}{\LDA} & = \eHc{2} + 2\eJ{12} + \ldots. -\end{align} -\end{subequations} +%The KS orbital energies are given by +%\begin{subequations} +%\begin{align} +% \eps{1}{\LDA} +% & = \eHc{1} + 2\eJ{11} +% + \frac{1}{2} \int \left. \fdv{\E{\xc}{\LDA}[\n{}{}]}{\n{}{}} \right|_{\n{}{} = \n{}{(0)}(\br{})} \n{}{(0)}(\br{}) d\br{}, +% \\ +% \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{align} - \eps{1}{\eLDA} & = \eHc{1} + 2\eJ{11} + \ldots, - \\ - \eps{2}{\eLDA} & = \eHc{2} + 2\eJ{12} + \ldots. -\end{align} -\end{subequations} +%\begin{subequations} +%\begin{align} +% \eps{1}{\eLDA} & = \eHc{1} + 2\eJ{11} + \ldots, +% \\ +% \eps{2}{\eLDA} & = \eHc{2} + 2\eJ{12} + \ldots. +%\end{align} +%\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{1}{\ew{}} & = 2 \eHc{2} + \ldots. + \Eps{0}{\ew{}} & = 2(1-\ew{}) \eps{1}{\ew{}}, + & + \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 %%%