\documentclass[aip,jcp,reprint,noshowkeys]{revtex4-1} \usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,mhchem,longtable,wrapfig} \usepackage{natbib} \usepackage[extra]{tipa} \bibliographystyle{achemso} \AtBeginDocument{\nocite{achemso-control}} \usepackage{mathpazo,libertine} \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=blue, filecolor=blue, urlcolor=blue, citecolor=blue } \newcommand{\alert}[1]{\textcolor{red}{#1}} \definecolor{darkgreen}{HTML}{009900} \usepackage[normalem]{ulem} \newcommand{\titou}[1]{\textcolor{red}{#1}} \newcommand{\jt}[1]{\textcolor{purple}{#1}} \newcommand{\manu}[1]{\textcolor{darkgreen}{#1}} \newcommand{\toto}[1]{\textcolor{brown}{#1}} \newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}} \newcommand{\trashJT}[1]{\textcolor{purple}{\sout{#1}}} \newcommand{\trashMG}[1]{\textcolor{darkgreen}{\sout{#1}}} \newcommand{\trashAS}[1]{\textcolor{brown}{\sout{#1}}} \newcommand{\MG}[1]{\manu{(\underline{\bf MG}: #1)}} \newcommand{\JT}[1]{\juju{(\underline{\bf JT}: #1)}} \newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}} \newcommand{\AS}[1]{\toto{(\underline{\bf TOTO}: #1)}} \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=blue, filecolor=blue, urlcolor=blue, citecolor=blue } \newcommand{\mc}{\multicolumn} \newcommand{\fnm}{\footnotemark} \newcommand{\fnt}{\footnotetext} \newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}} \newcommand{\SI}{\textcolor{blue}{supporting information}} \newcommand{\QP}{\textsc{quantum package}} % methods \newcommand{\evGW}{evGW} \newcommand{\qsGW}{qsGW} \newcommand{\GOWO}{G$_0$W$_0$} \newcommand{\GW}{GW} \newcommand{\GnWn}[1]{G$_{#1}$W$_{#1}$} % operators \newcommand{\hH}{\Hat{H}} % energies \newcommand{\Ec}{E_\text{c}} \newcommand{\EHF}{E_\text{HF}} \newcommand{\EKS}{E_\text{KS}} \newcommand{\EcK}{E_\text{c}^\text{Klein}} \newcommand{\EcRPA}{E_\text{c}^\text{RPA}} \newcommand{\EcGM}{E_\text{c}^\text{GM}} \newcommand{\EcMP}{E_c^\text{MP2}} \newcommand{\Egap}{E_\text{gap}} \newcommand{\IP}{\text{IP}} \newcommand{\EA}{\text{EA}} \newcommand{\RH}{R_{\ce{H2}}} \newcommand{\RF}{R_{\ce{F2}}} \newcommand{\RBeO}{R_{\ce{BeO}}} % orbital energies \newcommand{\nDIIS}{N^\text{DIIS}} \newcommand{\maxDIIS}{N_\text{max}^\text{DIIS}} \newcommand{\nSat}[1]{N_{#1}^\text{sat}} \newcommand{\eSat}[2]{\epsilon_{#1,#2}} \newcommand{\e}[1]{\epsilon_{#1}} \newcommand{\eHF}[1]{\epsilon^\text{HF}_{#1}} \newcommand{\teHF}[1]{\Tilde{\epsilon}^\text{HF}_{#1}} \newcommand{\eKS}[1]{\epsilon^\text{KS}_{#1}} \newcommand{\eQP}[1]{\epsilon^\text{QP}_{#1}} \newcommand{\eGOWO}[1]{\epsilon^\text{\GOWO}_{#1}} \newcommand{\eGW}[1]{\epsilon^\text{\GW}_{#1}} \newcommand{\eGnWn}[2]{\epsilon^\text{\GnWn{#2}}_{#1}} \newcommand{\de}[1]{\Delta\epsilon_{#1}} \newcommand{\deHF}[1]{\Delta\epsilon^\text{HF}_{#1}} \newcommand{\Om}[1]{\Omega_{#1}} \newcommand{\eHOMO}{\epsilon_\text{HOMO}} \newcommand{\eLUMO}{\epsilon_\text{LUMO}} \newcommand{\HOMO}{\text{HOMO}} \newcommand{\LUMO}{\text{LUMO}} % Matrix elements \newcommand{\A}[1]{A_{#1}} \newcommand{\B}[1]{B_{#1}} \newcommand{\tA}{\Tilde{A}} \newcommand{\tB}{\Tilde{B}} \renewcommand{\S}[1]{S_{#1}} \newcommand{\G}[1]{G_{#1}} \newcommand{\Po}[1]{P_{#1}} \newcommand{\W}[1]{W_{#1}} \newcommand{\Wc}[1]{W^\text{c}_{#1}} \newcommand{\vc}[1]{v_{#1}} \newcommand{\SigX}[1]{\Sigma^\text{x}_{#1}} \newcommand{\SigC}[1]{\Sigma^\text{c}_{#1}} \newcommand{\tSigC}[1]{\Tilde{\Sigma}^\text{c}_{#1}} \newcommand{\SigCp}[1]{\Sigma^\text{p}_{#1}} \newcommand{\SigCh}[1]{\Sigma^\text{h}_{#1}} \newcommand{\SigGW}[1]{\Sigma^\text{\GW}_{#1}} \newcommand{\Z}[1]{Z_{#1}} % Matrices \newcommand{\bG}{\boldsymbol{G}} \newcommand{\bW}{\boldsymbol{W}} \newcommand{\bvc}{\boldsymbol{v}} \newcommand{\bSig}{\boldsymbol{\Sigma}} \newcommand{\bSigX}{\boldsymbol{\Sigma}^\text{x}} \newcommand{\bSigC}{\boldsymbol{\Sigma}^\text{c}} \newcommand{\bSigGW}{\boldsymbol{\Sigma}^\text{\GW}} \newcommand{\be}{\boldsymbol{\epsilon}} \newcommand{\bDelta}{\boldsymbol{\Delta}} \newcommand{\beHF}{\boldsymbol{\epsilon}^\text{HF}} \newcommand{\beGW}{\boldsymbol{\epsilon}^\text{\GW}} \newcommand{\beGnWn}[1]{\boldsymbol{\epsilon}^\text{\GnWn{#1}}} \newcommand{\bdeGnWn}[1]{\Delta\boldsymbol{\epsilon}^\text{\GnWn{#1}}} \newcommand{\bde}{\boldsymbol{\Delta\epsilon}} \newcommand{\bdeHF}{\boldsymbol{\Delta\epsilon}^\text{HF}} \newcommand{\bdeGW}{\boldsymbol{\Delta\epsilon}^\text{GW}} \newcommand{\bOm}{\boldsymbol{\Omega}} \newcommand{\bA}{\boldsymbol{A}} \newcommand{\bB}{\boldsymbol{B}} \newcommand{\bX}{\boldsymbol{X}} \newcommand{\bY}{\boldsymbol{Y}} \newcommand{\bZ}{\boldsymbol{Z}} \newcommand{\fc}{f_\text{c}} \newcommand{\Vc}{V_\text{c}} \newcommand{\MO}[1]{\phi_{#1}} % coordinates \newcommand{\br}[1]{\mathbf{r}_{#1}} \renewcommand{\b}[1]{\mathbf{#1}} \renewcommand{\d}{\text{d}} \newcommand{\dbr}[1]{d\br{#1}} \renewcommand{\bra}[1]{\ensuremath{\langle #1 \vert}} \renewcommand{\ket}[1]{\ensuremath{\vert #1 \rangle}} \renewcommand{\braket}[2]{\ensuremath{\langle #1 \vert #2 \rangle}} \newcommand{\ISCD}{Institut des Sciences du Calcul et des Donn\'ees, Sorbonne Universit\'e, Paris, France} \newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France} \newcommand{\LCT}{Laboratoire de Chimie Th\'eorique (UMR 7616), Sorbonne Universit\'e, CNRS, Paris, France} \begin{document} \title{A Density-Based Basis Set Correction for GW Methods} \author{Bath\'elemy Pradines} \affiliation{\LCT} \affiliation{\ISCD} \author{Emmanuel Giner} \affiliation{\LCT} \author{Anthony Scemama} \affiliation{\LCPQ} \author{Julien Toulouse} \affiliation{\LCT} \author{Pierre-Fran\c{c}ois Loos} \email[Corresponding author: ]{loos@irsamc.ups-tlse.fr} \affiliation{\LCPQ} \begin{abstract} %\begin{wrapfigure}[13]{o}[-1.25cm]{0.5\linewidth} % \centering % \includegraphics[width=\linewidth]{TOC} %\end{wrapfigure} \end{abstract} \maketitle %%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} \label{sec:intro} %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Theory} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Many-body Green-function theory with DFT basis-set correction} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Following Ref.~\onlinecite{GinPraFerAssSavTou-JCP-18}, we start by defining, for a $N$-electron system with nuclei-electron potential $v_\text{ne}(\b{r})$, the approximate ground-state energy for one-electron densities $n$ which are ``representable'' in a finite basis set ${\cal B}$ \begin{equation} E_0^{\cal B} = \min_{n \in {\cal D}^{\cal B}} \left\{ F[n] + \int v_\text{ne}(\b{r}) n(\b{r}) \d\b{r}\right\}, \label{E0B} \end{equation} where ${\cal D}^{\cal B}$ is the set of $N$-representable densities which can be extracted from a wave function $\Psi^{\cal B}$ expandable in the Hilbert space generated by ${\cal B}$. In this expression, $F[n]=\min_{\Psi\to n} \bra{\Psi} \hat{T} + \hat{W}_\text{ee}\ket{\Psi}$ is the exact Levy-Lieb universal density functional, where $\hat{T}$ and $\hat{W}_\text{ee}$ are the kinetic and electron-electron interaction operators, which is then decomposed as \begin{equation} F[n] = F^{\cal B}[n] + \bar{E}^{\cal B}[n], \label{Fn} \end{equation} where $F^{\cal B}[n]$ is the Levy-Lieb density functional with wave functions $\Psi^{\cal B}$ expandable in the Hilbert space generated by ${\cal B}$ \begin{equation} F^{\cal B}[n] = \min_{\Psi^{\cal B}\to n} \bra{\Psi^{\cal B}} \hat{T} + \hat{W}_\text{ee}\ket{\Psi^{\cal B}}, \end{equation} and $\bar{E}^{\cal B}[n]$ is the complementary basis-correction density functional. In the present work, instead of using wave-function methods for calculating $F^{\cal B}[n]$, we reexpress it with a contrained search over $N$-representable one-electron Green functions $G^{\cal B}(\b{r},\b{r}',\omega)$ representable in the basis set ${\cal B}$ \begin{equation} F^{\cal B}[n] = \min_{G^{\cal B}\to n} \Omega^{\cal B}[G^{\cal B}], \label{FBn} \end{equation} where $\Omega^{\cal B}[G]$ is chosen to be a Klein-like energy functional of the Green function (see, e.g., Refs.~\onlinecite{SteLee-BOOK-13,MarReiCep-BOOK-16,DahLee-JCP-05,DahLeeBar-IJQC-05,DahLeeBar-PRA-06}) \begin{equation} \Omega^{\cal B}[G] = \Tr \left[\ln ( - G ) \right] - \Tr \left[ (G_\text{s}^{\cal B})^{-1} G -1 \right] + \Phi_\text{Hxc}^{\cal B}[G], \label{OmegaB} \end{equation} where $(G_\text{s}^{\cal B})^{-1}$ is the projection into ${\cal B}$ of the inverse free-particle Green function $(G_\text{s})^{-1}(\b{r},\b{r}',\omega)= (\omega + (1/2) \nabla_\b{r}^2 )\delta(\b{r}-\b{r}')$ and we have used the notation $\Tr [A B] = 1/(2\pi i) \int_{-\infty}^{+\infty} \! \d \omega \, e^{i \omega 0^+} \! \iint \! \d \b{r} \d \b{r}' A(\b{r},\b{r}',\omega) B(\b{r}',\b{r},\omega)$. In Eq.~(\ref{OmegaB}), $\Phi_\text{Hxc}^{\cal B}[G]$ is a Hartree-exchange-correlation (Hxc) functional of the Green functional such as its functional derivatives yields the Hxc self-energy in the basis: $\delta \Phi_\text{Hxc}^{\cal B}[G]/\delta G(\b{r},\b{r}',\omega) = \Sigma_\text{Hxc}^{\cal B}[G](\b{r},\b{r}',\omega)$. Inserting Eqs.~(\ref{Fn}) and~(\ref{FBn}) into Eq.~(\ref{E0B}), we finally arrive at \begin{equation} E_0^{\cal B} = \min_{G^{\cal B}} \left\{ \Omega^{\cal B}[G^{\cal B}] + \int v_\text{ne}(\b{r}) n_{G^{\cal B}}(\b{r}) \d\b{r} + \bar{E}^{\cal B}[n_{G^{\cal B}}] \right\}, \label{E0BGB} \end{equation} where the minimization is over $N$-representable one-electron Green functions $G^{\cal B}(\b{r},\b{r}',\omega)$ representable in the basis set ${\cal B}$. The stationary condition from Eq.~(\ref{E0BGB}) gives the following Dyson equation \begin{equation} (G^{\cal B})^{-1} = (G_\text{0}^{\cal B})^{-1}- \Sigma_\text{Hxc}^{\cal B}[G^{\cal B}]- \bar{\Sigma}^{\cal B}[n_{G^{\cal B}}], \label{Dyson} \end{equation} where $(G_\text{0}^{\cal B})^{-1}$ is the basis projection of the inverse non-interacting Green function with potential $v_\text{ne}(\b{r})$, $(G_\text{0})^{-1}(\b{r},\b{r}',\omega)= (\omega + (1/2) \nabla_\b{r}^2 + v_\text{ne}(\b{r}) + \lambda)\delta(\b{r}-\b{r}')$ with the chemical potential $\lambda$, and $\bar{\Sigma}^{\cal B}$ is a frequency-independent local self-energy coming from functional derivative of the complementary basis-correction density functional \begin{equation} \bar{\Sigma}^{\cal B}[n](\b{r},\b{r}') = \bar{v}^{\cal B}[n](\b{r}) \delta(\b{r}-\b{r}'), \end{equation} with $\bar{v}^{\cal B}[n](\b{r}) = \delta \bar{E}^{\cal B}[n] / \delta n(\b{r})$. The solution of the Dyson equation~(\ref{Dyson}) gives the Green function $G^{\cal B}(\b{r},\b{r}',\omega)$ which is not exact but should converge more rapidly with the basis set thanks to the presence of the basis-set correction $\bar{\Sigma}^{\cal B}$. Of course, in the complete-basis-set (CBS) limit, the basis-set correction vanishes, $\bar{\Sigma}^{{\cal B}\to \text{CBS}}=0$, and the Green function becomes exact, $G^{{\cal B}\to \text{CBS}}=G$. %From Julien: %\begin{equation} %\fdv{E[n_G]}{G(r,r',\omega)} = \int \fdv{E[n_G]}{n(r'')}] \fdv{n_G(r'')}{G(r,r',w)} dr'' %\end{equation} % %\begin{equation} %n_G(r'') = i \int G(r'',r'',w) d\omega %\end{equation} % % %\begin{equation} %\fdv{n_G(r'')}{G(r,r',w)} = \delta(r -r') \delta (r'-r'') %\end{equation} % % %\begin{equation} %\begin{split} % \fdv{E[n_G]}{G(r,r',w)} % & = \int \fdv{E[n_G]}{n(r'')} \delta(r -r') \delta (r'-r'') dr'' % \\ % & = \fdv{E[n_G]}{n(r)} \delta(r -r') % \\ % & = v[n_G](r) \delta(r -r') %\end{split} %\end{equation} \subsection{The GW Approximation} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Here, we provide self-contained summary of the main equations and quantities behind {\GOWO} and {\evGW}. More details can be found, for example, in Refs.~\citenum{vanSetten_2013, Kaplan_2016, Bruneval_2016}. For a given (occupied or virtual) orbital $p$, the correlation part of the self-energy is conveniently split in its hole (h) and particle (p) contributions \begin{equation} \label{eq:SigC} \SigC{p}(\omega) = \SigCp{p}(\omega) + \SigCh{p}(\omega), \end{equation} which, within the GW approximation, read \begin{subequations} \begin{align} \label{eq:SigCh} \SigCh{p}(\omega) & = 2 \sum_{i}^\text{occ} \sum_{x} \frac{[pi|x]^2}{\omega - \e{i} + \Om{x} - i \eta}, \\ \label{eq:SigCp} \SigCp{p}(\omega) & = 2 \sum_{a}^\text{virt} \sum_{x} \frac{[pa|x]^2}{\omega - \e{a} - \Om{x} + i \eta}, \end{align} \end{subequations} where $\eta$ is a positive infinitesimal. The screened two-electron integrals \begin{equation} [pq|x] = \sum_{ia} (pq|ia) (\bX+\bY)_{ia}^{x} \end{equation} are obtained via the contraction of the bare two-electron integrals \cite{Gill_1994} $(pq|rs)$ and the transition densities $(\bX+\bY)_{ia}^{x}$ originating from a random phase approximation (RPA) calculation \cite{Casida_1995, Dreuw_2005} \begin{equation} \label{eq:LR} \begin{pmatrix} \bA & \bB \\ \bB & \bA \\ \end{pmatrix} \begin{pmatrix} \bX \\ \bY \\ \end{pmatrix} = \bOm \begin{pmatrix} \boldsymbol{1} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{-1} \\ \end{pmatrix} \begin{pmatrix} \bX \\ \bY \\ \end{pmatrix}, \end{equation} with \begin{align} \label{eq:RPA} A_{ia,jb} & = \delta_{ij} \delta_{ab} (\epsilon_a - \epsilon_i) + 2 (ia|jb), & B_{ia,jb} & = 2 (ia|bj), \end{align} and $\delta_{pq}$ is the Kronecker delta. \cite{NISTbook} The one-electron energies $\epsilon_p$ in \eqref{eq:SigCh}, \eqref{eq:SigCp} and \eqref{eq:RPA} are either the HF or the GW quasiparticle energies. Equation \eqref{eq:LR} also provides the neutral excitation energies $\Om{x}$. In practice, there exist two ways of determining the {\GOWO} QP energies. \cite{Hybertsen_1985a, vanSetten_2013} In its ``graphical'' version, they are provided by one of the many solutions of the (non-linear) QP equation \begin{equation} \label{eq:QP-G0W0} \omega = \eHF{p} + \Re[\SigC{p}(\omega)]. \end{equation} In this case, special care has to be taken in order to select the ``right'' solution, known as the QP solution. In particular, it is usually worth calculating its renormalization weight (or factor), $\Z{p}(\eHF{p})$, where \begin{equation} \label{eq:Z} \Z{p}(\omega) = \qty[ 1 - \pdv{\Re[\SigC{p}(\omega)]}{\omega} ]^{-1}. \end{equation} Because of sum rules, \cite{Martin_1959, Baym_1961, Baym_1962, vonBarth_1996} the other solutions, known as satellites, share the remaining weight. In a well-behaved case (belonging to the weakly correlated regime), the QP weight is much larger than the sum of the satellite weights, and of the order of $0.7$-$0.9$. Within the linearized version of {\GOWO}, one assumes that \begin{equation} \label{eq:SigC-lin} \SigC{p}(\omega) \approx \SigC{p}(\eHF{p}) + (\omega - \eHF{p}) \left. \pdv{\SigC{p}(\omega)}{\omega} \right|_{\omega = \eHF{p}}, \end{equation} that is, the self-energy behaves linearly in the vicinity of $\omega = \eHF{p}$. Substituting \eqref{eq:SigC-lin} into \eqref{eq:QP-G0W0} yields \begin{equation} \label{eq:QP-G0W0-lin} \eGOWO{p} = \eHF{p} + \Z{p}(\eHF{p}) \Re[\SigC{p}(\eHF{p})]. \end{equation} Unless otherwise stated, in the remaining of this paper, the {\GOWO} QP energies are determined via the linearized method. In the case of {\evGW}, the QP energy, $\eGW{p}$, are obtained via Eq.~\eqref{eq:QP-G0W0}, which has to be solved self-consistently due to the QP energy dependence of the self-energy [see Eq.~\eqref{eq:SigC}]. \cite{Hybertsen_1986, Shishkin_2007, Blase_2011, Faber_2011} At least in the weakly correlated regime where a clear QP solution exists, we believe that, within {\evGW}, the self-consistent algorithm should select the solution of the QP equation \eqref{eq:QP-G0W0} with the largest renormalization weight $\Z{p}(\eGW{p})$. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Basis Set Correction} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The present basis set correction is a two-level correction. First, one has to correct the neutral excitations $\Om{x}$ from the RPA calculation. The corrected matrix elements read \begin{align} \label{eq:RPA} \tA_{ia,jb} & = \A{ia,jb} + (ia|\fc|jb), & \tB_{ia,jb} & = \B{ia,jb} + (ia|\fc|bj), \end{align} where the elements $\A{ia,jb}$ and $\B{ia,jb}$ are given by Eq.~\eqref{eq:RPA}. \begin{equation} \fc(\br{1},\br{2})= \frac{\delta^2 \Ec}{\delta n(\br{1})\delta n(\br{2})} \end{equation} In a second time, we correct the GW energy \begin{equation} \tSigC{p} = \SigC{p} + (p|\Vc|p) \end{equation} with \begin{equation} \Vc(\br{}) = \fdv{\Ec}{n(\br{})} \end{equation} %%%%%%%%%%%%%%%%%%%%%%%% \section{Computational details} \label{sec:compdetails} %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% \section{Results and Discussion} \label{sec:res} %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% \section{Conclusion} \label{sec:ccl} %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %\section*{Supporting Information Available} %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% \begin{acknowledgements} This work was performed using HPC resources from GENCI-TGCC (Grant No.~2018-A0040801738), CALMIP (Toulouse) under allocation 2019-18005 and the Jarvis-Alpha cluster from the \textit{Institut Parisien de Chimie Physique et Th\'eorique}. \end{acknowledgements} %%%%%%%%%%%%%%%%%%%%%%%% \bibliography{GW-srDFT,GW-srDFT-control,biblio} \end{document}