802 lines
49 KiB
TeX
802 lines
49 KiB
TeX
\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
|
|
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,bbold,siunitx,xspace}
|
|
\usepackage[version=4]{mhchem}
|
|
|
|
\usepackage[utf8]{inputenc}
|
|
\usepackage[T1]{fontenc}
|
|
|
|
\usepackage{hyperref}
|
|
\hypersetup{
|
|
colorlinks,
|
|
linkcolor={red!50!black},
|
|
citecolor={red!70!black},
|
|
urlcolor={red!80!black}
|
|
}
|
|
|
|
\usepackage{listings}
|
|
\definecolor{codegreen}{rgb}{0.58,0.4,0.2}
|
|
\definecolor{codegray}{rgb}{0.5,0.5,0.5}
|
|
\definecolor{codepurple}{rgb}{0.25,0.35,0.55}
|
|
\definecolor{codeblue}{rgb}{0.30,0.60,0.8}
|
|
\definecolor{backcolour}{rgb}{0.98,0.98,0.98}
|
|
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
|
|
|
|
\definecolor{sqred}{rgb}{0.85,0.1,0.1}
|
|
\definecolor{sqgreen}{rgb}{0.25,0.65,0.15}
|
|
\definecolor{sqorange}{rgb}{0.90,0.50,0.15}
|
|
\definecolor{sqblue}{rgb}{0.10,0.3,0.60}
|
|
|
|
\lstdefinestyle{mystyle}{
|
|
backgroundcolor=\color{backcolour},
|
|
commentstyle=\color{codegreen},
|
|
keywordstyle=\color{codeblue},
|
|
numberstyle=\tiny\color{codegray},
|
|
stringstyle=\color{codepurple},
|
|
basicstyle=\ttfamily\footnotesize,
|
|
breakatwhitespace=false,
|
|
breaklines=true,
|
|
captionpos=b,
|
|
keepspaces=true,
|
|
numbers=left,
|
|
numbersep=5pt,
|
|
numberstyle=\ttfamily\tiny\color{mygray},
|
|
showspaces=false,
|
|
showstringspaces=false,
|
|
showtabs=false,
|
|
tabsize=2
|
|
}
|
|
|
|
\lstset{style=mystyle}
|
|
|
|
\newcommand{\titou}[1]{\textcolor{red}{#1}}
|
|
\newcommand{\trashPFL}[1]{\textcolor{\red}{\sout{#1}}}
|
|
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
|
|
\newcommand{\ant}[1]{\textcolor{teal}{#1}}
|
|
\newcommand{\ANT}[1]{\ant{(\underline{\bf ANT}: #1)}}
|
|
|
|
% addresses
|
|
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
|
|
|
|
\begin{document}
|
|
|
|
\title{A Similarity Renormalization Group Approach To Many-Body Perturbation Theory}
|
|
|
|
\author{Antoine \surname{Marie}}
|
|
\email{amarie@irsamc.ups-tlse.fr}
|
|
\affiliation{\LCPQ}
|
|
|
|
\author{Pierre-Fran\c{c}ois \surname{Loos}}
|
|
\email{loos@irsamc.ups-tlse.fr}
|
|
\affiliation{\LCPQ}
|
|
|
|
\begin{abstract}
|
|
Here comes the abstract.
|
|
%\bigskip
|
|
%\begin{center}
|
|
% \boxed{\includegraphics[width=0.5\linewidth]{TOC}}
|
|
%\end{center}
|
|
%\bigskip
|
|
\end{abstract}
|
|
|
|
\maketitle
|
|
|
|
%=================================================================%
|
|
\section{Introduction}
|
|
\label{sec:intro}
|
|
%=================================================================%
|
|
|
|
One-body Green's functions provide a natural and elegant way to access the charged excitation energies of a physical system. \cite{CsanakBook,FetterBook,Martin_2016,Golze_2019}
|
|
The non-linear Hedin equations consist of a closed set of equations leading to the exact interacting one-body Green's function and, therefore, to a wealth of properties such as the total energy, density, ionization potentials, electron affinities, as well as spectral functions, without the explicit knowledge of the wave functions associated with the neutral and charged states of the system. \cite{Hedin_1965}
|
|
Unfortunately, solving exactly Hedin's equations is usually out of reach and one must resort to approximations.
|
|
In particular, the $GW$ approximation, \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019,Bruneval_2021} which has been first introduced in the context of solids \cite{Strinati_1980,Strinati_1982a,Strinati_1982b,Hybertsen_1985,Hybertsen_1986,Godby_1986,Godby_1987,Godby_1987a,Godby_1988,Blase_1995} and is now widely applied to molecular systems, \cite{Rohlfing_1999a,Horst_1999,Puschnig_2002,Tiago_2003,Rocca_2010,Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018,Blase_2018,Liu_2020,Li_2017,Li_2019,Li_2020,Li_2021,Blase_2020,Holzer_2018a,Holzer_2018b,Loos_2020e,Loos_2021,McKeon_2022} yields accurate charged excitation energies for weakly correlated systems \cite{Hung_2017,vanSetten_2015,vanSetten_2018,Caruso_2016,Korbel_2014,Bruneval_2021} at a relatively low computational cost. \cite{Foerster_2011,Liu_2016,Wilhelm_2018,Forster_2021,Duchemin_2019,Duchemin_2020,Duchemin_2021}
|
|
|
|
The $GW$ method approximates the self-energy $\Sigma$ which relates the exact interacting Green's function $G$ to a non-interacting reference version $G_0$ through a Dyson equation of the form
|
|
\begin{equation}
|
|
\label{eq:dyson}
|
|
G(1,2) = G_0(1,2) + \int d(34) G_0(1,3)\Sigma(3,4) G(4,2),
|
|
\end{equation}
|
|
where $1 = (\bx_1, t_1)$ is a composite coordinate gathering spin-space and time variables.
|
|
The self-energy encapsulates all the Hartree-exchange-correlation effects which are not taken into account in the reference system.
|
|
%Throughout this manuscript the references are chosen to be the Hartree-Fock (HF) ones so that the self-energy only account for the missing correlation.
|
|
Approximating $\Sigma$ as the first-order term of its perturbative expansion with respect to the screened Coulomb potential $W$ yields the so-called $GW$ approximation \cite{Hedin_1965,Martin_2016}
|
|
\begin{equation}
|
|
\label{eq:gw_selfenergy}
|
|
\Sigma^{\GW}(1,2) = \ii G(1,2) W(1,2).
|
|
\end{equation}
|
|
Diagrammatically, $GW$ corresponds to a resummation of the (time-dependent) direct ring diagrams via the computation of the random-phase approximation (RPA) polarizability \cite{Ren_2012,Chen_2017} and is thus particularly well suited for weak correlation.
|
|
|
|
\titou{Alternatively, one can choose to define $\Sigma$ as the $n$th-order expansion in terms of the bare Coulomb interaction $v$ leading to the GF($n$) class of approximations. \cite{SzaboBook,Ortiz_2013,Hirata_2015,Hirata_2017}
|
|
The GF(2) approximation \cite{Casida_1989,Casida_1991,Phillips_2014,Phillips_2015,Rusakov_2014,Rusakov_2016,Backhouse_2021,Backhouse_2020b,Backhouse_2020a,Pokhilko_2021a,Pokhilko_2021b,Pokhilko_2022} is also known as the second Born approximation in condensed matter physics. \cite{Stefanucci_2013}}
|
|
|
|
Despite a wide range of successes, many-body perturbation theory has well-documented limitations. \cite{Kozik_2014,Stan_2015,Rossi_2015,Tarantino_2017,Schaefer_2013,Schaefer_2016,Gunnarsson_2017,vanSetten_2015,Maggio_2017a,Duchemin_2020}
|
|
For example, modeling core electron spectroscopy requires core ionization energies which have been proven to be challenging for routine $GW$ calculations. \cite{Golze_2018,Golze_2020,Li_2022}
|
|
Many-body perturbation theory can also be used to access optical excitation energies through the Bethe-Salpeter equation. \cite{Salpeter_1951,Strinati_1988,Blase_2018,Blase_2020} However, the accuracy is not yet satisfying for triplet excited states, where instabilities often occur. \cite{Bruneval_2015,Jacquemin_2017a,Jacquemin_2017b,Holzer_2018a}
|
|
Therefore, even if $GW$ offers a good trade-off between accuracy and computational cost, some situations might require higher precision.
|
|
Unfortunately, defining a systematic way to go beyond $GW$ via the inclusion of vertex corrections has been demonstrated to be a tricky task. \cite{Baym_1961,Baym_1962,DeDominicis_1964a,DeDominicis_1964b,Bickers_1989a,Bickers_1989b,Bickers_1991,Hedin_1999,Bickers_2004,Shirley_1996,DelSol_1994,Schindlmayr_1998,Morris_2007,Shishkin_2007b,Romaniello_2009a,Romaniello_2012,Gruneis_2014,Hung_2017,Maggio_2017b,Mejuto-Zaera_2022}
|
|
For example, Lewis and Berkelbach have shown that naive vertex corrections can even worsen the quasiparticle energies with respect to $GW$. \cite{Lewis_2019}
|
|
We refer the reader to the recent review by Golze and co-workers (see Ref.~\onlinecite{Golze_2019}) for an extensive list of current challenges in many-body perturbation theory.
|
|
|
|
Recently, it has been shown that a variety of physical quantities, such as charged and neutral excitations energies as well as correlation and total energies, computed within many-body perturbation theory exhibit unphysical discontinuities. \cite{Loos_2018b,Veril_2018,Loos_2020e,Berger_2021,DiSabatino_2021,Monino_2022,Scott_2023}
|
|
Even more worrying, these discontinuities can happen in the weakly correlated regime where the $GW$ approximation is supposed to be valid.
|
|
These discontinuities have been traced back to a transfer of spectral weight between two solutions of the quasi-particle equation, \cite{Monino_2022} and is another occurrence of the infamous intruder-state problem.\cite{Andersson_1994,Andersson_1995a,Roos_1995,Forsberg_1997,Olsen_2000,Choe_2001}
|
|
In addition, systems, where the quasiparticle equation admits two solutions with similar spectral weights, are known to be particularly difficult to converge for partially self-consistent $GW$ schemes. \cite{Veril_2018,Forster_2021,Monino_2022}
|
|
|
|
In a recent study, Monino and Loos showed that these discontinuities could be removed by the introduction of a regularizer inspired by the similarity renormalization group (SRG) in the quasiparticle equation. \cite{Monino_2022}
|
|
Encouraged by the recent successes of regularization schemes in many-body quantum chemistry methods, as in single- and multi-reference perturbation theory, \cite{Lee_2018a,Shee_2021,Evangelista_2014b,ChenyangLi_2019a,Battaglia_2022} the present work investigates the application of the SRG formalism to many-body perturbation theory in its $GW$ \titou{and GF(2) variants}.
|
|
The SRG has been developed independently by Wegner \cite{Wegner_1994} and Glazek and Wilson \cite{Glazek_1993,Glazek_1994} in the context of condensed matter systems and light-front quantum field theories, respectively.
|
|
This formalism has been introduced in quantum chemistry by White \cite{White_2002} before being explored in more detail by Evangelista and coworkers in the context of multi-reference electron correlation theories. \cite{Evangelista_2014b,ChenyangLi_2015, ChenyangLi_2016,ChenyangLi_2017,ChenyangLi_2018,ChenyangLi_2019a,Zhang_2019,ChenyangLi_2021,Wang_2021,Wang_2023}
|
|
The SRG has also been successful in the context of nuclear structure theory, where it was first developed as a mature computational tool thanks to the work of several research groups.
|
|
\cite{Bogner_2007,Tsukiyama_2011,Tsukiyama_2012,Hergert_2013,Hergert_2016,Frosini_2022a,Frosini_2022b,Frosini_2022c}
|
|
See Ref.~\onlinecite{Hergert_2016a} for a recent review in this field.
|
|
|
|
The SRG transformation aims at decoupling an internal (or reference) space from an external space while incorporating information about their coupling in the reference space.
|
|
This process often results in the appearance of intruder states. \cite{Evangelista_2014b,ChenyangLi_2019a}
|
|
However, SRG is particularly well-suited to avoid them because the decoupling of each external configuration is inversely proportional to its energy difference with the reference space.
|
|
By definition, intruder states have energies that are close to the reference energy, and therefore are the last to be decoupled.
|
|
By stopping the SRG transformation once all external configurations except the intruder states have been decoupled,
|
|
correlation effects between the internal and external spaces can be incorporated (or folded) without the presence of intruder states.
|
|
|
|
The goal of this manuscript is to determine if the SRG formalism can effectively address the issue of intruder states in many-body perturbation theory, as it has in other areas of electronic and nuclear structure theory.
|
|
This open question will lead us to an intruder-state-free static approximation of the self-energy derived from first-principles that can be employed in \titou{qs$GW$} calculations.
|
|
\titou{Here, we focus on the $GW$ approximation but the subsequent derivations can be straightforwardly applied to other approximations such as GF(2) or $T$-matrix.}
|
|
|
|
The manuscript is organized as follows.
|
|
We begin by reviewing the $GW$ approximation in Sec.~\ref{sec:gw} and then briefly review the SRG formalism in Sec.~\ref{sec:srg}.
|
|
A perturbative analysis of SRG applied to $GW$ is presented in Sec.~\ref{sec:srggw}.
|
|
The computational details are provided in Sec.~\ref{sec:comp_det} before turning to the results section.
|
|
This section starts by
|
|
|
|
%=================================================================%
|
|
%\section{Theoretical background}
|
|
% \label{sec:theoretical}
|
|
%=================================================================%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{The $GW$ approximation}
|
|
\label{sec:gw}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
The central equation of many-body perturbation theory based on Hedin's equations is the so-called quasi-particle equation which, within the $GW$ approximation, reads
|
|
\begin{equation}
|
|
\label{eq:quasipart_eq}
|
|
\qty[ \bF + \bSig(\omega = \epsilon_p) ] \psi_p(\bx) = \epsilon_p \psi_p(\bx),
|
|
\end{equation}
|
|
where $\bF$ is the Fock matrix \cite{SzaboBook} and $\bSig(\omega)$ is (the correlation part of) the $GW$ self-energy.
|
|
Both are $K \times K$ matrices with $K$ the number of one-electron orbitals.
|
|
The self-energy can be physically understood as a correction to the Hartree-Fock (HF) problem (represented by $\bF$) accounting for dynamical screening effects.
|
|
Similarly to the HF case, Eq.~\eqref{eq:quasipart_eq} has to be solved self-consistently.
|
|
\titou{Note that $\bSig(\omega)$ is dynamical, \ie it depends on the one-electron orbitals $\psi_p(\bx)$ and their corresponding energies $\epsilon_p$, while $\bF$ depends only on the orbitals.}
|
|
|
|
The matrix elements of $\bSig(\omega)$ have the following analytic expression \cite{Hedin_1999,Tiago_2006,Bruneval_2012,vanSetten_2013,Bruneval_2016}
|
|
\begin{equation}
|
|
\label{eq:GW_selfenergy}
|
|
\Sigma_{pq}(\omega)
|
|
= \sum_{i\nu} \frac{W_{pi\nu} W_{qi\nu}}{\omega - \epsilon_i + \Omega_{\nu} - \ii \eta}
|
|
+ \sum_{a\nu} \frac{W_{pa\nu}W_{qa\nu}}{\omega - \epsilon_a - \Omega_{\nu} + \ii \eta},
|
|
\end{equation}
|
|
where $\eta$ is a positive infinitesimal and the screened two-electron integrals are
|
|
\begin{equation}
|
|
\label{eq:GW_sERI}
|
|
W_{pq\nu} = \sum_{ia}\eri{pi}{qa}\qty( \bX_{\nu}+\bY_{\nu})_{ia},
|
|
\end{equation}
|
|
where $\bX$ and $\bY$ are the components of the eigenvectors of the particle-hole direct RPA problem defined as
|
|
\begin{equation}
|
|
\mqty( \bA & \bB \\ -\bA & -\bB ) \mqty( \bX \\ \bY ) = \mqty( \bX \\ \bY ) \boldsymbol{\Omega},
|
|
\end{equation}
|
|
with
|
|
\begin{subequations}
|
|
\begin{align}
|
|
A_{ia,jb} & = (\epsilon_i - \epsilon_a) \delta_{ij}\delta_{ab} + \eri{ib}{aj}.
|
|
\\
|
|
B_{ia,jb} & = \eri{ij}{ab}.
|
|
\end{align}
|
|
\end{subequations}
|
|
The diagonal matrix $\boldsymbol{\Omega}$ contains the eigenvalues and its elements $\Omega_\nu$ appear in Eq.~\eqref{eq:GW_selfenergy}.
|
|
\titou{The TDA case is discussed in Appendix \ref{sec:nonTDA}.}
|
|
Throughout the manuscript, the indices $p,q,r,s$ are general orbitals while $i,j,k,l$ and $a,b,c,d$ refers to occupied and virtual orbitals, respectively.
|
|
The indices $\mu$ and $\nu$ are composite indices, \eg $\nu=(ia)$, referring to neutral excitations.
|
|
|
|
Because of the frequency dependence, fully solving the quasi-particle equation is a rather complicated task.
|
|
Hence, several approximate schemes have been developed to bypass self-consistency.
|
|
The most popular one is the one-shot (perturbative) scheme, known as $G_0W_0$, where the self-consistency is completely abandoned, and the off-diagonal elements of Eq.~\eqref{eq:quasipart_eq} are neglected.
|
|
\titou{Assuming a HF starting point,} this results in $K$ quasi-particle equations that read
|
|
\begin{equation}
|
|
\label{eq:G0W0}
|
|
\epsilon_p^{\HF} + \Sigma_{pp}(\omega) - \omega = 0,
|
|
\end{equation}
|
|
where $\Sigma_{pp}(\omega)$ are the diagonal elements of $\bSig$ and $\epsilon_p^{\HF}$ are the HF orbital energies.
|
|
The previous equations are non-linear with respect to $\omega$ and therefore can have multiple solutions $\epsilon_{p,s}$ for a given $p$ (where the index $s$ is numbering solutions).
|
|
These solutions can be characterized by their spectral weight given by the renormalisation factor $Z_{p,s}$
|
|
\begin{equation}
|
|
\label{eq:renorm_factor}
|
|
0 \leq Z_{p,s} = \qty[ 1 - \eval{\pdv{\Sigma_{pp}(\omega)}{\omega}}_{\omega=\epsilon_{p,s}} ]^{-1} \leq 1.
|
|
\end{equation}
|
|
The solution with the largest weight is referred to as the quasi-particle while the others are known as satellites (or shake-up transitions).
|
|
However, in some cases, Eq.~\eqref{eq:G0W0} can have two (or more) solutions with similar weights, hence the quasi-particle is not well-defined.
|
|
These additional solutions with large weights are the previously mentioned intruder states.
|
|
|
|
One obvious flaw of the one-shot scheme mentioned above is its starting point dependence.
|
|
Indeed, in Eq.~\eqref{eq:G0W0} we chose to use the HF orbital energies but this is arbitrary and one could have chosen Kohn-Sham orbitals for example.
|
|
Therefore, one can \titou{optimize} the starting point to obtain the best one-shot energies possible, which is commonly done. \cite{Korzdorfer_2012,Marom_2012,Bruneval_2013,Gallandi_2015,Caruso_2016, Gallandi_2016}
|
|
\PFL{Maybe it is worth mentioning here that is is a fairly heuristic approach that is obviously system dependent?}
|
|
Alternatively, one could solve this set of quasi-particle equations self-consistently leading to the eigenvalue-only self-consistent scheme (ev$GW$). \cite{Shishkin_2007a,Blase_2011b,Marom_2012,Kaplan_2016,Wilhelm_2016}
|
|
The solutions $\epsilon_p$ are used to build Eq.~\eqref{eq:G0W0} instead of the HF ones and then these equation are solved for $\omega$ again.
|
|
This procedure is iterated until convergence for $\epsilon_p$ is reached.
|
|
\PFL{This is not quite right. It is probably going to be easier to explain when you're going to introduce the explicit expressions of these quantities.}
|
|
\ANT{Is it better now ?}
|
|
However, if one of the quasi-particle equations does not have a well-defined quasi-particle solution, reaching self-consistency can be quite difficult, if not impossible.
|
|
Even at convergence, the starting point dependence is not totally removed as the results still depend on the initial molecular orbitals. \cite{Marom_2012}
|
|
|
|
In order to update both the orbital energies and coefficients, one must consider the off-diagonal elements in $\bSig(\omega)$.
|
|
To take into account the off-diagonal elements without solving the dynamic quasi-particle equation [Eq.~\eqref{eq:quasipart_eq}], one can resort to the quasi-particle self-consistent (qs) $GW$ scheme in which $\bSig(\omega)$ is replaced by a static approximation $\bSig^{\qs}$.
|
|
Then the qs$GW$ problem is solved using the usual HF algorithm with $\bF$ replaced by $\bF + \bSig^{\qs}$.
|
|
Various choices for $\bSig^\qs$ are possible but the most popular one is the following Hermitian approximation
|
|
\begin{equation}
|
|
\label{eq:sym_qsgw}
|
|
\Sigma_{pq}^\qs = \frac{1}{2}\Re[\Sigma_{pq}(\epsilon_p) + \Sigma_{pq}(\epsilon_q) ].
|
|
\end{equation}
|
|
This form has first been introduced by Faleev and co-workers \cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007} before being derived as the effective Hamiltonian that minimizes the length of the gradient of the Klein functional for non-interacting Green's function. \cite{Ismail-Beigi_2017}
|
|
One of the main results of this manuscript is the derivation from first principles of an alternative static Hermitian form.
|
|
This will be done in the next sections.
|
|
|
|
Once again, in cases where multiple solutions have large spectral weights, qs$GW$ self-consistency can be difficult to reach.
|
|
Multiple solutions of Eq.~\eqref{eq:G0W0} arise due to the $\omega$ dependence of the self-energy.
|
|
Therefore, by suppressing this dependence the static approximation relies on the fact that there is one well-defined quasi-particle solution.
|
|
If it is not the case, the qs$GW$ self-consistent scheme will oscillate between the solutions with large weights. \cite{Forster_2021}
|
|
|
|
The satellites causing convergence problems are the so-called intruder states.
|
|
The intruder state problem can be dealt with by introducing \textit{ad hoc} regularisers.
|
|
The $\ii \eta$ term that is usually added in the denominators of the self-energy [see Eq.~(\ref{eq:GW_selfenergy})] is the usual imaginary-shift regulariser used in various other theories flawed by intruder states. \cite{Battaglia_2022} \ant{more ref...}
|
|
Various other regularisers are possible and in particular one of us has shown that a regulariser inspired by the SRG had some advantages over the imaginary shift. \cite{Monino_2022}
|
|
But it would be more rigorous, and more instructive, to obtain this regulariser from first principles by applying the SRG formalism to many-body perturbation theory.
|
|
This is the aim of the rest of this work.
|
|
|
|
Applying the SRG to $GW$ could gradually remove the coupling between the quasi-particle and the satellites resulting in a renormalized quasi-particle.
|
|
However, to do so one needs to identify the coupling terms in Eq.~\eqref{eq:quasipart_eq}, which is not straightforward.
|
|
The way around this problem is to transform Eq.~\eqref{eq:quasipart_eq} to its upfolded version and the coupling terms will elegantly appear in the process.
|
|
The upfolded $GW$ quasi-particle equation is
|
|
\begin{equation}
|
|
\label{eq:GWlin}
|
|
\begin{pmatrix}
|
|
\bF & \bV^{\text{2h1p}} & \bV^{\text{2p1h}} \\
|
|
\T{(\bV^{\text{2h1p}})} & \bC^{\text{2h1p}} & \bO \\
|
|
\T{(\bV^{\text{2p1h}})} & \bO & \bC^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
\bX \\
|
|
\bY^{\text{2h1p}} \\
|
|
\bY^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
=
|
|
\begin{pmatrix}
|
|
\bX \\
|
|
\bY^{\text{2h1p}} \\
|
|
\bY^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
\boldsymbol{\epsilon},
|
|
\end{equation}
|
|
where $\boldsymbol{\epsilon}$ is a diagonal matrix collecting the quasi-particle and satellite energies, the 2h1p and 2p1h matrix elements are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
C^\text{2h1p}_{ija,klc} & = \qty[ \qty( \epsilon_{i} + \epsilon_{j} - \epsilon_{a} ) \delta_{jl} \delta_{ac} - \eri{jc}{al} ] \delta_{ik} ,
|
|
\\
|
|
C^\text{2p1h}_{iab,kcd} & = \qty[ \qty( \epsilon_{a} + \epsilon_{b} - \epsilon_{i} ) \delta_{ik} \delta_{ac} + \eri{ak}{ic} ] \delta_{bd},
|
|
\end{align}
|
|
\end{subequations}
|
|
and the corresponding coupling blocks read
|
|
\begin{align}
|
|
V^\text{2h1p}_{p,klc} & = \eri{pc}{kl},
|
|
&
|
|
V^\text{2p1h}_{p,kcd} & = \eri{pk}{dc}.
|
|
\end{align}
|
|
The usual $GW$ non-linear equation can be obtained by applying L\"odwin partitioning technique \cite{Lowdin_1963} to Eq.~\eqref{eq:GWlin} which gives the following expression for the self-energy \cite{Bintrim_2021}
|
|
\begin{equation}
|
|
\begin{split}
|
|
\bSig(\omega)
|
|
& = \bV^{\hhp} \qty(\omega \mathbb{1} - \bC^{\hhp})^{-1} \T{(\bV^{\hhp})}
|
|
\\
|
|
& + \bV^{\pph} \qty(\omega \mathbb{1} - \bC^{\pph})^{-1} \T{(\bV^{\pph})},
|
|
\end{split}
|
|
\end{equation}
|
|
which can be further developed to give exactly Eq.~(\ref{eq:GW_selfenergy}).
|
|
|
|
Equations \eqref{eq:GWlin} and \eqref{eq:quasipart_eq} have exactly the same solutions but one is linear and the other not.
|
|
The price to pay for this linearity is that the size of the matrix in the former is $\order{K^3}$ while it is $\order{K}$ in the latter.
|
|
We refer to Ref.~\onlinecite{Bintrim_2021} for a detailed discussion of the up/downfolding processes of the $GW$ equations (see also Ref.~\cite{Tolle_2022}).
|
|
|
|
As can be readily seen in Eq.~\eqref{eq:GWlin}, the blocks $V^\text{2h1p}$ and $ V^\text{2p1h}$ are coupling the 1h and 1p configuration to the dressed 2h1p and 2p1h configurations.
|
|
Therefore, these blocks will be the target of the SRG transformation but before going into more detail we will review the SRG formalism.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{The similarity renormalization group}
|
|
\label{sec:srg}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
The similarity renormalization group method aims at continuously transforming a Hamiltonian to a diagonal form, or more often to a block-diagonal form.
|
|
Therefore, the transformed Hamiltonian
|
|
\begin{equation}
|
|
\label{eq:SRG_Ham}
|
|
\bH(s) = \bU(s) \, \bH \, \bU^\dag(s)
|
|
\end{equation}
|
|
depends on a flow parameter $s$, such that $\bH(s=0)$ is the initial untransformed Hamiltonian and $\bH(s=\infty)$ is the (block-)diagonal Hamiltonian.
|
|
An evolution equation for $\bH(s)$ can be easily obtained by differentiating Eq.~\eqref{eq:SRG_Ham} with respect to $s$, yielding the flow equation
|
|
\begin{equation}
|
|
\label{eq:flowEquation}
|
|
\dv{\bH(s)}{s} = \comm{\boldsymbol{\eta}(s)}{\bH(s)},
|
|
\end{equation}
|
|
where $\boldsymbol{\eta}(s)$, the flow generator, is defined as
|
|
\begin{equation}
|
|
\boldsymbol{\eta}(s) = \dv{\bU(s)}{s} \bU^\dag(s) = - \boldsymbol{\eta}^\dag(s).
|
|
\end{equation}
|
|
To solve this equation at a lower cost than the one of diagonalizing the initial Hamiltonian, one must introduce an approximate form for $\boldsymbol{\eta}(s)$.
|
|
\titou{Before defining such an approximation, we need to define what are the blocks to suppress to obtain a block-diagonal Hamiltonian.
|
|
The Hamiltonian is separated into two parts as
|
|
\begin{equation}
|
|
\bH(s) = \underbrace{\bH^\text{d}(s)}_{\text{diagonal}} + \underbrace{\bH^\text{od}(s)}_{\text{off-diagonal}},
|
|
\end{equation}
|
|
and, by definition, we have the following condition on $\bH^\text{od}$
|
|
\begin{equation}
|
|
\bH^\text{od}(s=\infty) = \boldsymbol{0}.
|
|
\end{equation}}
|
|
\PFL{Move this part at the start of the section.}
|
|
|
|
In this work, we consider Wegner's canonical generator which is defined as \cite{Wegner_1994}
|
|
\begin{equation}
|
|
\boldsymbol{\eta}^\text{W}(s) = \comm{\bH^\text{d}(s)}{\bH(s)} = \comm{\bH^\text{d}(s)}{\bH^\text{od}(s)}.
|
|
\end{equation}
|
|
If this generator is used, the following condition is verified \cite{Kehrein_2006}
|
|
\begin{equation}
|
|
\label{eq:derivative_trace}
|
|
\dv{s}\text{Tr}\left[ \bH^\text{od}(s)^2 \right] \leq 0
|
|
\end{equation}
|
|
which implies that the matrix elements of the off-diagonal part decrease in a monotonic way throughout the transformation.
|
|
Even more, the coupling coefficients associated with the highest-energy determinants are removed first as we shall evidence in the perturbative analysis below.
|
|
The main drawback of this generator is that it generates a stiff set of ODE which is therefore difficult to solve numerically. \ant{ref}
|
|
However, here we will not tackle the full SRG problem but only consider analytical low-order perturbative expressions so we will not be affected by this problem. \cite{Evangelista_2014,Hergert_2016}
|
|
|
|
Let us now perform the perturbative analysis of the SRG equations.
|
|
For $s=0$, the initial problem is
|
|
\begin{equation}
|
|
\bH(0) = \bH^\text{d}(0) + \lambda \bH^\text{od}(0)
|
|
\end{equation}
|
|
where $\lambda$ is the usual perturbation parameter and the off-diagonal part of the Hamiltonian has been defined as the perturber.
|
|
For finite values of $s$, we have the following perturbation expansion of the Hamiltonian
|
|
\begin{equation}
|
|
\label{eq:perturbation_expansionH}
|
|
\bH(s) = \bH^{(0)}(s) + \lambda ~ \bH^{(1)}(s) + \lambda^2 \bH^{(2)}(s) + \cdots
|
|
\end{equation}
|
|
Hence, the generator $\boldsymbol{\eta}(s)$ admits a perturbation expansion as well.
|
|
Then, one can collect order by order the terms in Eq.~\eqref{eq:flowEquation} and solve analytically the low-order differential equations.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Regularized $GW$ approximation}
|
|
\label{sec:srggw}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
Finally, the SRG formalism exposed above will be applied to $GW$.
|
|
The first step is to define the diagonal and off-diagonal parts of the $GW$ effective Hamiltonian.
|
|
As hinted at the end of Sec.~\ref{sec:gw}, the diagonal and off-diagonal parts are defined as
|
|
\begin{align}
|
|
\label{eq:diag_and_offdiag}
|
|
\bH^\text{d}(s) &=
|
|
\begin{pmatrix}
|
|
\bF & \bO & \bO \\
|
|
\bO & \bC^{\text{2h1p}} & \bO \\
|
|
\bO & \bO & \bC^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
& \\
|
|
\bH^\text{od}(s) &=
|
|
\begin{pmatrix}
|
|
\bO & \bV^{\text{2h1p}} & \bV^{\text{2p1h}} \\
|
|
(\bV^{\text{2h1p}})^{\mathrm{T}} & \bO & \bO \\
|
|
(\bV^{\text{2p1h}})^{\mathrm{T}} & \bO & \bO \\
|
|
\end{pmatrix}
|
|
\end{align}
|
|
where we have omitted the $s$ dependence of the matrix elements for the sake of brevity.
|
|
Then, the aim is to solve order by order the flow equation [see Eq.~\eqref{eq:flowEquation}] knowing that the initial conditions are
|
|
\begin{align}
|
|
\bHd{0}(0) &= \begin{pmatrix}
|
|
\bF{}{} & \bO \\
|
|
\bO & \bC{}{}
|
|
\end{pmatrix} &
|
|
\bHod{0}(0) &=\begin{pmatrix}
|
|
\bO & \bO \\
|
|
\bO & \bO
|
|
\end{pmatrix} \notag \\
|
|
\bHd{1}(0) &=\begin{pmatrix}
|
|
\bO & \bO \\
|
|
\bO & \bO
|
|
\end{pmatrix} &
|
|
\bHod{1}(0) &= \begin{pmatrix}
|
|
\bO & \bV{}{} \\
|
|
\bV{}{\dagger} & \bO \notag
|
|
\end{pmatrix} \notag
|
|
\end{align}
|
|
where we have defined the matrices $\bC$ and $\bV$ that collects the 2h1p and 2p1h channels for the sake of conciseness.
|
|
Once the analytical low-order perturbative expansions are known they can be inserted in Eq.~\eqref{eq:GWlin} before downfolding to obtain a renormalized quasi-particle equation.
|
|
In particular, in this manuscript, the focus will be on the second-order renormalized quasi-particle equation.
|
|
|
|
%///////////////////////////%
|
|
\subsection{Zeroth-order matrix elements}
|
|
%///////////////////////////%
|
|
|
|
There is only one zeroth order term in the right-hand side of the flow equation
|
|
\begin{equation}
|
|
\dv{\bH^{(0)}}{s} = \comm{\comm{\bHd{0}}{\bHod{0}}}{\bH^{(0)}},
|
|
\end{equation}
|
|
and performing the block matrix products gives the following system of equations
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\dv{\bF^{(0)}}{s} &= \bO \\
|
|
\dv{\bC^{(0)}}{s} &= \bO \\
|
|
\dv{\bV^{(0),\dagger}}{s} &= 2 \bC^{(0)}\bV^{(0),\dagger}\bF^{(0)} - \bV^{(0),\dagger}(\bF^{(0)})^2 - (\bC^{(0)})^2\bV^{(0),\dagger} \\
|
|
\dv{\bV^{(0)}}{s} &= 2 \bF^{(0)}\bV^{(0)}\bC^{(0)} - (\bF^{(0)})^2\bV^{(0)} - \bV^{(0)}(\bC^{(0)})^2
|
|
\end{align}
|
|
\end{subequations}
|
|
where the $s$ dependence of $\bV^{(0)}$ and $\bV^{(0),\dagger}$ has been dropped in the last two equations.
|
|
$\bF^{(0)}$ and $\bC^{(0)}$ do not depend on $s$ as a consequence of the first two equations.
|
|
The last equation can be solved by introducing $\bU$ the matrix that diagonalizes $\bC^{(0)} = \bU \bD^{(0)} \bU^{-1}$ such that the differential equation for $\bV^{(0)}$ becomes
|
|
\begin{equation}
|
|
\label{eq:eqdiffW0}
|
|
\dv{\bW^{(0)}}{s} = 2 \bF^{(0)}\bW^{(0)} \bD^{(0)} - (\bF^{(0)})^2\bW^{(0)} - \bW^{(0)} (\bD^{(0)})^2
|
|
\end{equation}
|
|
where $\bW^{(0)}(s)= \bV^{(0)}(s) \bU$.
|
|
The matrix elements of $\bU$ and $\bD^{(0)}$ are
|
|
\begin{align}
|
|
U_{(p,v),(q,w)} &= \delta_{pq} \bX_{v,w} \\
|
|
D_{(p,v),(q,w)}^{(0)} &= \left(\epsilon_p + \text{sign}(\epsilon_p-\epsilon_\text{F})\Omega_v\right)\delta_{pq}\delta_{vw}
|
|
\end{align}
|
|
where $\epsilon_\text{F}$ is the Fermi level.
|
|
Note that the matrix $\bU$ is also used in the downfolding process of Eq.~\eqref{eq:GWlin}. \cite{Bintrim_2021}
|
|
|
|
Thanks to the diagonal structure of $\bF^{(0)}$ and $\bD^{(0)}$, Eq.~\eqref{eq:eqdiffW0} can be easily solved and give
|
|
\begin{equation}
|
|
W_{p,(q,v)}^{(0)}(s) = W_{p,(q,v)}^{(0)}(0) e^{- (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 s}
|
|
\end{equation}
|
|
Due to the initial conditions $\bV^{(0)}(0) = \bO$, we have $\bW^{(0)}(s)=\bO$ and therefore $\bV^{(0)}(s)=\bO=\bV^{(0)}(0) $.
|
|
Therefore, the zeroth order Hamiltonian is
|
|
\begin{equation}
|
|
\bH^{(0)}(s) = \bH^{(0)}(0),
|
|
\end{equation}
|
|
\ie it is independent of $s$.
|
|
|
|
%///////////////////////////%
|
|
\subsection{First-order matrix elements}
|
|
%///////////////////////////%
|
|
|
|
Knowing that $\bHod{0}(s)=\bO$, the first order flow equation is
|
|
\begin{equation}
|
|
\dv{\bH^{(1)}}{s} = \comm{\comm{\bHd{0}}{\bHod{1}}}{\bHd{0}}
|
|
\end{equation}
|
|
which gives the same system of equations as in the previous subsection except that $\bV^{(0)}$ and $\bV^{(0),\dagger}$ should be replaced by $\bV^{(1)}$ and $\bV^{(1),\dagger}$.
|
|
|
|
Once again the two first equations are easily solved
|
|
\begin{align}
|
|
\bF^{(1)}(s) &= \bF^{(1)}(0) = \bO & \bC^{(1)}(s) &= \bC^{(1)}(0) = \bO.
|
|
\end{align}
|
|
and the first order coupling elements are given by (up to a multiplication by $\bU^{-1}$)
|
|
\begin{align}
|
|
W_{p,(q,v)}^{(1)}(s) &= W_{p,(q,v)}^{(1)}(0) e^{- (F_{pp}^{(0)} - D_{(q,v),(q,v)}^{(0)})^2 s} \\
|
|
&= W_{p,(q,v)}^{(1)}(0) e^{- (\epsilon_p - \epsilon_q - \text{sign}(\epsilon_q-\epsilon_F)\Omega_v)^2 s} \notag
|
|
\end{align}
|
|
At $s=0$ the elements $W_{p,(q,v)}^{(1)}(0)$ are equal to the two-electron screened integrals defined in Eq.~\eqref{eq:GW_sERI} while for $s\to\infty$ they go to zero.
|
|
Therefore, $W_{p,(q,v)}^{(1)}(s)$ are renormalized two-electrons screened integrals.
|
|
Note the close similarity of the first-order element expressions with the ones of Evangelista in Ref.~\onlinecite{Evangelista_2014b} obtained in a second quantization formalism (see also Ref.~\onlinecite{Hergert_2016}).
|
|
|
|
%///////////////////////////%
|
|
\subsection{Second-order matrix elements}
|
|
% ///////////////////////////%
|
|
|
|
The second-order renormalized quasi-particle equation is given by
|
|
\begin{equation}
|
|
\label{eq:GW_renorm}
|
|
\left( \widetilde{\bF}(s) + \widetilde{\bSig}(\omega; s) \right) \bX = \omega \bX
|
|
\end{equation}
|
|
with
|
|
\begin{align}
|
|
\widetilde{\bF}(s) &= \bF^{(0)}+\bF^{(2)}(s)\\
|
|
\label{eq:srg_sigma}
|
|
\widetilde{\bSig}(\omega; s) &= \bV^{(1)}(s) \left(\omega \mathbb{1} - \bC^{(0)}\right)^{-1} (\bV^{(1)}(s))^{\dagger}
|
|
\end{align}
|
|
|
|
As can be readily seen above, $\bF^{(2)}$ is the only second-order block of the effective Hamiltonian contributing to the second-order SRG quasi-particle equation.
|
|
Collecting every second-order terms in the flow equation and performing the block matrix products results in the following differential equation for $\bF^{(2)}$
|
|
\begin{equation}
|
|
\label{eq:diffeqF2}
|
|
\dv{\bF^{(2)}}{s} = \bF^{(0)}\bV^{(1)}\bV^{(1),\dagger} + \bV^{(1)}\bV^{(1),\dagger}\bF^{(0)} - 2 \bV^{(1)}\bC^{(0)}\bV^{(1),\dagger} .
|
|
\end{equation}
|
|
This can be solved by simple integration along with the initial condition $\bF^{(2)}=\bO$ to give
|
|
\begin{align}
|
|
&F_{pq}^{(2)}(s) = \\
|
|
&\sum_{r,v} \frac{\Delta_{prv}+ \Delta_{qrv}}{\Delta_{prv}^2 + \Delta_{qrv}^2} W_{p,(r,v)} W^{\dagger}_{(r,v),q}\left(1 - e^{-(\Delta_{prv}^2 + \Delta_{qrv}^2) s}\right), \notag
|
|
\end{align}
|
|
with $\Delta_{prv} = \epsilon_p - \epsilon_r - \text{sign}(\epsilon_r-\epsilon_F)\Omega_v$.
|
|
|
|
At $s=0$, this second-order correction is null while for $s\to\infty$ it tends towards the following static limit
|
|
\begin{equation}
|
|
\label{eq:static_F2}
|
|
F_{pq}^{(2)}(\infty) = \sum_{r,v} \frac{\Delta_{prv}+ \Delta_{qrv}}{\Delta_{prv}^2 + \Delta_{qrv}^2} W_{p,(r,v)} W^{\dagger}_{(r,v),q}.
|
|
\end{equation}
|
|
Note that in the $s\to\infty$ limit the dynamic part of the self-energy [see Eq.~\eqref{eq:srg_sigma}] tends to zero.
|
|
Therefore, the SRG flow transforms the dynamic part of $\bSig(\omega)$ into a static correction.
|
|
This transformation is done gradually starting from the states that have the largest denominators in Eq.~\eqref{eq:static_F2}.
|
|
|
|
Interestingly, the static limit, \ie $s\to\infty$ limit, of Eq.~\eqref{eq:GW_renorm} defines an alternative qs$GW$ approximation to the one defined by Eq.~\eqref{eq:sym_qsgw} which matrix elements read as
|
|
\begin{align}
|
|
\label{eq:sym_qsGW}
|
|
&\Sigma_{pq}^{\text{qs}GW}(\eta) = \\
|
|
&\frac{1}{2} \sum_{r,v} \left( \frac{\Delta_{prv}}{\Delta_{prv}^2 + \eta^2} +\frac{\Delta_{qrv}}{\Delta_{qrv}^2 + \eta^2} \right) W_{p,(r,v)} W^{\dagger}_{(r,v),q}. \notag
|
|
\end{align}
|
|
Yet, both approximation are closely related as they share the same diagonal terms when $\eta=0$.
|
|
Also, note that the SRG static approximation is naturally Hermitian as opposed to the symmetrized case where it is enforced by symmetrization.
|
|
|
|
However, as will be discussed in more detail in the results section, the convergence of the qs$GW$ scheme using $\widetilde{\bF}(\infty)$ is very poor.
|
|
This is similar to the symmetric case when the imaginary shift $\ii \eta$ is set to zero.
|
|
Indeed, in qs$GW$ calculation using the symmetrized static form, increasing $\eta$ to ensure convergence in difficult cases is most often unavoidable.
|
|
Therefore, we will define the SRG-qs$GW$ static effective Hamiltonian as
|
|
\begin{align}
|
|
\label{eq:SRG_qsGW}
|
|
&\Sigma_{pq}^{\text{SRG}}(s) = \\
|
|
&\sum_{r,v} \frac{\Delta_{prv}+ \Delta_{qrv}}{\Delta_{prv}^2 + \Delta_{qrv}^2} W_{p,(r,v)} W_{q,(r,v)}\left(1 - e^{-(\Delta_{prv}^2 + \Delta_{qrv}^2) s}\right) \notag
|
|
\end{align}
|
|
which depends on one regularising parameter $s$ analogously to $\eta$ in the usual case.
|
|
The fact that the $s\to\infty$ static limit does not always converge when used in a qs$GW$ calculation could have been predicted because in this limit even the intruder states have been included in $\tilde{\bF}$.
|
|
Therefore, we should use a value of $s$ large enough to include almost every states but small enough to avoid intruder states.
|
|
|
|
To conclude this section, we will discuss the case of discontinuities.
|
|
Indeed, previously we mentioned that intruder states are responsible for both the poor convergence of qs$GW$ and discontinuities in physical quantities at the $\GOWO$ level.
|
|
So is it possible to use the SRG machinery developed above to remove discontinuities?
|
|
In fact, not directly because discontinuities are due to intruder states in the dynamic part while we have seen just above that a finite value of $s$ is well-designed to avoid the intruder states in the static part.
|
|
However, doing a change of variable such that
|
|
\begin{align}
|
|
e^{-s} &= 1-e^{-t} & 1 - e^{-s} &= e^{-t}
|
|
\end{align}
|
|
hence choosing a finite value of $t$ is well-designed to avoid discontinuities in the dynamic.
|
|
In fact, the dynamic part after the change of variable is closely related to the SRG-inspired regulariser introduced by Monino and Loos in Ref.~\onlinecite{Monino_2022}.
|
|
|
|
%=================================================================%
|
|
\section{Computational details}
|
|
\label{sec:comp_det}
|
|
%=================================================================%
|
|
|
|
The two qs$GW$ variants considered in this work have been implemented in an in-house program.
|
|
The $GW$ implementation closely follows the one of mol$GW$. \cite{Bruneval_2016}
|
|
The geometries have been optimized at the CC3 level in the aug-cc-pVTZ basis set without frozen core using the CFOUR program.
|
|
The reference CCSD(T) IP energies have been obtained using default parameters of Gaussian 16.
|
|
This means that the cations used an unrestricted HF reference while the neutral ground-state energies have been obtained in a restricted formalism.
|
|
|
|
%=================================================================%
|
|
\section{Results}
|
|
\label{sec:results}
|
|
%=================================================================%
|
|
|
|
The results section is divided in two parts.
|
|
The first step will be to analyse in depth the behavior of the two static self-energy approximations in a few test cases.
|
|
Then the accuracy of the IP yielded by the Sym and SRG schemes will be statistically gauged over a test set of molecules.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\subsection{Flow parameter dependence of the SRG/Sym-qs$GW$/GF(2) schemes}
|
|
\label{sec:flow_param_dep}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%% FIG 1 %%%
|
|
\begin{figure*}
|
|
\centering
|
|
\includegraphics[width=\linewidth]{fig1.pdf}
|
|
\caption{
|
|
Add caption \ANT{Should we add $G_0W_0$?}
|
|
\label{fig:fig1}}
|
|
\end{figure*}
|
|
%%% %%% %%% %%%
|
|
|
|
This section starts by considering the Neon atom and the water molecule in the aug-cc-pVTZ cartesian basis set in Fig.~\ref{fig:fig1}.
|
|
Figure~\ref{fig:fig1} shows the error of various methods for the principal IP with respect to CCSD(T).
|
|
The HF IPs (cyan lines) are overestimated, this is a consequence of the missing correlation, a result which is now well-understood. \cite{Lewis_2019} \ANT{I should maybe search for another ref as well.}
|
|
The usual Sym-qs$GW$ scheme (green lines) brings a quantitative improvement as both IP energies are now within \SI{0.5}{\electronvolt} of the reference.
|
|
The Neon atom is a well-behaved system and could be converged without regularisation parameter while for water $\eta$ was set to 0.01 to help convergence.
|
|
|
|
Figure~\ref{fig:fig1} also displays the SRG-qs$GW$ IPs as a function of the flow parameter (blue curves).
|
|
At $s=0$, the IPs are equal to their HF counterparts as expected from the discussion of Sec.~\ref{sec:srggw}.
|
|
For $s\to\infty$ both IPs reach a plateau at an error that is significantly smaller than their $s=0$ starting point.
|
|
Even more, the values associated with these plateau are more accurate than their Sym-qs$GW$ counterparts.
|
|
However, the SRG-qs$GW$ error do not decrease smoothly between the HF values and their limits as for small $s$ values they are actually worst than the HF IPs.
|
|
|
|
In addition, we also considered the qs$GW$ and SRG-qs$GW$ methods based on a TDA screening.
|
|
The TDA IPs are now understimated unlike their RPA counterparts.
|
|
The SRG-qs$GW$ absolute error for the IPs clearly deteriorates when going from RPA to TDA.
|
|
On the other hand, for qs$GW$ the TDA-based IP is better than the RPA one for Neon while it is the other way around for water.
|
|
This trend will be investigated in more details in the next subsection.
|
|
|
|
\ANT{Maybe we should add GF(2) because it allows us to explain the behavior of the SRG curve using perturbation theory.}
|
|
The behavior of the SRG-qsGF2 IPS is similar to the SRG-qs$GW$ one.
|
|
Add sentence about $GW$ better than GF2 when the results will be here.
|
|
The decrease and then increase behavior of the IPs can be rationalised using results from perturbation theory for GF(2).
|
|
We refer the reader to the chapter 8 of Ref.~\onlinecite{Schirmer_2018} for more details about this analysis.
|
|
The GF(2) IP admits the following perturbation expansion...
|
|
|
|
Because $GW$ relies on an infinite resummation of diagram such a perturbation analysis is difficult to make in this case.
|
|
But the mechanism causing the increase/decrease of the $GW$ IPs as a function of $s$ should be closely related to the GF(2) one exposed above.
|
|
|
|
\begin{itemize}
|
|
\item Li2 interesting because HF underestimate IP
|
|
\item LiH interesting because Sym-qs$GW$ worst than HF
|
|
\item BeO interesting because usually difficult to converge in qs$GW$ (cf Forster 2021)
|
|
\end{itemize}
|
|
|
|
%%% FIG 2 %%%
|
|
\begin{figure*}
|
|
\centering
|
|
\includegraphics[width=\linewidth]{fig2.pdf}
|
|
\caption{
|
|
Add caption
|
|
\label{fig:fig2}}
|
|
\end{figure*}
|
|
%%% %%% %%% %%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\subsection{Comparison of the Sym-qs and SRG-qs approximations}
|
|
\label{sec:SRG_vs_Sym}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
The test set considered in this study is composed of the GW20 set of molecules introduced by Lewis and Berkelbach. \cite{Lewis_2019}
|
|
This set is made of the 20 smallest atoms and molecules of the GW100 benchmark set.
|
|
We also added the MgO and O3 molecules which are part of GW100 and are known to be difficult to converged for qs$GW$.
|
|
In addition, we considered the Quest 1 and 2 sets which is made of small and medium size organic molecules.
|
|
|
|
|
|
%=================================================================%
|
|
\section{Conclusion}
|
|
\label{sec:conclusion}
|
|
%=================================================================%
|
|
|
|
Here comes the conclusion.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\acknowledgements{This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
|
|
The authors thank Francesco Evangelista for inspiring discussions.
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section*{Data availability statement}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
The data that supports the findings of this study are available within the article.% and its supplementary material.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\bibliography{SRGGW}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
\appendix
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Non-TDA $GW$ and $GW_{\text{TDHF}}$ equations}
|
|
\label{sec:nonTDA}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
The matrix elements of the $GW$ self-energy without TDA are the same as in Eq.~\eqref{eq:GW_selfenergy} but the screened integrals are now defined as
|
|
\begin{equation}
|
|
\label{eq:GWnonTDA_sERI}
|
|
W_{p,(q,v)} = \sum_{ia}\eri{pi}{qa}\qty( \bX_{v} + \bY_{v})_{ia},
|
|
\end{equation}
|
|
where $\bX$ and $\bY$ are the eigenvector matrices of the full particle-hole dRPA problem defined as
|
|
\begin{equation}
|
|
\label{eq:full_dRPA}
|
|
\begin{pmatrix}
|
|
\bA & \bB \\
|
|
- \bB & \bA \\
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
\bX \\
|
|
\bY \\
|
|
\end{pmatrix} = \boldsymbol{\Omega}
|
|
\begin{pmatrix}
|
|
\bX \\
|
|
\bY \\
|
|
\end{pmatrix},
|
|
\end{equation}
|
|
with
|
|
\begin{align}
|
|
A_{ij,ab} &= (\epsilon_i - \epsilon_a) \delta_{ij}\delta_{ab} + \eri{ib}{aj}, \\
|
|
B_{ij,ab} &= \eri{ij}{ab}.
|
|
\end{align}
|
|
$\boldsymbol{\Omega}$ is the diagonal matrix of eigenvalues. Note that $\boldsymbol{\Omega}$ in this case has the same size as in the TDA because we consider only the positive excitations of the full dRPA problem.
|
|
|
|
Defining an unfold version of this equation that does not require a diagonalization of the RPA problem before unfolding is a tricky task (see supplementary material of Ref.~\onlinecite{Bintrim_2021}).
|
|
However, because we will eventually downfold again the upfolded matrix, we can use the following matrix \cite{Tolle_2022}
|
|
\begin{equation}
|
|
\label{eq:nonTDA_upfold}
|
|
\begin{pmatrix}
|
|
\bF & \bW^{\text{2h1p}} & \bW^{\text{2p1h}} \\
|
|
(\bW^{\text{2h1p}})^{\mathrm{T}} & \bD^{\text{2h1p}} & \bO \\
|
|
(\bW^{\text{2p1h}})^{\mathrm{T}} & \bO & \bD^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX \\
|
|
\bY^{\text{2h1p}} \\
|
|
\bY^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
=
|
|
\begin{pmatrix}
|
|
\bX \\
|
|
\bY^{\text{2h1p}} \\
|
|
\bY^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\boldsymbol{\epsilon},
|
|
\end{equation}
|
|
which already depends on the screened integrals and therefore require the knowledge of the eigenvectors of the dRPA problem defined in Eq.~\eqref{eq:full_dRPA}.
|
|
|
|
|
|
where $\boldsymbol{\epsilon}$ is a diagonal matrix collecting the quasi-particle and satellite energies, the 2h1p and 2p1h matrix elements are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
D^\text{2h1p}_{ija,klc} & = \delta_{ik}\delta_{jl} \delta_{ac} \qty[\epsilon_i - \Omega_{ja}] ,
|
|
\\
|
|
D^\text{2p1h}_{iab,kcd} & = \delta_{ik}\delta_{ac} \delta_{bd} \qty[\epsilon_a + \Omega_{ib}] ,
|
|
\end{align}
|
|
\end{subequations}
|
|
and the corresponding coupling blocks read
|
|
\begin{align}
|
|
W^\text{2h1p}_{p,klc} & = \sum_{ia}\eri{pi}{ka} \qty( \bX_{lc} + \bY_{lc})_{ia} \\
|
|
W^\text{2p1h}_{p,kcd} & = \sum_{ia}\eri{pi}{ca} \qty( \bX_{kd} + \bY_{kd})_{ia}
|
|
\end{align}
|
|
|
|
Using the SRG on this matrix instead of Eq.~\eqref{eq:GWlin} gives the same expression for $\bW^{(1)}$, $\bF^{(2)}$ and $\bSig^{\text{SRG}}$ but now the screened integrals are the one of Eq.~\eqref{eq:GWnonTDA_sERI} and the eigenvalues $\Omega$ and eigenvectors $\bX$ and $\bY$ are the ones of the full RPA problem defined in Eq.~\eqref{eq:full_dRPA}.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{GF(2) equations \ant{NOT SURE THAT WE KEEP IT}}
|
|
\label{sec:GF2}
|
|
%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
The GF($n$) formalism is defined such that the self-energy includes every diagram up to $n$-th order of M\"oller-Plesset perturbation theory.
|
|
The matrix elements of its second-order version read as
|
|
\begin{align}
|
|
\label{eq:GF2_selfenergy}
|
|
\Sigma_{pq}^{\text{GF(2)}}(\omega)
|
|
&= \frac{1}{\sqrt{2}} \sum_{ija} \frac{\aeri{pa}{ij}\aeri{qa}{ij}}{\omega + \epsilon _a -\epsilon_i -\epsilon_j - \ii \eta} \\
|
|
&+ \frac{1}{\sqrt{2}} \sum_{iab} \frac{\aeri{pi}{ab}\aeri{qi}{ab}}{\omega + \epsilon _i -\epsilon_a -\epsilon_b + \ii \eta}
|
|
\end{align}
|
|
This self-energy can be upfolded similarly to the $GW$ case and one obtain the following ``super-matrix''
|
|
\begin{equation}
|
|
\label{eq:unfolded_matrice}
|
|
\bH =
|
|
\begin{pmatrix}
|
|
\bF & \bV^{\text{2h1p}} & \bV^{\text{2p1h}} \\
|
|
(\bV^{\text{2h1p}})^{\mathsf{T}} & \bC^{\text{2h1p}} & \bO \\
|
|
(\bV^{\text{2p1h}})^{\mathsf{T}} & \bO & \bC^{\text{2p1h}} \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
The expression of the coupling blocks $\bV{}{}$ and the diagonal blocks $\bC{}{}$ is given below.
|
|
\begin{align}
|
|
\label{eq:GF2_unfolded}
|
|
V^\text{2h1p}_{p,ija} & = \frac{1}{\sqrt{2}}\aeri{pa}{ij}
|
|
\\
|
|
V^\text{2p1h}_{p,iab} & = \frac{1}{\sqrt{2}}\aeri{pi}{ab}
|
|
\\
|
|
C^\text{2h1p}_{ija,klc} & = \qty( \epsilon_i + \epsilon_j - \epsilon_a) \delta_{jl} \delta_{ac} \delta_{ik}
|
|
\\
|
|
C^\text{2p1h}_{iab,kcd} & = \qty( \epsilon_a + \epsilon_b - \epsilon_i) \delta_{ik} \delta_{ac} \delta_{bd}
|
|
\end{align}
|
|
Note that this matrix is exactly the ADC(2) matrix for charged excitations.
|
|
The fact that the integrals are not screened in GF(2) manifests itself in the fact that the $\bC$ matrices are already diagonal.
|
|
|
|
Applying the SRG formalism to this matrix is completely analog to the derivation exposed in the main text.
|
|
We only give the analytical expressions of the matrix elements needed for the second-order SRG-GF(2) quasiparticle equations.
|
|
|
|
\begin{equation}
|
|
(V^\text{2h1p}_{p,ija})^{(1)}(s) = \frac{1}{\sqrt{2}}\aeri{pa}{ij} e^{- (\epsilon_p + \epsilon_a - \epsilon_i - \epsilon_j)^2 s}
|
|
\end{equation}
|
|
\begin{equation}
|
|
(V^\text{2h1p}_{p,iab})^{(1)}(s) = \frac{1}{\sqrt{2}}\aeri{pi}{ab} e^{- (\epsilon_p + \epsilon_i - \epsilon_a - \epsilon_b)^2 s}
|
|
\end{equation}
|
|
|
|
We define $ \Delta_{pq,rs} = \epsilon_p + \epsilon_q - \epsilon_r - \epsilon_s $
|
|
|
|
\begin{align}
|
|
F_{pq}^{(2)}(s) &= \sum_{ria} \frac{\epsilon_{p} + \epsilon_{q} - 2 (\epsilon_r \pm \Omega_v)}{(\epsilon_p - \epsilon_r \pm \Omega_v)^2 + (\epsilon_q - \epsilon_r \pm \Omega_v)^2} W_{p,(r,v)} \notag \\
|
|
&\times W^{\dagger}_{(r,v),q}\left(1 - e^{-(\epsilon_p - \epsilon_r \pm \Omega_v)^2s} e^{-(\epsilon_q - \epsilon_r \pm \Omega_v)^2s}\right).
|
|
\end{align}
|
|
|
|
|
|
\end{document}
|