625 lines
40 KiB
TeX
625 lines
40 KiB
TeX
\documentclass[aps,prb,reprint,noshowkeys,superscriptaddress]{revtex4-1}
|
|
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
|
|
\usepackage[version=4]{mhchem}
|
|
|
|
\usepackage{natbib}
|
|
\usepackage[extra]{tipa}
|
|
\bibliographystyle{achemso}
|
|
\AtBeginDocument{\nocite{achemso-control}}
|
|
|
|
|
|
\usepackage[utf8]{inputenc}
|
|
\usepackage[T1]{fontenc}
|
|
\usepackage{txfonts}
|
|
|
|
\usepackage[
|
|
colorlinks=true,
|
|
citecolor=blue,
|
|
breaklinks=true
|
|
]{hyperref}
|
|
\urlstyle{same}
|
|
|
|
\newcommand{\ie}{\textit{i.e.}}
|
|
\newcommand{\eg}{\textit{e.g.}}
|
|
\newcommand{\alert}[1]{\textcolor{red}{#1}}
|
|
\definecolor{darkgreen}{HTML}{009900}
|
|
\usepackage[normalem]{ulem}
|
|
\newcommand{\titou}[1]{\textcolor{red}{#1}}
|
|
\newcommand{\denis}[1]{\textcolor{purple}{#1}}
|
|
\newcommand{\xavier}[1]{\textcolor{darkgreen}{#1}}
|
|
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
|
|
\newcommand{\trashDJ}[1]{\textcolor{purple}{\sout{#1}}}
|
|
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
|
|
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
|
|
\renewcommand{\DJ}[1]{\denis{(\underline{\bf DJ}: #1)}}
|
|
\newcommand{\XB}[1]{\xavier{(\underline{\bf XB}: #1)}}
|
|
|
|
\newcommand{\mc}{\multicolumn}
|
|
\newcommand{\fnm}{\footnotemark}
|
|
\newcommand{\fnt}{\footnotetext}
|
|
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
|
|
\newcommand{\SI}{\textcolor{blue}{supporting information}}
|
|
\newcommand{\QP}{\textsc{quantum package}}
|
|
\newcommand{\T}[1]{#1^{\intercal}}
|
|
|
|
% coordinates
|
|
\newcommand{\br}[1]{\mathbf{r}_{#1}}
|
|
\newcommand{\dbr}[1]{d\br{#1}}
|
|
|
|
% methods
|
|
\newcommand{\evGW}{ev$GW$}
|
|
\newcommand{\qsGW}{qs$GW$}
|
|
\newcommand{\GOWO}{$G_0W_0$}
|
|
\newcommand{\Hxc}{\text{Hxc}}
|
|
\newcommand{\xc}{\text{xc}}
|
|
\newcommand{\Ha}{\text{H}}
|
|
\newcommand{\co}{\text{x}}
|
|
|
|
%
|
|
\newcommand{\Norb}{N}
|
|
\newcommand{\Nocc}{O}
|
|
\newcommand{\Nvir}{V}
|
|
\newcommand{\IS}{\lambda}
|
|
|
|
% operators
|
|
\newcommand{\hH}{\Hat{H}}
|
|
|
|
% energies
|
|
\newcommand{\Enuc}{E^\text{nuc}}
|
|
\newcommand{\Ec}{E_\text{c}}
|
|
\newcommand{\EHF}{E^\text{HF}}
|
|
\newcommand{\EBSE}{E^\text{BSE}}
|
|
\newcommand{\EcRPA}{E_\text{c}^\text{RPA}}
|
|
\newcommand{\EcRPAx}{E_\text{c}^\text{RPAx}}
|
|
\newcommand{\EcBSE}{E_\text{c}^\text{BSE}}
|
|
\newcommand{\IP}{\text{IP}}
|
|
\newcommand{\EA}{\text{EA}}
|
|
|
|
% orbital energies
|
|
\newcommand{\e}[1]{\epsilon_{#1}}
|
|
\newcommand{\eHF}[1]{\epsilon^\text{HF}_{#1}}
|
|
\newcommand{\eKS}[1]{\epsilon^\text{KS}_{#1}}
|
|
\newcommand{\eQP}[1]{\epsilon^\text{QP}_{#1}}
|
|
\newcommand{\eGOWO}[1]{\epsilon^\text{\GOWO}_{#1}}
|
|
\newcommand{\eGW}[1]{\epsilon^{GW}_{#1}}
|
|
\newcommand{\eevGW}[1]{\epsilon^\text{\evGW}_{#1}}
|
|
\newcommand{\eGnWn}[2]{\epsilon^\text{\GnWn{#2}}_{#1}}
|
|
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
|
|
|
|
% Matrix elements
|
|
\newcommand{\A}[2]{A_{#1}^{#2}}
|
|
\newcommand{\tA}[2]{\Tilde{A}_{#1}^{#2}}
|
|
\newcommand{\B}[2]{B_{#1}^{#2}}
|
|
\renewcommand{\S}[1]{S_{#1}}
|
|
\newcommand{\ABSE}[2]{A_{#1}^{#2,\text{BSE}}}
|
|
\newcommand{\BBSE}[2]{B_{#1}^{#2,\text{BSE}}}
|
|
\newcommand{\ARPA}[2]{A_{#1}^{#2,\text{RPA}}}
|
|
\newcommand{\BRPA}[2]{B_{#1}^{#2,\text{RPA}}}
|
|
\newcommand{\ARPAx}[2]{A_{#1}^{#2,\text{RPAx}}}
|
|
\newcommand{\BRPAx}[2]{B_{#1}^{#2,\text{RPAx}}}
|
|
\newcommand{\G}[1]{G_{#1}}
|
|
\newcommand{\LBSE}[1]{L_{#1}}
|
|
\newcommand{\XiBSE}[1]{\Xi_{#1}}
|
|
\newcommand{\Po}[1]{P_{#1}}
|
|
\newcommand{\W}[2]{W_{#1}^{#2}}
|
|
\newcommand{\Wc}[1]{W^\text{c}_{#1}}
|
|
\newcommand{\vc}[1]{v_{#1}}
|
|
\newcommand{\Sig}[1]{\Sigma_{#1}}
|
|
\newcommand{\SigGW}[1]{\Sigma^{GW}_{#1}}
|
|
\newcommand{\Z}[1]{Z_{#1}}
|
|
\newcommand{\MO}[1]{\phi_{#1}}
|
|
\newcommand{\ERI}[2]{(#1|#2)}
|
|
\newcommand{\sERI}[2]{[#1|#2]}
|
|
|
|
%% bold in Table
|
|
\newcommand{\bb}[1]{\textbf{#1}}
|
|
\newcommand{\rb}[1]{\textbf{\textcolor{red}{#1}}}
|
|
\newcommand{\gb}[1]{\textbf{\textcolor{darkgreen}{#1}}}
|
|
|
|
% excitation energies
|
|
\newcommand{\OmRPA}[2]{\Omega_{#1}^{#2,\text{RPA}}}
|
|
\newcommand{\OmRPAx}[2]{\Omega_{#1}^{#2,\text{RPAx}}}
|
|
\newcommand{\OmBSE}[2]{\Omega_{#1}^{#2,\text{BSE}}}
|
|
|
|
\newcommand{\spinup}{\downarrow}
|
|
\newcommand{\spindw}{\uparrow}
|
|
\newcommand{\singlet}{\uparrow\downarrow}
|
|
\newcommand{\triplet}{\uparrow\uparrow}
|
|
|
|
% Matrices
|
|
\newcommand{\bO}{\mathbf{0}}
|
|
\newcommand{\bI}{\mathbf{1}}
|
|
\newcommand{\bvc}{\mathbf{v}}
|
|
\newcommand{\bSig}{\mathbf{\Sigma}}
|
|
\newcommand{\bSigX}{\mathbf{\Sigma}^\text{x}}
|
|
\newcommand{\bSigC}{\mathbf{\Sigma}^\text{c}}
|
|
\newcommand{\bSigGW}{\mathbf{\Sigma}^{GW}}
|
|
\newcommand{\be}{\mathbf{\epsilon}}
|
|
\newcommand{\beGW}{\mathbf{\epsilon}^{GW}}
|
|
\newcommand{\beGnWn}[1]{\mathbf{\epsilon}^\text{\GnWn{#1}}}
|
|
\newcommand{\bde}{\mathbf{\Delta\epsilon}}
|
|
\newcommand{\bdeHF}{\mathbf{\Delta\epsilon}^\text{HF}}
|
|
\newcommand{\bdeGW}{\mathbf{\Delta\epsilon}^{GW}}
|
|
\newcommand{\bOm}[1]{\mathbf{\Omega}^{#1}}
|
|
\newcommand{\bA}[1]{\mathbf{A}^{#1}}
|
|
\newcommand{\btA}[1]{\Tilde{\mathbf{A}}^{#1}}
|
|
\newcommand{\bB}[1]{\mathbf{B}^{#1}}
|
|
\newcommand{\bX}[1]{\mathbf{X}^{#1}}
|
|
\newcommand{\bY}[1]{\mathbf{Y}^{#1}}
|
|
\newcommand{\bZ}[1]{\mathbf{Z}^{#1}}
|
|
\newcommand{\bK}{\mathbf{K}}
|
|
\newcommand{\bP}[1]{\mathbf{P}^{#1}}
|
|
|
|
% units
|
|
\newcommand{\IneV}[1]{#1 eV}
|
|
\newcommand{\InAU}[1]{#1 a.u.}
|
|
\newcommand{\InAA}[1]{#1 \AA}
|
|
\newcommand{\kcal}{kcal/mol}
|
|
|
|
\newcommand{\NEEL}{Universit\'e Grenoble Alpes, CNRS, Institut NEEL, F-38042 Grenoble, France}
|
|
\newcommand{\CEISAM}{Laboratoire CEISAM - UMR CNRS 6230, Universit\'e de Nantes, 2 Rue de la Houssini\`ere, BP 92208, 44322 Nantes Cedex 3, France}
|
|
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
|
|
\newcommand{\CEA}{Universit\'e Grenoble Alpes, CEA, IRIG-MEM-L Sim, 38054 Grenoble, France}
|
|
|
|
\begin{document}
|
|
|
|
\title{Pros and Cons of the Bethe-Salpeter Formalism for Ground-State Potential Energy Surfaces}
|
|
|
|
\author{Pierre-Fran\c{c}ois \surname{Loos}}
|
|
\email{loos@irsamc.ups-tlse.fr}
|
|
\affiliation{\LCPQ}
|
|
\author{Anthony \surname{Scemama}}
|
|
\email{scemama@irsamc.ups-tlse.fr}
|
|
\affiliation{\LCPQ}
|
|
\author{Ivan \surname{Duchemin}}
|
|
\email{ivan.duchemin@cea.fr}
|
|
\affiliation{\CEA}
|
|
\author{Denis \surname{Jacquemin}}
|
|
\email{denis.jacquemin@univ-nantes.fr}
|
|
\affiliation{\CEISAM}
|
|
\author{Xavier \surname{Blase}}
|
|
\email{xavier.blase@neel.cnrs.fr }
|
|
\affiliation{\NEEL}
|
|
|
|
\begin{abstract}
|
|
%\begin{wrapfigure}[12]{o}[-1.25cm]{0.4\linewidth}
|
|
% \centering
|
|
% \includegraphics[width=\linewidth]{TOC}
|
|
%\end{wrapfigure}
|
|
The combined many-body Green's function $GW$ and Bethe-Salpeter equation (BSE) formalism has 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 bond length.
|
|
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 and equilibrium bond distances.
|
|
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.
|
|
\end{abstract}
|
|
|
|
\maketitle
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\section{Introduction}
|
|
%\label{sec:intro}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
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 studying 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,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 exchange-correlation functional selected for the underlying DFT calculation). \cite{Jacquemin_2015,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,Loos_2019}
|
|
|
|
A significant limitation of the BSE formalism, as compared to TD-DFT, lies in the lack of analytical 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}
|
|
While calculations of the $GW$ quasiparticle energy ionic gradients is becoming popular,
|
|
\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 Kohn-Sham (KS) LDA forces as its ground-state analog.
|
|
|
|
Contrary to TD-DFT which relies on 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 theorem (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}
|
|
Such a modified BSE polarization propagator was inspired by a previous study on the homogeneous electron gas. \cite{Maggio_2016}
|
|
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 (see below) 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.
|
|
The location of the minima on the ground-state PES is of particular interest.
|
|
This study is a first preliminary step towards the development of analytical nuclear gradients within the BSE@$GW$ formalism.
|
|
Thanks to comparisons 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 and equilibrium distances.
|
|
However, we also observe that, 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}
|
|
|
|
%The paper is organized as follows.
|
|
%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}.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\section{Theory}
|
|
%\label{sec:theo}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\subsection{The Bethe-Salpeter equation}
|
|
%\label{sec:BSE}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
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}
|
|
\begin{multline}
|
|
\label{eq:BSE}
|
|
\LBSE{}(1,2,1',2') = \LBSE{0}(1,2,1',2')
|
|
\\
|
|
+ \int d3 d4 d5 d6 \LBSE{0}(1,4,1',3) \XiBSE{}(3,5,4,6) \LBSE{}(6,2,5,2')
|
|
\end{multline}
|
|
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.
|
|
In the $GW$ approximation, \cite{Hedin_1965,Aryasetiawan_1998,Onida_2002,Martin_2016,Reining_2017} we have
|
|
\begin{equation}
|
|
\SigGW{\xc}(1,2) = i \G{}(1,2) \W{}{}(1^+,2),
|
|
\end{equation}
|
|
where $\W{}{}$ is the screened Coulomb operator, and hence the BSE reduces to
|
|
\begin{equation}
|
|
\XiBSE{}(3,5,4,6) = \delta(3,4) \delta(5,6) \vc{}(3,6) - \delta(3,6) \delta(4,5) \W{}{}(3,4),
|
|
\end{equation}
|
|
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)$.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\subsection{BSE in a finite basis}
|
|
%\label{sec:BSE_basis}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
For a closed-shell system in a finite basis, to compute the singlet BSE excitation energies (within the static approximation) of the physical system (\ie, $\IS = 1$), one must solve the following linear response problem \cite{Casida,Dreuw_2005,Martin_2016}
|
|
\begin{equation}
|
|
\label{eq:LR}
|
|
\begin{pmatrix}
|
|
\bA{\IS} & \bB{\IS} \\
|
|
-\bB{\IS} & -\bA{\IS} \\
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
\bX{\IS}_m \\
|
|
\bY{\IS}_m \\
|
|
\end{pmatrix}
|
|
=
|
|
\Om{m}{\IS}
|
|
\begin{pmatrix}
|
|
\bX{\IS}_m \\
|
|
\bY{\IS}_m \\
|
|
\end{pmatrix},
|
|
\end{equation}
|
|
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
|
|
\begin{equation}
|
|
\label{eq:small-LR}
|
|
(\bA{\IS} - \bB{\IS})^{1/2} (\bA{\IS} + \bB{\IS}) (\bA{\IS} - \bB{\IS})^{1/2} \bZ{\IS} = (\bOm{\IS})^2 \bZ{\IS},
|
|
\end{equation}
|
|
where the excitation amplitudes are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\bX{\IS} + \bY{\IS} = (\bOm{\IS})^{-1/2} (\bA{\IS} - \bB{\IS})^{+1/2} \bZ{\IS},
|
|
\\
|
|
\bX{\IS} - \bY{\IS} = (\bOm{\IS})^{+1/2} (\bA{\IS} - \bB{\IS})^{-1/2} \bZ{\IS}.
|
|
\end{align}
|
|
\end{subequations}
|
|
Introducing the so-called Mulliken notation for the bare two-electron integrals
|
|
\begin{equation}
|
|
\ERI{pq}{rs} = \iint \frac{\MO{p}(\br{}) \MO{q}(\br{}) \MO{r}(\br{}') \MO{s}(\br{}')}{\abs*{\br{} - \br{}'}} \dbr{} \dbr{}',
|
|
\end{equation}
|
|
and the corresponding (static) screened Coulomb potential matrix elements at coupling strength $\IS$
|
|
\begin{equation}
|
|
\W{pq,rs}{\IS} = \iint \MO{p}(\br{}) \MO{q}(\br{}) \W{}{\IS}(\br{},\br{}') \MO{r}(\br{}') \MO{s}(\br{}') \dbr{} \dbr{}',
|
|
\end{equation}
|
|
the BSE matrix elements read
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_BSE-A}
|
|
\ABSE{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \IS \qty[ 2 \ERI{ia}{jb} - \W{ij,ab}{\IS} ],
|
|
\\
|
|
\label{eq:LR_BSE-B}
|
|
\BBSE{ia,jb}{\IS} & = \IS \qty[ 2 \ERI{ia}{bj} - \W{ib,aj}{\IS} ],
|
|
\end{align}
|
|
\end{subequations}
|
|
where $\eGW{p}$ are the $GW$ quasiparticle energies.
|
|
%In the standard BSE implementation, the screened Coulomb potential $\W{}{\IS}$ is taken to be static $(\omega \rightarrow 0)$.
|
|
In the standard BSE approach, $\W{}{\IS}$ is built within the direct RPA scheme, \ie,
|
|
\begin{subequations}
|
|
\label{eq:wrpa}
|
|
\begin{align}
|
|
\W{}{\IS}(\br{},\br{}')
|
|
& = \int \frac{\epsilon_{\IS}^{-1}(\br{},\br{}''; \omega=0)}{\abs*{\br{}' - \br{}''}} \dbr{}'' ,
|
|
\\
|
|
\epsilon_{\IS}(\br{},\br{}'; \omega)
|
|
& = \delta(\br{}-\br{}') - \IS \int \frac{\chi_{0}(\br{},\br{}''; \omega)}{\abs*{\br{}' - \br{}''}} \dbr{}'' ,
|
|
\end{align}
|
|
\end{subequations}
|
|
with $\epsilon_{\IS}$ the dielectric function at coupling constant $\IS$ and $\chi_{0}$ the non-interacting polarizability. In the occupied-to-virtual orbital product basis, the spectral representation of $\W{}{\IS}$ can be written as follows in the case of real spatial orbitals
|
|
\begin{multline}
|
|
\label{eq:W}
|
|
\W{ij,ab}{\IS}(\omega) = \ERI{ij}{ab} + 2 \sum_m^{\Nocc \Nvir} \sERI{ij}{m} \sERI{ab}{m}
|
|
\\
|
|
\times \qty(\frac{1}{\omega - \OmRPA{m}{\IS} + i \eta} - \frac{1}{\omega + \OmRPA{m}{\IS} - i \eta}),
|
|
\end{multline}
|
|
where the spectral weights at coupling strength $\IS$ read
|
|
\begin{equation}
|
|
\sERI{pq}{m} = \sum_i^{\Nocc} \sum_a^{\Nvir} \ERI{pq}{ia} (\bX{\IS}_m + \bY{\IS}_m)_{ia}.
|
|
\end{equation}
|
|
In the case of complex molecular orbitals, see Ref.~\onlinecite{Holzer_2019} for a correct use of complex conjugation in the spectral representation of $\W{}{}$.
|
|
|
|
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
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-A}
|
|
\ARPA{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eHF{a} - \eHF{i}) + 2 \IS \ERI{ia}{jb},
|
|
\\
|
|
\label{eq:LR_RPA-B}
|
|
\BRPA{ia,jb}{\IS} & = 2 \IS \ERI{ia}{bj},
|
|
\end{align}
|
|
\end{subequations}
|
|
where $\eHF{p}$ are the HF orbital energies.
|
|
|
|
The relationship between the BSE formalism and the well-known RPAx (\ie, RPA with exchange) approach can be obtained by switching off the screening
|
|
%namely setting $\epsilon_{\IS}({\bf r},{\bf r}'; \omega) = \delta({\bf r}-{\bf r}')$
|
|
so that $\W{}{\IS}$ reduces to the bare Coulomb potential $\vc{}$.
|
|
In this limit, the $GW$ quasiparticle energies reduce to the Hartree-Fock (HF) eigenvalues, and Eqs.~\eqref{eq:LR_BSE-A} and \eqref{eq:LR_BSE-B} to the RPAx equations:
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPAx-A}
|
|
\ARPAx{ia,jb}{\IS} & = \delta_{ij} \delta_{ab} (\eHF{a} - \eHF{i}) + \IS \qty[ 2 \ERI{ia}{jb} - \ERI{ij}{ab} ],
|
|
\\
|
|
\label{eq:LR_RPAx-B}
|
|
\BRPAx{ia,jb}{\IS} & = \IS \qty[ 2 \ERI{ia}{bj} - \ERI{ib}{aj} ].
|
|
\end{align}
|
|
\end{subequations}
|
|
%This allows to understand that the strength parameter $\IS$ enters twice in the $\IS W^{\IS}$ contribution, one time to renormalize the screening efficiency, and a second time to renormalize the direct electron-hole interaction.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\subsection{Ground-state BSE energy}
|
|
%\label{sec:BSE_energy}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
The key quantity to define in the present context is the total BSE ground-state energy $\EBSE$.
|
|
Although this choice is not unique, \cite{Holzer_2018} we propose here to define it as
|
|
\begin{equation}
|
|
\label{eq:EtotBSE}
|
|
\EBSE = \Enuc + \EHF + \EcBSE,
|
|
\end{equation}
|
|
where $\Enuc$ and $\EHF$ are the nuclear repulsion energy and ground-state HF energy (respectively), and
|
|
\begin{equation}
|
|
\label{eq:EcBSE}
|
|
\EcBSE = \frac{1}{2} \int_0^1 \Tr(\bK \bP{\IS}) d\IS
|
|
\end{equation}
|
|
is the ground-state BSE correlation energy computed in the adiabatic connection framework, where
|
|
\begin{equation}
|
|
\label{eq:K}
|
|
\bK =
|
|
\begin{pmatrix}
|
|
\btA{\IS=1} & \bB{\IS=1} \\
|
|
\bB{\IS=1} & \btA{\IS=1} \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
is the interaction kernel \cite{Angyan_2011, Holzer_2018} [with $\tA{ia,jb}{\IS} = \IS \ERI{ia}{bj}$],
|
|
\begin{equation}
|
|
\label{eq:2DM}
|
|
\bP{\IS} =
|
|
\begin{pmatrix}
|
|
\bY{\IS} \T{(\bY{\IS})} & \bY{\IS} \T{(\bX{\IS})} \\
|
|
\bX{\IS} \T{(\bY{\IS})} & \bX{\IS} \T{(\bX{\IS})} \\
|
|
\end{pmatrix}
|
|
-
|
|
\begin{pmatrix}
|
|
\bO & \bO \\
|
|
\bO & \bI \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
is the correlation part of the two-electron density matrix at interaction strength $\IS$, and $\Tr$ denotes the matrix trace.
|
|
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 been named ``XBS'' for ``extended Bethe Salpeter'' by Holzer \textit{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 maintained as $\IS$ varies.
|
|
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, \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 Eqs.~\eqref{eq:LR_RPA-A} and \eqref{eq:LR_RPA-B}, and their RPAx analogs in Eqs.~\eqref{eq:LR_RPAx-A} and \eqref{eq:LR_RPAx-B}.
|
|
In the following, we will refer to these two types of calculations as RPA@HF and RPAx@HF, respectively.
|
|
Finally, we will also consider the RPA@$GW$@HF scheme which consists in replacing the HF orbital energies in Eq.~\eqref{eq:LR_RPA-A} by the $GW$ quasiparticles energies.
|
|
|
|
Several important comments are in order here.
|
|
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. }
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\section{Computational details}
|
|
%\label{sec:comp_details}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
All the $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies have been done using a (restricted) 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.
|
|
In the case of {\GOWO}, the quasiparticle energies have been obtained by linearizing the non-linear, frequency-dependent quasiparticle equation.
|
|
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018, Veril_2018}.
|
|
Finally, the infinitesimal $\eta$ has been set to zero for all calculations.
|
|
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}
|
|
The computational cost of these methods, in their usual implementation, scale as $\order*{N^5}$, $\order*{N^5}$, $\order*{N^6}$, and $\order*{N^7}$, respectively.
|
|
All the other calculations have been performed with our locally developed $GW$ software. \cite{Loos_2018,Veril_2018}
|
|
As one-electron basis sets, we employ the Dunning family (cc-pVXZ) defined with cartesian gaussian functions.
|
|
Unless, otherwise stated, the frozen-core approximation has not been enforced in our calculations in order to provide a fair comparison between methods.
|
|
However, we have found that the conclusions drawn in the present study hold within the frozen-core approximation (see {\SI}).
|
|
|
|
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.
|
|
This step is, by far, the computational bottleneck in our current implementation.
|
|
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}
|
|
|
|
%%% FIG 1 %%%
|
|
\begin{figure*}
|
|
\includegraphics[width=0.49\linewidth]{H2_GS_VQZ}
|
|
\includegraphics[width=0.49\linewidth]{LiH_GS_VQZ}
|
|
\caption{
|
|
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-pVQZ 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*}
|
|
%%% %%% %%%
|
|
|
|
%%% FIG 2 %%%
|
|
\begin{figure*}
|
|
\includegraphics[height=0.35\linewidth]{LiF_GS_VQZ}
|
|
\includegraphics[height=0.35\linewidth]{HCl_GS_VQZ}
|
|
\caption{
|
|
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-pVQZ 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*}
|
|
%%% %%% %%%
|
|
|
|
%%% FIG 3 %%%
|
|
\begin{figure*}
|
|
\includegraphics[height=0.26\linewidth]{N2_GS_VQZ}
|
|
\includegraphics[height=0.26\linewidth]{CO_GS_VQZ}
|
|
\includegraphics[height=0.26\linewidth]{BF_GS_VQZ}
|
|
\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-pVQZ 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*}
|
|
%%% %%% %%%
|
|
|
|
%%% FIG 4 %%%
|
|
\begin{figure}
|
|
\includegraphics[width=\linewidth]{F2_GS_VQZ}
|
|
\caption{
|
|
Ground-state PES of \ce{F2} around its equilibrium geometry obtained at various levels of theory with the cc-pVQZ basis set.
|
|
%Additional graphs for other basis sets and within the frozen-core approximation can be found in the {\SI}.
|
|
\label{fig:PES-F2}
|
|
}
|
|
\end{figure}
|
|
%%% %%% %%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\section{Potential energy surfaces}
|
|
%\label{sec:PES}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
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{HCl}, \ce{N2}, \ce{CO}, \ce{BF}, and \ce{F2}.
|
|
The PES of these molecules for various methods 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 for various basis sets are gathered in Table \ref{tab:Req}.
|
|
Additional graphs for other basis sets can be found in the {\SI}.
|
|
|
|
%%% TABLE I %%%
|
|
\begin{table*}
|
|
\caption{
|
|
Equilibrium bond length (in bohr) of the ground state of diatomic molecules obtained at various levels of theory and basis sets.
|
|
The reference CC3 and corresponding BSE@{\GOWO}@HF data are highlighted in bold black and bold red for visual convenience, respectively.
|
|
The values in parenthesis have been obtained by fitting a Morse potential to the PES.
|
|
}
|
|
\label{tab:Req}
|
|
\begin{ruledtabular}
|
|
\begin{tabular}{llcccccccc}
|
|
& & \mc{8}{c}{Molecules} \\
|
|
\cline{3-10}
|
|
Method & Basis & \ce{H2} & \ce{LiH} & \ce{LiF} & \ce{HCl} & \ce{N2} & \ce{CO} & \ce{BF} & \ce{F2} \\
|
|
\hline
|
|
CC3 & cc-pVDZ & 1.438 & 3.043 & 3.012 & 2.435 & 2.114 & 2.166 & 2.444 & 2.740 \\
|
|
& cc-pVTZ & 1.403 & 3.011 & 2.961 & 2.413 & 2.079 & 2.143 & 2.392 & 2.669 \\
|
|
& cc-pVQZ &\bb{1.402} &\bb{3.019} &\bb{2.963} &\bb{2.403} &\bb{2.075} &\bb{2.136} &\bb{2.390} &\bb{2.663} \\
|
|
CCSD & cc-pVDZ & 1.438 & 3.044 & 3.006 & 2.433 & 2.101 & 2.149 & 2.435 & 2.695 \\
|
|
& cc-pVTZ & 1.403 & 3.012 & 2.954 & 2.409 & 2.064 & 2.126 & 2.382 & 2.629 \\
|
|
& cc-pVQZ & 1.402 & 3.020 & 2.953 & 2.398 & 2.059 & 2.118 & 2.118 & 2.621 \\
|
|
CC2 & cc-pVDZ & 1.426 & 3.046 & 3.026 & 2.427 & 2.146 & 2.187 & 2.444 & 2.710 \\
|
|
& cc-pVTZ & 1.393 & 3.008 & 2.995 & 2.406 & 2.109 & 2.163 & 2.394 & 2.664 \\
|
|
& cc-pVQZ & 1.391 & 2.989 & 2.982 & 2.396 & 2.106 & 2.156 & 2.393 & 2.665 \\
|
|
MP2 & cc-pVDZ & 1.426 & 3.041 & 3.010 & 2.426 & 2.133 & 2.166 & 2.431 & 2.681 \\
|
|
& cc-pVTZ & 1.393 & 3.004 & 2.968 & 2.405 & 2.095 & 2.144 & 2.383 & 2.636 \\
|
|
& cc-pVQZ & 1.391 & 3.008 & 2.970 & 2.395 & 2.091 & 2.137 & 2.382 & 2.634 \\
|
|
BSE@{\GOWO}@HF & cc-pVDZ & 1.437 & 3.042 & 3.000 & 2.454 & 2.107 & 2.153 & 2.407 & (2.698) \\
|
|
& cc-pVTZ & 1.404 & 3.023 & (2.982) & 2.410 & 2.068 & 2.116 & (2.389) & (2.647) \\
|
|
& cc-pVQZ &\rb{1.399} &\rb{3.017} &\rb{(2.974)} &\gb{(2.408)} &\gb{(2.070)} &\gb{(2.130)} &\gb{(2.383)} &\rb{(2.640)}\\
|
|
RPA@{\GOWO}@HF & cc-pVDZ & 1.426 & 3.019 & 2.994 & 2.436 & 2.083 & 2.144 & 2.403 & (2.629) \\
|
|
& cc-pVTZ & 1.388 & 2.988 & (2.965) & 2.408 & 2.055 & 2.114 & (2.370) & (2.584) \\
|
|
& cc-pVQZ & 1.382 & 2.997 & (2.965) &\gb{(2.389)} &\gb{(2.045)} &\gb{(2.110)} &\gb{(2.367)} & (2.571) \\
|
|
RPAx@HF & cc-pVDZ & 1.428 & 3.040 & 2.998 & 2.424 & 2.077 & 2.130 & 2.417 & 2.611 \\
|
|
& cc-pVTZ & 1.395 & 3.003 & 2.943 & 2.400 & 2.046 & 2.110 & 2.368 & 2.568 \\
|
|
& cc-pVQZ & 1.394 & 3.011 & 2.944 &\gb{(2.393)} &\gb{(2.041)} &\gb{(2.105)} &\gb{(2.367)} &\gb{(2.563)} \\
|
|
RPA@HF & cc-pVDZ & 1.431 & 3.021 & 2.999 & 2.424 & 2.083 & 2.134 & 2.416 & 2.623 \\
|
|
& cc-pVTZ & 1.388 & 2.978 & 2.939 & 2.396 & 2.045 & 2.110 & 2.362 & 2.579 \\
|
|
& cc-pVQZ & 1.386 & 2.994 & 2.946 &\gb{(2.385)} &\gb{(2.042)} &\gb{(2.104)} &\gb{(2.365)} &\gb{(2.571)} \\
|
|
\end{tabular}
|
|
\end{ruledtabular}
|
|
\end{table*}
|
|
|
|
Let us start with the two smallest molecules, \ce{H2} and \ce{LiH}, which are both held by covalent bonds.
|
|
Their corresponding PES computed with the cc-pVQZ basis are reported in Fig.~\ref{fig:PES-H2-LiH}.
|
|
For \ce{H2}, we take as reference the full configuration interaction (FCI) energies \cite{QP2} 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 underestimate the FCI energy, while RPAx@HF and BSE@{\GOWO}@HF slightly over and undershoot 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.
|
|
|
|
Albeit the shallow nature of the \ce{LiH} PES, the scenario is almost identical for \ce{LiH} for which we report the CC2, CCSD and CC3 energies in addition to MP2.
|
|
In this case, RPAx@HF and BSE@{\GOWO}@HF nestle the CCSD and CC3 energy curves, and they are almost perfectly parallel.
|
|
Here again, the BSE@{\GOWO}@HF equilibrium bond length (obtained with cc-pVQZ) is extremely accurate ($3.017$ bohr) as compared to FCI ($3.019$ bohr).
|
|
|
|
The cases of \ce{LiF} and \ce{HCl} (see Fig.~\ref{fig:PES-LiF-HCl}) are interesting as they corresponds to strongly polarized bonds towards the halogen atoms which are much more electronegative than the first row elements.
|
|
For these ionic bonds, 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 bound for much longer in the case of \ce{HCl}.
|
|
Maybe surprisingly, BSE@{\GOWO}@HF is on par with both CC2 and CCSD, and outperforms 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.
|
|
As observed in Refs.~\onlinecite{vanSetten_2015,Maggio_2017,Loos_2018} and explained in details in Refs.~\onlinecite{Veril_2018,Duchemin_2020}, these irregularities, which makes particularly tricky the location of the minima, are due to ``jumps'' between distinct solutions of the $GW$ quasiparticle equation.
|
|
Including a broadening via the increasing the value of $\eta$ in the $GW$ self-energy and the screened Coulomb operator soften the problem, but does not remove it completely.
|
|
Note that these irregularities would be genuine discontinuities in the case of {\evGW}. \cite{Veril_2018}
|
|
|
|
|
|
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}, and systematically outperforms both CC2 and CCSD.
|
|
|
|
The \ce{F2} molecule is a notoriously difficult case to treat due to the weakness of its covalent bond (see Fig.~\ref{fig:PES-F2}), hence its relatively long equilibrium bond length ($2.663$ bohr at the CC3/cc-pVQZ level).
|
|
Similarly to what we have observed for \ce{LiF} and \ce{BF}, there is an irregularities near the minimum of the {\GOWO}-based curves.
|
|
However, BSE@{\GOWO}@HF is the closest to the CC3 curve
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\section{Conclusion}
|
|
%\label{sec:conclusion}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
In this Letter, we hope to have illustrated that the ACFDT@BSE formalism is a promising methodology for the computation of accurate ground-state PES and their corresponding equilibrium structures.
|
|
To do so, we have shown that calculating the BSE correlation energy computed within the ACFDT framework yield extremely accurate PES around equilibrium.
|
|
(Their accuracy near the dissociation limit remains an open question.)
|
|
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}
|
|
We believe that this central issue must be resolved if one wants to expand the applicability of the present method.
|
|
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.
|
|
We hope to be able to report on this in the near future.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\begin{acknowledgements}
|
|
%PFL would like to thank Julien Toulouse for enlightening discussions about RPA, and
|
|
XB is indebted to Valerio Olevano for numerous discussions.
|
|
This work was performed using HPC resources from GENCI-TGCC (Grant No.~2018-A0040801738) and CALMIP (Toulouse) under allocation 2019-18005.
|
|
Funding from the \textit{``Centre National de la Recherche Scientifique''} is acknowledged.
|
|
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''.}
|
|
\end{acknowledgements}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section*{Supporting Information}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
See {\SI} for additional potential energy curves with other basis sets and within the frozen-core approximation.
|
|
|
|
\bibliography{BSE-PES,BSE-PES-control}
|
|
|
|
\end{document}
|