fig4

Manuscript/ufGW.tex

@@ -127,7 +127,7 @@
 \begin{abstract}
 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.
 Considering the $GW$ self-energy as an effective Hamiltonian, it is shown that these issues are key signatures of strong correlation in the $(N\pm1)$-electron states and can be directly related to the intruder state problem.
-A simple and efficient regularization procedure inspired by the similarity renormalization group is proposed to avoid such issues.
+A simple and efficient regularization procedure inspired by the similarity renormalization group is proposed to avoid such issues and speed up convergence of partially self-consistent $GW$ calculations.
 %\bigskip
 %\begin{center}
 %	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
@@ -140,7 +140,7 @@ A simple and efficient regularization procedure inspired by the similarity renor
 \section{Introduction}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-The $GW$ approximation of many-body perturbation theory, \cite{Hedin_1965,Martin_2016} allows to compute accurate charged excitation (\ie, ionization potentials, electron affinities and fundamental gaps) in solids and molecules. \cite{Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019}
+The $GW$ approximation of many-body perturbation theory \cite{Hedin_1965,Martin_2016} allows to compute accurate charged excitation (\ie, ionization potentials, electron affinities and fundamental gaps) in solids and molecules. \cite{Aryasetiawan_1998,Onida_2002,Reining_2017,Golze_2019}
 Its popularity in the molecular electronic 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 computational cost \cite{Foerster_2011,Liu_2016,Wilhelm_2018,Forster_2021,Duchemin_2021} and somehow surprising accuracy for weakly-correlated systems. \cite{Korbel_2014,vanSetten_2015,Caruso_2016,Hung_2017,vanSetten_2018,Bruneval_2021}
 
 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:
@@ -172,21 +172,21 @@ Electron correlation is then explicitly incorporated into one-body quantities vi
 %\end{subequations}
 %where $v$ is the bare Coulomb interaction, $\delta(12)$ is Dirac's delta function and $(1)$ is a composite coordinate gathering spin, space and time variables $(\sigma_1,\boldsymbol{r}_1,t_1)$.
 
-In recent studies, \cite{Loos_2018b,Veril_2018,Loos_2020e,Berger_2021,DiSabatino_2021} we discovered that one can observe (unphysical) irregularities and/or discontinuities in the energy surfaces of several key quantities (ionization potential, electron affinity, fundamental gap, total and correlation energies, as well as vertical excitation energies) even in the weakly-correlated regime. 
+In recent studies, \cite{Loos_2018b,Veril_2018,Loos_2020e,Berger_2021,DiSabatino_2021} we discovered that one can observe (unphysical) irregularities and/or discontinuities in the energy surfaces of several key quantities (ionization potential, electron affinity, fundamental and optical gaps, total and correlation energies, as well as excitation energies) even in the weakly-correlated regime. 
 These issues were discovered in Ref.~\onlinecite{Loos_2018b} while studying a model two-electron system \cite{Seidl_2007,Loos_2009a,Loos_2009c} and they were further investigated in Ref.~\onlinecite{Veril_2018}, where we provided additional evidences and explanations of these undesirable features in real molecular systems.
 In particular, we showed that each branch of the self-energy $\Sigma$ is associated with a distinct quasiparticle solution, and that each switch between solutions implies a significant discontinuity in the quasiparticle energy due to the transfer of weight between two solutions of the quasiparticle equation. \cite{Veril_2018}
 Multiple solution issues in $GW$ appears frequently, \cite{vanSetten_2015,Maggio_2017,Duchemin_2020} especially for orbitals that are energetically far from the Fermi level, such as in core ionized states. \cite{Golze_2018,Golze_2020}
 
-In addition to obvious irregularities on potential energy surfaces that hampers the accurate determination of properties such as equilibrium bond lengths and harmonic vibrational frequencies, \cite{Loos_2020e,Berger_2021} one direct consequence of these discontinuities is the difficulty to converge (partially) self-consistent $GW$ calculations as the self-consistent procedure jumps erratically from one solution to the other even if convergence accelerator techniques such as DIIS are employed. \cite{Pulay_1980,Pulay_1982,Veril_2018}
+In addition to obvious irregularities in potential energy surfaces that hampers the accurate determination of properties such as equilibrium bond lengths and harmonic vibrational frequencies, \cite{Loos_2020e,Berger_2021} one direct consequence of these discontinuities is the difficulty to converge (partially) self-consistent $GW$ calculations as the self-consistent procedure jumps erratically from one solution to the other even if convergence accelerator techniques such as DIIS are employed. \cite{Pulay_1980,Pulay_1982,Veril_2018}
 Note in passing that the present issues do not only appear in $GW$ as the $T$-matrix \cite{Romaniello_2012,Zhang_2017,Li_2021b,Loos_2022} and second-order Green's function (or second Born) formalisms \cite{SzaboBook,Casida_1989,Casida_1991,Stefanucci_2013,Ortiz_2013,Phillips_2014,Phillips_2015,Rusakov_2014,Rusakov_2016,Hirata_2015,Hirata_2017} exhibit the same drawbacks.
 
 It was shown that these problems can be tamed by using a static Coulomb-hole plus screened-exchange (COHSEX) \cite{Hedin_1965,Hybertsen_1986,Hedin_1999,Bruneval_2006} self-energy \cite{Berger_2021} or by considering a fully self-consistent $GW$ scheme, \cite{Stan_2006,Stan_2009,Rostgaard_2010,Caruso_2012,Caruso_2013,Caruso_2013a,Caruso_2013b,Koval_2014,Wilhelm_2018} where one considers not only the quasiparticle solution but also the satellites at each iteration. \cite{DiSabatino_2021}
 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 propose 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 and speed up convergence of partially self-consistent $GW$ calculations.
 
-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 such as ev$GW$, \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016} where one updates only the quasiparticle energies or qs$GW$, \cite{Gui_2018,Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016} where both quasiparticle energies and orbitals are updated at each iteration. 
+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 such as ev$GW$, \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016} where one updates only the quasiparticle energies, and qs$GW$, \cite{Gui_2018,Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016} where both quasiparticle energies and orbitals are updated at each iteration. 
 Moreover, we consider a Hartree-Fock (HF) starting point but it can be straightforwardly extended to a Kohn-Sham 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 vacant orbitals, while $m$ labels single excitations $i \to a$.
 Atomic units are used throughout.
@@ -211,7 +211,7 @@ Within the Tamm-Dancoff approximation (that we enforce here for the sake of simp
 \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 random-phase approximation (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, \ie, 
 \begin{equation}
 	\bA{}{\RPA} \cdot \bX{m}{\RPA} = \Om{m}{\RPA} \bX{m}{\RPA}
 \end{equation}
@@ -277,7 +277,7 @@ The size of this eigenvalue problem is $1 + \Nocc^2 \Nvir + \Nocc \Nvir^2 = \ord
 Thus, this step scales as $\order*{\Norb^9}$ with conventional diagonalization algorithms.
 Note, however, that the blocks $\bC{}{\text{2h1p}}$ and $\bC{}{\text{2p1h}}$ do not need to be recomputed for each orbital.
 
-It is paramount to understand that diagonalizing $\bH^{(p)}$ [see Eq.~\eqref{eq:Hp}] is completely equivalent to solving the quasiparticle equation \eqref{eq:qp_eq}.
+It is crucial to understand that diagonalizing $\bH^{(p)}$ [see Eq.~\eqref{eq:Hp}] is completely equivalent to solving the quasiparticle equation \eqref{eq:qp_eq}.
 This can be further illustrated by expanding the secular equation associated with Eq.~\eqref{eq:Hp}
 \begin{equation}
 	\det[ \bH^{(p)} - \omega \bI ] = 0 
@@ -340,7 +340,7 @@ This example was already considered in our previous work \cite{Veril_2018} but h
 The downfolded and upfolded {\GOWO} schemes have been implemented in the electronic structure package QuAcK \cite{QuAcK} which is freely available at \url{https://github.com/pfloos/QuAcK}.
 These calculations are based on restricted HF eigenvalues and orbitals.
 We denote as $\ket*{1\Bar{1}}$ the $\Ne$-electron ground-state Slater determinant where the orbital 1 is occupied by one spin-up and one spin-down electron.
-Similarly notations will be employed for the $(\Ne\pm1)$-electron configurations.
+Similar notations will be employed for the $(\Ne\pm1)$-electron configurations.
 
 In Fig.~\ref{fig:H2}, we report the variation of the quasiparticle energies of the four orbitals as functions of the internuclear distance $\RHH$. 
 One can easily diagnose two problematic regions showing obvious discontinuities around $\RHH = \SI{1.2}{\angstrom}$ for the LUMO$+1$ ($p = 3$) and $\RHH = \SI{0.5}{\angstrom}$ for the LUMO$+2$ ($p = 4$).
@@ -353,7 +353,7 @@ Inspection of their corresponding eigenvectors reveals that the $(\Ne+1)$-electr
 By construction, the quasiparticle solution diabatically follows the reference determinant $\ket*{1\Bar{1}3}$ through the avoided crossing (thick lines in Fig.~\ref{fig:H2_zoom}) which is precisely the origin of the energetic discontinuity.
 
 A similar scenario is at play in the region around $\RHH = \SI{0.5}{\angstrom}$ for the LUMO$+2$ (right panel of Fig.~\ref{fig:H2_zoom}) but it now involves three solutions ($s = 5$, $s = 6$, and $s = 7$).
-The electronic configurations of the Slater determinant involved are the $\ket*{1\Bar{1}4}$ reference determinant as well as two external determinants of configuration $\ket*{1\Bar{?}?}$ and $\ket*{1\Bar{?}?}$.
+The electronic configurations of the Slater determinant involved are the $\ket*{1\Bar{1}4}$ reference determinant as well as two external determinants of configuration \titou{$\ket*{1\Bar{?}?}$} and \titou{$\ket*{1\Bar{?}?}$}.
 These forms two avoided crossings in rapid successions, which create two discontinuities in the energy surface (see Fig.~\ref{fig:H2}).
 In this region, although the ground-state wave function is well described by the $\Ne$-electron HF determinant, a situation that can be safely labeled as single-reference, one can see that the $(\Ne+1)$-electron state involves three Slater determinants and can then be labeled as a multi-reference (or strongly-correlated) situation with near-degenerate electronic configurations. 
 Therefore, one can conclude that this downfall of $GW$ is a key signature of strong correlation in the $(\Ne\pm1)$-electron states that yields a significant redistribution of weights amongst electronic configurations.
@@ -380,7 +380,8 @@ Therefore, one can conclude that this downfall of $GW$ is a key signature of str
 %	FIGURE 4
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{figure}
-	\includegraphics[width=\linewidth]{fig4}
+%	\includegraphics[width=\linewidth]{fig4a}
+	\includegraphics[width=\linewidth]{fig4b}
 	\caption{
 	\label{fig:H2reg}
 	Difference between regularized and non-regularized quasiparticle energies $\reps{p}{\GW} - \eps{p}{\GW}$ as functions of the internuclear distance $\RHH$ (in \si{\angstrom}) of \ce{H2} at the {\GOWO}@HF/6-31G level.
@@ -439,6 +440,7 @@ However, this diabatization is more or less accurate depending on the value of $
 For $\eta = 10$, the value is clearly too large inducing a large difference between the two sets of quasiparticle energies (purple curves).
 For $\eta = 0.1$, we have the opposite scenario where the value is too small and some irregularities remain (green curves).
 We have found that $\eta = 1$ is a good compromise that does not alter significantly the quasiparticle energies while providing a smooth transition between the two solutions.
+\titou{This values could be further refined for specfici applications.}
 
 To further evidence this, Fig.~\ref{fig:H2reg} reports the difference between regularized and non-regularized quasiparticle energies as functions of $\RHH$.