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. 

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. 

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). 

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. 

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{
}

\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{
}

\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{
}

\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 ``normal'' $G_0W_0$ scheme.
}

\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. The authors thank the reviewer for this valuable comment.
}

\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{
}

\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. 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.
%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  
}

\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{
}

\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{
}

\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{
}

\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{
}

\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{
}

\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{
}

\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{
}

\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{
}

\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{
}

\item 
{For all plots, the authors should include the units for $\eta$.}
\\
\alert{
}

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

\end{enumerate}

\end{letter}
\end{document}