saving work

Manuscript/ufGW.tex

@@ -49,7 +49,8 @@
 \newcommand{\RPA}{\text{RPA}}
 
 % 
-\newcommand{\Norb}{N}
+\newcommand{\Ne}{N}
+\newcommand{\Norb}{K}
 \newcommand{\Nocc}{O}
 \newcommand{\Nvir}{V}
 
@@ -123,8 +124,8 @@
 	\affiliation{\LCPQ}
 
 \begin{abstract}
-By upfolding the frequency-dependent $GW$ quasiparticle equation, we explain the appearance of multiple solutions and unphysical discontinuities in various physical quantities computed within the $GW$ approximation.
+By recasting the non-linear frequency-dependent $GW$ quasiparticle equation into a linear eigenvalue problem, we explain the appearance of multiple solutions and unphysical discontinuities in various physical quantities computed within the $GW$ approximation.
-By considering the $GW$ self-energy as an effective Hamiltonian, the appearance of these multiple solutions and discontinuities can be directly related to the intruder state problem.
+Considering the $GW$ self-energy as an effective Hamiltonian, these issues can be directly related to the intruder state problem.
 A simple and efficient regularization procedure is proposed to avoid such issues.
 %\bigskip
 %\begin{center}
@@ -141,7 +142,7 @@ A simple and efficient regularization procedure is proposed to avoid such issues
 The $GW$ approximation of many-body perturbation theory \cite{Hedin_1965,Martin_2016} allows to compute accurate charged excitation (\ie, ionization potentials and electron affinities) in solids and molecules. \cite{Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019}
 Its popularity in the molecular structure community is rapidly growing \cite{Ke_2011,Bruneval_2012,Bruneval_2013,Bruneval_2015,Blase_2016,Bruneval_2016, Bruneval_2016a,Koval_2014,Hung_2016,Blase_2018,Boulanger_2014,Li_2017,Hung_2016,Hung_2017,vanSetten_2015,vanSetten_2018,vanSetten_2015, Maggio_2017,vanSetten_2018,Richard_2016,Gallandi_2016,Knight_2016,Dolgounitcheva_2016,Bruneval_2015,Krause_2015,Govoni_2018,Caruso_2016} thanks to its relatively low cost and somehow surprising accuracy for weakly-correlated systems. \cite{Bruneval_2021}
 
-The idea behind the $GW$ approximation (or Green's function-based methods in general) is to recast the many-body problem into a set of non-linear one-body equations. The introduction of the self-energy $\Sigma$ links the non-interacting Green's function $G_0$ to its fully-interacting version $G$ via the following Dyson equation:
+The idea behind the $GW$ approximation is to recast the many-body problem into a set of non-linear one-body equations. The introduction of the self-energy $\Sigma$ links the non-interacting Green's function $G_0$ to its fully-interacting version $G$ via the following Dyson equation:
 \begin{equation}
 	G = G_0 + G_0 \Sigma G
 \end{equation}
@@ -177,10 +178,11 @@ It was shown that these problems could be alleviated by using a static Coulomb-h
 However, none of these solutions is completely satisfying as a static approximation of the self-energy can induce significant loss in accuracy and fully self-consistent calculations can be quite challenging in terms of implementation and cost.
 
 In the present article, via an upfolding process of the non-linear $GW$ equation, \cite{Bintrim_2021a} we provide further physical insights into the origin of these discontinuities by highlighting, in particular, the role of intruder states.
-Inspired by regularized electronic structure theories, \cite{Lee_2018a,Evangelista_2014b} these new insights allow us to provide a cheap and efficient regularization scheme in order to avoid these issues.
+Inspired by regularized electronic structure theories, \cite{Lee_2018a,Evangelista_2014b} these new insights allow us to propose a cheap and efficient regularization scheme in order to avoid these issues.
 
 Here, we consider the one-shot {\GOWO} \cite{Strinati_1980,Hybertsen_1985a,Hybertsen_1986,Godby_1988,Linden_1988,Northrup_1991,Blase_1994,Rohlfing_1995,Shishkin_2007} for the sake of simplicity but the same analysis can be performed in the case of (partially) self-consistent schemes.\cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016,Gui_2018,Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016}
 Moreover, we consider a restricted Hartree-Fock (HF) starting point but it can be straightforwardly extended to a Kohn-Sham (KS) starting point.
+Throughout this article, $p$ and $q$ are general (spatial) orbitals, $i$, $j$, $k$, and $l$ denotes occupied orbitals, $a$, $b$, $c$, and $d$ are (unoccupied) virtual orbitals, while $m$ labels single excitations $i \to a$.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Downfold: The non-linear $GW$ problem}
@@ -191,17 +193,18 @@ Within the {\GOWO} approximation, in order to obtain the quasiparticle energies
 \label{eq:qp_eq}
 	\eps{p}{\HF} + \SigC{p}(\omega) - \omega = 0
 \end{equation}
-where $\eps{p}{\HF}$ is the $p$th HF orbital energy and the correlation part of the {\GOWO} self-energy reads
+where $\eps{p}{\HF}$ is the $p$th HF orbital energy and the correlation part of the {\GOWO} self-energy
 \begin{equation}
 	\SigC{p}(\omega) 
 	= \sum_{im} \frac{2\ERI{pi}{m}^2}{\omega - \eps{i}{\HF} + \Om{m}{\RPA}}
 	+ \sum_{am} \frac{2\ERI{pa}{m}^2}{\omega - \eps{a}{\HF} - \Om{m}{\RPA}}
 \end{equation}
+is constituted by a hole (h) and a particle (p) term.
 Within the Tamm-Dancoff approximation, the screened two-electron integrals are given by
 \begin{equation}
 	\ERI{pq}{m} = \sum_{ia} \ERI{pq}{ia} X_{ia,m}^\RPA
 \end{equation}
-where $\Om{m}{\RPA}$ and $\bX{m}{\RPA}$ are respectively the $m$th eigenvalue and eigenvector of the RPA problem in the Tamm-Dancoff approximation, \ie, 
+where $\Om{m}{\RPA}$ and $\bX{m}{\RPA}$ are respectively the $m$th eigenvalue and eigenvector of the random-phase approximation (RPA) problem in the Tamm-Dancoff approximation, \ie, 
 \begin{equation}
 	\bA{}{\RPA} \cdot \bX{m}{\RPA} = \Om{m}{\RPA} \bX{m}{\RPA}
 \end{equation}
@@ -211,10 +214,10 @@ with
 \end{equation}
 and 
 \begin{equation}
-	\ERI{pq}{rs} = \iint \MO{p}(\br_1) \MO{q}(\br_1) \frac{1}{\abs{\br_1 - \br_2}} \MO{r}(\br_2) \MO{s}(\br_2) d\br_1 \dbr_2
+	\ERI{pq}{ia} = \iint \MO{p}(\br_1) \MO{q}(\br_1) \frac{1}{\abs{\br_1 - \br_2}} \MO{i}(\br_2) \MO{a}(\br_2) d\br_1 \dbr_2
 \end{equation}
 are two-electron integrals over the HF (spatial) orbitals $\MO{p}(\br)$. 
-As a non-linear equation, Eq.~\eqref{eq:qp_eq} has many solutions $\eps{p,s}{\GW}$ and their corresponding weight is given by the value of the so-called renormalization factor
+As a non-linear equation, Eq.~\eqref{eq:qp_eq} has many solutions $\eps{p,s}{\GW}$ and their corresponding weight are given by the value of the following renormalization factor
 \begin{equation}
 	0 \le Z_{p,s} = \qty[ 1 - \eval{\pdv{\SigC{p}(\omega)}{\omega}}_{\omega = \eps{p,s}{\GW}} ]^{-1} \le 1
 \end{equation} 
@@ -259,8 +262,8 @@ and the corresponding coupling blocks read
 	&
 	V^\text{2p1h}_{p,kcd} & = \sqrt{2} \ERI{pd}{kc}
 \end{align}
-The size of this eigenvalue problem is $N = 1 + N^\text{2h1p} + N^\text{2p1h} = 1 + O^2 V + O V^2$, and this eigenvalue problem has to be solved for each orbital that one wants to correct.
+The size of this eigenvalue problem is $1 + O^2 V + O V^2$ (where $O$ and $V$ are the number of occupied and virtual orbitals, respectively), and this eigenvalue problem has to be solved for each orbital that one wishes to correct.
-Note, however, that the block $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$ do not need to be recomputed for each orbital.
+Note, however, that the blocks $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$ do not need to be recomputed for each orbital.
 Because the renormalization factor corresponds to the projection of the vector $\bc{}{(p,s)}$ onto the reference space, the weight of a solution $(p,s)$ is given by the  the first coefficient of their corresponding eigenvector $\bc{}{(p,s)}$, \ie, 
 \begin{equation}
 	Z_{p,s} = \qty[ c_{1}^{(p,s)} ]^{2}
@@ -279,10 +282,14 @@ and comparing it with Eq.~\eqref{eq:qp_eq} by setting
 	+ \bV{p}{\text{2p1h}} \cdot \qty(\omega \bI - \bC{}{\text{2p1h}} )^{-1} \cdot \T{\qty(\bV{p}{\text{2p1h}})} 
 \end{multline}
 where $\bI$ is the identity matrix.
-The main difference between the two approaches is that, by diagonalizing Eq.~\eqref{eq:Hp}, one has directly access to the eigenvectors associated with each quasiparticle and satellites.
+The main mathematical difference between the two approaches is that, by diagonalizing Eq.~\eqref{eq:Hp}, one has directly access to the eigenvectors associated with each quasiparticle and satellites.
 
 One can see this downfolding process as the construction of a frequency-dependent effective Hamiltonian where the reference (zeroth-order) space is composed by a single determinant of the 1h or 1p sector and the external (first-order) space by all the 2h1p and 2p1h configurations.  
 
+The element $\eps{p}{\HF}$ of $\bH^{(p)}$ [see Eq.~\eqref{eq:Hp}] corresponds to the relative energy of the $(\Ne\pm1)$-electron reference determinant (compared to the $\Ne$-electron HF determinant) while the diagonal elements of the blocks $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$ provide an estimate of the relative energy of the 2h1p and 2p1h determinants.
+In some situations, one of these determinant from the outer space may become of similar energy than the reference determinant, a situation that one could label as intruder state problem.
+Hence, the two diabatic electronic configurations may cross and form an avoided crossing.
+As we shall see below, this is when discontinuities occur and is ubiquitous in molecular systems.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{An illustrative example}
@@ -293,6 +300,9 @@ Therefore, such issues are ubiquitous when one wants to compute core ionized sta
 In order to illustrate the appearance and the origin of these multiple solutions, we consider the hydrogen molecule in the 6-31G basis set which corresponds to a system with 2 electrons and 4 spatial orbitals (one occupied and three virtual).
 This example was already considered in our previous work \cite{Veril_2018} but here we provide further insights on the origin of the appearances of these multiple solutions.
 
+The algorithm diabatically follows the solution, while it should be adiabatic.
+This is not the $\Ne$-electron situation where one has to check if it is multireference, but for the $\Ne\pm1$-electron situations.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %	FIGURE 1
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -314,23 +324,23 @@ As one can see there are two problematic regions showing obvious discontinuities
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 One way to hamper such issues is to resort to regularization of the $GW$ self-energy.
-Of course, the way of regularizing the self-energy is not unique but here we consider 3 different ways directly imported from MP2 theory.
+Of course, the way of regularizing the self-energy is not unique but here we consider three different ways directly imported from MP2 theory.
 
 This helps greatly convergence for (partially) self-consistent $GW$ methods.
 
 The regularized self-energy is 
 \begin{equation}
 \begin{split}
-	\rSigC{p}(\omega;s) 
+	\rSigC{p}(\omega;\kappa) 
-	& = \sum_{im} 2\ERI{pi}{m}^2 f_s(\omega - \eps{i}{\HF} + \Om{m}{\RPA})
+	& = \sum_{im} 2\ERI{pi}{m}^2 f_\kappa(\omega - \eps{i}{\HF} + \Om{m}{\RPA})
 	\\
-	& + \sum_{am} 2\ERI{pa}{m}^2 f_s(\omega - \eps{a}{\HF} - \Om{m}{\RPA})
+	& + \sum_{am} 2\ERI{pa}{m}^2 f_\kappa(\omega - \eps{a}{\HF} - \Om{m}{\RPA})
 \end{split}
 \end{equation}
-where various choices for the ``regularizer'' $f_s$ are possible.
+where various choices for the ``regularizer'' $f_\kappa$ are possible.
 Our investigation have shown that
 \begin{equation}
-	f_s(\Delta) = \frac{1-e^{-2s\Delta^2}}{\Delta}
+	f_\kappa(\Delta) = \frac{1-e^{-2\kappa\Delta^2}}{\Delta}
 \end{equation}
 is a very convenient form and has been derived from the flow equation within driven-similarity renormalization group \cite{Evangelista_2014}
 Of course, by construction, we have