diff --git a/BSE-PES.tex b/BSE-PES.tex index c9a9ecf..9f0dcb3 100644 --- a/BSE-PES.tex +++ b/BSE-PES.tex @@ -150,18 +150,18 @@ \newcommand{\InAA}[1]{#1 \AA} \newcommand{\kcal}{kcal/mol} -\newcommand{\NEEL}{Univ. Grenoble Alpes, CNRS, Institut NEEL, F-38042 Grenoble, France} +\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}{ Univ. Grenoble Alpes, CEA, IRIG-MEM-L Sim, 38054 Grenoble, France } +\newcommand{\CEA}{Universit\'e Grenoble Alpes, CEA, IRIG-MEM-L Sim, 38054 Grenoble, France} \begin{document} \title{Ground-State Potential Energy Surfaces Within the Bethe-Salpeter Formalism: Pros and Cons} -\author{Xavier \surname{Blase}} -\email{xavier.blase@neel.cnrs.fr } -\affiliation{\NEEL} +\author{Pierre-Fran\c{c}ois \surname{Loos}} +\email{loos@irsamc.ups-tlse.fr} +\affiliation{\LCPQ} \author{Ivan \surname{Duchemin}} \email{ivan.duchemin@cea.fr} \affiliation{\CEA} @@ -171,9 +171,9 @@ \author{Denis \surname{Jacquemin}} \email{denis.jacquemin@univ-nantes.fr} \affiliation{\CEISAM} -\author{Pierre-Fran\c{c}ois \surname{Loos}} -\email{loos@irsamc.ups-tlse.fr} -\affiliation{\LCPQ} +\author{Xavier \surname{Blase}} +\email{xavier.blase@neel.cnrs.fr } +\affiliation{\NEEL} \begin{abstract} %\begin{wrapfigure}[12]{o}[-1.25cm]{0.4\linewidth} @@ -182,7 +182,7 @@ %\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 distance. +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} @@ -195,28 +195,28 @@ However, we sometimes observe unphysical irregularities on the ground-state PES %%%%%%%%%%%%%%%%%%%%%%%% 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} +\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,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 exchange-correlation 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} +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 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} +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 KS-DFT (LDA) forces as its ground-state analog. +\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 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. +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 (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} +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 leads to two different versions of the BSE correlation energy. \cite{Holzer_2018} +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 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 and equilibrium distances. +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. @@ -299,12 +299,12 @@ In the absence of instabilities (\ie, $\bA{\IS} - \bB{\IS}$ is positive-definite 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}, \\ - \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} -With Mulliken's notation of the bare two-electron integrals +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} @@ -324,7 +324,7 @@ the BSE matrix elements read \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, the screened Coulomb potential $\W{}{\IS}$ is built within the direct RPA scheme: +In the standard BSE approach, $\W{}{\IS}$ is built within the direct RPA scheme, \ie, \begin{subequations} \label{eq:wrpa} \begin{align} @@ -335,7 +335,7 @@ In the standard BSE approach, the screened Coulomb potential $\W{}{\IS}$ is buil & = \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 molecular orbitals product basis, the spectral representation of $\W{}{\IS}$ can be written as follows in the case of real spatial orbitals \footnote{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$.} +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 \footnote{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$.} \begin{multline} \label{eq:W} \W{ij,ab}{\IS}(\omega) = \ERI{ij}{ab} + 2 \sum_m^{\Nocc \Nvir} \sERI{ij}{m} \sERI{ab}{m} @@ -362,7 +362,7 @@ 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 eigenvalues, and Eqs.~\eqref{eq:LR_BSE-A} and \eqref{eq:LR_BSE-B} to the RPAx equations: +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} @@ -437,7 +437,7 @@ Even for weakly correlated systems, triplet instabilities are much more common, %\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) Hartree-Fock (HF) starting point, which is a very adequate choice in the case of the (small) systems that we have considered here. +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}. @@ -520,7 +520,6 @@ The reference CC3 and corresponding BSE@{\GOWO}@HF data are highlighted in bold 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} \\ @@ -540,45 +539,17 @@ The values in parenthesis have been obtained by fitting a Morse potential to the & 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.074) & 2.116 & (2.389) & (2.647) \\ + & 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.065(2.048) & 2.114 & (2.370) & (2.584) \\ - & cc-pVQZ & 1.382 &3.013(2.998) & (2.965) &\gb{(2.389)} &\gb{(2.045)} &\gb{(2.110)} &\gb{(2.367)} & (2.571) \\ + & 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 &\gb{(2.943)} &\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 &\gb{(2.946)} &\gb{(2.385)} &\gb{(2.042)} &\gb{(2.104)} &\gb{(2.365)} &\gb{(2.571)} \\ -% FROZEN CORE VERSION -% Method & Basis & \ce{H2} & \ce{LiH}& \ce{LiF}& \ce{N2} & \ce{CO} & \ce{BF} & \ce{F2} & \ce{HCl}\\ -% \hline -% CC3 & cc-pVDZ & 1.438 & 3.052 & 3.014 & 2.115 & 2.167 & 2.447 & 2.741 & 2.438 \\ -% & cc-pVTZ & 1.403 & 3.036 & 2.985 & 2.087 & 2.150 & 2.405 & 2.672 & 2.414 \\ -% & cc-pVQZ & 1.402 & 3.037 & 2.985 & 2.080 & 2.142 & 2.398 & 2.667 & 2.413 \\ -% CCSD & cc-pVDZ & 1.438 & 3.044 & 3.006 & 2.101 & 2.149 & 2.435 & 2.695 & 2.433 \\ -% & cc-pVTZ & 1.403 & 3.012 & 2.954 & 2.064 & 2.126 & 2.382 & 2.629 & 2.409 \\ -% & cc-pVQZ & 1.402 & 3.020 & 2.953 & 2.059 & 2.118 & 2.380 & 2.621 & 2.398 \\ -% CC2 & cc-pVDZ & 1.426 & & & & & & & \\ -% & cc-pVTZ & 1.393 & & & & & & & \\ -% & cc-pVQZ & 1.391 & & & & & & & \\ -% MP2 & cc-pVDZ & 1.426 & 3.049 & 3.012 & 2.134 & 2.167 & 2.433 & 2.681 & 2.429 \\ -% & cc-pVTZ & 1.393 & 3.026 & 2.990 & 2.104 & 2.151 & 2.395 & 2.640 & 2.407 \\ -% & cc-pVQZ & 1.391 & 3.026 & 2.990 & 2.098 & 2.144 & 2.389 & 2.638 & 2.405 \\ -% BSE@{\GOWO}@HF & cc-pVDZ & 1.437 & & & & & & & \\ -% & cc-pVTZ & 1.404 & & & & & & & \\ -% & cc-pVQZ & 1.399 & & & & & & & \\ -% RPA@{\GOWO}@HF & cc-pVDZ & 1.426 & & & & & & & \\ -% & cc-pVTZ & 1.388 & & & & & & & \\ -% & cc-pVQZ & 1.382 & & & & & & & \\ -% RPAx@HF & cc-pVDZ & 1.428 & & & & & & & \\ -% & cc-pVTZ & 1.395 & & & & & & & \\ -% & cc-pVQZ & 1.394 & & & & & & & \\ -% RPA@HF & cc-pVDZ & 1.431 & & & & & & & \\ -% & cc-pVTZ & 1.388 & & & & & & & \\ -% & cc-pVQZ & 1.386 & & & & & & & \\ - \end{tabular} \end{ruledtabular} \end{table*} @@ -632,7 +603,8 @@ See {\SI} for additional potential energy curves with other basis sets and withi %%%%%%%%%%%%%%%%%%%%%%%% \begin{acknowledgements} -PFL would like to thank Julien Toulouse for enlightening discussions about RPA, and XB is indebted to Valerio Olevano for numerous discussions. +%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''.}