1055 lines
68 KiB
TeX
1055 lines
68 KiB
TeX
\documentclass[aip,jcp,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}
|
|
\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}}
|
|
\usepackage[normalem]{ulem}
|
|
\newcommand{\titou}[1]{\textcolor{red}{#1}}
|
|
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
|
|
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
|
|
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #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}{\mathbf{r}}
|
|
\newcommand{\dbr}{d\br}
|
|
|
|
% 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{c}}
|
|
\newcommand{\x}{\text{x}}
|
|
|
|
%
|
|
\newcommand{\Norb}{N_\text{orb}}
|
|
\newcommand{\Nocc}{O}
|
|
\newcommand{\Nvir}{V}
|
|
|
|
% operators
|
|
\newcommand{\hH}{\Hat{H}}
|
|
\newcommand{\hS}{\Hat{S}}
|
|
|
|
% methods
|
|
\newcommand{\KS}{\text{KS}}
|
|
\newcommand{\HF}{\text{HF}}
|
|
\newcommand{\RPA}{\text{RPA}}
|
|
\newcommand{\BSE}{\text{BSE}}
|
|
\newcommand{\dBSE}{\text{dBSE}}
|
|
\newcommand{\GW}{GW}
|
|
\newcommand{\stat}{\text{stat}}
|
|
\newcommand{\dyn}{\text{dyn}}
|
|
\newcommand{\TDA}{\text{TDA}}
|
|
|
|
% 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{\EcBSE}{E_\text{c}^\text{BSE}}
|
|
|
|
% orbital energies
|
|
\newcommand{\e}[1]{\eps_{#1}}
|
|
\newcommand{\eHF}[1]{\eps^\text{HF}_{#1}}
|
|
\newcommand{\eKS}[1]{\eps^\text{KS}_{#1}}
|
|
\newcommand{\eQP}[1]{\eps^\text{QP}_{#1}}
|
|
\newcommand{\eGOWO}[1]{\eps^\text{\GOWO}_{#1}}
|
|
\newcommand{\eGW}[1]{\eps^{GW}_{#1}}
|
|
\newcommand{\eevGW}[1]{\eps^\text{\evGW}_{#1}}
|
|
\newcommand{\eGnWn}[2]{\eps^\text{\GnWn{#2}}_{#1}}
|
|
\newcommand{\Om}[2]{\Omega_{#1}^{#2}}
|
|
\newcommand{\tOm}[2]{\Tilde{\Omega}_{#1}^{#2}}
|
|
|
|
\newcommand{\homu}{\frac{{\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{\tW}[2]{\widetilde{W}_{#1}^{#2}}
|
|
\newcommand{\Wc}[1]{W^\text{c}_{#1}}
|
|
\newcommand{\vc}[1]{v_{#1}}
|
|
\newcommand{\Sig}[2]{\Sigma_{#1}^{#2}}
|
|
\newcommand{\SigC}[1]{\Sigma^\text{c}_{#1}}
|
|
\newcommand{\SigX}[1]{\Sigma^\text{x}_{#1}}
|
|
\newcommand{\SigXC}[1]{\Sigma^\text{xc}_{#1}}
|
|
\newcommand{\Z}[1]{Z_{#1}}
|
|
\newcommand{\MO}[1]{\phi_{#1}}
|
|
\newcommand{\ERI}[2]{(#1|#2)}
|
|
\newcommand{\rbra}[1]{(#1|}
|
|
\newcommand{\rket}[1]{|#1)}
|
|
\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}[1]{\Omega_{#1}^{\text{RPA}}}
|
|
\newcommand{\OmRPAx}[1]{\Omega_{#1}^{\text{RPAx}}}
|
|
\newcommand{\OmBSE}[1]{\Omega_{#1}^{\text{BSE}}}
|
|
|
|
|
|
% 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}[2]{\mathbf{A}_{#1}^{#2}}
|
|
\newcommand{\bB}[2]{\mathbf{B}_{#1}^{#2}}
|
|
\newcommand{\bX}[2]{\mathbf{X}_{#1}^{#2}}
|
|
\newcommand{\bY}[2]{\mathbf{Y}_{#1}^{#2}}
|
|
\newcommand{\bZ}[2]{\mathbf{Z}_{#1}^{#2}}
|
|
\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}
|
|
|
|
% orbitals, gaps, etc
|
|
\newcommand{\eps}{\varepsilon}
|
|
\newcommand{\IP}{I}
|
|
\newcommand{\EA}{A}
|
|
\newcommand{\HOMO}{\text{HOMO}}
|
|
\newcommand{\LUMO}{\text{LUMO}}
|
|
\newcommand{\Eg}{E_\text{g}}
|
|
\newcommand{\EgFun}{\Eg^\text{fund}}
|
|
\newcommand{\EgOpt}{\Eg^\text{opt}}
|
|
\newcommand{\EB}{E_B}
|
|
|
|
\newcommand{\sig}{\sigma}
|
|
\newcommand{\bsig}{{\Bar{\sigma}}}
|
|
\newcommand{\sigp}{{\sigma'}}
|
|
\newcommand{\bsigp}{{\Bar{\sigma}'}}
|
|
\newcommand{\taup}{{\tau'}}
|
|
|
|
\newcommand{\up}{\uparrow}
|
|
\newcommand{\dw}{\downarrow}
|
|
\newcommand{\upup}{\uparrow\uparrow}
|
|
\newcommand{\updw}{\uparrow\downarrow}
|
|
\newcommand{\dwup}{\downarrow\uparrow}
|
|
\newcommand{\dwdw}{\downarrow\downarrow}
|
|
\newcommand{\spc}{\text{sc}}
|
|
\newcommand{\spf}{\text{sf}}
|
|
|
|
% addresses
|
|
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
|
|
|
|
\begin{document}
|
|
|
|
\title{Spin-Conserved and Spin-Flip Optical Excitations From the Bethe-Salpeter Equation Formalism}
|
|
|
|
\author{Enzo \surname{Monino}}
|
|
\affiliation{\LCPQ}
|
|
\author{Pierre-Fran\c{c}ois \surname{Loos}}
|
|
\email{loos@irsamc.ups-tlse.fr}
|
|
\affiliation{\LCPQ}
|
|
|
|
\begin{abstract}
|
|
Like adiabatic time-dependent density-functional theory (TD-DFT), the Bethe-Salpeter equation (BSE) formalism of many-body perturbation theory, in its static approximation, is ``blind'' to double (and higher) excitations, which are, for example, ubiquitous in conjugated molecules like polyenes.
|
|
Here, we apply the spin-flip \textit{ansatz} (which considers the lowest triplet state as the reference configuration instead of the singlet ground state) to the BSE formalism in order to access double excitations.
|
|
The present scheme is based on a spin-unrestricted version of the $GW$ approximation employed to compute the charged excitations and screened Coulomb potential required for the BSE calculations.
|
|
Dynamical corrections to the static BSE optical excitations are taken into account via an unrestricted generalization of our recently developed (renormalized) perturbative treatment.
|
|
The performance of the present spin-flip BSE formalism is illustrated by computing the vertical excitation energies of the beryllium atom, the hydrogen molecule at various bond lengths, and cyclobutadiene in its rectangular and square-planar geometries.
|
|
%\bigskip
|
|
%\begin{center}
|
|
% \boxed{\includegraphics[width=0.5\linewidth]{TOC}}
|
|
%\end{center}
|
|
%\bigskip
|
|
\end{abstract}
|
|
|
|
\maketitle
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Introduction}
|
|
\label{sec:intro}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
Due to the ubiquitous influence of processes involving electronic excited states in physics, chemistry, and biology, their faithful description from first principles has been one of the grand challenges faced by theoretical chemists since the dawn of computational chemistry.
|
|
Accurately predicting ground- and excited-state energies (hence excitation energies) is particularly valuable in this context, and it has concentrated most of the efforts within the community.
|
|
An armada of theoretical and computational methods have been developed to this end, each of them being plagued by its own flaws. \cite{Roos_1996,Piecuch_2002,Dreuw_2005,Krylov_2006,Sneskov_2012,Gonzales_2012,Laurent_2013,Adamo_2013,Ghosh_2018,Blase_2020,Loos_2020d,Casanova_2020}
|
|
The fact that none of these methods is successful in every chemical scenario has encouraged chemists to carry on the development of new excited-state methodologies, their main goal being to get the most accurate excitation energies (and properties) at the lowest possible computational cost in the most general context. \cite{Loos_2020d}
|
|
|
|
Originally developed in the framework of nuclear physics, \cite{Salpeter_1951} and popularized in condensed-matter physics, \cite{Sham_1966,Strinati_1984,Delerue_2000} one of the new emerging method in the computational chemistry landscape is the Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988,Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999,Blase_2018,Blase_2020} from many-body perturbation theory \cite{Onida_2002,Martin_2016} which, based on an underlying $GW$ calculation to compute accurate charged excitations and the dynamically-screened Coulomb potential, \cite{Hedin_1965,Golze_2019} is able to provide accurate optical (\ie, neutral) excitations for molecular systems at a rather modest computational cost.\cite{Rohlfing_1999a,Horst_1999,Puschnig_2002,Tiago_2003,Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018,Blase_2018,Liu_2020,Blase_2020,Holzer_2018a,Holzer_2018b,Loos_2020e}
|
|
Most of BSE implementations rely on the so-called static approximation, \cite{Blase_2018,Bruneval_2016,Krause_2017,Liu_2020} which approximates the dynamical (\ie, frequency-dependent) BSE kernel by its static limit.
|
|
Like adiabatic time-dependent density-functional theory (TD-DFT), \cite{Runge_1984,Casida_1995,Petersilka_1996,UlrichBook} the static BSE formalism is plagued by the lack of double (and higher) excitations, which are, for example, ubiquitous in conjugated molecules like polyenes \cite{Maitra_2004,Cave_2004,Saha_2006,Watson_2012,Shu_2017,Barca_2018a,Barca_2018b,Loos_2019} or the ground state of open-shell molecules. \cite{Casida_2005,Huix-Rotllant_2011,Loos_2020f}
|
|
Indeed, both adiabatic TD-DFT \cite{Levine_2006,Tozer_2000,Elliott_2011,Maitra_2012,Maitra_2016} and static BSE \cite{ReiningBook,Romaniello_2009b,Sangalli_2011,Loos_2020h,Authier_2020} can only access (singlet and triplet) single excitations with respect to the reference determinant usually taken as the closed-shell singlet ground state.
|
|
Double excitations are even challenging for state-of-the-art methods, \cite{Loos_2018a,Loos_2019,Loos_2020c,Loos_2020d,Veril_2020} like the approximate third-order coupled-cluster (CC3) method \cite{Christiansen_1995b,Koch_1997} or equation-of-motion coupled-cluster with singles, doubles and triples (EOM-CCSDT). \cite{Kucharski_1991,Kallay_2004,Hirata_2000,Hirata_2004}
|
|
|
|
One way to access double excitations is via the spin-flip formalism established by Krylov in 2001, \cite{Krylov_2001a,Krylov_2001b,Krylov_2002} with earlier attempts by Bethe, \cite{Bethe_1931} as well as Shibuya and McKoy. \cite{Shibuya_1970}
|
|
The idea behind the spin-flip \textit{ansatz} is rather simple: instead of considering the singlet ground state as reference, the reference is taken as the lowest triplet state.
|
|
In such a way, one can access the singlet ground state and the singlet doubly-excited state via a spin-flip deexcitation and excitation (respectively), the difference of these two excitation energies providing an estimate of the double excitation.
|
|
We refer the interested reader to Refs.~\onlinecite{Krylov_2006,Krylov_2008,Casanova_2020} for detailed reviews on spin-flip methods.
|
|
Note that a similar idea has been exploited by the group of Weito Yang to access double excitations in the context of the particle-particle random-phase approximation. \cite{Peng_2013,Yang_2013b,Yang_2014a,Peng_2014,Zhang_2016,Sutton_2018}
|
|
|
|
One obvious issue of spin-flip methods is that not all double excitations are accessible in such a way.
|
|
Moreover, spin-flip methods are usually hampered by spin contamination \cite{Casanova_2020} (\ie, artificial mixing with configurations of different spin multiplicities) due to spin incompleteness of the configuration interaction expansion as well as the possible spin contamination of the reference configuration. \cite{Krylov_2000b}
|
|
This issue can be alleviated by increasing the excitation order at a significant cost or by selectively complementing the spin-incomplete configuration set with the missing configurations. \cite{Sears_2003,Casanova_2008,Huix-Rotllant_2010,Li_2010,Li_2011a,Li_2011b,Zhang_2015,Lee_2018}
|
|
|
|
Nowadays, spin-flip techniques are widely available for many types of methods such as equation-of-motion coupled cluster (EOM-CC), \cite{Krylov_2001a,Levchenko_2004,Manohar_2008,Casanova_2009a,Dutta_2013} configuration interaction (CI), \cite{Krylov_2001b,Krylov_2002,Mato_2018,Casanova_2008,Casanova_2009b} TD-DFT, \cite{Shao_2003,Wang_2004,Li_2011a,Bernard_2012,Zhang_2015} the algebraic-diagrammatic construction (ADC) scheme,\cite{Lefrancois_2015,Lefrancois_2016} and others \cite{Mayhall_2014a,Mayhall_2014b,Bell_2013,Mayhall_2014c} with successful applications in bond breaking processes, \cite{Golubeva_2007} radical chemistry, \cite{Slipchenko_2002,Wang_2005,Slipchenko_2003,Rinkevicius_2010,Ibeji_2015,Hossain_2017,Orms_2018,Luxon_2018} and photochemistry in general \cite{Casanova_2012,Gozem_2013,Nikiforov_2014,Lefrancois_2016} to mention a few.
|
|
|
|
Here we apply the spin-flip technique to the BSE formalism in order to access, in particular, double excitations, \cite{Authier_2020} but not only.
|
|
The present BSE calculations are based on the spin-unrestricted version of both $GW$ (Sec.~\ref{sec:UGW}) and BSE (Sec.~\ref{sec:UBSE}).
|
|
To the best of our knowledge, the present study is the first to apply the spin-flip formalism to the BSE method.
|
|
Moreover, we also go beyond the static approximation by taking into account dynamical effects (Sec.~\ref{sec:dBSE}) via an unrestricted generalization of our recently developed (renormalized) perturbative correction which builds on the seminal work of Strinati, \cite{Strinati_1982,Strinati_1984,Strinati_1988} Romaniello and collaborators, \cite{Romaniello_2009b,Sangalli_2011} and Rohlfing and coworkers. \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b,Lettmann_2019}
|
|
We also discuss the computation of oscillator strengths (Sec.~\ref{sec:os}) and the expectation value of the spin operator $\expval*{\hS^2}$ as a diagnostic of the spin contamination for both ground and excited states (Sec.~\ref{sec:spin}).
|
|
Computational details are reported in Sec.~\ref{sec:compdet} and our results for the beryllium atom \ce{Be} (Subsec.~\ref{sec:Be}), the hydrogen molecule \ce{H2} (Subsec.~\ref{sec:H2}), and cyclobutadiene \ce{C4H4} (Subsec.~\ref{sec:CBD}) are discussed in Sec.~\ref{sec:res}.
|
|
Finally, we draw our conclusions in Sec.~\ref{sec:ccl}.
|
|
Unless otherwise stated, atomic units are used.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Unrestricted $GW$ formalism}
|
|
\label{sec:UGW}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
Let us consider an electronic system consisting of $n = n_\up + n_\dw$ electrons (where $n_\up$ and $n_\dw$ are the number of spin-up and spin-down electrons, respectively) and $N$ one-electron basis functions.
|
|
The number of spin-up and spin-down occupied orbitals are $O_\up = n_\up$ and $O_\dw = n_\dw$, respectively, and, assuming the absence of linear dependencies in the one-electron basis set, there is $V_\up = N - O_\up$ and $V_\dw = N - O_\dw$ spin-up and spin-down virtual (\ie, unoccupied) orbitals.
|
|
The number of spin-conserved (sc) single excitations is then $S^\spc = S_{\up\up}^\spc + S_{\dw\dw}^\spc = O_\up V_\up + O_\dw V_\dw$, while the number of spin-flip (sf) excitations is $S^\spf = S_{\up\dw}^\spf + S_{\dw\up}^\spf = O_\up V_\dw + O_\dw V_\up$.
|
|
Let us denote as $\MO{p_\sig}(\br)$ the $p$th (spin)orbital with spin $\sig$ (where $\sig =$ $\up$ or $\dw$) and $\e{p_\sig}{}$ its one-electron energy.
|
|
It is important to understand that, in a spin-conserved excitation the hole orbital $\MO{i_\sig}$ and particle orbital $\MO{a_\sig}$ have the same spin $\sig$.
|
|
In a spin-flip excitation, the hole and particle states, $\MO{i_\sig}$ and $\MO{a_\bsig}$, have opposite spins, $\sig$ and $\bsig$.
|
|
We assume real quantities throughout this manuscript, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, $p$, $q$, $r$, and $s$ indicate arbitrary orbitals, and $m$ labels single excitations.
|
|
Moreover, we consider systems with collinear spins and a spin-independent Hamiltonian without contributions such as spin-orbit interaction.
|
|
|
|
%================================
|
|
\subsection{The dynamical screening}
|
|
%================================
|
|
The pillar of Green's function many-body perturbation theory is the (time-ordered) one-body Green's function, which has poles at the charged excitations (i.e., ionization potentials and electron affinities) of the system. \cite{ReiningBook}
|
|
The spin-$\sig$ component of the one-body Green's function reads \cite{ReiningBook,Bruneval_2016a}
|
|
\begin{equation}
|
|
\label{eq:G}
|
|
G^{\sig}(\br_1,\br_2;\omega)
|
|
= \sum_i \frac{\MO{i_\sig}(\br_1) \MO{i_\sig}(\br_2)}{\omega - \e{i_\sig}{} - i\eta}
|
|
+ \sum_a \frac{\MO{a_\sig}(\br_1) \MO{a_\sig}(\br_2)}{\omega - \e{a_\sig}{} + i\eta}
|
|
\end{equation}
|
|
where $\eta$ is a positive infinitesimal.
|
|
As readily seen in Eq.~\eqref{eq:G}, the Green's function can be evaluated at different levels of theory depending on the choice of orbitals and energies, $\MO{p_\sig}$ and $\e{p_\sig}{}$.
|
|
For example, $G_{\KS}^{\sig}$ is the independent-particle Green's function built with Kohn-Sham (KS) orbitals $\MO{p_\sig}^{\KS}(\br)$ and one-electron energies $\e{p_\sig}^{\KS}$. \cite{Hohenberg_1964,Kohn_1965,ParrBook}
|
|
Within self-consistent schemes, these quantities can be replaced by quasiparticle energies and orbitals evaluated within the $GW$ approximation (see below). \cite{Hedin_1965,Golze_2019}
|
|
|
|
Based on the spin-up and spin-down components of $G$ defined in Eq.~\eqref{eq:G}, one can easily compute the non-interacting polarizability (which is a sum over spins)
|
|
\begin{equation}
|
|
\label{eq:chi0}
|
|
\chi_0(\br_1,\br_2;\omega) = - \frac{i}{2\pi} \sum_\sig \int G^{\sig}(\br_1,\br_2;\omega+\omega') G^{\sig}(\br_1,\br_2;\omega') d\omega'
|
|
\end{equation}
|
|
and subsequently the dielectric function
|
|
\begin{equation}
|
|
\label{eq:eps}
|
|
\epsilon(\br_1,\br_2;\omega) = \delta(\br_1 - \br_2) - \int \frac{\chi_0(\br_1,\br_3;\omega) }{\abs{\br_2 - \br_3}} d\br_3
|
|
\end{equation}
|
|
where $\delta(\br)$ is the Dirac delta function.
|
|
Based on this latter ingredient, one can access the dynamically-screened Coulomb potential
|
|
\begin{equation}
|
|
\label{eq:W}
|
|
W(\br_1,\br_2;\omega) = \int \frac{\epsilon^{-1}(\br_1,\br_3;\omega) }{\abs{\br_2 - \br_3}} d\br_3
|
|
\end{equation}
|
|
which is naturally spin independent as the bare Coulomb interaction $\abs{\br_1 - \br_2}^{-1}$ does not depend on spin coordinates.
|
|
|
|
Within the $GW$ formalism, \cite{Hedin_1965,Onida_2002,Golze_2019} the dynamical screening is computed at the random-phase approximation (RPA) level by considering only the manifold of the spin-conserved neutral excitations.
|
|
In the orbital basis, the spectral representation of $W$ is
|
|
\begin{multline}
|
|
\label{eq:W_spectral}
|
|
W_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega) = \ERI{p_\sig q_\sig}{r_\sigp s_\sigp}
|
|
+ \sum_m \ERI{p_\sig q_\sig}{m}\ERI{r_\sigp s_\sigp}{m}
|
|
\\
|
|
\times \qty[ \frac{1}{\omega - \Om{m}{\spc,\RPA} + i \eta} - \frac{1}{\omega + \Om{m}{\spc,\RPA} - i \eta} ]
|
|
\end{multline}
|
|
where the bare two-electron integrals are \cite{Gill_1994}
|
|
\begin{equation}
|
|
\ERI{p_\sig q_\tau}{r_\sigp s_\taup} = \int \frac{\MO{p_\sig}(\br_1) \MO{q_\tau}(\br_1) \MO{r_\sigp}(\br_2) \MO{s_\taup}(\br_2)}{\abs{\br_1 - \br_2}} d\br_1 d\br_2
|
|
\end{equation}
|
|
and the screened two-electron integrals (or spectral weights) are explicitly given by
|
|
\begin{equation}
|
|
\label{eq:sERI}
|
|
\ERI{p_\sig q_\sig}{m} = \sum_{ia\sigp} \ERI{p_\sig q_\sig}{r_\sigp s_\sigp} (\bX{m}{\spc,\RPA}+\bY{m}{\spc,\RPA})_{i_\sigp a_\sigp}
|
|
\end{equation}
|
|
In Eqs.~\eqref{eq:W_spectral} and \eqref{eq:sERI}, the spin-conserved RPA neutral excitations $\Om{m}{\spc,\RPA}$ and their corresponding eigenvectors, $\bX{m}{\spc,\RPA}$ and $\bY{m}{\spc,\RPA}$, are obtained by solving a linear response system of the form
|
|
\begin{equation}
|
|
\label{eq:LR-RPA}
|
|
\begin{pmatrix}
|
|
\bA{}{} & \bB{}{} \\
|
|
-\bB{}{} & -\bA{}{} \\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX{m}{} \\
|
|
\bY{m}{} \\
|
|
\end{pmatrix}
|
|
=
|
|
\Om{m}{}
|
|
\begin{pmatrix}
|
|
\bX{m}{} \\
|
|
\bY{m}{} \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
where the expressions of the matrix elements of $\bA{}{}$ and $\bB{}{}$ are specific of the method and of the spin manifold.
|
|
The spin structure of these matrices, though, is general
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR-RPA-AB-sc}
|
|
\bA{}{\spc} & = \begin{pmatrix}
|
|
\bA{}{\upup,\upup} & \bA{}{\upup,\dwdw} \\
|
|
\bA{}{\dwdw,\upup} & \bA{}{\dwdw,\dwdw} \\
|
|
\end{pmatrix}
|
|
&
|
|
\bB{}{\spc} & = \begin{pmatrix}
|
|
\bB{}{\upup,\upup} & \bB{}{\upup,\dwdw} \\
|
|
\bB{}{\dwdw,\upup} & \bB{}{\dwdw,\dwdw} \\
|
|
\end{pmatrix}
|
|
\\
|
|
\label{eq:LR-RPA-AB-sf}
|
|
\bA{}{\spf} & = \begin{pmatrix}
|
|
\bA{}{\updw,\updw} & \bO \\
|
|
\bO & \bA{}{\dwup,\dwup} \\
|
|
\end{pmatrix}
|
|
&
|
|
\bB{}{\spf} & = \begin{pmatrix}
|
|
\bO & \bB{}{\updw,\dwup} \\
|
|
\bB{}{\dwup,\updw} & \bO \\
|
|
\end{pmatrix}
|
|
\end{align}
|
|
\end{subequations}
|
|
In the absence of instabilities, the linear eigenvalue problem \eqref{eq:LR-RPA} has particle-hole symmetry which means that the eigenvalues are obtained by pairs $\pm \Om{m}{}$.
|
|
In such a case, $(\bA{}{}-\bB{}{})^{1/2}$ is positive definite, and Eq.~\eqref{eq:LR-RPA} can be recast as a Hermitian problem of half its original dimension
|
|
\begin{equation}
|
|
\label{eq:small-LR}
|
|
(\bA{}{} - \bB{}{})^{1/2} \cdot (\bA{}{} + \bB{}{}) \cdot (\bA{}{} - \bB{}{})^{1/2} \cdot \bZ{}{} = \bOm{2} \cdot \bZ{}{}
|
|
\end{equation}
|
|
where the excitation amplitudes are
|
|
\begin{equation}
|
|
\bX{}{} + \bY{}{} = \bOm{-1/2} \cdot (\bA{}{} - \bB{}{})^{1/2} \cdot \bZ{}{}
|
|
\end{equation}
|
|
Within the Tamm-Dancoff approximation (TDA), the coupling terms between the resonant and anti-resonant parts, $\bA{}{}$ and $-\bA{}{}$, are neglected, which consists in setting $\bB{}{} = \bO$.
|
|
In such a case, Eq.~\eqref{eq:LR-RPA} reduces to a straightforward Hermitian problem of the form:
|
|
\begin{equation}
|
|
\bA{}{} \cdot \bX{m}{} = \Om{m}{} \bX{m}{}
|
|
\end{equation}
|
|
Note that, for spin-flip excitations, it is quite common to enforce the TDA especially when one considers a triplet reference as the first ``excited-state'' is usually the ground state of the closed-shell system (hence, corresponding to a negative excitation energy).
|
|
|
|
At the RPA level, the matrix elements of $\bA{}{}$ and $\bB{}{}$ are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-A}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} & = \delta_{ij} \delta_{ab} \delta_{\sig \sigp} \delta_{\tau \taup} (\e{a_\tau} - \e{i_\sig}) + \ERI{i_\sig a_\tau}{b_\sigp j_\taup}
|
|
\\
|
|
\label{eq:LR_RPA-B}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} & = \ERI{i_\sig a_\tau}{j_\sigp b_\taup}
|
|
\end{align}
|
|
\end{subequations}
|
|
from which we obtain the following expressions
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-Asc}
|
|
\A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} & = \delta_{ij} \delta_{ab} \delta_{\sig \sigp} (\e{a_\sig} - \e{i_\sig}) + \ERI{i_\sig a_\sig}{b_\sigp j_\sigp}
|
|
\\
|
|
\label{eq:LR_RPA-Bsc}
|
|
\B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} & = \ERI{i_\sig a_\sig}{j_\sigp b_\sigp}
|
|
\end{align}
|
|
\end{subequations}
|
|
for the spin-conserved excitations and
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_RPA-Asf}
|
|
\A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\RPA} & = \delta_{ij} \delta_{ab} (\e{a_\bsig} - \e{i_\sig})
|
|
\\
|
|
\label{eq:LR_RPA-Bsf}
|
|
\B{i_\sig a_\bsig,j_\bsig b_\sig}{\spf,\RPA} & = 0
|
|
\end{align}
|
|
\end{subequations}
|
|
for the spin-flip excitations.
|
|
|
|
%================================
|
|
\subsection{The $GW$ self-energy}
|
|
%================================
|
|
Within the acclaimed $GW$ approximation, \cite{Hedin_1965,Golze_2019} the exchange-correlation (xc) part of the self-energy
|
|
\begin{equation}
|
|
\label{eq:Sig}
|
|
\begin{split}
|
|
\Sig{}{\xc,\sig}(\br_1,\br_2;\omega)
|
|
& = \Sig{}{\x,\sig}(\br_1,\br_2) + \Sig{}{\co,\sig}(\br_1,\br_2;\omega)
|
|
\\
|
|
& = \frac{i}{2\pi} \int G^{\sig}(\br_1,\br_2;\omega+\omega') W(\br_1,\br_2;\omega') e^{i \eta \omega'} d\omega'
|
|
\end{split}
|
|
\end{equation}
|
|
is, like the one-body Green's function, spin-diagonal, and its spectral representation reads
|
|
\begin{subequations}
|
|
\begin{gather}
|
|
\Sig{p_\sig q_\sig}{\x}
|
|
= - \sum_{i} \ERI{p_\sig i_\sig}{i_\sig q_\sig}
|
|
\\
|
|
\begin{split}
|
|
\Sig{p_\sig q_\sig}{\co}(\omega)
|
|
& = \sum_{im} \frac{\ERI{p_\sig i_\sig}{m} \ERI{q_\sig i_\sig}{m}}{\omega - \e{i_\sig} + \Om{m}{\spc,\RPA} - i \eta}
|
|
\\
|
|
& + \sum_{am} \frac{\ERI{p_\sig a_\sig}{m} \ERI{q_\sig a_\sig}{m}}{\omega - \e{a_\sig} - \Om{m}{\spc,\RPA} + i \eta}
|
|
\end{split}
|
|
\end{gather}
|
|
\end{subequations}
|
|
where the self-energy has been split in its exchange (x) and correlation (c) contributions.
|
|
The Dyson equation linking the Green's function and the self-energy holds separately for each spin component
|
|
\begin{equation}
|
|
\label{eq:Dyson_G}
|
|
\begin{split}
|
|
\qty[ G^{\sig}(\br_1,\br_2;\omega) ]^{-1}
|
|
& = \qty[ G_{\KS}^{\sig}(\br_1,\br_2;\omega) ]^{-1}
|
|
\\
|
|
& + \Sig{}{\xc,\sig}(\br_1,\br_2;\omega) - v^{\xc}(\br_1) \delta(\br_1 - \br_2)
|
|
\end{split}
|
|
\end{equation}
|
|
where $v^{\xc}(\br)$ is the KS (local) exchange-correlation potential.
|
|
The target quantities here are the quasiparticle energies $\eGW{p_\sig}$, \ie, the poles of $G$ [see Eq.~\eqref{eq:G}], which correspond to well-defined addition/removal energies (unlike the KS orbital energies).
|
|
Because the exchange-correlation part of the self-energy is, itself, constructed with the Green's function [see Eq.~\eqref{eq:Sig}], the present process is, by nature, self-consistent.
|
|
The same comment applies to the dynamically-screened Coulomb potential $W$ entering the definition of $\Sig{}{\xc}$ [see Eq.~\eqref{eq:Sig}] which is also constructed from $G$ [see Eqs.~\eqref{eq:chi0}, \eqref{eq:eps}, and \eqref{eq:W}].
|
|
|
|
%================================
|
|
\subsection{Level of self-consistency}
|
|
%================================
|
|
This is where $GW$ schemes differ.
|
|
In its simplest perturbative (\ie, one-shot) version, known as {\GOWO}, \cite{Strinati_1980,Hybertsen_1985a,Hybertsen_1986,Godby_1988,Linden_1988,Northrup_1991,Blase_1994,Rohlfing_1995,Shishkin_2007} a single iteration is performed, and the quasiparticle energies $\eGW{p_\sig}$ are obtained by solving the frequency-dependent quasiparticle equation
|
|
\begin{equation}
|
|
\label{eq:QP-eq}
|
|
\omega = \e{p_\sig}{} + \Sig{p_\sig}{\xc}(\omega) - V_{p_\sig}^{\xc}
|
|
\end{equation}
|
|
where $\Sig{p_\sig}{\xc}(\omega) \equiv \Sig{p_\sig p_\sig}{\xc}(\omega)$ and its offspring quantities have been constructed at the KS level, and
|
|
\begin{equation}
|
|
V_{p_\sigma}^{\xc} = \int \MO{p_\sig}(\br) v^{\xc}(\br) \MO{p_\sig}(\br) d\br
|
|
\end{equation}
|
|
Because, from a practical point of view, one is usually interested by the so-called quasiparticle solution (or peak), the quasiparticle equation \eqref{eq:QP-eq} is often linearized around $\omega = \e{p_\sig}^{\KS}$, yielding
|
|
\begin{equation}
|
|
\label{eq:G0W0_lin}
|
|
\eGW{p_\sig}
|
|
= \e{p_\sig}^{\KS} + Z_{p_\sig} [\Sig{p_\sig}{\xc}(\e{p_\sig}^{\KS}) - V_{p_\sig}^{\xc} ]
|
|
\end{equation}
|
|
where
|
|
\begin{equation}
|
|
\label{eq:Z_GW}
|
|
Z_{p_\sig} = \qty[ 1 - \left. \pdv{\Sig{p_\sig}{\xc}(\omega)}{\omega} \right|_{\omega = \e{p_\sig}^{\KS}} ]^{-1}
|
|
\end{equation}
|
|
is a renormalization factor (with $0 \le Z_{p_\sig} \le 1$) which also represents the spectral weight of the quasiparticle solution.
|
|
In addition to the principal quasiparticle peak which, in a well-behaved case, contains most of the spectral weight, the frequency-dependent quasiparticle equation \eqref{eq:QP-eq} generates a finite number of satellite resonances with smaller weights. \cite{Loos_2018b}
|
|
|
|
Within the ``eigenvalue'' self-consistent $GW$ scheme (known as ev$GW$), \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016,Gui_2018} several iterations are performed during which only the one-electron energies entering the definition of the Green's function [see Eq.~\eqref{eq:G}] are updated by the quasiparticle energies obtained at the previous iteration (the corresponding orbitals remain evaluated at the KS level).
|
|
|
|
Finally, within the quasiparticle self-consistent $GW$ (qs$GW$) scheme, \cite{Faleev_2004,vanSchilfgaarde_2006,Kotani_2007,Ke_2011,Kaplan_2016} both the one-electron energies and the orbitals are updated until convergence is reached.
|
|
These are obtained via the diagonalization of an effective Fock matrix which includes explicitly a frequency-independent and Hermitian self-energy defined as
|
|
\begin{equation}
|
|
\Tilde{\Sigma}_{p_\sig q_\sig}^{\xc} = \frac{1}{2} \qty[ \Sig{p_\sig q_\sig}{\xc}(\e{p_\sig}{}) + \Sig{q_\sig p_\sig}{\xc}(\e{p_\sig}{}) ]
|
|
\end{equation}
|
|
|
|
%================================
|
|
\section{Unrestricted Bethe-Salpeter equation formalism}
|
|
\label{sec:UBSE}
|
|
%================================
|
|
Like its TD-DFT cousin, \cite{Runge_1984,Casida_1995,Petersilka_1996,Dreuw_2005} the BSE formalism \cite{Salpeter_1951,Strinati_1988,Albrecht_1998,Rohlfing_1998,Benedict_1998,vanderHorst_1999} deals with the calculation of (neutral) optical excitations as measured by absorption spectroscopy. \cite{Rohlfing_1999a,Horst_1999,Puschnig_2002,Tiago_2003,Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018}
|
|
Using the BSE formalism, one can access the spin-conserved and spin-flip excitations.
|
|
In a nutshell, BSE builds on top of a $GW$ calculation by adding up excitonic effects (\ie, the electron-hole binding energy) to the $GW$ fundamental gap which is itself a corrected version of the KS gap.
|
|
The purpose of the underlying $GW$ calculation is to provide quasiparticle energies and a dynamically-screened Coulomb potential that are used to build the BSE Hamiltonian from which the vertical excitations of the system are extracted.
|
|
|
|
%================================
|
|
\subsection{Static approximation}
|
|
\label{sec:BSE}
|
|
%================================
|
|
|
|
Within the so-called static approximation of BSE, the Dyson equation that links the generalized four-point susceptibility $L^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)$ and the BSE kernel $\Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6)$ is \cite{ReiningBook,Bruneval_2016a}
|
|
\begin{multline}
|
|
L^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
|
|
= L_{0}^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
|
|
\\
|
|
+ \int L_{0}^{\sig\sigp}(\br_1,\br_4;\br_1',\br_3;\omega)
|
|
\Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6)
|
|
\\
|
|
\times L^{\sig\sigp}(\br_6,\br_2;\br_5,\br_2';\omega)
|
|
d\br_3 d\br_4 d\br_5 d\br_6
|
|
\end{multline}
|
|
where
|
|
\begin{multline}
|
|
L_{0}^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
|
|
\\
|
|
= \frac{1}{2\pi} \int G^{\sig}(\br_1,\br_2';\omega+\omega') G^{\sig}(\br_1',\br_2;\omega') d\omega'
|
|
\end{multline}
|
|
is the non-interacting analog of the two-particle correlation function $L$.
|
|
|
|
Within the $GW$ approximation, the static BSE kernel is
|
|
\begin{multline}
|
|
\label{eq:kernel}
|
|
i \Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6)
|
|
= \frac{\delta(\br_3 - \br_4) \delta(\br_5 - \br_6) }{\abs{\br_3-\br_6}}
|
|
\\
|
|
- \delta_{\sig\sigp} W(\br_3,\br_4;\omega = 0) \delta(\br_3 - \br_6) \delta(\br_4 - \br_6)
|
|
\end{multline}
|
|
where, as usual, we have not considered the higher-order terms in $W$ by neglecting the derivative $\partial W/\partial G$. \cite{Hanke_1980,Strinati_1982,Strinati_1984,Strinati_1988}
|
|
|
|
As readily seen in Eq.~\eqref{eq:kernel}, the static approximation consists in neglecting the frequency dependence of the dynamically-screened Coulomb potential.
|
|
In this case, the spin-conserved and spin-flip BSE optical excitations are obtained by solving the usual Casida-like linear response (eigen)problem:
|
|
\begin{equation}
|
|
\label{eq:LR-BSE}
|
|
\begin{pmatrix}
|
|
\bA{}{\BSE} & \bB{}{\BSE} \\
|
|
-\bB{}{\BSE} & -\bA{}{\BSE} \\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX{m}{\BSE} \\
|
|
\bY{m}{\BSE} \\
|
|
\end{pmatrix}
|
|
=
|
|
\Om{m}{\BSE}
|
|
\begin{pmatrix}
|
|
\bX{m}{\BSE} \\
|
|
\bY{m}{\BSE} \\
|
|
\end{pmatrix}
|
|
\end{equation}
|
|
Defining the elements of the static screening as $W^{\stat}_{p_\sig q_\sig,r_\sigp s_\sigp} = W_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega = 0)$, the general expressions of the BSE matrix elements are
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_BSE-A}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{\BSE} & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\taup a_\tau}
|
|
\\
|
|
\label{eq:LR_BSE-B}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{\BSE} & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\taup,j_\sigp a_\tau}
|
|
\end{align}
|
|
\end{subequations}
|
|
from which we obtain the following expressions for the spin-conserved and spin-flip BSE excitations:
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_BSE-Asc}
|
|
\A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\BSE} & = \A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\sigp a_\sig}
|
|
\\
|
|
\label{eq:LR_BSE-Bsc}
|
|
\B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\BSE} & = \B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\sigp,j_\sigp a_\sig}
|
|
\\
|
|
\label{eq:LR_BSE-Asf}
|
|
\A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\BSE} & = \A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\RPA} - W^{\stat}_{i_\sig j_\sig,b_\bsig a_\bsig}
|
|
\\
|
|
\label{eq:LR_BSE-Bsf}
|
|
\B{i_\sig a_\bsig,j_\bsig b_\sig}{\spf,\BSE} & = - W^{\stat}_{i_\sig b_\sig,j_\bsig a_\bsig}
|
|
\end{align}
|
|
\end{subequations}
|
|
|
|
At this stage, it is of particular interest to discuss the form of the spin-flip matrix elements defined in Eqs.~\eqref{eq:LR_BSE-Asf} and \eqref{eq:LR_BSE-Bsf}.
|
|
As readily seen from Eq.~\eqref{eq:LR_RPA-Asf}, at the RPA level, the spin-flip excitations are given by the difference of one-electron energies, hence missing out on key exchange and correlation effects.
|
|
This is also the case at the TD-DFT level when one relies on (semi-)local functionals.
|
|
This explains why most of spin-flip TD-DFT calculations are performed with hybrid functionals containing a substantial amount of Hartree-Fock exchange as only the exact exchange integral of the form $\ERI{i_\sig j_\sig}{b_\bsig a_\bsig}$ survive spin-symmetry requirements.
|
|
At the BSE level, these matrix elements are, of course, also present thanks to the contribution of $W^{\stat}_{i_\sig j_\sig,b_\bsig a_\bsig}$ as evidenced in Eq.~\eqref{eq:W_spectral} but it also includes correlation effects.
|
|
|
|
%================================
|
|
\subsection{Dynamical correction}
|
|
\label{sec:dBSE}
|
|
%================================
|
|
In order to go beyond the ubiquitous static approximation of BSE \cite{Strinati_1988,Rohlfing_2000,Sottile_2003,Myohanen_2008,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Huix-Rotllant_2011,Sakkinen_2012,Zhang_2013,Rebolini_2016,Olevano_2019,Lettmann_2019} (which is somehow similar to the adiabatic approximation of TD-DFT \cite{Casida_2005,Huix-Rotllant_2011,Casida_2016,Maitra_2004,Cave_2004,Elliott_2011,Maitra_2012}), we have recently implemented, following Strinati's seminal work \cite{Strinati_1982,Strinati_1984,Strinati_1988} (see also the work of Romaniello \textit{et al.} \cite{Romaniello_2009b} and Sangalli \textit{et al.} \cite{Sangalli_2011}), a renormalized first-order perturbative correction in order to take into consideration the dynamical nature of the screened Coulomb potential $W$. \cite{Loos_2020h,Authier_2020}
|
|
This dynamical correction to the static BSE kernel (dubbed as dBSE in the following) does permit to recover additional relaxation effects coming from higher excitations.
|
|
|
|
Our implementation follows closely the work of Rohlfing and co-workers \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b} in which they computed the dynamical correction in the TDA and plasmon-pole approximation.
|
|
However, our scheme goes beyond the plasmon-pole approximation as the spectral representation of the dynamically-screened Coulomb potential is computed exactly at the RPA level consistently with the underlying $GW$ calculation:
|
|
\begin{multline}
|
|
\widetilde{W}_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega) = \ERI{p_\sig q_\sig}{r_\sigp s_\sigp}
|
|
+ \sum_m \ERI{p_\sig q_\sig}{m}\ERI{r_\sigp s_\sigp}{m}
|
|
\\
|
|
\times \Bigg[ \frac{1}{\omega - (\e{s_\sigp}{} - \e{q_\sig}{}) - \Om{m}{\spc,\RPA} + i \eta}
|
|
\\
|
|
+ \frac{1}{\omega - (\e{r_\sigp}{} - \e{p_\sig}{}) - \Om{m}{\spc,\RPA} + i \eta} \Bigg]
|
|
\end{multline}
|
|
The dBSE non-linear response problem is
|
|
\begin{multline}
|
|
\label{eq:LR-dyn}
|
|
\begin{pmatrix}
|
|
\bA{}{\dBSE}(\Om{m}{\dBSE}) & \bB{}{\dBSE}(\Om{m}{\dBSE})
|
|
\\
|
|
-\bB{}{\dBSE}(-\Om{m}{\dBSE}) & -\bA{}{\dBSE}(-\Om{m}{\dBSE})
|
|
\\
|
|
\end{pmatrix}
|
|
\cdot
|
|
\begin{pmatrix}
|
|
\bX{m}{\dBSE} \\
|
|
\bY{m}{\dBSE} \\
|
|
\end{pmatrix}
|
|
\\
|
|
=
|
|
\Om{m}{\dBSE}
|
|
\begin{pmatrix}
|
|
\bX{m}{\dBSE} \\
|
|
\bY{m}{\dBSE} \\
|
|
\end{pmatrix}
|
|
\end{multline}
|
|
where the dynamical matrices are generally defined as
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:LR_dBSE-A}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{\dBSE}(\omega) & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} \widetilde{W}_{i_\sig j_\sigp,b_\taup a_\tau}(\omega)
|
|
\\
|
|
\label{eq:LR_dBSE-B}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{\dBSE}(\omega) & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} \widetilde{W}_{i_\sig b_\taup,j_\sigp a_\tau}(\omega)
|
|
\end{align}
|
|
\end{subequations}
|
|
from which one can easily obtained the matrix elements for the spin-conserved and spin-flip manifolds similarly to Eqs.~\eqref{eq:LR_BSE-Asc}, \eqref{eq:LR_BSE-Bsc}, \eqref{eq:LR_BSE-Asf}, and \eqref{eq:LR_BSE-Bsf}.
|
|
Following Rayleigh-Schr\"odinger perturbation theory, we then decompose the non-linear eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static (\ie, linear) reference and a first-order dynamic (\ie, non-linear) perturbation such that
|
|
\begin{multline}
|
|
\label{eq:LR-PT}
|
|
\begin{pmatrix}
|
|
\bA{}{\dBSE}(\omega) & \bB{}{\dBSE}(\omega) \\
|
|
-\bB{}{\dBSE}(-\omega) & -\bA{}{\dBSE}(-\omega) \\
|
|
\end{pmatrix}
|
|
\\
|
|
=
|
|
\begin{pmatrix}
|
|
\bA{}{(0)} & \bB{}{(0)}
|
|
\\
|
|
-\bB{}{(0)} & -\bA{}{(0)}
|
|
\\
|
|
\end{pmatrix}
|
|
+
|
|
\begin{pmatrix}
|
|
\bA{}{(1)}(\omega) & \bB{}{(1)}(\omega) \\
|
|
-\bB{}{(1)}(-\omega) & -\bA{}{(1)}(-\omega) \\
|
|
\end{pmatrix}
|
|
\end{multline}
|
|
with
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:BSE-A0}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{(0)} & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\BSE}
|
|
\\
|
|
\label{eq:BSE-B0}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{(0)} & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\BSE}
|
|
\end{align}
|
|
\end{subequations}
|
|
and
|
|
\begin{subequations}
|
|
\begin{align}
|
|
\label{eq:BSE-A1}
|
|
\A{i_\sig a_\tau,j_\sigp b_\taup}{(1)}(\omega) & = - \delta_{\sig \sigp} \widetilde{W}_{i_\sig j_\sigp,b_\taup a_\tau}(\omega) + \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\taup a_\tau}
|
|
\\
|
|
\label{eq:BSE-B1}
|
|
\B{i_\sig a_\tau,j_\sigp b_\taup}{(1)}(\omega) & = - \delta_{\sig \sigp} \widetilde{W}_{i_\sig b_\taup,j_\sigp a_\tau}(\omega) + \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\taup,j_\sigp a_\tau}
|
|
\end{align}
|
|
\end{subequations}
|
|
The dBSE excitation energies are then obtained via
|
|
\begin{equation}
|
|
\Om{m}{\dBSE} = \Om{m}{\BSE} + \zeta_{m} \Om{m}{(1)}
|
|
\end{equation}
|
|
where $\Om{m}{\BSE} \equiv \Om{m}{(0)}$ are the static (zeroth-order) BSE excitation energies obtained by solving Eq.~\eqref{eq:LR-BSE}, and
|
|
\begin{equation}
|
|
\label{eq:Om1-TDA}
|
|
\Om{m}{(1)} = \T{(\bX{m}{\BSE})} \cdot \bA{}{(1)}(\Om{m}{\BSE}) \cdot \bX{m}{\BSE}
|
|
\end{equation}
|
|
are first-order corrections (with $\bX{m}{\BSE} \equiv \bX{m}{(0)}$) obtained within the dynamical TDA (dTDA) with the renormalization factor
|
|
\begin{equation}
|
|
\label{eq:Z}
|
|
\zeta_{m} = \qty[ 1 - \T{(\bX{m}{\BSE})} \cdot \left. \pdv{\bA{}{(1)}(\omega)}{\omega} \right|_{\omega = \Om{m}{\BSE}} \cdot \bX{m}{\BSE} ]^{-1}
|
|
\end{equation}
|
|
which, unlike the $GW$ case [see Eq.~\eqref{eq:Z_GW}], is not restricted to be between $0$ and $1$.
|
|
In most cases, the value of $\zeta_{m}$ is close to unity which indicates that the perturbative expansion behaves nicely.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\subsection{Oscillator strengths}
|
|
\label{sec:os}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
Oscillator strengths, \ie, transition dipole moments from the ground to the corresponding excited state, are key quantities that are linked to experimental intensities and are usually used to probe the quality of excited-state calculations. \cite{Harbach_2014,Kannar_2014,Chrayteh_2021,Sarkar_2021}
|
|
|
|
For the spin-conserved transitions, the $x$ component of the transition dipole moment is
|
|
\begin{equation}
|
|
\mu_{x,m}^{\spc} = \sum_{ia\sig} (i_\sig|x|a_\sig)(\bX{m}{\spc}+\bY{m}{\spc})_{i_\sig a_\sig}
|
|
\end{equation}
|
|
where
|
|
\begin{equation}
|
|
(p_\sig|x|q_\sigp) = \int \MO{p_\sig}(\br) \, x \, \MO{q_\sigp}(\br) d\br
|
|
\end{equation}
|
|
are one-electron integrals in the orbital basis.
|
|
The total oscillator strength in the so-called length gauge \cite{Sarkar_2021} is given by
|
|
\begin{equation}
|
|
f_{m}^{\spc} = \frac{2}{3} \Om{m}{\spc} \qty[ \qty(\mu_{x,m}^{\spc})^2 + \qty(\mu_{x,m}^{\spc})^2 + \qty(\mu_{x,m}^{\spc})^2 ]
|
|
\end{equation}
|
|
For spin-flip transitions, we have $f_{m}^{\spf} = 0$ as the transition matrix elements $(i_\sig|x|a_\bsig)$ vanish via integration over the spin coordinate.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\subsection{Spin contamination}
|
|
\label{sec:spin}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
One of the key issues of linear response formalism based on unrestricted references is spin contamination or the artificial mixing with configurations of different spin multiplicities.
|
|
As nicely explained in Ref.~\onlinecite{Casanova_2020}, there are two sources of spin contamination: i) spin contamination of the reference configuration for which, for example, $\expval*{\hS^2} > 2$ for high-spin triplets, and ii) spin contamination of the excited states due to spin incompleteness of the CI expansion.
|
|
The latter issue is an important source of spin contamination in the present context as BSE is limited to single excitations with respect to the reference configuration.
|
|
Specific schemes have been developed to palliate these shortcomings and we refer the interested reader to Ref.~\onlinecite{Casanova_2020} for a detailed discussion on this matter.
|
|
|
|
In order to monitor closely how contaminated are these states, we compute
|
|
\begin{equation}
|
|
\expval*{\hS^2}_m = \expval*{\hS^2}_0 + \Delta \expval*{\hS^2}_m
|
|
\end{equation}
|
|
where
|
|
\begin{equation}
|
|
\expval*{\hS^2}_{0}
|
|
= \frac{n_{\up} - n_{\dw}}{2} \qty( \frac{n_{\up} - n_{\dw}}{2} + 1 )
|
|
+ n_{\dw} - \sum_p (p_{\up}|p_{\dw})^2
|
|
\end{equation}
|
|
is the expectation value of $\hS^2$ for the reference configuration, the first term corresponding to the exact value of $\expval*{\hS^2}$, and
|
|
\begin{equation}
|
|
\label{eq:OV}
|
|
(p_\sig|q_\sigp) = \int \MO{p_\sig}(\br) \MO{q_\sigp}(\br) d\br
|
|
\end{equation}
|
|
are overlap integrals between spin-up and spin-down orbitals.
|
|
|
|
For a given single excitation $m$, the explicit expressions of $\Delta \expval*{\hS^2}_m^{\spc}$ and $\Delta \expval*{\hS^2}_m^{\spf}$ can be found in the Appendix of Ref.~\onlinecite{Li_2011a} for spin-conserved and spin-flip excitations, and are functions of the vectors $\bX{m}{}$ and $\bY{m}{}$ as well as the orbital overlaps defined in Eq.~\eqref{eq:OV}.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Computational details}
|
|
\label{sec:compdet}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
All the systems under investigation have a closed-shell electronic structure and we consider the lowest triplet state as reference for the spin-flip calculations adopting the unrestricted formalism throughout this work.
|
|
The {\GOWO} calculations performed to obtain the screened Coulomb potential and the quasiparticle energies required to compute the BSE neutral excitations are performed using an unrestricted HF (UHF) starting point, and the {\GOWO} quasiparticle energies are obtained by linearizing the frequency-dependent quasiparticle equation [see Eq.~\eqref{eq:G0W0_lin}].
|
|
Note that the entire set of orbitals and energies is corrected.
|
|
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018b,Veril_2018,Loos_2020e,Loos_2020h,Berger_2021}.
|
|
|
|
Here, we do not investigate how the starting orbitals affect the BSE@{\GOWO} excitation energies.
|
|
This is left for future work.
|
|
However, it is worth mentioning that, for the present (small) molecular systems, HF is usually a good starting point, \cite{Loos_2020a,Loos_2020e,Loos_2020h} although improvements could certainly be obtained with starting orbitals and energies computed with, for example, optimally-tuned range-separated hybrid functionals. \cite{Stein_2009,Stein_2010,Refaely-Abramson_2012,Kronik_2012}
|
|
Besides, {\GOWO}@UHF and ev$GW$@UHF yield similar quasiparticle energies, while {\GOWO} allows us to avoid rather laborious iterations as well as the significant additional computational effort of ev$GW$. \cite{Loos_2020e,Loos_2020h,Berger_2021}
|
|
In the following, all linear response calculations are performed within the TDA to ensure consistency between the spin-conserved and spin-flip results.
|
|
Finally, the infinitesimal $\eta$ is set to $100$ meV for all calculations.
|
|
|
|
All the static and dynamic BSE calculations (labeled in the following as SF-BSE and SF-dBSE respectively) are performed with the software \texttt{QuAcK}, \cite{QuAcK} developed in our group and freely available on \texttt{github}.
|
|
The standard and extended spin-flip ADC(2) calculations [SF-ADC(2)-s and SF-ADC(2)-x, respectively] as well as the SF-ADC(3) \cite{Lefrancois_2015} are performed with Q-CHEM 5.2.1. \cite{qchem4}
|
|
Spin-flip TD-DFT calculations \cite{Shao_2003} considering the BLYP, \cite{Becke_1988,Lee_1988} B3LYP, \cite{Becke_1988,Lee_1988,Becke_1993a} and BH\&HLYP \cite{Lee_1988,Becke_1993b} functionals with contains $0\%$, $20\%$, and $50\%$ of exact exchange are labeled as SF-TD-BLYP, SF-TD-B3LYP, and SF-TD-BH\&HLYP, respectively, and are also performed with Q-CHEM 5.2.1.
|
|
EOM-CCSD excitation energies \cite{Koch_1990,Stanton_1993,Koch_1994} are computed with Gaussian 09. \cite{g09}
|
|
As a consistency check, we systematically perform SF-CIS calculations \cite{Krylov_2001a} with both \texttt{QuAcK} and Q-CHEM, and make sure that they yield identical excitation energies.
|
|
Throughout this work, all spin-flip and spin-conserved calculations are performed with a UHF reference.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Results}
|
|
\label{sec:res}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%===============================
|
|
\subsection{Beryllium atom}
|
|
\label{sec:Be}
|
|
%===============================
|
|
|
|
As a first example, we consider the simple case of the beryllium atom in a small basis (6-31G) which was considered by Krylov in two of her very first papers on spin-flip methods. \cite{Krylov_2001a,Krylov_2001b}
|
|
It was also considered in later studies thanks to its pedagogical value. \cite{Sears_2003,Casanova_2020}
|
|
Beryllium has a $^1S$ ground state with $1s^2 2s^2$ configuration.
|
|
The excitation energies corresponding to the first singlet and triplet single excitations $2s \to 2p$ with $P$ spatial symmetries as well as the first singlet and triplet double excitations $2s^2 \to 2p^2$ with $D$ and $P$ spatial symmetries (respectively) are reported in Table \ref{tab:Be} and depicted in Fig.~\ref{fig:Be}.
|
|
|
|
On the left side of Fig.~\ref{fig:Be}, we report SF-TD-DFT excitation energies (red lines) obtained with the BLYP, B3LYP, and BH\&HLYP functionals, which correspond to an increase of exact exchange from 0\% to 50\%.
|
|
As mentioned in Ref.~\onlinecite{Casanova_2020}, the $^3P(1s^2 2s^1 2p^1)$ and the $^1P(1s^2 2s^1 2p^1)$ states are degenerate at the SF-TD-BLYP level.
|
|
Indeed, due to the lack of coupling terms in the spin-flip block of the SD-TD-DFT equations (see Subsec.~\ref{sec:BSE}), their excitation energies are given by the energy difference between the $2s$ and $2p$ orbitals and both states are strongly spin contaminated.
|
|
Including exact exchange, like in SF-TD-B3LYP and SF-TD-BH\&HLYP, lifts this degeneracy and improves the description of both states.
|
|
However, the SF-TD-BH\&HLYP excitation energy of the $^1P(1s^2 2s^1 2p^1)$ state is still off by $1.6$ eV as compared to the FCI reference.
|
|
For the other states, the agreement between SF-TD-BH\&HLYP and FCI is significantly improved.
|
|
|
|
The center part of Fig.~\ref{fig:Be} shows the SF-(d)BSE results (blue lines) alongside the SF-CIS excitation energies (purple lines).
|
|
All of these are computed with 100\% of exact exchange with the additional inclusion of correlation in the case of SF-BSE and SF-dBSE thanks to the introduction of the static and dynamical screening, respectively.
|
|
Overall, the SF-CIS and SF-BSE excitation energies are closer to FCI than the SF-TD-DFT ones, except for the lowest triplet state where the SF-TD-BH\&HLYP excitation energy is more accurate probably due to error compensation.
|
|
At the exception of the $^1D$ state, SF-BSE improves over SF-CIS with a rather small contribution from the additional dynamical effect included in the SF-dBSE scheme.
|
|
Note that the exact exchange seems to spin purified the $^3P(1s^2 2s^1 2p^1)$ state while the singlet states at the SF-BSE level are slightly more spin contaminated than their SF-CIS counterparts.
|
|
|
|
The right side of Fig.~\ref{fig:Be} illustrates the performance of the SF-ADC methods.
|
|
Interestingly, SF-BSE and SF-ADC(2)-s have rather similar accuracies, except again for the $^1D$ state where SF-ADC(2)-s has clearly the edge over SF-BSE.
|
|
Finally, both SF-ADC(2)-x and SF-ADC(3) yield excitation energies very close to FCI for this simple system with significant improvements for the lowest $^3P$ state and the $^1D$ doubly-excited state.
|
|
|
|
%%% TABLE I %%%
|
|
%\begin{squeezetable}
|
|
\begin{table*}
|
|
\caption{
|
|
Excitation energies (in eV) with respect to the $^1S(1s^2 2s^2)$ singlet ground state of \ce{Be} obtained for various methods with the 6-31G basis set.
|
|
All the spin-flip calculations have been performed with a UHF reference.
|
|
The $\expval*{\hS^2}$ value associated with each state is reported in parenthesis (when available).
|
|
\label{tab:Be}}
|
|
\begin{ruledtabular}
|
|
\begin{tabular}{lcccccccccc}
|
|
& \mc{5}{c}{Excitation energies (eV)} \\
|
|
% & \mc{5}{c}{6-31G} & \mc{5}{c}{aug-cc-pVQZ} \\
|
|
\cline{2-6} %\cline{7-11}
|
|
Method & $^1S(1s^2 2s^2)$ & $^3P(1s^2 2s^1 2p^1)$ & $^1P(1s^2 2s^1 2p^1)$
|
|
& $^3P(1s^22 p^2)$ & $^1D(1s^22p^2)$ \\
|
|
% & $^1S(1s^2 2s^2)$ & $^3P(1s^2 2s^1 2p^1)$ & $^1P(1s^2 2s^1 2p^1)$
|
|
% & $^3P(1s^22 p^2)$ & $^1D(1s^22p^2)$ \\
|
|
\hline
|
|
SF-TD-BLYP\fnm[1] & (0.002) & 3.210(1.000) & 3.210(1.000) & 6.691(1.000) & 7.598(0.013) \\
|
|
SF-TD-B3LYP\fnm[1] & (0.001) & 3.332(1.839) & 4.275(0.164) & 6.864(1.000) & 7.762(0.006) \\
|
|
SF-TD-BH\&HLYP\fnm[1] & (0.000) & 2.874(1.981) & 4.922(0.023) & 7.112(1.000) & 8.188(0.002) \\
|
|
SF-CIS\fnm[2] & (0.002) & 2.111(2.000) & 6.036(0.014) & 7.480(1.000) & 8.945(0.006) \\
|
|
SF-BSE@{\GOWO} & (0.004) & 2.399(1.999) & 6.191(0.023) & 7.792(1.000) & 9.373(0.013) \\
|
|
% SF-BSE@ev$GW$ & (0.004) & 2.407(1.999) & 6.199(0.023) & 7.788(1.000) & 9.388(0.013) \\
|
|
SF-dBSE@{\GOWO} & & 2.363 & 6.263 & 7.824 & 9.424 \\
|
|
% SF-dBSE@ev$GW$ & & 2.369 & 6.273 & 7.820 & 9.441 \\
|
|
SF-ADC(2)-s & & 2.433 & 6.255 & 7.745 & 9.047 \\
|
|
SF-ADC(2)-x & & 2.866 & 6.581 & 7.664 & 8.612 \\
|
|
SF-ADC(3) & & 2.863 & 6.579 & 7.658 & 8.618 \\
|
|
FCI\fnm[2] & (0.000) & 2.862(2.000) & 6.577(0.000) & 7.669(2.000) & 8.624(0.000) \\
|
|
% & (0.000) & 2.718(2.000) & 5.277(0.000) & 6.450(2.000) & 7.114(0.000) \\
|
|
\end{tabular}
|
|
\end{ruledtabular}
|
|
\fnt[1]{Excitation energies taken from Ref.~\onlinecite{Casanova_2020}.}
|
|
\fnt[2]{Excitation energies taken from Ref.~\onlinecite{Krylov_2001a}.}
|
|
\end{table*}
|
|
%\end{squeezetable}
|
|
%%% %%% %%% %%%
|
|
|
|
%%% FIG 1 %%%
|
|
\begin{figure}
|
|
\includegraphics[width=\linewidth]{Be}
|
|
\caption{
|
|
Excitation energies (in eV) with respect to the $^1S(1s^2 2s^2)$ singlet ground state of \ce{Be} obtained with the 6-31G basis for various levels of theory:
|
|
SF-TD-DFT \cite{Casanova_2020} (red), SF-CIS \cite{Krylov_2001a} (purple), SF-BSE (blue), SF-ADC (orange), and FCI \cite{Krylov_2001a} (black).
|
|
All these spin-flip calculations have been performed with a UHF reference.
|
|
\label{fig:Be}}
|
|
\end{figure}
|
|
%%% %%% %%% %%%
|
|
|
|
%===============================
|
|
\subsection{Hydrogen molecule}
|
|
\label{sec:H2}
|
|
%===============================
|
|
|
|
Our second example deals with the dissociation of the \ce{H2} molecule, which is a prototypical system for testing new electronic structure methods and, specifically, their accuracy in the presence of strong correlation (see, for example, Refs.~\onlinecite{Caruso_2013,Barca_2014,Vuckovic_2017,Li_2021}, and references therein).
|
|
The $\text{X}\,{}^1 \Sigma_g^+$ ground state of \ce{H2} has an electronic configuration $(1\sigma_g)^2$ configuration.
|
|
The variation of the excitation energies associated with the three lowest singlet excited states with respect to the elongation of the \ce{H-H} bond are of particular interest here.
|
|
The lowest singly excited state $\text{B}\,{}^1 \Sigma_u^+$ has a $(1\sigma_g )(1\sigma_u)$ configuration, while the singly excited state $\text{E}\,{}^1 \Sigma_g^+$ and the doubly excited state $\text{F}\,{}^1 \Sigma_g^+$ have $(1\sigma_g ) (2\sigma_g)$ and $(1\sigma_u )(1\sigma_u)$ configurations, respectively.
|
|
Because these latter two excited states interact strongly and form an avoided crossing around $R(\ce{H-H}) = 1.4$ \AA, they are usually labeled as the $\text{EF}\,{}^1 \Sigma_g^+$ state.
|
|
Note that this avoided crossing is not visible with non-spin-flip methods restricted to single excitations (such as CIS, TD-DFT, and BSE) as these are ``blind'' to double excitations.
|
|
Three methods, in their standard and spin-flip versions, are studied here (CIS, TD-BH\&HLYP and BSE) and are compared to the reference EOM-CCSD excitation energies (that is equivalent to FCI in the case of \ce{H2}).
|
|
All these calculations are performed with the cc-pVQZ basis.
|
|
|
|
The top panel of Fig.~\ref{fig:H2} shows the CIS (dotted lines) and SF-CIS (dashed lines) excitation energies as a function of $R(\ce{H-H})$.
|
|
The EOM-CCSD reference energies are represented by solid lines.
|
|
We observe that both CIS and SF-CIS poorly describe the $\text{B}\,{}^1\Sigma_u^+$ state in the dissociation limit with an error greater than $1$ eV, while CIS, unlike SF-CIS, is much more accurate around the equilibrium geometry.
|
|
Similar observations can be made for the $\text{E}\,{}^1\Sigma_g^+$ state with a good description at the CIS level for all bond lengths.
|
|
SF-CIS does not model accurately the $\text{E}\,{}^1\Sigma_g^+$ state before the avoided crossing, but the agreement between SF-CIS and EOM-CCSD is much satisfactory for bond length greater than $1.6$ \AA.
|
|
Oppositely, SF-CIS describes better the $\text{F}\,{}^1\Sigma_g^+$ state before the avoided crossing than after, while this state is completely absent at the CIS level.
|
|
Indeed, as mentioned earlier, CIS is unable to locate any avoided crossing as it cannot access double excitations.
|
|
At the SF-CIS level, the avoided crossing between the $\text{E}$ and $\text{F}$ states is qualitatively reproduced and placed at a slightly larger bond length than at the EOM-CCSD level.
|
|
|
|
In the central panel of Fig.~\ref{fig:H2}, we report the (SF-)TD-BH\&HLYP results.
|
|
Similar graphs for (SF-)TD-BLYP and (SF-)TD-B3LYP are reported in the {\SI}.
|
|
SF-TD-BH\&HLYP shows, at best, qualitative agreement with EOM-CCSD, while the TD-BH\&HLYP excitation energies of the $\text{B}$ and $\text{E}$ states are only trustworthy around equilibrium but inaccurate at dissociation.
|
|
Note that \ce{H2} is a rather challenging system for (SF)-TD-DFT from a general point of view. \cite{Vuckovic_2017,Cohen_2008a,Cohen_2008c,Cohen_2012}
|
|
|
|
In the bottom panel of Fig.~\ref{fig:H2}, (SF-)BSE excitation energies for the same three singlet states are represented.
|
|
SF-BSE provides surprisingly accurate excitation energies for the $\text{B}\,{}^1\Sigma_u^+$ state with errors between $0.05$ and $0.3$ eV, outperforming in the process the standard BSE formalism.
|
|
However SF-BSE does not describe well the $\text{E}\,{}^1\Sigma_g^+$ state with error ranging from half an eV to $1.6$ eV.
|
|
Similar performances are observed at the BSE level around equilibrium with a clear improvement in the dissociation limit.
|
|
Remarkably, SF-BSE shows a good agreement with EOM-CCSD for the $\text{F}\,{}^1\Sigma_g^+$ doubly-excited state, resulting in an avoided crossing around $R(\ce{H-H}) = 1.6$ \AA.
|
|
A similar graph comparing (SF-)dBSE and EOM-CCSD excitation energies can be found in the {\SI} where it is shown that dynamical effects do not affect the present conclusions.
|
|
|
|
The right side of Fig.~\ref{fig:H2} shows the amount of spin contamination as a function of the bond length for SF-CIS (top), SF-TD-BH\&HLYP (center), and SF-BSE (bottom).
|
|
Overall, one can see that $\expval*{\hS^2}$ behaves similarly for SF-CIS and SF-BSE with a small spin contamination of the $\text{B}\,{}^1\Sigma_u^+$ at short bond length. In contrast, the $\text{B}$ state is much more spin contaminated at the SF-TD-BH\&HLYP level.
|
|
For all spin-flip methods, the $\text{E}$ is strongly spin contaminated as expected, while the $\expval*{\hS^2}$ values associated with the $\text{F}$ state
|
|
only deviate significantly from zero for short bond length and around the avoided crossing where it strongly couples with the spin contaminated $\text{E}$ state.
|
|
|
|
%%% FIG 2 %%%
|
|
\begin{figure*}
|
|
\includegraphics[width=0.4\linewidth]{H2_CIS}
|
|
\hspace{0.05\linewidth}
|
|
\includegraphics[width=0.4\linewidth]{H2_CIS_S2}
|
|
\vspace{0.025\linewidth}
|
|
\\
|
|
\includegraphics[width=0.4\linewidth]{H2_BHHLYP}
|
|
\hspace{0.05\linewidth}
|
|
\includegraphics[width=0.4\linewidth]{H2_BHHLYP_S2}
|
|
\vspace{0.025\linewidth}
|
|
\\
|
|
\includegraphics[width=0.4\linewidth]{H2_BSE}
|
|
\hspace{0.05\linewidth}
|
|
\includegraphics[width=0.4\linewidth]{H2_BSE_S2}
|
|
\caption{
|
|
Excitation energies with respect to the $\text{X}\,{}^1 \Sigma_g^+$ ground state (left) and expectation value of the spin operator $\expval*{\hS^2}$ (right) of the $\text{B}\,{}^1\Sigma_u^+$ (red), $\text{E}\,{}^1\Sigma_g^+$ (black), and $\text{F}\,{}^1\Sigma_g^+$ (blue) states of \ce{H2} obtained with the cc-pVQZ basis at the (SF-)CIS (top), (SF-)TD-BH\&HLYP (middle), and (SF-)BSE (bottom) levels of theory.
|
|
The reference EOM-CCSD excitation energies are represented as solid lines, while the results obtained with and without spin-flip are represented as dashed and dotted lines, respectively.
|
|
All the spin-conserved and spin-flip calculations have been performed with an unrestricted reference.
|
|
The raw data are reported in the {\SI}.
|
|
\label{fig:H2}}
|
|
\end{figure*}
|
|
%%% %%% %%%
|
|
|
|
%===============================
|
|
\subsection{Cyclobutadiene}
|
|
\label{sec:CBD}
|
|
%===============================
|
|
|
|
Cyclobutadiene (CBD) is an interesting example as the electronic character of its ground state can be tuned via geometrical deformation. \cite{Balkova_1994,Levchenko_2004,Manohar_2008,Karadakov_2008,Li_2009,Shen_2012,Lefrancois_2015,Casanova_2020,Vitale_2020}
|
|
%with potential large spin contamination.
|
|
In the $D_{2h}$ rectangular geometry of the $A_g$ singlet ground state, the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) are non-degenerate, and the singlet ground state can be safely labeled as single-reference with well-defined doubly-occupied orbitals.
|
|
However, in the $D_{4h}$ square-planar geometry of the $A_{2g}$ triplet state, the HOMO and LUMO are strictly degenerate, and the electronic ground state (which is still of singlet nature with $B_{1g}$ spatial symmetry, hence violating Hund's rule) is strongly multi-reference with singly occupied orbitals (\ie, singlet open-shell state).
|
|
In this case, single-reference methods notoriously fail.
|
|
Nonetheless, the lowest triplet state of symmetry $^3 A_{2g}$ remains of single-reference character and is then a perfect starting point for spin-flip calculations.
|
|
The $D_{2h}$ and $D_{4h}$ optimized geometries of the $^1 A_g$ and $^3 A_{2g}$ states of CBD have been extracted from Ref.~\onlinecite{Manohar_2008} and have been obtained at the CCSD(T)/cc-pVTZ level.
|
|
For comparison purposes, EOM-SF-CCSD and SF-ADC excitation energies have been extracted from Ref.~\onlinecite{Manohar_2008} and Ref.~\onlinecite{Lefrancois_2015}, respectively.
|
|
All of them have been obtained with a UHF reference like the SF-BSE calculations performed here.
|
|
|
|
Tables~\ref{tab:CBD_D2h} and \ref{tab:CBD_D4h} report excitation energies (with respect to the singlet ground state) obtained at the $D_{2h}$ and $D_{4h}$ geometries, respectively, for several methods using the spin-flip \textit{ansatz}.
|
|
All these results are represented in Fig.~\ref{fig:CBD}.
|
|
For each geometry, three excited states are under investigation:
|
|
i) the $1\,{}^3B_{1g}$, $1\,{}^1B_{1g}$, and $2\,{}^1A_{g}$ states of the $D_{2h}$ geometry;
|
|
ii) the $1\,{}^3 A_{2g}$, $2\,{}^1 A_{1g}$, and $1\,{}^1 B_{2g}$ states of the $D_{4h}$ geometry.
|
|
It is important to mention that the $2\,{}^1A_{1g}$ state of the rectangular geometry has a significant double excitation character, \cite{Loos_2019} and is then hardly described by second-order methods [such as CIS(D), \cite{Head-Gordon_1994,Head-Gordon_1995} ADC(2), \cite{Trofimov_1997,Dreuw_2015} CC2, \cite{Christiansen_1995a} or EOM-CCSD \cite{Koch_1990,Stanton_1993,Koch_1994}] and remains a real challenge for third-order methods [as, for example, ADC(3), \cite{Trofimov_2002,Harbach_2014,Dreuw_2015} CC3, \cite{Christiansen_1995b} or EOM-CCSDT \cite{Kucharski_1991,Kallay_2004,Hirata_2000,Hirata_2004}].
|
|
|
|
Comparing the present SF-BSE@{\GOWO} results for the rectangular geometry (see Table \ref{tab:CBD_D2h}) to the most accurate ADC level, \ie, SF-ADC(3), we have a difference in excitation energy of $0.017$ eV for the $1\,^3B_{1g}$ state.
|
|
This difference grows to $0.572$ eV for the $1\,^1B_{1g}$ state and then shrinks to $0.212$ eV for the $2\,^1A_{g}$ state.
|
|
Overall, adding dynamical corrections via the SF-dBSE@{\GOWO} scheme does not improve the accuracy of the excitation energies [as compared to SF-ADC(3)] with errors of $0.052$, $0.393$, and $0.293$ eV for the $1\,^3B_{1g}$, $1\,^1B_{1g}$, and $2\,^1 A_{g}$ states, respectively.
|
|
|
|
Now, looking at Table \ref{tab:CBD_D4h} which gathers the results for the square-planar geometry, we see that, at the SF-BSE@{\GOWO} level, the first two states are wrongly ordered with the triplet $1\,^3B_{1g}$ state lower than the singlet $1\,^1A_g$ state.
|
|
(The same observation can be made at the SF-TD-B3LYP level.)
|
|
This is certainly due to the poor Hartree-Fock reference which lacks opposite-spin correlation and this issue could be potentially alleviated by using a better starting point for the $GW$ calculation, as discussed in Sec.~\ref{sec:compdet}.
|
|
Nonetheless, it is pleasing to see that adding the dynamical correction in SF-dBSE@{\GOWO} not only improves the agreement with SF-ADC(3) but also retrieves the right state ordering.
|
|
Then, CBD stands as an excellent example for which dynamical corrections are necessary to get the right chemistry at the SF-BSE level.
|
|
|
|
%%% FIG 3 %%%
|
|
\begin{figure*}
|
|
\includegraphics[width=0.4\linewidth]{CBD_D2h}
|
|
\hspace{0.05\linewidth}
|
|
\includegraphics[width=0.4\linewidth]{CBD_D4h}
|
|
\caption{
|
|
Vertical excitation energies of CBD.
|
|
Left: $1\,{}^3B_{1g}$, $1\,{}^1B_{1g}$, and $2\,{}^1A_{1g}$ states at the $D_{2h}$ rectangular equilibrium geometry of the $\text{X}\,{}^1 A_{g}$ ground state (see Table \ref{tab:CBD_D2h} for the raw data).
|
|
Right: $1\,{}^3A_{2g}$, $2\,{}^1A_{1g}$, and $1\,{}^1B_{2g}$ states at the $D_{4h}$ square-planar equilibrium geometry of the $1\,{}^3 A_{2g}$ state (see Table \ref{tab:CBD_D4h} for the raw data).
|
|
All the spin-flip calculations have been performed with a UHF reference and the cc-pVTZ basis set.
|
|
\label{fig:CBD}}
|
|
\end{figure*}
|
|
%%% %%% %%%
|
|
|
|
%%% TABLE II %%%
|
|
\begin{table}
|
|
\caption{
|
|
Vertical excitation energies (with respect to the singlet $\text{X}\,{}^1A_{g}$ ground state) of the $1\,{}^3B_{1g}$, $1\,{}^1B_{1g}$, and $2\,{}^1A_{g}$ states of CBD at the $D_{2h}$ rectangular equilibrium geometry of the $\text{X}\,{}^1 A_{g}$ ground state.
|
|
All the spin-flip calculations have been performed with a UHF reference and the cc-pVTZ basis set.
|
|
\label{tab:CBD_D2h}}
|
|
\begin{ruledtabular}
|
|
\begin{tabular}{lrrr}
|
|
& \mc{3}{c}{Excitation energies (eV)} \\
|
|
\cline{2-4}
|
|
Method & $1\,{}^3B_{1g}$ & $1\,{}^1B_{1g}$ & $2\,{}^1A_{g}$ \\
|
|
\hline
|
|
SF-TD-B3LYP\fnm[3] & $1.750$ & $2.260$ & $4.094$ \\
|
|
SF-TD-BH\&HLYP\fnm[3] & $1.583$ & $2.813$ & $4.528$ \\
|
|
SF-CIS\fnm[1] & $1.521$ & $3.836$ & $5.499$ \\
|
|
EOM-SF-CCSD\fnm[1] & $1.654$ & $3.416$ & $4.360$ \\
|
|
EOM-SF-CCSD(fT)\fnm[1] & $1.516$ & $3.260$ & $4.205$ \\
|
|
EOM-SF-CCSD(dT)\fnm[1] & $1.475$ & $3.215$ & $4.176$ \\
|
|
SF-ADC(2)-s\fnm[2] & $1.573$ & $3.208$ & $4.247$ \\
|
|
SF-ADC(2)-x\fnm[2] & $1.576$ & $3.141$ & $3.796$ \\
|
|
SF-ADC(3)\fnm[2] & $1.456$ & $3.285$ & $4.334$ \\
|
|
SF-BSE@{\GOWO}\fnm[3] & $1.438$ & $2.704$ & $4.540$ \\
|
|
SF-dBSE@{\GOWO}\fnm[3] & $1.403$ & $2.883$ & $4.621$ \\
|
|
\end{tabular}
|
|
\end{ruledtabular}
|
|
\fnt[1]{Spin-flip EOM-CC values from Ref.~\onlinecite{Manohar_2008}.}
|
|
\fnt[2]{Values from Ref.~\onlinecite{Lefrancois_2015}.}
|
|
\fnt[3]{This work.}
|
|
\end{table}
|
|
%%% %%% %%% %%%
|
|
|
|
%%% TABLE III %%%
|
|
\begin{table}
|
|
\caption{
|
|
Vertical excitation energies (with respect to the singlet $\text{X}\,{}^1B_{1g}$ ground state) of the $1\,{}^3A_{2g}$, $2\,{}^1A_{1g}$, and $1\,{}^1B_{2g}$ states of CBD at the $D_{4h}$ square-planar equilibrium geometry of the $1\,{}^3A_{2g}$ state.
|
|
All the spin-flip calculations have been performed with a UHF reference and the cc-pVTZ basis set.
|
|
\label{tab:CBD_D4h}}
|
|
\begin{ruledtabular}
|
|
\begin{tabular}{lrrr}
|
|
& \mc{3}{c}{Excitation energies (eV)} \\
|
|
\cline{2-4}
|
|
Method & $1\,{}^3A_{2g}$ & $2\,{}^1A_{1g}$ & $1\,{}^1B_{2g}$ \\
|
|
\hline
|
|
SF-TD-B3LYP\fnm[3] & $-0.020$ & $0.486$ & $0.547$ \\
|
|
SF-TD-BH\&HLYP\fnm[3] & $0.048$ & $1.282$ & $1.465$ \\
|
|
SF-CIS\fnm[1] & $0.317$ & $3.125$ & $2.650$ \\
|
|
EOM-SF-CCSD\fnm[1] & $0.369$ & $1.824$ & $2.143$ \\
|
|
EOM-SF-CCSD(fT)\fnm[1] & $0.163$ & $1.530$ & $1.921$ \\
|
|
EOM-SF-CCSD(dT)\fnm[1] & $0.098$ & $1.456$ & $1.853$ \\
|
|
SF-ADC(2)-s\fnm[2] & $0.266$ & $1.664$ & $1.910$ \\
|
|
SF-ADC(2)-x\fnm[2] & $0.217$ & $1.123$ & $1.799$ \\
|
|
SF-ADC(3)\fnm[2] & $0.083$ & $1.621$ & $1.930$ \\
|
|
SF-BSE@{\GOWO}\fnm[3] & $-0.049$ & $1.189$ & $1.480$ \\
|
|
SF-dBSE@{\GOWO}\fnm[3] & $0.012$ & $1.507$ & $1.841$ \\
|
|
\end{tabular}
|
|
\end{ruledtabular}
|
|
\fnt[1]{Spin-flip EOM-CC values from Ref.~\onlinecite{Manohar_2008}.}
|
|
\fnt[2]{Values from Ref.~\onlinecite{Lefrancois_2015}.}
|
|
\fnt[3]{This work.}
|
|
\end{table}
|
|
%%% %%% %%% %%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Conclusion}
|
|
\label{sec:ccl}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
In this article, we have presented the extension of the BSE approach of many-body perturbation theory to the spin-flip formalism in order to access double excitations in realistic molecular systems.
|
|
The present spin-flip calculations rely on a spin-unrestricted version of the $GW$ approximation and the BSE formalism with, on top of this, a dynamical correction to the static BSE optical excitations via an unrestricted generalization of our recently developed renormalized perturbative treatment.
|
|
Taking the beryllium atom, the dissociation of the hydrogen molecule, and cyclobutadiene in two different geometries as examples, we have shown that the spin-flip BSE formalism can accurately model double excitations and seems to surpass systematically its spin-flip TD-DFT parent.
|
|
Further improvements could be obtained thanks to a better choice of the starting orbitals and their energies and we hope to investigate this in a forthcoming paper.
|
|
Techniques to alleviate the spin contamination in spin-flip calculations will also be explored in the near future.
|
|
We hope to these new encouraging results will stimulate new developments around the BSE formalism to further establish it as a valuable \text{ab inito} alternative to TD-DFT for the study of molecular excited states.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\acknowledgements{
|
|
We would like to thank Pina Romaniello, Xavier Blase, and Denis Jacquemin for insightful discussions.
|
|
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section*{Supporting information available}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
Additional graphs comparing (SF-)TD-BLYP, (SF-)TD-B3LYP, and (SF-)dBSE with EOM-CCSD for the \ce{H2} molecule, raw data associated with Fig.~\ref{fig:H2}, output files associated with all the calculations performed in the present article.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%\section*{Data availability statement}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%The data that supports the findings of this study are available within the article and its supplementary material.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\bibliography{sfBSE,sfBSE-control}
|
|
%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
\end{document}
|