ccl and discussion

BSEdyn.tex

@@ -226,7 +226,7 @@ Our calculations are benchmarked against high-level (coupled-cluster) calculatio
 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. 
 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 \titou{(\ie, the electron-hole binding energy) $\EB$} to the $GW$ HOMO-LUMO gap
+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
 \begin{equation}
 	\Eg^{\GW} = \eps_{\LUMO}^{\GW} - \eps_{\HOMO}^{\GW},
 \end{equation}
@@ -238,7 +238,7 @@ in order to approximate the optical gap
 \begin{equation}
 	\EgOpt = E_1^{N} - E_0^{N} = \EgFun + \EB,
 \end{equation}
-where \trashPFL{$\EB$ is the electron-hole binding energy and}
+where 
 \begin{equation} \label{eq:Egfun}
 	\EgFun = I^N - A^N
 \end{equation}
@@ -266,7 +266,7 @@ Zhang \textit{et al.}~have studied the frequency-dependent second-order Bethe-Sa
 Rebolini and Toulouse have performed a similar investigation in a range-separated context, and they have reported a modest improvement over its static counterpart. \cite{Rebolini_2016,Rebolini_PhD} 
 In these two latter studies, they also followed a (non-self-consistent) perturbative approach within the TDA with a renormalization of the first-order perturbative correction.
 
-It is important to note that, although all the studies mentioned above are clearly going beyond the static approximation of BSE, they are not able to recover additional excitations as the perturbative treatment \titou{accounts for dynamical effects only on excitations already present in the static limit.} However, it does permit to recover, for transitions with a dominant single-excitation character, additional relaxation effects coming from higher excitations.
+It is important to note that, although all the studies mentioned above are clearly going beyond the static approximation of BSE, they are not able to recover additional excitations as the perturbative treatment accounts for dynamical effects only on excitations already present in the static limit. However, it does permit to recover, for transitions with a dominant single-excitation character, additional relaxation effects coming from higher excitations.
 These higher excitations would be explicitly present in the BSE Hamiltonian by ``unfolding'' the dynamical BSE kernel, and one would recover a linear eigenvalue problem with, nonetheless, a much larger dimension. \cite{Loos_2020f}
 
 Based on a simple two-level model which permits to analytically solve the dynamical equations, Romaniello and coworkers \cite{Romaniello_2009b,Sangalli_2011} evidenced that one can genuinely access additional excitations by solving the non-linear, frequency-dependent eigenvalue problem. 
@@ -281,8 +281,10 @@ Following a similar strategy, Romaniello \textit{et al.} \cite{Romaniello_2009b}
 In this regard, MBPT provides key insights about what is missing in adiabatic TD-DFT, as discussed in details by Casida and Huix-Rotllant in Ref.~\onlinecite{Casida_2016}.
 
 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. \xavier{By comparison with higher-level calculations coupled-cluster (CC) we show that ... RESULTS AS IN CONCLUSION ... }
-%Moreover, we investigate quantitatively the effect of the TDA by computing both the resonant and anti-resonant dynamical corrections to the 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.
+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.
 
 %%%%%%%%%%%%%%%%%%%%%%%%
@@ -1033,6 +1035,7 @@ From this figure, it is quite clear that the dynamically-corrected excitation en
 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.
 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}.
 The magnitude of the dynamical correction for $n \ra \pis$ and $\pi \ra \pis$ transitions is much more important: $0.3$--$0.5$ eV for singlets and $0.3$--$0.7$ eV for triplets.
 
 Dynamical BSE does not quite reach the accuracy of second-order methods [CIS(D), ADC(2), CC2, and CCSD] for the singlet and triplet optical excitations of these small molecules.
@@ -1108,29 +1111,31 @@ As expected the static BSE excitation energies are much more accurate for these
 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.
 Rydberg states are again very slightly affected by dynamical effects, while the dynamical corrections associated with the $n \ra \pis$ and $\pi \ra \pis$ transitions are much larger and of the same magnitude ($0.3$--$0.6$ eV) for both types of transitions. 
 This latter observation is somehow different from the outcomes reached by Rohlfing and coworkers in previous works \cite{Ma_2009a,Ma_2009b} (see Sec.~\ref{sec:intro}) where they observed i) smaller corrections, and ii) that $n \ra \pis$ transitions are more affected by the dynamical screening than $\pi \ra \pis$ transitions. 
-\titou{The larger size of the molecules considered in Refs.~\onlinecite{Ma_2009a,Ma_2009b} may play a role on the magnitude of the corrections, even though we do not observe here a significant reduction going from small systems (\ce{N2}, \ce{CO}, \ldots) to larger ones (acrolein, butadiene, \ldots). 
+The larger size of the molecules considered in Refs.~\onlinecite{Ma_2009a,Ma_2009b} may play a role on the magnitude of the corrections, even though we do not observe here a significant reduction going from small systems (\ce{N2}, \ce{CO}, \ldots) to larger ones (acrolein, butadiene, \ldots). 
 We emphasize further that previous calculations \cite{Ma_2009a,Ma_2009b} were performed within the plasmon-pole approximation for modeling the dynamical behaviour of the screened Coulomb potential, while we go beyond this approximation in the present study [see Eq.~\eqref{eq:wtilde}]. 
-Finally, while errors were defined with respect to experimental data in Refs.~\onlinecite{Ma_2009a,Ma_2009b}, we consider here as reference high-level CC calculations performed with the very same geometries and basis sets. 
+Finally, while errors were defined with respect to experimental data in Refs.~\onlinecite{Ma_2009a,Ma_2009b}, we consider here as reference high-level CC calculations performed with the very same geometries and basis sets than our BSE calculations. 
 As pointed out in previous works, \cite{Loos_2018,Loos_2019b,Loos_2020g} a direct comparison between theoretical transition energies and experimental data is a delicate task, as many factors (such as zero-point vibrational energies and geometrical relaxation) must be taken into account for fair comparisons.
-Further investigations are required to better evaluate the impact of these considerations on the influence of dynamical screening.}
+Further investigations are required to better evaluate the impact of these considerations on the influence of dynamical screening.
+
+To provide further insight into the magnitude of the dynamical correction to valence, Rydberg, and CT excitations, let us consider a simple two-level systems with $i = j = h$ and $a = b = l$, where $(h,l)$ stand for HOMO and LUMO. 
+The dynamical correction associated with the HOMO-LUMO transition reads
+\begin{equation*}
+	\W{hh,ll}{\text{stat}} - \widetilde{W}_{hh,ll}( \Om{1}{} )
+	=     4 \sERI{hh}{hl} \sERI{ll}{hl} \frac{ \Om{hl}{1} - 2 \OmRPA{hl}{}}{\OmRPA{hl}{} ( \OmRPA{hl}{} - \Om{hl}{1} ) },
+\end{equation*}
+where the only RPA excitation energy, $\Om{1}{}$, is again the HOMO-LUMO transition, \ie, $m=hl$ [see Eq.~\eqref{eq:sERI}]. 
+For CT excitations with vanishing HOMO-LUMO overlap [\ie, $\ERI{h}{l} \approx 0$], $\sERI{hh}{hl} \approx 0$ and $\sERI{ll}{hl} \approx 0$, so that one can expect the dynamical correction to be weak. 
+Likewise, Rydberg transitions which are characterized by a delocalized LUMO state, that is, a small HOMO-LUMO overlap, are expected to undergo weak dynamical corrections.  
+The discussion for $\pi \ra \pis$ and $n \ra \pis$ transitions is certainly more complex and molecule-specific symmetry arguments must be invoked to understand the magnitude of the $\sERI{hh}{hl}$ and $\sERI{ll}{hl}$ terms. 
+
+%[ Taking the simple case of a $\pi$-conjugated planar  molecule located in the (xy) plane, then $\pi-\pi^*$ transition will yield $[hh|hl]$ and $[ll|hl]$ terms built from $(p_z^2 | p_z^2)$ contributions on each atom. On the contrary, $n-\pi^*$ transitions will involve $(sp_z|ss)$ or $(sp_z|p_z
+%^2)$ local contribution that are small by planar symmetry. Bah bon ben ... should lead to small corrections to $n-\pi^*$ ...]
 
 As a final comment, let us discuss the two singlet states of butadiene reported in Table \ref{tab:BigMol}.\cite{Maitra_2004,Cave_2004,Saha_2006,Watson_2012,Shu_2017,Barca_2018a,Barca_2018b,Loos_2019}
 As discussed in Sec.~\ref{sec:intro}, these corresponds to a bright state of $^1B_u$ symmetry with a clear single-excitation character, and a dark $^1A_g$ state including a substantial fraction of double excitation character (roughly $30\%$).
 Although they are both of $\pi \ra \pis$ nature, they are very slightly altered by dynamical screening with corrections of $-0.12$ and $-0.03$ eV for the $^1B_u$ and $^1A_g$ states, respectively. 
 The small correction on the $^1A_g$ state might be explained by its rather diffuse nature (similar to a Rydberg states). \cite{Boggio-Pasqua_2004}
 
-\xavier{DISCUSSION: To provide further insight into the magnitude of the dynamical correction to valence, Rydberg or CT excitations, let's consider a simple 2-level systems with (ij=hh) and (ab=ll) where (h,l) stand for HOMO and LUMO. The dynamical correction to the $HOMO \rightarrow LUMO$ transition reads:}
-\begin{equation*}
-	\W{hh,ll}{\text{stat}} - \widetilde{W}_{hh,ll}( \Om{1}{} )
-	=     4 \sERI{hh}{hl} \sERI{ll}{hl} \left( \frac{ \Om{hl}{1} - 2 \OmRPA{hl}{}}{\OmRPA{hl}{} ( \OmRPA{hl}{} - \Om{hl}{1} ) } \right),
-\end{equation*}
-where the only RPA excitation  energy is again the $HOMO-LUMO$  transition with m=hl (see equation 26). For CT excitations with vanishing (hl) overlap, $[hh|hl] \simeq 0$ and similarly  $[ll|hl] \simeq 0$ so that one can expect the dynamical correction to be weak. Similarly, Rydberg transitions characterized by a delocalized LUMO state, namely a small HOMO-LUMO (hl) overlap, are expected to undergo weak dynamical corrections.  The discussion on $\pi-\pi^*$ and $n-\pi^*$ transitions is certainly more complex and molecule-specific symmetry arguments must be invoked to understand the magnitude of the $[hh|hl]$ and $[ll|hl]$ terms. 
-
-[ Taking the simple case of a $\pi$-conjugated planar  molecule located in the (xy) plane, then $\pi-\pi^*$ transition will yield $[hh|hl]$ and $[ll|hl]$ terms built from $(p_z^2 | p_z^2)$ contributions on each atom. On the contrary, $n-\pi^*$ transitions will involve $(sp_z|ss)$ or $(sp_z|p_z
-^2)$ local contribution that are small by planar symmetry. Bah bon ben ... should lead to small corrections to $n-\pi^*$ ...]
-
-
-
 %%%%%%%%%%%%%%%%%%%%%%%%
 \section{Conclusion}
 \label{sec:conclusion}
@@ -1139,11 +1144,12 @@ The BSE formalism is quickly gaining momentum in the electronic structure commun
 It now stands as a genuine cost-effective excited-state method and is regarded as a valuable alternative to the popular TD-DFT method.
 However, the vast majority of the BSE calculations are performed within the static approximation in which, in complete analogy with the ubiquitous adiabatic approximation in TD-DFT, the dynamical BSE kernel is replaced by its static limit.
 One key consequence of this static approximation is the absence of higher excitations from the BSE optical spectrum.
-Following Strinati's footsteps, \cite{Strinati_1982,Strinati_1984,Strinati_1988} several groups have explored the BSE formalism beyond the static approximation by retaining (or reviving) the dynamical nature of the screened Coulomb potential  \cite{Sottile_2003,Romaniello_2009b,Sangalli_2011} or via a perturbative approach coupled with the plasmon-pole approximation. \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b}
+Following Strinati's footsteps, \cite{Strinati_1982,Strinati_1984,Strinati_1988} several groups have explored the BSE formalism beyond the static approximation by retaining (or reviving) the dynamical nature of the screened Coulomb potential \cite{Sottile_2003,Romaniello_2009b,Sangalli_2011} or via a perturbative approach coupled with the plasmon-pole approximation. \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b}
 
-In the present study, we compute 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 report a significant number of calculations for various molecular systems.
-Our calculations are 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.
+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.
+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.
 
 %%%%%%%%%%%%%%%%%%%%%%%%
 %\section*{Supplementary material}
@@ -1256,7 +1262,7 @@ one gets
 \end{split}
 \end{equation}
 %We now act on the $N$-electron ground-state wave function with
-Assuming now that $\lbrace \e{p} , \e{q} \rbrace$ are proper addition/removal energies, such as the $GW$ quasiparticle energies, one can use the following relations:  
+\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 relations
 \begin{subequations}
 \begin{align}
   e^{+i\hH \tau_{65} } \ha^{\dagger}_p \ket{N} &=