PTEROSOR/BSEdyn
10
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

ready for submission

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-07-27 14:19:48 +02:00
parent 59ee48009b
commit d994ade3b4
3 changed files with 41 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
BSEdyn.tex
View File

@ -205,10 +205,11 @@
\affiliation{\NEEL}
\begin{abstract}
Similar to the ubiquitous adiabatic approximation in time-dependent density-functional theory, the static approximation, which substitutes a dynamical (\ie, frequency-dependent) kernel by its static limit, is usually enforced in most implementations of the Bethe-Salpeter equation (BSE) formalism.
The Bethe-Salpeter equation (BSE) formalism is a computationally affordable method for the calculation of accurate optical excitation energies in molecular systems.
Similar to the ubiquitous adiabatic approximation of time-dependent density-functional theory, the static approximation, which substitutes a dynamical (\ie, frequency-dependent) kernel by its static limit, is usually enforced in most implementations of the BSE formalism.
Here, going beyond the static approximation, we compute the dynamical correction of the electron-hole screening for molecular excitation energies thanks to a renormalized first-order perturbative correction to the static BSE excitation energies.
The present dynamical correction goes beyond the plasmon-pole approximation as the dynamical screening of the Coulomb interaction is computed exactly within the random-phase approximation.
Our calculations are benchmarked against high-level (coupled-cluster) calculations, allowing to assess the clear improvements brought by the dynamical correction.
Our calculations are benchmarked against high-level (coupled-cluster) calculations, allowing to assess the clear improvement brought by the dynamical correction for both singlet and triplet optical transitions.
%\\
%\bigskip
%\begin{center}
@ -223,7 +224,7 @@ Our calculations are benchmarked against high-level (coupled-cluster) calculatio
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%
The Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} is to the $GW$ approximation \cite{Hedin_1965,Golze_2019} of many-body perturbation theory (MBPT) \cite{Martin_2016} what time-dependent density-functional theory (TD-DFT) \cite{Runge_1984,Casida_1995} is to Kohn-Sham density-functional theory (KS-DFT), \cite{Hohenberg_1964,Kohn_1965} an affordable way of computing the neutral (or optical) excitations of a given electronic system.
The Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} is to the $GW$ approximation \cite{Hedin_1965,Golze_2019} of many-body perturbation theory (MBPT) \cite{Onida_2002,Martin_2016} what time-dependent density-functional theory (TD-DFT) \cite{Runge_1984,Casida_1995} is to Kohn-Sham density-functional theory (KS-DFT), \cite{Hohenberg_1964,Kohn_1965} an affordable way of computing the neutral (or optical) excitations of a given electronic system.
In recent years, it has been shown to be a valuable tool for computational chemists with a large number of systematic benchmark studies on large families of molecular systems appearing in the literature \cite{Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018} (see Ref.~\onlinecite{Blase_2018} for a recent review).
Qualitatively, taking the optical gap (\ie, the lowest optical excitation energy) as an example, BSE builds on top of a $GW$ calculation by adding up excitonic effects (\ie, the electron-hole binding energy $\EB$) to the $GW$ HOMO-LUMO gap
@ -283,7 +284,7 @@ In this regard, MBPT provides key insights about what is missing in adiabatic TD
In the present study, we extend the work of Rohlfing and coworkers \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b} by proposing a renormalized first-order perturbative correction to the static BSE excitation energies.
Importantly, our correction goes beyond the plasmon-pole approximation as the dynamical screening of the Coulomb interaction is computed exactly.
In order to assess the accuracy of the present scheme, we report singlet and triplet excitation energies of various natures for small- and medium-size molecules.
Our calculations are benchmarked against high-level (coupled-cluster) calculations, allowing to clearly evidence the systematic improvements brought by the dynamical correction.
Our calculations are benchmarked against high-level coupled-cluster (CC) calculations, allowing to clearly evidence the systematic improvement brought by the dynamical correction.
In particular, we found that, although $n \ra \pis$ and $\pi \ra \pis$ transitions are systematically red-shifted by $0.3$--$0.6$ eV, dynamical effects have a much smaller magnitude for charge transfer (CT) and Rydberg states.
Unless otherwise stated, atomic units are used.
@ -292,7 +293,7 @@ Unless otherwise stated, atomic units are used.
\label{sec:theory}
%%%%%%%%%%%%%%%%%%%%%%%%
In this Section, following Strinati's seminal work, \cite{Strinati_1988} we first derive in some details the theoretical foundations leading to the dynamical BSE.
In this Section, following Strinati's seminal work, \cite{Strinati_1988} we first discuss in some details the theoretical foundations leading to the dynamical BSE.
%Additional details about this derivation are provided as {\SI}.
We present, in a second step, the perturbative implementation of the dynamical correction as compared to the standard static approximation.
@ -514,7 +515,7 @@ with
\B{ia,jb}{\RPA} & = 2 \ERI{ia}{bj},
\end{align}
\end{subequations}
where the $\e{p}$'s are taken as the HF orbital energies in the case of $G_0W_0$ \cite{Hybertsen_1985a, Hybertsen_1986} or as the $GW$ quasiparticle energies in the case of self-consistent schemes such as ev$GW$. \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011}
where the $\e{p}$'s are taken as the HF orbital energies in the case of $G_0W_0$ \cite{Hybertsen_1985a, Hybertsen_1986} or as the $GW$ quasiparticle energies in the case of self-consistent schemes such as ev$GW$. \cite{Hybertsen_1986,Shishkin_2007,Blase_2011,Faber_2011,Rangel_2016,Kaplan_2016,Gui_2018}
The RPA matrices $\bA{\RPA}$ and $\bB{\RPA}$ in Eq.~\eqref{eq:LR-RPA} are of size $\Nocc \Nvir \times \Nocc \Nvir$, where $\Nocc$ and $\Nvir$ are the number of occupied and virtual orbitals (\ie, $\Norb = \Nocc + \Nvir$ is the total number of spatial orbitals), respectively, and $\bX{m}{\RPA}$, and $\bY{m}{\RPA}$ are (eigen)vectors of length $\Nocc \Nvir$.
The analysis of the poles of the integrand in Eq.~\eqref{eq:wtilde} yields
@ -551,7 +552,7 @@ The analysis of the (off-diagonal) screened Coulomb potential matrix elements mu
\times \qty[ \frac{1}{\Om{ij}{S} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ba}{S} - \Om{m}{\RPA} + i\eta} ],
\end{multline}
reveals strong divergences even for low-lying excitations when, for example, $\Om{ba}{S} - \Om{m}{\RPA} = \Om{S}{} - \Om{m}{\RPA} - ( \eGW{a} - \eGW{b} ) \approx 0$.
Such divergences may explain that in previous studies dynamical effects were only accounted for at the TDA level. \cite{Strinati_1988,Rohlfing_2000,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Zhang_2013,Rebolini_2016}
Such divergences may explain that, in previous studies, dynamical effects were only accounted for at the TDA level. \cite{Strinati_1988,Rohlfing_2000,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Zhang_2013,Rebolini_2016}
To avoid confusions here, enforcing the TDA for the dynamical correction (which corresponds to neglecting the dynamical correction originating from the anti-resonant part of the BSE Hamiltonian) will be labeled as dTDA in the following.
Going beyond the dTDA is outside the scope of the present study but shall be addressed eventually.
@ -739,7 +740,7 @@ Although it might be reduced to $\order*{\Norb^4}$ operations with standard reso
All systems under investigation have a closed-shell singlet ground state.
We then adopt a restricted formalism throughout this work.
The $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies are done using a (restricted) HF starting point.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986} quasiparticle energies are employed as starting points to compute the BSE neutral excitations.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986,vanSetten_2013} quasiparticle energies are employed as starting points to compute the BSE neutral excitations.
These quasiparticle energies are obtained by linearizing the frequency-dependent quasiparticle equation, and the entire set of orbitals is corrected.
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018b,Veril_2018}.
Note that, for the present (small) molecular systems, {\GOWO}@HF and ev$GW$@HF yield similar quasiparticle energies and fundamental gap.
@ -748,7 +749,7 @@ In the present study, the zeroth-order Hamiltonian [see Eq.~\eqref{eq:LR-PT}] is
The dynamical correction, however, is computed in the dTDA throughout.
As one-electron basis sets, we employ the Dunning families cc-pVXZ and aug-cc-pVXZ (X = D, T, and Q) defined with cartesian Gaussian functions.
Finally, the infinitesimal $\eta$ is set to $100$ meV for all calculations.
It is important to mention that the small molecular systems considered here are particularly challenging for the BSE formalism which is known to work best for larger systems where the amount of screening is more important. \cite{Jacquemin_2017b,Rangel_2017}
It is important to mention that the small molecular systems considered here are particularly challenging for the BSE formalism, \cite{Hirose_2015,Loos_2018b} which is known to work best for larger systems where the amount of screening is more important. \cite{Jacquemin_2017b,Rangel_2017}
For comparison purposes, we employ the theoretical best estimates (TBEs) and geometries of Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} from which CIS(D), \cite{Head-Gordon_1994,Head-Gordon_1995} ADC(2), \cite{Trofimov_1997,Dreuw_2015} CC2, \cite{Christiansen_1995a} CCSD, \cite{Purvis_1982} and CC3 \cite{Christiansen_1995b} excitation energies are also extracted.
Various statistical quantities are reported in the following: the mean signed error (MSE), mean absolute error (MAE), root-mean-square error (RMSE), and the maximum positive [Max($+$)] and maximum negative [Max($-$)] errors.
@ -822,7 +823,7 @@ The \ce{N2} molecule is a very convenient example for this kind of study as it c
As we shall further illustrate below, the magnitude of the dynamical correction is characteristic of the type of transitions.
One key result of the present investigation is that the dynamical correction is quite basis set insensitive with a maximum variation of $0.03$ eV between aug-cc-pVDZ and aug-cc-pVQZ.
It is only for the smallest basis set (cc-pVDZ) that one can observe significant differences.
We can then safely conclude that the dynamical correction converges rapidly with respect to the size of the one-electron basis set, a triple-$\zeta$ or an augmented double-$\zeta$ basis set being enough to obtain near complete basis set limit values.
We can then safely conclude that the dynamical correction converges rapidly with respect to the size of the one-electron basis set, a triple-$\zeta$ or an augmented double-$\zeta$ basis being enough to obtain near complete basis set limit values.
This is quite a nice feature as it means that one does not need to compute the dynamical correction in a very large basis to get a meaningful estimate of its magnitude.
%The second important observation extracted from the results gathered in Table \ref{tab:N2} is that the dTDA is a rather satisfactory approximation, especially for the singlet states where one observes a maximum discrepancy of $0.03$ eV between the ``full'' and dTDA excitation energies.
%The story is different for the triplet states for which deviations of the order of $0.3$ eV is observed between the two sets, the dTDA of excitation energies being, as we shall see later on, more accurate.
@ -882,11 +883,11 @@ This is quite a nice feature as it means that one does not need to compute the d
& $^1\Pi$ & Ryd. & & 12.37 & 12.32 & -0.05 & 1.004
& 12.06 & 12.03 & 11.83 & 11.96 & 11.72 \\
\\
\ce{HNO} & $^1A''(n \ra \pis)$ & Val. & 11.71 & 2.46 & 1.98 & -0.48 & 1.035
& 1.80 & 1.68 & 1.74 & 1.76 & 1.74 \\
& $^1A'$ & Ryd. & & 7.05 & 7.01 & -0.04 & 1.003
& 5.81 & 5.73 & 5.72 & 6.30 & 6.27 \\
\\
% \ce{HNO} & $^1A''(n \ra \pis)$ & Val. & 11.71 & 2.46 & 1.98 & -0.48 & 1.035
% & 1.80 & 1.68 & 1.74 & 1.76 & 1.74 \\
% & $^1A'$ & Ryd. & & 7.05 & 7.01 & -0.04 & 1.003
% & 5.81 & 5.73 & 5.72 & 6.30 & 6.27 \\
% \\
\ce{C2H2} & $^1\Sigma_{u}^-(\pi \ra \pis)$ & Val. & 12.28 & 7.37 & 7.05 & -0.32 & 1.026
& 7.28 & 7.24 & 7.26 & 7.15 & 7.10 \\
& $^1\Delta_{u}(\pi \ra \pis)$ & Val. & & 7.74 & 7.46 & -0.29 & 1.025
@ -914,12 +915,12 @@ This is quite a nice feature as it means that one does not need to compute the d
& $^1A_1(\pi \ra \pis)$ & Val. & & 10.05 & 9.81 & -0.24 & 1.026
& 9.08 & 9.46 & 9.54 & 9.67 & 9.43 \\
\hline
MAE & & & & 0.65 & 0.50 & &
& 0.41 & 0.24 & 0.25 & 0.14 & 0.00 \\
MSE & & & & 0.65 & 0.48 & &
& 0.12 & 0.00 & 0.00 & 0.13 & 0.00 \\
RMSE & & & & 0.71 & 0.58 & &
& 0.54 & 0.34 & 0.33 & 0.19 & 0.00 \\
MAE & & & & 0.64 & 0.50 & &
& 0.43 & 0.24 & 0.25 & 0.15 & 0.00 \\
MSE & & & & 0.64 & 0.48 & &
& 0.14 & 0.02 & 0.03 & 0.14 & 0.00 \\
RMSE & & & & 0.70 & 0.58 & &
& 0.55 & 0.33 & 0.33 & 0.20 & 0.00 \\
Max($+$) & & & & 1.08 & 0.91 & &
& 1.06 & 0.54 & 0.57 & 0.44 & 0.00 \\
Max($-$) & & & & 0.20 & -0.22 & &
@ -970,11 +971,11 @@ This is quite a nice feature as it means that one does not need to compute the d
& $^3\Sigma_u^+$ & Ryd. & & 11.17 & 11.07 & -0.10 & 1.008
& 10.98 & 10.83 & 10.60 & 10.71 & 10.47 \\
\\
\ce{HNO} & $^3A''(n \ra \pis)$ & Val. & 11.71 & 1.27 & 0.67 & -0.60 & 1.036
& 0.91 & 0.78 & 0.84 & 0.85 & 0.88 \\
& $^3A'(\pi \ra \pis)$ & Val. & & 5.55 & 4.87 & -0.69 & 1.037
& 5.72 & 5.46 & 5.44 & 5.49 & 5.61 \\
\\
% \ce{HNO} & $^3A''(n \ra \pis)$ & Val. & 11.71 & 1.27 & 0.67 & -0.60 & 1.036
% & 0.91 & 0.78 & 0.84 & 0.85 & 0.88 \\
% & $^3A'(\pi \ra \pis)$ & Val. & & 5.55 & 4.87 & -0.69 & 1.037
% & 5.72 & 5.46 & 5.44 & 5.49 & 5.61 \\
% \\
\ce{C2H2} & $^3\Sigma_{u}^+(\pi \ra \pis)$ & Val. & 12.28 & 5.83 & 5.32 & -0.51 & 1.031
& 5.79 & 5.75 & 5.76 & 5.45 & 5.53 \\
& $^3\Delta_{u}(\pi \ra \pis)$ & Val. & & 6.64 & 6.23 & -0.41 & 1.028
@ -996,16 +997,16 @@ This is quite a nice feature as it means that one does not need to compute the d
& $^3B_2(n \ra 3s)$ & Ryd. & & 7.60 & 7.56 & -0.05 & 1.002
& 6.66 & 6.39 & 6.44 & 7.08 & 7.06 \\
\hline
MAE & & & & 0.39 & 0.29 & &
& 0.25 & 0.21 & 0.22 & 0.09 & 0.00 \\
MSE & & & & 0.39 & 0.01 & &
& 0.21 & 0.08 & 0.12 & 0.04 & 0.00 \\
RMSE & & & & 0.44 & 0.35 & &
& 0.30 & 0.27 & 0.29 & 0.13 & 0.00 \\
MAE & & & & 0.41 & 0.27 & &
& 0.27 & 0.21 & 0.24 & 0.10 & 0.00 \\
MSE & & & & 0.41 & 0.06 & &
& 0.23 & 0.10 & 0.14 & 0.05 & 0.00 \\
RMSE & & & & 0.45 & 0.33 & &
& 0.31 & 0.27 & 0.30 & 0.13 & 0.00 \\
Max($+$) & & & & 0.70 & 0.60 & &
& 0.63 & 0.57 & 0.63 & 0.29 & 0.00 \\
Max($-$) & & & & -0.06 & -0.74 & &
& -0.40 & -0.67 & -0.62 & -0.12 & 0.00 \\
Max($-$) & & & & 0.11 & -0.39 & &
& -0.40 & -0.67 & -0.62 & -0.11 & 0.00 \\
\end{tabular}
\end{ruledtabular}
\end{table*}
@ -1017,6 +1018,7 @@ The highly-accurate TBEs of Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} (
For each excitation, we report the static and dynamic excitation energies, $\Om{S}{\stat}$ and $\Om{S}{\dyn}$, as well as the value of the renormalization factor $Z_S$ defined in Eq.~\eqref{eq:Z}.
As one can see in Tables \ref{tab:BigTabSi} and \ref{tab:BigTabTr}, the value of $Z_S$ is always quite close to unity which shows that the perturbative expansion behaves nicely, and that a first-order correction is probably quite a good estimate of the non-perturbative result.
Moreover, we have observed that an iterative, self-consistent resolution [where the dynamically-corrected excitation energies are re-injected in Eq.~\eqref{eq:Om1}] yields basically the same results as its (cheaper) renormalized version.
Note that, unlike in $GW$ where the renormalization factor lies in between $0$ and $1$, the dynamical BSE renormalization factor $Z_S$ defined in Eq.~\eqref{eq:Z} can be smaller or greater than unity.
A clear general trend is the consistent red shift of the static BSE excitation energies induced by the dynamical correction, as anticipated in Sec.~\ref{sec:dynW}.
%%% FIG I %%%
@ -1032,7 +1034,7 @@ A clear general trend is the consistent red shift of the static BSE excitation e
The results gathered in Tables \ref{tab:BigTabSi} and \ref{tab:BigTabTr} are depicted in Fig.~\ref{fig:SiTr-SmallMol}, where we report the error (with respect to the TBEs) for the singlet and triplet excitation energies computed within the static and dynamic BSE formalism.
From this figure, it is quite clear that the dynamically-corrected excitation energies are systematically improved upon their static analogs, especially for singlet states.
(In the case of triplets, one would notice a few cases where the excitation energies is underestimated.)
In particular, the MAE is reduced from $0.65$ to $0.50$ eV for singlets, and from $0.39$ to $0.29$ eV for triplets.
In particular, the MAE is reduced from $0.64$ to $0.50$ eV for singlets, and from $0.41$ to $0.27$ eV for triplets.
The MSE and RMSE are also systematically improved when one takes into account dynamical effects.
The second important observation extracted from Fig.~\ref{fig:SiTr-SmallMol} is that the (singlet and triplet) Rydberg states are rather unaltered by the dynamical effects with a correction of few hundredths of eV in most cases.
The same comment applies to the CT excited state of \ce{HCl}.
@ -1105,7 +1107,7 @@ However, it is definitely an improvement in terms of performances as compared to
Table \ref{tab:BigMol} reports singlet and triplet excitation energies for larger molecules (acrolein \ce{H2C=CH-CH=O}, butadiene \ce{H2C=CH-CH=CH2}, diacetylene \ce{HC#C-C#CH}, glyoxal \ce{O=CH-CH=O}, and streptocyanine-C1 \ce{H2N-CH=NH2+}) at the static and dynamic BSE levels with the aug-cc-pVDZ basis set.
We also report the CC3 excitation energies computed in Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} with the same basis set.
These will be our reference as they are known to be extremely accurate ($0.03$--$0.04$ eV from the TBEs). \cite{Loos_2020g}
These will be our reference as they are known to be extremely accurate ($0.03$--$0.04$ eV from the TBEs). \cite{Loos_2018a,Loos_2019,Loos_2020b,Loos_2020g}
Errors associated with these excitation energies (with respect to CC3) are represented in Fig.~\ref{fig:SiTr-BigMol}.
As expected the static BSE excitation energies are much more accurate for these larger molecules with a MAE of $0.32$ eV, a MSE of $0.30$ eV, and a RMSE of $0.38$ eV.
Here again, the dynamical correction improves the accuracy of BSE by lowering the MAE, MSE, and RMSE to $0.23$, $0.00$, and $0.29$ eV, respectively.
@ -1149,7 +1151,7 @@ Following Strinati's footsteps, \cite{Strinati_1982,Strinati_1984,Strinati_1988}
In the present study, we have computed exactly the dynamical screening of the Coulomb interaction within the random-phase approximation, going effectively beyond both the usual static approximation and the plasmon-pole approximation.
In order to assess the accuracy of the present scheme, we have reported a significant number of calculations for various molecular systems.
Our calculations have been benchmarked against high-level (coupled-cluster) calculations, allowing to clearly evidence the systematic improvements brought by the dynamical correction for both singlet and triplet excited states.
Our calculations have been benchmarked against high-level CC calculations, allowing to clearly evidence the systematic improvements brought by the dynamical correction for both singlet and triplet excited states.
We have found that, although $n \ra \pis$ and $\pi \ra \pis$ transitions are systematically red-shifted by $0.3$--$0.6$ eV thanks to dynamical effects, their magnitude is much smaller for CT and Rydberg states.
%%%%%%%%%%%%%%%%%%%%%%%%
@ -1183,7 +1185,7 @@ Combining the Fourier transform (with respect to $t_1$) of $L_0(1,4;1',3)$
[L_0](\bx_1,4;\bx_{1'},3 \; | \; \omega_1 )
= -i \int dt_1 e^{i \omega_1 t_1 } G(1,3)G(4,1'),
\end{align}
(where $t_{1'} = t_1^{+}$) with the \trashPFL{Fourier expansion} \titou{inverse Fourier transform} of the Green's function, \eg,
(where $t_{1'} = t_1^{+}$) with the inverse Fourier transform of the Green's function, \eg,
\begin{align}
G(1,3) = \int \frac{ d\omega }{ 2\pi } G(\bx_1,\bx_3;\omega) e^{-i \omega \tau_{13} },
\end{align}
@ -1263,7 +1265,7 @@ one gets
\end{split}
\end{equation}
%We now act on the $N$-electron ground-state wave function with
\titou{Assuming now that the $\e{a}$'s and $\e{i}$'s are proper addition and removal energies (respectively)}, such as the $GW$ quasiparticle energies, one can use the following relationships
Assuming now that the $\lbrace \e{p} \rbrace$'s are proper addition/removal energies, such as the $GW$ quasiparticle energies, one can use the following relationships
\begin{subequations}
\begin{align}
e^{+i\hH \tau_{65} } \ha^{\dagger}_p \ket{N} &=

BIN
fig1a.pdf
View File

Binary file not shown.

BIN
fig1b.pdf
View File

Binary file not shown.