ufGW/Response_Letter/Response_Letter.tex

\documentclass[10pt]{letter}
\usepackage{UPS_letterhead,xcolor,mhchem,ragged2e,hyperref}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\definecolor{darkgreen}{HTML}{009900}


\begin{document}

\begin{letter}%
{To the Editors of the Journal of Chemical Physics,}

\opening{Dear Editors,}

\justifying
Please find attached a revised version of the manuscript entitled 
\begin{quote}
	\textit{``Unphysical Discontinuities, Intruder States and Regularization in $GW$ Methods''}.
\end{quote}
We thank the reviewers for their constructive comments.
Our detailed responses to their comments can be found below.
For convenience, changes are highlighted in red in the revised version of the manuscript. 

We look forward to hearing from you.

\closing{Sincerely, the authors.}

\newpage

%%% REVIEWER 1 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#1}
 
{This in an interesting article and definitely deserves publication in JCP. 
However, I think that there are few issues the authors need to address before publication of this article. 
In my opinion addressing of these issues would be helpful to both quantum chemistry audience and even Green's function practitioners. 
}
\\
\alert{
Thank you for supporting publication of the present manuscript.
As detailed below, we have taken into account the comments and suggestions of the reviewers that we believe have overall improved the quality of the present paper.}

{In my understanding methods that are solved through non-linear equations display multiple solutions due to the dependence on the starting point and presence of multiple local minima. 
While the presence of these multiple solutions may be connected to the disappearance of denominators, the problem is much more general. 
The problem of multiple solutions exists even if fully self-consistent GW is solved in a finite temperature formalism (on imaginary axis - which effectively removes the problem with vanishing denominators). 
Different solutions can be then illustrated by considering their physical properties for example $S^2$. 
In such cases, the potential energy curves are smooth within a region where one solution is dominant and one can show properties of these solutions as a function of $S^2$, see J. Chem. Phys. 155, 024101, and J. Chem. Phys. 155, 024119.}
\\
\alert{
The reviewer is correct; these issues are more general and do not originate solely from the vanishing denominator.
For example, we have recently studied multiple solutions in another non-linear method, coupled-cluster theory [J. Chem. Theory Comput. 17, 4756 (2021)], where their origin stems from the non-linear nature of the equations, only.
However, in the case of $GW$ at least, we have clearly observed that these multiple solution issues appear when denominators vanish.
More precisely, in most $GW$ methods, because one only considers the quasiparticle solution at each iteration while discarding the other solutions (known as satellites, discontinuities are observed. 
This has been studied in details in a previous study [J. Chem. Theory Comput. 14, 5220 (2018)].
In a fully self-consistent scheme where one takes into account all these solutions at each iteration, these issues do not appear, as we have recently discussed in a recent paper [Front. Chem. (9, 751054 (2021)]. 
The references provided by the reviewer has been added in due place and this point has been clarified in the revised version of the manuscript.}

{Convergence accelerators such as DIIS, KAIN, LCIS, can be successfully used and smooth PES can be produced when one is relatively close to a local minima, one should not expect that DIIS will help with convergence where multiple close lying solutions exists and when two solutions are competing, see J. Chem. Phys. 156, 094101 (2022).}
\\
\alert{
Indeed, convergence accelerators such as DIIS can be used to ease convergence but they will not make these discontinuities disappear as their origin is more profond.
This particular case is discussed in a recent paper [J. Chem. Theory Comput. 14, 5220 (2018)] where we have provided our implementation of DIIS within $GW$ methods (see Appendix).
This point has been further stressed in the revised version of the manuscript.}


{It is my understanding that the calculation on the illustrative example H2 in 6-31G basis illustrates the existence of multiple such solutions. 
The issue of introducing a regularizer to avoid vanishing denominators is more of an issue for GW on real axis. 
Of course these issues of vanishing denominators and multiple solutions due to the non-linear character of the Dyson equation are connected but I think it is worth mentioning explicitly how different versions of GW are affected by these issues. 
Frankly, I find it infuriating that in the GW community many times it is not clearly spelled out which flavor of the so-called GW is used and how the equations are solved. 
}
\\
\alert{We agree with the reviewer that, historically, the literature has been unclear on how to solve the $GW$ equations, but we would like to stress that our group has made a clear effort to provide all the working equations and details necessary to solve these equations in the different cases.
We have added additional details and references in order to guide the readers to the relevant papers when one can find all the necessary details regarding the implementation of $GW$ methods.}

{Here, I will describe more specific issues that I noticed during the reading of the manuscript:}
\\
\alert{
}

\begin{enumerate}

\item 
{The authors should specify on which axis (omega) they perform their evaluation since this will not be clear to a regular chemistry audience. 
The convergence characteristics on imaginary or real axis is very different. 
Please, specify it. It may seem trivial to a regular GW practitioner, it is not trivial for a regular theoretical chemist not very familiar with Green's functions.}
\\
\alert{
The evaluation of all the quantities of interest ares performed on the real axis, i.e., for real values of $\omega$.
This is now stated clearly in the revised version of the manuscript at the very beginning of Section II.
}

\item 
{It would be reassuring if the authors mentioned that they use $O(K^9)$ scheme just to illustrate multiple solutions and this should not be done in regular calculations. 
A reader not familiar with GW algorithms may conclude that this what they advocate that should be done. 
On the same note, it would be good if the authors mentioned that the scheme that scales as $O(K^6)$ is again not what one does in most GW schemes that are in use for calculations of larger systems.}
\\
\alert{
Thank you for pointing this out. 
We now clearly state that our $O(K^9)$ scheme is illustrative and that state-of-the-art $GW$ calculations scales as $O(K^3)$ instead of $O(K^6)$.
}

\item 
{In expressions such as $\epsilon^{GW}_{p,s}$ $s$ is numbering solutions. 
I do not think that it is clearly mentioned in the article when the notation appears first time.}
\\
\alert{
We now clearly state that $s$ numbers solutions just above Eq.~(8).
}

\item 
{I assume that in Fig.~1 authors plot the lowest solution after the diagonalization for each orbital $p$.}
\\
\alert{In Fig.~1 we plot the quasiparticle solution for each orbital, i.e, the solution with the largest spectral weight which is obtained using the $G_0W_0$ scheme.
This is, in particular, mentioned in the caption of Fig.~1.
In general, for each orbital, it exists additional (satellites) solutions with lower and higher energies than the quasiparticle solution.
}

\item 
{In Fig.~2, the thick and thin solid lines are hardly visible. 
A better way of displaying is necessary. 
Otherwise it all is a bit incomprehensible.}
\\
\alert{Visibility of the different solid lines has been improved. 
We thank the reviewer for this valuable comment.
Figure 2 illustrates one of the key results of our work and must be as clear as possible.
}

\item 
{In Fig.~3, it is unclear to me that the traditional shift and the shift introduced by Evangelista result in demonstrating that Evangelista shift is superior. 
It seems to me that the difference between the values of qp energies before and post shifting are of the same order of magnitude for both regularizers. 
Could authors elaborate what they see differently and what my untrained eyes could not see?}
\\
\alert{We agree that this part was incomplete and, therefore, we have made extensive efforts in order to improve it.
In particular, we have changed the notations regarding the various regularizers that we have studied.
We now use $\eta$ for the traditional regularizer and $\kappa$ for Evangelista's regularizer.
We have also included additional graphs for different values of $\eta$ and $\kappa$ which shows how the quasiparticle energies are altered by the choice of the regularizing function and the values of $\eta$ and $\kappa$.
}

\item 
{In Fig.~4, again I do not understand author's conclusions. 
I understand that the regularization seems to introduce only 10 meV correction to the Homo and Lumo which is good. 
However, I cannot understand how the authors can claim that the correction introduced to the $p=3$ and 4 is viable. 
These qp energies are different by 2-3 eV. How are the smooth curves advantageous if the results are so incorrect? 
Could authors elaborate?}
\\
\alert{Indeed the HOMO and LUMO orbitals do not show discontinuities along the dissociation coordinate so no need for a correction. Thus, it is an important feature that the regularization introduces only a small correction for these orbitals. It is also true that the regularization introduces a correction of few eVs for the LUMO+1 ($p=3$) and LUMO+2 ($p=4$) orbitals but we have to note that the quasiparticle solutions of Eq.~2 for these orbitals appear at the poles of the self-energy. So the regularized self-energy has to do a large correction which leads to large error on the quasiparticle energies. Moreover, it is also essential to notice that we talk here about the $G_0W_0$ scheme but in case of a partially self-consistent scheme then the use of regularization seems critical.
}

\item 
{How the values of Fig.~4 depend on different choices of $\eta$ magnitude? This is crucial for assessing if a regularizer scheme is viable. }
\\
\alert{Several graphs have been added in the supporting information for different values of $\eta$ and $\kappa$.
The discussion has been expanded accordingly.
}

\item 
{In Fig.~5, it would be interesting to know how the total electronic energy depends on the choice of the regularizer and its magnitude. 
Without showing such a graph it is hard to know. 
Also would the regularizer help when the HF eigenvalues of Homo and Lumo are known to become degenerate at the stretched distance?}
\\
\alert{We add few graphs in the supporting information where we use different values for the $\eta$ and the $\kappa$ regularizer. 
We can see that the value $\kappa = 1$ gives a much better result than the $\eta = 1$ one.
}

\item 
{If I wanted to perform an actual calculation how I should choose a right regularizer? 
It is reasonable to expect that the value of $\eta$ is system dependent. 
How could I recognize which value is right?}
\\
\alert{
Like in other regularized methods, $\eta$ is an empirical parameter that must be chosen.
There is no definite answer but, depending on the type of properties studied, the value of $\eta$ must be chosen carefully.
This point is mentioned in the manuscript (Section V).
We hope to report further on this in a forthcoming paper but this requires extensive calculations. 
}

\item 
{In Fig.~5, I find very surprising that authors find that the curve is so "bumpy" already at the equilibrium distance. 
Usually in the quantum chemistry community, it is a standard practice to pass a solution from the previous geometry to the subsequent one as a starting point. 
This removes "bumpiness" in regions where denominators do not vanish. 
Are indeed the denominators already vanishing here at the equilibrium? 
If they do, then I find it highly unusual. 
Or are these competing solutions that could be removed by passing the starting point between neighboring geometries without a regularizer?}
\\
\alert{
}


\item 
{In my opinion, the novelty of this work is in showing how the unfolding scheme introduced by Booth group can be used to access different quasiparticle solutions. 
I am less convinced about the goodness of the regularizer since I suspect that the solutions depend quite a bit on the magnitude of $\eta$. 

Can this paper be published as a rapid communication? 
Yes, but a plot showing that the total energy is relatively insensitive as a function of regularizer $\eta$ would make the argument more convincing to broaden the applicability of the partially self-consistent Green's function on the real axis.}
\\
\alert{
As stated above, the regularizer section has been improved and we are happy if our manuscript is published as a Regular Article (as recommended by the editor).
}

\end{enumerate}

%%% REVIEWER 2 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#2}
 
{The article by Monino and Loos presents a way to avoid numerical difficulties arising from GW calculations in moderately correlated electronic systems, where a small perturbation in a particular parameter of the problem, such as the interatomic distance, leads to a sudden jump in the main quasiparticle peak, as computed within the GW approximation. The article is written clearly and includes new, original insights that could be useful for an audience of specialists interested in electronic-structure calculations. I support the article for publication, but ask the authors to address my points below: }
\\
\alert{Thank you for these positive comments and for supporting publication of our manuscript. 
Below, we address the points raised by Reviewer \#2.
}

\begin{enumerate}

\item 
{In standard GW calculations in solids, one assigns a quasiparticle peak to the solution of Dyson's equation that is associated with the largest value of the spectral function. 
This is the only general method that works in solids as, rigorously, there is no strict solution for Eq.~2 in a generic system - unless one lets $\omega$ be complex. 
Making this connection, especially around or before Eq.~8, would clarify the method to a readership that is not only interested in molecular systems.}
\\
\alert{
We have modified the manuscript below Eq.~(8) to clarify this point.
}

\item 
{I am not very happy with the term "intruder state" and with calling the procedure a "regularized GW method". 
These terms give the impression that there is a spurious solution that doesn't belong to a particular part of the Hilbert space and which leads to a sudden jump of the quasiparticle energy. 
Instead, my understanding is that the physics behind the jumps has nothing to do with spurious states, but one's insistence on making a one-to-one correspondence of the mean-field energy with a peak in the interacting spectral function. 
For strongly correlated systems, the spectral function may have pronounced satellite peaks, whose intensities may vary significantly with a small change in an external parameter. 
This is a point that the authors explain well in the text, but the terminology of intruder state confuses this otherwise simple physical picture. 
In addition, the term "regularized GW method" gives the impression that the GW approach has unphysical features. 
Instead, this procedure is nothing but a small broadening parameter that smooths out sudden jumps between neighboring peaks in the spectral function. }
\\
\alert{
We understand the point of the reviewer but "intruder state" and "regularization" are well-defined terms in the electronic structure community which are not linked with the appearance of spurious solutions.
The intruder state problem is well documented in multireference perturbation theory and comes usually from a poor choice of the active space.
By definition, an intruder state has a similar energy than the zeroth-order wave function and should be then moved in the model space; this is exactly what is happening in the case of $GW$.
}

\item 
{How does the result depend on the details of the "regularization" function? 
The authors should be aware that traditional GW calculations performed in solids often use a small broadening, and hence such a "regularization" is naturally captured by many codes - whether or not on purpose.}
\\
\alert{
We have specified this point in Section V of the revised manuscript.
Moreover, as detailed below (see the answers to Reviewer \#1), we have thoroughly modified and expanded this section to test the effect of the regularization function and its parameter.
In particular, we have changed the notations regarding the various regularizers that we have studied.
We now use $\eta$ for the traditional regularizer and $\kappa$ for Evangelista's regularizer.
We have also included additional graphs for different values of $\eta$ and $\kappa$ which shows how the quasiparticle energies are altered by the choice of the regularizing function and the values of $\eta$ and $\kappa$.
}

\item 
{I'm confused with the top row in Figure 3. 
Shouldn't the brown dashed curve with $\eta=0.01$ go to the solution with the main quasiparticle peak for the left and right limits of the plots? 
My naive intuition is that a small $\eta$ means that we are doing less smoothing to the self-energy/spectral function (or, in the language of the authors, performing a weaker "renormalization"). 
So, why doesn't the curve with $\eta=0.01$ jump abruptly to the solutions computed without the smoothing/"regularization" procedure?}
\\
\alert{
This was an unfortunate mistake that has been corrected in the revised version of the manuscript.
Thank you for spotting it.
}

\item 
{For all plots, the authors should include the units for $\eta$.}
\\
\alert{We specified the units for the $\eta$ and $\kappa$ for each graph.
}

\item 
{The authors should specify the meaning of the various acronyms in Fig.~5. 
}
\\
\alert{
The meaning of the acronyms in Fig.~5 have been clarified accordingly.
}

\end{enumerate}

\end{letter}
\end{document}