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


\begin{document}

\begin{letter}%
{To the Editors of the Journal of Physical Chemistry Letters}

\opening{Dear Editors,}

\justifying
Please find attached a revised version of the manuscript entitled 
\begin{quote}
\textit{``Pros and Cons of the Bethe-Salpeter Formalism for Ground-State Energies''}.
\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. 
Concerning Reviewer \#3's comment that suggests to contact the TURBOMOLE authors, this is certainly fine with us. 
However spotting out potential differences may require a large amount of time and we thus seek the Editor's advice on this issue. 
We describe in our reply all the checks we performed that leads us to conclude that our number accurately match the method we describe.

We look forward to hearing from you.

\closing{Sincerely, the authors.}

%%% REVIEWER 1 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#1}

\begin{itemize}

	\item 
	{This manuscript is very well written and presents an important contribution to the field of physical chemistry as a whole. 
	I do have a few minor remarks I would like the authors to consider.}
	\\
	\alert{We thank the reviewer for his/her kind comments.
	We have taken all his/her comments into account and our response to these comments can be found below.}

	\item 
	{The GW and BSE theory, approximations, and computational details are all very well described. 
	To my mind it would be good if the authors also add a bit more details on the CC calculations. 
	Details on whether Unrestricted Hartree-Fock or Restricted Hartree-Fock is used, if spatial symmetry is fixed may be useful. 
	I have recently seen examples where these details start to become visible, especially in the case where, like here, the results of the different types of calculations start to be very close.}
	\\
	\alert{All our calculations have been performed within the restricted Hartree-Fock formalism. 
	This is now clearly stated in the Computational Details section (see bottom of the right row in page 3).}

	\item 
	{For \ce{H2} figure 1 reports FCI as a reference but table 1 reports CC3 as a reference. 
	If both are used as reported it may be good to explicitly mention this to avoid confusion.}
	\\
	\alert{We now mentioned in Table 1 that both CCSD and CC3 are equivalent to FCI for the two-electron \ce{H2} molecule.}
	
	\item 
	{Finally in most, except \ce{LiH} molecules the CC3 and CCSD curves are considerably different, as compared to the difference between CC3 and BSE. 
	For instance for \ce{LiF} and \ce{HCl} I could easily imagine the difference between for example CCSDT(Q) and BSE may be even smaller. 
	For \ce{BF} and \ce{CO} higher levels of CC could even pass beyond the BSE curve. 
	Is there a clear reason not to go beyond CC3? 
	Could the authors provide an error bar on the CC3 results?}	
	\\
	\alert{Calculations with the cc-pVQZ beyond CC3 are extremely expensive.
	Nonetheless, in order to quantity the error in CC3, we have performed, at the CC3 equilibrium geometries, CCSDT  and CCSDT(Q) calculations.
	For the smallest systems (\ce{H2} and \ce{LiH}), we were also able to perform CCSDTQ calculations.
	The results, which are reported in the Supplementary information, show clearly that CC3 is very close to CCSDT, and more generally to the FCI limit.
	This statement is also mentioned in the revised version of the manuscript in the Computational Details section.}	

\end{itemize}

%%% REVIEWER 2 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#2}

\begin{itemize}
	\item 
	{This paper entitled ``Pros and Cons of the Bethe-Salpeter Formalism for Ground-State Energies'' is a theoretical/numerical study that calculates the correlation energies for 8 diatomic molecules around their equilibrium bond length.
	The reported results are in remarkable agreement with high-level wavefunctions methods, such as coupled-cluster (CC3).
	The calculations all are presented for diatomic molecules. 
	The generalization to more complex and less symmetric molecules would require an expression for the energy gradients, which is unfortunately missing in the literature.
	Though the figures in the paper are good, the BSE correlation energies is unlikely to have a significant impact in the community without a formula for the gradients.
	The results are limited to too simple cases to have a significant impact in a foreseeable future. 
	Therefore it may be not within the scope of J.~Phys.~Chem.~Lett.}
	\\
	\alert{We fully agree with the reviewer that it is unfortunate that we have not access to the analytical nuclear gradients.
	However, let us mention that we are currently working on this very actively.
	Moreover, as mentioned in the Introduction of our paper, \textit{``this study is a first necessary step towards the development of analytical nuclear gradients within the BSE@GW formalism''}, and thanks to the \textit{``remarkable agreement with high-level wavefunctions methods''}, we are eager to further develop the present methods, which we believe could have a significant impact on the entire electronic structure community in a foreseeable future.
	In short, we believe that we must first have a proof of concept of the accuracy of the method before embarking on developing the complex machinery of analytic gradients. 
	}
	
	\item 
	{The ``pros and cons'' mentioned in the title is misleading. 
	The ``cons'' are not really discussed in the text.}
	\\
	\alert{The paper is entitled ``pros and cons'' because (i) we demonstrate that ACFDT@BSE is extremely accurate for absolute energies around their equilibrium geometry (as pointed out by the reviewer), but (ii) irregularities are observed in the potential energy curves due to issues in the underlying GW calculation.
	We mentioned several times this fact in the original paper (in the abstract, end of the introduction, several times in the discussion, and in the conclusion) and we believe that these irregularities are very serious ``cons'', which hamper the applicability of $GW$ and BSE.
	We thus believe that the title of our manuscript is justified and appropriate.
	We emphasize that these irregularities do not condemn the usefulness of the present approach, but just question the way the community extract the ``quasiparticle energy'' from the full GW spectral function. }
	
	\item 
	{Furthermore, more stringent tests to the BSE correlation energies would be welcome.
	Correlation energies are fine at equilibrium distance. 
	But what about atomization energies?
	In the past, Kurth and Perdew proposed the ``RPA+'' method [Phys. Rev. B 59, 10461 (1999)] that was terrific to fix the absolute correlation energies, but was never used in practice because it deteriorated the atomization energies.
	Also in the previous studies, the full dissociation curves were chosen as test cases for ACFD correlation energies [Gonze et al. J. Chem. Phys. 122, 094116 (2005); Ren, et al. J. Mat. Science 47, 7447 (2012)].
	I think the manuscript would benefit from such an addition.}
	\\
	\alert{As mentioned several times in our manuscript, we do consider only diatomics around their equilibrium geometry, and we did not consider atomization energies, which have been considered in previous studies by several groups [see Refs. 66,67,93,104,105 in the revised version of the manuscript].
	As a matter of fact, it appears from our work that the difficulties related to static correlations in the dissociation limit have unfairly hidden the ``pros'' of the BSE correlation energy close to the equilibrium distance.  
	Along the line of Referee \#2's comment on the importance of gradients, the possibility to use BSE to relax structures in the ground-state and/or close to the minimum of the S1 surface (far away from the dissociation limit in most cases) is an important goal that should not be dismissed by the rather general dissociation limit problem.}
		
	\item 
	{As the authors concentrate on the quadratic region around the equilibrium, it would be informative to have a table presenting the vibrational frequencies.}	
	\\
	\alert{Again, vibrational frequencies and other related quantities of interest will be thoroughly studied when analytical gradients will be available.
	The present paper, which is already quite long for a letter, focused specifically on energetic and comparison with high-level CC methods.}
	\item 
	{Please remove the claim in the conclusion about getting 99\% of the correlation energy when \ce{H2} and \ce{LiH} (2 molecules out of 8) fail much beyond this error bar.}	
	\\
	\alert{The reviewer is right; this claim has been removed.}	
	
\end{itemize}

%%% REVIEWER 3 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#3}

\begin{itemize}
	\item 
	{This is a very nice letter on an interesting subject. 
	It is presented well and I recommend publication in The Journal of Physical Chemistry Letters. 
	The authors have investigated potential energy curves of small molecules such as \ce{H2}, \ce{LiH}, \ce{LiF}, \ce{HCl}, \ce{N2}, \ce{CO}, \ce{BF}, and \ce{F2} to test a new computational approach with respect to computing the Bethe-Salpeter correlation energy. 
	One interesting aspect is the fact that they start from the Hartree-Fock approximation to the electronic wave function (in place of Kohn-Sham theory). 
	The results are amazingly close to those obtain from high-level coupled-cluster calculations. 
	Publication is recommended, but major revision is needed. 
	In particular, the authors should address the following points:}
	\\
	\alert{We thank the referee for supporting publication of the present letter.
	His/her comments are addressed below.}
	
	\item 
	{The computation of the screened Coulomb matrix elements according to Eq.~(13) requires the RPA energies $\Omega_m^{\lambda,\text{RPA}}$.  
	Is it true that these RPA energies are computed by solving the (direct) RPA equations with the matrix elements of Eq.~(15), that is, with Hartree-Fock orbital energies on the diagonal? 
	Usually, quasiparticle energies are used to construct the screened Coulomb matrix. 
	This should be discussed.}
	\\
	\alert{Indeed, the screened Coulomb matrix elements are computed with the (direct) RPA neutral excitations, starting with HF orbital energies.
	The use of mean-field (DFT or HF) orbital energies, instead of self-consistent quasiparticle energies, is the standard ``single'' shot $G_0W_0$ approach.
	We have clarified this point in due place.
	In the case of ev$GW$ and qs$GW$, the quasiparticle energies are used to compute $W$ instead of the HF orbital energies.
	Our implementation was using the $G_0W_0$ quasiparticle energies to compute $W$ in the BSE calculation.
	For the sake of consistency (and accordingly to the equations reported in our manuscript), we have modified this and we now use the HF orbital energies to compute $W$ during the BSE calculation.
	The correlation energies in Table 1 have been updated accordingly (this does not alter the equilibrium geometry).
	As the reviewer will notice, the alteration on the correlation energies is very small (less than $1\%$), and does not explain the differences between our results and the TURBOMOLE numbers (see below).
	As a matter of fact, for these small systems, the close agreement between RPA@HF and RPA@G0W0@HF shows that there is no point in performing more expensive self-consistent $GW$ calculations.}
		
	\item 
	{Just below Eq.~(19), the matrix element $\Tilde{A}$ is defined. 
	The matrix $B$ of Eq.~(19) is not defined, but if this matrix is the same as in Eq.~(15b), then $B=2A$. 
	This factor of two does not seem to be correct.}
	\\
	\alert{We thank the reviewer for spotting this unfortunate typo. 
	We have indeed forgotten a factor 2 in the definition of $\Tilde{A}$.
	This has been corrected.
	Note that this factor 2 was present in our code.}
		
	\item 
	{The authors point out that the approach that they denote as ``BSE'' had been named ``XBS'' by the authors of Ref.~67. 
	It would be very helpful if results could be presented not only for this ``XBS'' approach but also for what the authors of Ref.~67 denote as ``BSE''.}
	\\
	\alert{We have added in the supplementary information the ``BSE'' correlation energies (as described in Ref.~67) computed at the XBS equilibrium geometries.
	The various PES obtained at this level of theory are not reported as they are very similar to the XBS ones.
	As one can see from these new values, we found the ``BSE'' and ``XBS'' schemes yield very similar results, and that the ``BSE'' scheme seems to provide slightly more accurate results due to compensation of errors.}	
	
	\item 
	{I have tried to reproduce the results presented in Table 1 by computing the correlation energy for the given bond length using the TURBOMOLE package of programs. 
	Unfortunately, I was unable to reproduce the BSE@G0W0@HF results of Table 1 (and neither two of the CC2 results). 
	This may be due to a somewhat different scheme for computing the linearized G0W0 quasiparticle energies. 
	I recommend to add all necessary computational details required to reproduce the G0W0 quasiparticle energies. 
	These could also be given explicitly in the Supporting Information. 
	Furthermore, the authors may wish to contact the authors of Ref.~67 to resolve the discrepancies by running TURBOMOLE computations.}	
	\\
	\alert{
	We thank the reviewer for his careful study of our manuscript.
	Indeed, one of our CC2 equilibrium bond length was wrongly reported (\ce{LiH}).
	This has been corrected ($3.010$ instead of $2.989$ bohr).
	For the CC2 equilibrium bond length and energies of LiF, we have double-checked our calculation and the value that we have originally reported is correct.
	Note that the error between our value and the one provided by the referee is only $0.2\%$.
	The difference between our XBS@G0W0@HF results and the ones provided by the reviewer is troublesome as they are significant (around 10\%).
	All the necessary computational details required to reproduce our results have been provided in the original manuscript.
	We don't believe any further information are required.
	\begin{enumerate}
	\item Let us mentioned that the discrepancies between our values and the TURBOMOLE numbers cannot originate from the linearization of the G0W0 quasiparticle equation as the referee was able to reproduce our RPA@G0W0@HF results which requires such quantities.
	Therefore, the GW part of the calculation is consistent between the two implementations.
	\item Because we cannot reproduce the BSE@G0W0@HF TURBOMOLE data, it means that the difference is not due to the $\lambda$-dependency of the screened Coulomb elements, as the  BSE@G0W0@HF does not require to recompute $W$ for each value of $\lambda$ (in BSE@G0W0@HF, $W$ is just scaled by $\lambda$ similarly to the Coulomb and exchange integrals in RPA and RPAx).
	\item Because the reviewer can reproduce all our RPA and RPAx numbers, this means that the linear response part of the calculation is correct.
	Moreover, because the plasmon and ACFDT scheme yield the same correlation energy, the numerical integration over the interaction strength $\lambda$ can be validated.
	\item We have performed extensive comparisons between our implementation and the MOLGW software, and we have been able to reproduce almost perfectly the G0W0 QP energies as well as the BSE excitation energies with our own code. Therefore, we believe that our BSE implementation is correct.
	\item We have studied the dependency of the XBS and BSE correlation energies with respect to various quantities. 
	For example, we have modified the level of approximations of $W$ or of the QP energies. 
	Our results are very robust and stable with variations of the order of $0.1$--$0.2$ mHa.
	\item We have also performed TURBOMOLE calculations with the latest commercial version (7.4.1), and we could perfectly reproduce the QP energies at the linearized G0W0 level.
	However, the BSE excitation energies (with the \texttt{noqpw} keyword to match our implementation) were slightly different from ours (which we recall are nearly identical to the MOLGW ones).
	Our guess is that there might be an inconsistency in the BSE calculation but we are not sure about its exact origin.
	For example, the BSE calculation is performed within the resolution of the identity (RI) in TURBOMOLE which is not the case in our code.
	However, we are surprised by the magnitude of the difference as the error associated with the RI approximation is usually much smaller than this.
	We will contact the authors of Ref.~67 to resolve the discrepancies between our implementation and TURBOMOLE.
	\end{enumerate}
}	
\end{itemize}

% H2
%                   R/a0   -Ec/mEh
% CCSD        1.402    40.4   v
% CC2           1.391    33.3   v
% MP2           1.391    33.2   v
% XBS@G0W0@HF   1.399    43.5
% BSE@G0W0@HF   1.399    45.6
% RPA@G0W0@HF   1.382    57.6   v
% RPAx@HF       1.394    37.9   v
% RPA@HF         1.386    57.3   v
%
% LiH
%                   R/a0   -Ec/mEh
% CCSD        3.020    69.8   v
% CC2           2.989    58.1   <---!
% MP2           3.008    57.9   v
% XBS@G0W0@HF   3.017    73.7
% BSE@G0W0@HF   3.017    76.9
% RPA@G0W0@HF   2.997   101.1   v
% RPAx@HF       3.011    65.2   v
% RPA@HF         2.994   100.2   v
%
% LiF
%                   R/a0   -Ec/mEh
% CCSD        2.953   372.6   v
% CC2           2.982   377.5   <---!
% MP2           2.970   373.0   v
% XBS@G0W0@HF   2.974   373.9   
% BSE@G0W0@HF   2.974   383.1
% RPA@G0W0@HF   2.965   473.0   v 
% RPAx@HF       2.944   343.6   v
% RPA@HF         2.946   465.9   v
%
% HCl
%                   R/a0   -Ec/mEh
% CCSD        2.398   370.2  (v)
% CC2           2.396   356.7  (v)
% MP2           2.395   355.5  (v)
% XBS@G0W0@HF   2.400   373.0
% BSE@G0W0@HF   2.400   381.6
% RPA@G0W0@HF   2.370   450.8  (v)
% RPAx@HF       2.391   344.1  (v)
% RPA@HF         2.382   442.3  (v)
%
% N2
%                   R/a0   -Ec/mEh
% CCSD        2.059   470.6   v
% CC2           2.106   488.0   v
% MP2           2.091   477.9   v 
% XBS@G0W0@HF   2.065   476.2
% BSE@G0W0@HF   2.065   488.3
% RPA@G0W0@HF   2.043   580.2   v 
% RPAx@HF       2.041   427.2   v
% RPA@HF         2.042   569.3   v 
%
% CO
%                   R/a0   -Ec/mEh
% CCSD        2.118   455.2   v
% CC2           2.156   465.5   v
% MP2           2.137   455.0   v
% XBS@G0W0@HF   2.134   460.3
% BSE@G0W0@HF   2.134   471.9
% RPA@G0W0@HF   2.132   566.5   v
% RPAx@HF       2.104   416.3   v
% RPA@HF         2.103   555.8   v
%
% BF
%                   R/a0   -Ec/mEh
% CCSD        2.380   432.9   v
% CC2           2.393   427.3   v
% MP2           2.382   421.6   v
% XBS@G0W0@HF   2.385   435.2
% BSE@G0W0@HF   2.385   446.2
% RPA@G0W0@HF   2.365   545.5   v
% RPAx@HF       2.366   399.0   v
% RPA@HF         2.364   537.6   v
%
% F2
%                   R/a0   -Ec/mEh
% CCSD        2.621   644.0   v
% CC2           2.665   654.9   v
% MP2           2.634   644.3   v
% XBS@G0W0@HF   2.640   647.2
% BSE@G0W0@HF   2.640   663.1
% RPA@G0W0@HF   2.571   794.2   v
% RPAx@HF       2.565   586.1   v
% RPA@HF         2.573   781.2   v
 
\end{letter}
\end{document}






