The combined many-body Green's function $GW$ and Bethe-Salpeter equation (BSE) formalisms have shown to be a promising alternative to time-dependent density-functional theory (TD-DFT) in order to compute vertical transition energies of molecular systems.
The BSE formalism can also be employed to compute ground-state correlation energies thanks to the adiabatic-connection fluctuation-dissipation theorem (ACFDT).
Here, we study the topological features of the ground-state potential energy surfaces (PES) of several diatomic molecules near their equilibrium distance.
Thanks to comparisons with state-of-art computational approaches, we show that ACFDT@BSE is surprisingly accurate, and can even compete with coupled cluster methods in terms of total energies, equilibrium distances or \titou{harmonic vibrational frequencies}.
However, we sometimes observe unphysical irregularities on the ground-state PES, in relation with the appearance of satellite resonances with a weight similar to that of the $GW$ quasiparticle peak.
With a similar computational cost to time-dependent density-functional theory (TD-DFT), \cite{Runge_1984,Casida} the many-body Green's function Bethe-Salpeter equation (BSE) formalism
\cite{Salpeter_1951,Strinati_1988,Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999} is a valuable alternative that has gained momentum in the past few years for the study of molecular systems \cite{Ma_2009,Pushchnig_2002,Tiago_2003,Palumno_2009,Rocca_2010,Sharifzadeh_2012,Cudazzo_2012,Boulanger_2014,Ljungberg_2015,Hirose_2015,Cocchi_2015,Ziaei_2017,Abramson_2017}
It now stands as a computationally inexpensive method that can effectively model excited states \cite{Gonzales_2012,Loos_2020a} with a typical error of $0.1$--$0.3$ eV according to large and systematic benchmark calculations. \cite{Jacquemin_2015,Bruneval_2015,Blase_2016,Jacquemin_2016,Hung_2016,Hung_2017,Krause_2017,Jacquemin_2017,Blase_2018}
One of the main advantages of BSE compared to TD-DFT is that it allows a faithful description of charge-transfer states. \cite{Lastra_2011,Blase_2011b,Baumeier_2012,Duchemin_2012,Cudazzo_2013,Ziaei_2016}
Moreover, when performed on top of a (partially) self-consistently {\evGW} calculation, \cite{Hybertsen_1986, Shishkin_2007, Blase_2011, Faber_2011,Rangel_2016,Kaplan_2016,Gui_2018} BSE@{\evGW} has been shown to be weakly dependent on its starting point (\ie, on the xc functional selected for the underlying DFT calculation). \cite{Jacquemin_2016,Gui_2018}
However, similar to adiabatic TD-DFT, \cite{Levine_2006,Tozer_2000,Huix-Rotllant_2010,Elliott_2011} the static version of BSE cannot describe multiple excitations. \cite{Romaniello_2009a,Sangalli_2011}
A significant limitation of the BSE formalism, as compared to TD-DFT, lies in the lack of analytic nuclear gradients (\ie, the first derivatives of the energy with respect to the nuclear displacements) for both the ground and excited states, \cite{Furche_2002} preventing efficient applications to the study of chemoluminescence, fluorescence and other related processes \cite{Navizet_2011} associated with geometric relaxation of ground and excited states, and structural changes upon electronic excitation. \cite{Bernardi_1996,Olivucci_2010,Robb_2007}
\cite{Lazzeri_2008,Faber_2011b,Yin_2013,Faber_2015,Montserrat_2016,Zhenglu_2019} only one pioneering study of the excited-state BSE gradients has been published so far. \cite{Beigi_2003} In this study devoted to small molecules (\ce{CO} and \ce{NH3}), only the BSE excitation energy gradients were calculated, while computing the KS-DFT (LDA) forces as its ground-state analog.
Contrary to TD-DFT which relies on Kohn-Sham density-functional theory (KS-DFT) \cite{Hohenberg_1964,Kohn_1965,ParrBook} as its ground-state counterpart, the BSE ground-state energy is not a well-defined quantity, and no clear consensus has been found regarding its formal definition.
It then remains in its infancy with very few available studies for atomic and molecular systems. \cite{Olsen_2014,Holzer_2018,Li_2019,Li_2020}
As a matter of fact, in the largest recent available benchmark study \cite{Holzer_2018} of the total energies of the atoms \ce{H}--\ce{Ne}, the atomization energies of the 26 small molecules forming the HEAT test set, \cite{Harding_2008} and the bond lengths and harmonic vibrational frequencies of $3d$ transition-metal monoxides, the BSE correlation energy, as evaluated within the adiabatic-connection fluctuation-dissipation (ACFDT) framework, \cite{Furche_2005} was mostly discarded from the set of tested techniques due to instabilities (negative frequency modes in the BSE polarization propagator) and replaced by an approximate (RPAsX) approach where the screened-Coulomb potential matrix elements was removed from the resonant electron-hole contribution. \cite{Maggio_2016,Holzer_2018}
With such an approximation, amounting to neglect excitonic effects in the electron-hole propagator, the question of using either KS-DFT or $GW$ eigenvalues in the construction of the propagator becomes further relevant, increasing accordingly the number of possible definitions for the ground-state correlation energy.
Finally, renormalizing or not the Coulomb interaction by the interaction strength $\IS$ in the Dyson equation for the interacting polarizability leads to two different versions of the BSE correlation energy. \cite{Holzer_2018}
Here, in analogy to the random-phase approximation (RPA)-type formalisms \cite{Furche_2008,Toulouse_2009,Toulouse_2010,Angyan_2011,Ren_2012} and similarly to Refs.~\onlinecite{Olsen_2014,Maggio_2016,Holzer_2018}, the ground-state BSE energy is calculated in the adiabatic connection framework.
Embracing this definition, the purpose of the present study is to investigate the quality of ground-state PES near equilibrium obtained within the BSE approach for several diatomic molecules.
Thanks to comparison with both similar and state-of-art computational approaches, we show that the ACFDT@BSE@$GW$ approach is surprisingly accurate, and can even compete with high-order coupled cluster (CC) methods in terms of absolute energies, equilibrium distances or \titou{harmonic vibrational frequencies}.
However, we also observe, in some cases, unphysical irregularities on the ground-state PES, which are due to the appearance of a satellite resonance with a weight similar to that of the $GW$ quasiparticle peak. \cite{vanSetten_2015,Maggio_2017,Loos_2018,Veril_2018,Duchemin_2020}
%In Sec.~\ref{sec:theo}, we introduce the equations behind the BSE formalism.
%In particular, the general equations are reported in Subsec.~\ref{sec:BSE}, and the corresponding equations obtained in a finite basis in Subsec.~\ref{sec:BSE_basis}.
%Subsection \ref{sec:BSE_energy} defines the BSE total energy, and various other quantities of interest for this study.
%Computational details needed to reproduce the results of the present work are reported in Sec.~\ref{sec:comp_details}.
%Section \ref{sec:PES} reports ground-state PES for various diatomic molecules.
%Finally, we draw our conclusion in Sec.~\ref{sec:conclusion}.
In order to compute the neutral (optical) excitations of the system and their associated oscillator strengths, the BSE expresses the two-body propagator \cite{Strinati_1988, Martin_2016}
as the linear response of the one-body Green's function $\G{}$ with respect to a general non-local external potential
\begin{equation}
\XiBSE{}(3,5,4,6) = i \fdv{[\vc{\Ha}(3) \delta(3,4) + \Sig{\xc}(3,4)]}{\G{}(6,5)},
\end{equation}
which takes into account the self-consistent variation of the Hartree potential
\begin{equation}
\vc{\Ha}(1) = - i \int d2 \vc{}(2) \G{}(2,2^+),
\end{equation}
(where $\vc{}$ is the bare Coulomb operator) and the exchange-correlation self-energy $\Sig{\xc}$.
In Eq.~\eqref{eq:BSE}, $\LBSE{0}(1,2,1',2')=- i \G{}(1,2')\G{}(2,1')$, and $(1)=(\br{1},\sigma_1,t_1)$ is a composite index gathering space, spin and time variables.
where, as commonly done, we have neglected the term $\delta\W{}{}/\delta\G{}$ in the functional derivative of the self-energy. \cite{Hanke_1980,Strinati_1984,Strinati_1982}
Finally, the static approximation is enforced, \ie, $\W{}{}(1,2)=\W{}{}(\{\br{1}, \sigma_1, t_1\},\{\br{2}, \sigma_2, t_2\})\delta(t_1-t_2)$, which corresponds to restricting $\W{}{}$ to its static limit, \ie, $\W{}{}(1,2)=\W{}{}(\{\br{1},\sigma_1\},\{\br{2},\sigma_2\}; \omega=0)$.
For a closed-shell system, to compute the singlet BSE excitation energies (within the static approximation) of the physical system (\ie, $\IS=1$) in a finite basis, one must solve the following linear response problem \cite{Casida,Dreuw_2005,Martin_2016}
where $\Om{m}{\IS}$ is the $m$th excitation energy with eigenvector $\T{(\bX{\IS}_m \,\bY{\IS}_m)}$ at interaction strength $\IS$, $\T{}$ is the matrix transpose, and we assume real-valued spatial orbitals $\{\MO{p}(\br{})\}_{1\le p \le\Norb}$.
The matrices $\bA{\IS}$, $\bB{\IS}$, $\bX{\IS}$, and $\bY{\IS}$ are all of size $\Nocc\Nvir\times\Nocc\Nvir$ where $\Nocc$ and $\Nvir$ are the number of occupied and virtual orbitals (\ie, $\Norb=\Nocc+\Nvir$ is the total number of spatial orbitals), respectively.
In the following, the index $m$ labels the $\Nocc\Nvir$ single excitations, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals.
In the absence of instabilities (\ie, $\bA{\IS}-\bB{\IS}$ is positive-definite), \cite{Dreuw_2005} Eq.~\eqref{eq:LR} is usually transformed into an Hermitian eigenvalue problem of smaller dimension
are the bare two-electron integrals, and $\delta_{pq}$ is the Kronecker delta. \cite{NISTbook}
In Eq.~\eqref{eq:W}, $\eta$ is a positive infinitesimal, and $\OmRPA{m}{\IS}$ are the direct (\ie, without exchange) RPA neutral excitation energies computed by solving the linear eigenvalue problem \eqref{eq:LR} with the following matrix elements
Note that the present definition of the BSE correlation energy [see Eq.~\eqref{eq:EcBSE}], which we refer to as BSE@$GW$@HF here, has be labeled ``XBS'' for ``extended Bethe Salpeter'' by Holzer et al. \cite{Holzer_2018}
Contrary to DFT where the electron density is fixed along the adiabatic path, in the present formalism, the density is not maintain with respect to $\IS$.
Therefore, an additional contribution to Eq.~\eqref{eq:EcBSE} originating from the variation of the Green's function along the adiabatic connection should be added.
However, as commonly done within RPA and RPAx (\ie, RPA with exchange), \cite{Toulouse_2009, Toulouse_2010, Holzer_2018} we shall neglect it in the present study.
Equation \eqref{eq:EcBSE} can also be straightforwardly applied to RPA and RPAx, the only difference being the expressions of $\bA{\IS}$ and $\bB{\IS}$ used to obtain the eigenvectors $\bX{\IS}$ and $\bY{\IS}$ entering the definition of $\bP{\IS}$ [see Eq.~\eqref{eq:2DM}].
For RPA, these expressions have been provided in Eq.~\eqref{eq:LR_RPA}, and their RPAx analogs read
These two types of calculations will be refer to as RPA@HF and RPAx@HF respectively in the following.
Finally, we will also consider the RPA@$GW$@HF scheme which consists in replacing the HF orbital energies in Eq.~\eqref{eq:LR_RPA} by the $GW$ quasiparticles energies.
For spin-restricted closed-shell molecular systems around their equilibrium geometry (such as the ones studied here), it is rare to encounter singlet instabilities as these systems can be classified as weakly correlated.
However, singlet instabilities may appear in the presence of strong correlation, \eg, when the bond is stretched, hampering in particular the calculation of atomization energies. \cite{Holzer_2018}
Even for weakly correlated systems, triplet instabilities are much more common, but triplet excitations do not contribute to the correlation energy in the ACFDT formulation. \cite{Toulouse_2009, Toulouse_2010, Angyan_2011}
%\xavier{ IS THIS NEEDED NOW ? However, contrary to the plasmon formulation (an alternative expression of the BSE correlation energy), \cite{Schuck_Book, Sawada_1957b, Rowe_1968} the triplet excitations do not contribute in the ACFDT formulation, which is an indisputable advantage of this approach.
%Indeed, although the plasmon and adiabatic connection formulations are equivalent for RPA, \cite{Sawada_1957b, Gell-Mann_1957, Fukuta_1964, Furche_2008} this is not the case at the BSE and RPAx levels. \cite{Toulouse_2009, Toulouse_2010, Angyan_2011, Li_2020}
%For RPAx, an alternative plasmon expression (equivalent to its ACFDT analog) can be found if exchange is also added to the interaction kernel [see Eq.~\eqref{eq:K}]. \cite{Angyan_2011}
%However, to the best of our knowledge, such alternative plasmon expression does not exist for BSE. }
All the $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies have been done using a (restricted) Hartree-Fock (HF) starting point, which is a very adequate choice in the case of the (small) systems that we have considered here.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986} calculations are employed as starting point to compute the BSE neutral excitations.
The numerical integration required to compute the correlation energy along the adiabatic path [see Eq.~\eqref{eq:EcBSE}] has been performed with a 21-point Gauss-Legendre quadrature.
Comparison with the so-called plasmon (or trace) formula \cite{Furche_2008} at the RPA level has confirmed the excellent accuracy of this quadrature scheme over $\IS$.
For comparison purposes, we have also computed the PES at the MP2, CC2 \cite{Christiansen_1995}, CCSD, \cite{Purvis_1982} and CC3 \cite{Christiansen_1995b} levels of theory using DALTON. \cite{dalton}
Unless, otherwise stated, the frozen-core approximation has not been enforced in our calculations in order to provide a fair comparison between methods.
Because Eq.~\eqref{eq:EcBSE} requires the entire BSE singlet excitation spectrum for each quadrature point, we perform several complete diagonalization of the $\Nocc\Nvir\times\Nocc\Nvir$ BSE linear response matrix [see Eq.~\eqref{eq:small-LR}], which corresponds to a $\order{\Nocc^3\Nvir^3}=\order{\Norb^6}$ computational cost.
However, we are currently pursuing different avenues to lower this cost by computing the two-electron density matrix of Eq.~\eqref{eq:2DM} via a quadrature in frequency space. \cite{Duchemin_2019,Duchemin_2020}
In order to illustrate the performance of the BSE-based adiabatic connection formulation, we have computed the ground-state PES of several closed-shell diatomic molecules around their equilibrium geometry: \ce{H2}, \ce{LiH}, \ce{LiF}, \ce{N2}, \ce{CO}, \ce{BF}, \ce{F2}, and \ce{HCl}.
The PES of these molecules for various methods and Dunning's triple-$\zeta$ basis cc-pVTZ are represented in Figs.~\ref{fig:PES-H2-LiH}, \ref{fig:PES-LiF-HCl}, \ref{fig:PES-N2-CO-BF}, and \ref{fig:PES-F2}, while the computed equilibrium distances are gathered in Table \ref{tab:Req}.
Let us first start with the two smallest molecules, \ce{H2} and \ce{LiH} which are both linked by covalent bonds (see Fig.~\ref{fig:PES-H2-LiH}).
For \ce{H2}, we take as reference the full configuration interaction (FCI) energies and we also report the MP2 curve and its third-order variant (MP3), which improves upon MP2 towards FCI.
RPA@HF and RPA@{\GOWO}@HF yield almost identical results, and significantly overestimate (in absolute value) the FCI energy, while RPAx@HF and BSE@{\GOWO}@HF slightly underestimate and overestimate the FCI energy, respectively, RPAx@HF being the best match in the case of \ce{H2}.
Interestingly though, the BSE@{\GOWO}@HF scheme yields a more accurate equilibrium bond length than any other method irrespectively of the basis set.
For example, with the cc-pVQZ basis set, BSE@{\GOWO}@HF is only off by $0.003$ bohr compared to FCI, while RPAx@HF, MP2, and CC2 underestimate the bond length by $0.008$, $0.011$, and $0.011$ bohr, respectively.
The RPA-based schemes are much less accurate, with even shorter equilibrium bond lengths.
Ground-state PES of \ce{H2} (left) and \ce{LiH} (right) around their respective equilibrium geometry obtained at various levels of theory with the cc-pVTZ basis set.
Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-H2-LiH}
}
\end{figure*}
%%% %%% %%%
The cases of \ce{LiF} and \ce{HCl} (see Fig.~\ref{fig:PES-LiF-HCl}) are interesting as they corresponds to a strongly polarized bond towards the halogen atoms which are much more electronegative than the first row elements.
For these ionic bond, the performance of BSE@{\GOWO}@HF are terrific with an almost perfect match to the CC3 curve.
For \ce{LiF}, the two curves starting to deviate a few tenths of bohr after the equilibrium geometry, but they remain tightly ... for much longer in the case of \ce{HCl}.
Maybe surprisingly, BSE@{\GOWO}@HF outperforms both CC2 and CCSD, as well as RPAx@HF by a big margin for these two molecules exhibiting charge transfer.
However, in the case of \ce{LiF}, the attentive reader would have observed a small glitch in the $GW$-based curves very close to their minimum.
Ground-state PES of \ce{LiF} (left) and \ce{HCl} (right) around their respective equilibrium geometry obtained at various levels of theory with the cc-pVTZ basis set.
Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-LiF-HCl}
}
\end{figure*}
%%% %%% %%%
Let us now look at the isoelectronic series \ce{N2}, \ce{CO}, and \ce{BF}, which have a decreasing bond order (from triple bond to single bond).
In that case again, the performance of BSE@{\GOWO}@HF are outstanding as shown in Fig.~\ref{fig:PES-N2-CO-BF}.
%%% FIG 3 %%%
\begin{figure*}
\includegraphics[width=0.33\linewidth]{N2_GS_VTZ}
\includegraphics[width=0.33\linewidth]{CO_GS_VTZ}
\includegraphics[width=0.33\linewidth]{BF_GS_VTZ}
\caption{
Ground-state PES of the isoelectronic series \ce{N2} (left), \ce{CO} (center), and \ce{BF} (right) around their respective equilibrium geometry obtained at various levels of theory with the cc-pVTZ basis set.
Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
\label{fig:PES-N2-CO-BF}
}
\end{figure*}
%%% %%% %%%
The \ce{F2} molecule is a notoriously difficult case to treat due to the relative weakness of its covalent bond (see Fig.~\ref{fig:PES-F2}).
%%% FIG 4 %%%
\begin{figure}
\includegraphics[width=\linewidth]{F2_GS_VTZ}
\caption{
Ground-state PES of \ce{F2} around its equilibrium geometry obtained at various levels of theory with the cc-pVTZ basis set.
Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
In this Letter, we have shown that calculating the BSE correlation energy in the ACFDT framework yield extremely accurate PES around equilibrium.
We have illustrated this for 8 diatomic molecules for which we have also computed reference ground-state energies using coupled cluster methods (CC2, CCSD, and CC3).
However, we have also observed that, in some cases, unphysical irregularities on the ground-state PES due to the appearance of a satellite resonance with a weight similar to that of the $GW$ quasiparticle peak.
This shortcoming, which is entirely due to the quasiparticle nature of the underlying $GW$ calculation, has been thoroughly described in several previous studies.\cite{vanSetten_2015,Maggio_2017,Loos_2018,Veril_2018,Duchemin_2020}
In the perspective of developing analytical nuclear gradients within the BSE@$GW$ formalism, we are currently investigating the accuracy of the ACFDT@BSE scheme for excited-state PES.
This work has been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''. }
