GW approximation for open-shell molecules: a first-principles study

The basic quantity in the GWA is the non-local frequency-dependent electron self-energy Σ_xc( r , r ', ω), which describes the exchange and correlation effects. It is defined by

$$\begin{equation}{G}^{-1}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )={G}_{0}^{-1}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )+\left({{\Sigma}}_{xc}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )-{{\Sigma}}_{0,xc}(\boldsymbol{r},{\boldsymbol{r}}^{\prime })\right),\end{equation} \tag{ 1 }$$

where G( r , r ', ω) and G₀( r , r ', ω) represent the time-ordered interacting and non-interacting Green's functions [3, 4, 35]. The latter is typically obtained from a single-particle effective Hamiltonian H₀ as follows [2]

$$\begin{equation}\begin{aligned}{G}_{0}^{-1}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )& =\omega -{H}_{0}(\boldsymbol{r},{\boldsymbol{r}}^{\prime })\\ & =\left\{\delta (\boldsymbol{r}-{\boldsymbol{r}}^{\prime })\left[\omega -{H}_{H}(\boldsymbol{r})\right]-{{\Sigma}}_{0,\text{xc}}(\boldsymbol{r},{\boldsymbol{r}}^{\prime })\right\},\end{aligned}\end{equation} \tag{ 2 }$$

where Σ_0,xc and H_H are an effective static xc operator and the Hartree Hamiltonian, respectively. For instance, if H₀ is provided by the Hartree–Fock (HF) method, the xc operator Σ_0,xc will be the customary Fock exchange operator Σ_x.

According to Hedin [1, 2], the electron self-energy Σ_xc(ω), including electron–electron interaction effects beyond the Hartree potential, can be approximated by

$$\begin{equation}{{\Sigma}}_{\text{xc}}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )=\underset{\eta \to {0}^{+}}{\mathrm{lim}}\enspace \frac{\mathrm{i}}{2\pi }{\int }_{-\infty }^{+\infty }G(\boldsymbol{r},{\boldsymbol{r}}^{\prime };\omega +{\omega }^{\prime })W(\boldsymbol{r},{\boldsymbol{r}}^{\prime },{\omega }^{\prime }){\text{e}}^{\text{i}\eta {\omega }^{\prime }}\enspace \mathrm{d}{\omega }^{\prime },\end{equation} \tag{ 3 }$$

where η is a positive infinitesimal shift away from the real frequency axis and W(ω) is the screened interaction, as will be discussed in section 2.4. Such a frequency (energy) dependence distinguishes the GWA from static mean-field methods, so that, the GW self-energy now contains a dynamical Coulomb-hole (COH) term plus a dynamically screened exchange (SEX) part. In the static approximation, the GW self-energy reduces to the static COHSEX, giving rise to an effective Hamiltonian containing the non-local screened-exchange plus a local term representing the Coulomb interaction of the quasiparticle with the screening charge distribution (hole) around it. For more details we refer the readers to reference [3–5, 10].

Equations (1)–(3) must be satisfied simultaneously giving rise to a self-consistent GW procedure that implies a high computational cost [8, 17, 21, 29]. Fortunately, many single-particle Hamiltonians feature rather accurate approximations Σ_0,xc to the xc operator, so that the calculation of the electron self-energy can be accurately approximated using the non-interacting Green's function (2) in equation (3). This perturbative procedure is known as the one-shot GWA or G₀ W₀. Moreover, it has been reported that the fully self-consistent GWA might not yield the best results, while a partial self-consistency or G₀ W₀ frequently results in a better agreement with the experiment [5, 22].

Besides the one-shot GWA, a diagonal approximation is often employed; so that, the interacting Green's function G^σ ( r , r ', ω) can be approximately expanded in terms of the eigenstates ψ^n,σ ( r ) of the single-particle Hamiltonian ${H}_{0}^{\sigma }$ (σ = ↑, ↓) as

$$\begin{equation}{G}^{\sigma }(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )=\sum\limits _{n}\frac{{\psi }^{n,\sigma }(\boldsymbol{r}){\psi }^{{\ast}n,\sigma }({\boldsymbol{r}}^{\prime })}{\omega -{E}^{n,\sigma }{\pm}\mathrm{i}\eta }.\end{equation} \tag{ 4 }$$

In this approximate ansatz, the quasiparticle energies E^n,σ do differ from the eigenenergies ${E}_{0}^{n,\sigma }$ of the single-particle Hamiltonian ${H}_{0}^{\sigma }$. Inserting the ansatz (4) and the exact Lehmann representation for the non-interacting Green's function ${G}_{0}^{\sigma }(\omega )$ into the Dyson equation (1), we obtain a fixed-point equation for the GW quasiparticle energies [3, 4]

$$\begin{equation}{E}^{n,\sigma }={E}_{0}^{n,\sigma }+\left\langle {\psi }^{n,\sigma }\left\vert {{\Sigma}}_{\text{xc}}^{\sigma ,0}({E}^{n,\sigma })-{{\Sigma}}_{0,\text{xc}}^{\sigma }\right\vert {\psi }^{n,\sigma }\right\rangle ,\end{equation} \tag{ 5 }$$

where ${{\Sigma}}_{\text{xc}}^{\sigma ,0}({E}^{n,\sigma })$ is the G₀ W₀ self-energy defined by replacing the interacting Green's function G^σ (ω) with the non-interacting Green's function ${G}_{0}^{\sigma }(\omega )$ in equation (3).

In this study, we used the graphical method to solve the nonlinear quasiparticle equation (5); details are given in the supporting information [36]. For the HOMO this method often leads to a unique solution, meaning that the intersection occurs in a pole-less region of the self-energy [19, 37]. As a result, the slight slope of the self-energy gives rise to a large quasiparticle weight, associated with a single pronounced peak in the diagonal elements of the spectral function, ${A}_{jj}^{\sigma }(\omega )=\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}\left\vert {G}_{jj}^{\sigma }(\omega )\right\vert$ [24, 38]. However, for some cases, e.g. when the initial mean-field energy takes place in the vicinity of the self-energy pole position, equation (5) might provide more than one solution [19]. Ionization of lower lying states can also result in multiple solutions due to the coupling between created hole state and the neighboring excitations, producing several final states [39].

For open-shell systems, in particular, the occurrence of multiple solutions must be carefully considered. As an example, Lischner et al [14] have shown that precise knowledge of the self-energy on a fine frequency grid can capture the multiplet splittings, in agreement with the photo-ionization spectroscopy. These authors also point toward the importance of partial self-consistent calculations to accurately obtain the position of the self-energy pole(s). Although this study is limited to the standard G₀ W₀ approach, we give an example of this multiplet splitting in section 4.2.2.

An optimal choice of the single-particle Hamiltonian is important for the accuracy of the one-shot GW results, because of the well-known starting point dependence of the method [3]. In this study, we use unrestricted HF (UHF) [40, 41] and unrestricted Kohn–Sham (KS) [42–44] methods. HF Hamiltonian is a natural choice for Hedin's GWA, because the exchange part of Hedin's self-energy (3) coincides with the Fock operator [2, 3]. In fact, the GW self-energy can be split into a static non-local Fock exchange term and a frequency-dependent correlation part as ${{\Sigma}}_{\text{xc}}^{\sigma }(\omega )={{\Sigma}}_{\text{x}}^{\sigma }+{{\Sigma}}_{\text{c}}^{\sigma }(\omega )$ [3, 5, 11, 19].

Among different types of functionals used within the KS method, hybrid functionals tend to be the most accurate for molecules [45]. Hybrid functionals make use of the non-local Fock exchange operator similar to the HF method, frequently in combination with an effective screened Coulomb interaction, plus a density functional to describe correlation effects. Therefore, it does not come as a surprise that hybrid functionals provide optimal starting points for the G₀ W₀ calculations [3, 18–21]. In this study, we employ both full-range and range-separated hybrid functionals. As representatives for the former, the fraction of exact exchange varies from 0.20 for B3LYP [46] to 0.25 for PBE0 [47]. Range-separated functionals are also represented by the CAM–B3LYP [48] and HSE06 [49], comprising different amounts of the HF exchange at short- and long-range.

In the random phase approximation, with purely Coulombic interactions, the screened interaction W( r , r ', ω) is a spin-independent operator [3, 9] and reads

$$\begin{equation}W(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )=v(\boldsymbol{r},{\boldsymbol{r}}^{\prime })+\iint v(\boldsymbol{r},{\boldsymbol{r}}^{\prime \prime }){\chi }_{0}({\boldsymbol{r}}^{\prime \prime },{\boldsymbol{r}}^{\prime \prime \prime },\omega )W({\boldsymbol{r}}^{\prime \prime \prime },{\boldsymbol{r}}^{\prime },\omega )\mathrm{d}{\boldsymbol{r}}^{\prime \prime }\enspace \mathrm{d}{\boldsymbol{r}}^{\prime \prime \prime },\end{equation} \tag{ 6 }$$

where χ₀( r , r ', ω) is the non-interacting density response function [50, 51], which has an explicit expression in terms of the single-particle eigenstates ψ^n,σ ( r ) and their eigenenergies ${E}_{0}^{n,\sigma }$ [52–54]

$$\begin{equation}{\chi }_{0}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )=\sum\limits _{n,m,\sigma }\frac{({f}^{n,\sigma }-{f}^{m,\sigma }){\psi }^{n,{\sigma }^{{\ast}}}(\boldsymbol{r}){\psi }^{m,\sigma }(\boldsymbol{r}){\psi }^{m,{\sigma }^{{\ast}}}({\boldsymbol{r}}^{\prime }){\psi }^{n,\sigma }({\boldsymbol{r}}^{\prime })}{\omega -{E}_{0}^{m,\sigma }+{E}_{0}^{n,\sigma }+\text{i}\eta },\end{equation} \tag{ 7 }$$

where f^n,σ and f^m,σ are the occupation factors, n and m represent the MO indices, and ${E}_{0}^{n,\sigma }$ and ${E}_{0}^{m,\sigma }$ are MO energies, respectively.

One finds that the screened interaction W( r , r ', ω) is obtained from the spin-averaged χ₀( r , r ', ω) as defined above. Consequently, the spin-dependence of the G₀ W₀ self-energy ${{\Sigma}}_{\text{xc}}^{\sigma ,0}(\boldsymbol{r},{\boldsymbol{r}}^{\prime },\omega )$ stems from that of the non-interacting Green's function G₀ [55], where G₀ in turn inherits the spin indices of the single-particle Hamiltonian H₀, as defined in equation (2). In other words, the spin structure of the G₀ W₀ self-energy is determined by that of the Hamiltonian of the starting mean-field, which in the cases considered here is spin-diagonal.

It is well-known that a comparison between theoretical IEs and experimental values is a difficult test for electronic structure methods [19, 20, 28]. In order to demonstrate the validity of the GWA in the calculation of IEs and to evaluate its accuracy with respect to other approaches, we also computed IEs of our test-set molecules within the Δ-framework (usually is referred to as ΔSCF in the context of mean-field schemes), that is, the energy difference between two separate calculations for the neutral and cationic species [5, 56]. To this aim, we use several unrestricted correlated methods such as CCSD, CCSD(T), second-order Møller–Plesset perturbation theory (MP2), configuration interaction limited to single and double excitations (CISD), and density functional theory (DFT) with the B3LYP functional. Exchange-only methods such as UHF and restricted open-shell HF (ROHF) [57] were also employed.

Among all, we use the ΔCCSD(T) results as a reference since CCSD(T) has proven to be an accurate method for small to medium-sized closed-shell molecules and it remains the gold standard method in quantum chemistry [17, 18, 20, 29, 58, 59]. Likewise, the CCSD(T) approach is currently employed in many occasions as the reference method for open-shell atoms and molecules [58, 60]. More importantly, the G₀ W₀ method discussed in this paper is in essence a single reference approach, that is, a perturbative correction of the single configuration HF or KS models. Therefore, it seems reasonable to test its performance against highly accurate single reference methods, such as CCSD(T).

In the case of small finite systems, Δ-methods often lead to the IEs similar to the GWA quasiparticle excitation energies [10, 29, 56]. Nonetheless, a straightforward application of the Δ-methods is limited to the lowest IE. Calculations of deeper levels require a well-defined and efficient scheme to constrain the hole, left in the system after the excitation, in a particular state. This might be difficult in some cases [56, 61]. Moreover, Δ-methods present conceptual problems to deal with infinite periodic systems [5, 10], limiting their range of applicability. In contrast, the physics behind the GWA in terms of dynamical screened interactions in both finite and extended phases provides a general framework [5].

We would like to comment on the differences in the computational cost between the G₀ W₀ method and other applied approaches. Numerical implementations of the G₀ W₀ approach typically scale as O(N³)–O(N⁴) [62], where N is the number of basis functions. The computational demands in wave function-based methods are typically much higher in terms of processor time and storage requirements and present a worse scaling with basis size, e.g. O(N⁷) for CCSD(T) [63], prohibiting calculations for molecules containing more than a few tens of atoms [64].

In order to test the validity of our G₀ W₀ approach, especially in those cases where mean-fields and G₀ W₀ provide different qualitative results, we compare MO energy diagrams with Dyson orbitals obtained within the coupled-cluster formalism. Dyson orbitals are typically used to compute Compton scattering profiles [65], electron momentum spectra [66], and represent a very powerful interpretation tool [67–69].

The Dyson orbital ${\psi }_{n}^{d}(\mathbf{r})$ is a one-electron function that characterizes the probability amplitude distribution of obtaining electrons with a given binding energy during the ionization of a molecular species. It is defined as the overlap between the ground-state N-electron wave function ${{\Psi}}_{0}^{N}$ and the wave functions of the different states of the system containing N − 1 electrons ${{\Psi}}_{n}^{N-1}$

$$\begin{equation}{\psi }_{n}^{d}({\mathbf{r}}_{1})=\sqrt{N}\int {{\Psi}}_{0}^{N}({\mathbf{r}}_{1},{\mathbf{r}}_{2},\dots ,{\mathbf{r}}_{n}){{\Psi}}_{n}^{N-1}({\mathbf{r}}_{2},\dots ,{\mathbf{r}}_{n})\mathrm{d}{\mathbf{r}}_{2},\dots ,\mathrm{d}{\mathbf{r}}_{n}.\end{equation} \tag{ 8 }$$

In order to go beyond the Koopman's approximation [70] and account for electron correlation effects on the Dyson orbitals, we employ coupled-cluster theory with single and double excitations for the ground state, while the ionized (N − 1 electron) state is obtained with the ionized version of the equation-of-motion formalism (EOM-IP-CCSD) [71].

In this subsection, we evaluate the vertical IE of 42 open-shell molecules. IEs were calculated within the Δ-method and G₀ W₀ on top of different unrestricted mean-field solutions. Additionally, the electron affinities of the cationic species (EA⁺) are computed. As discussed in section 2.5, the ΔCCSD(T) method was chosen as the reference. Experimental IEs as reported in the NIST database [80] are also collected in the supporting information [36]. We find that IEs calculated within the ΔCCSD(T) agree with the available experiments, measured the vertical excitation energies, with a mean absolute error (MAE) of 0.11 eV.

For brevity, a color map representation of the deviations with respect to the ΔCCSD(T) reference is presented in figure 1 and statistical deviations were outlined in table 1. The reader is referred to the supporting information [36] for the molecule-resolved numerical results. In the case of Δ-methods, while UHF and ROHF variants show the highest error rates due to the lack of explicit electron correlation, unrestricted MP2, CISD, and B3LYP on average provide results in agreement with the reference. Note that, although B3LYP functional in the Δ-scheme leads to reasonable IEs, the errors of B3LYP KS eigenvalues can be rather large, in the range of several eVs.

Zoom In Zoom Out Reset image size

On the other side, G₀ W₀@UHF (G₀ W₀ using the UHF starting point) presents an MAE of 0.55 eV, which is reduced to 0.12, 0.14, and 0.20 eV in @PBE0, @CAM–B3LYP, and @B3LYP, respectively. These deviations are comparable to those found for closed-shell molecules [11, 19–21, 29]. Not surprisingly, the G₀ W₀ results show an undesired dependence on the starting point as an accuracy-limiting factor. G₀ W₀@UHF and @CAM–B3LYP tend to overestimate the IEs with a mean signed error (MSE) of 0.55 and 0.12 eV while the computed IEs by G₀ W₀@B3LYP and @PBE0 are underestimated with an MSE of −0.18 and −0.09 eV, respectively.

Figure 2 shows the error distribution of the computed IEs using G₀ W₀ as a function of the mean-field starting point. While IEs calculated by G₀ W₀@UHF features a maximal absolute deviation (MAD) of 1.46 eV, G₀ W₀@B3LYP and @PBE0 and @CAM–B3LYP for the majority of the molecules show a slight deviation. Particularly, IEs computed within G₀ W₀@B3LYP and @PBE0 are underestimated with the exception of those obtained from effectively one-electron atoms, namely, Li and Na. This seems to point toward an overestimation of the correlation energy in such atoms, which can be tentatively ascribed to the self-screening of the single valence electron [81]. In the supporting information [36], moreover, we discussed the low quality of the quasiparticle solution for these two atoms. In the literature [31], it is indicated that for such atoms G₀ W₀@HF can provide IE comparable to those obtained through the self-consistent GW flavor.

Zoom In Zoom Out Reset image size

When carrying unrestricted calculations, especially UHF, it is necessary to check the degree of spin contamination. The spin contamination can be judged by the expectation value of total spin operator $\left\langle {S}^{2}\right\rangle$ [24, 27, 41, 82]. Hence, we explored the difference between calculated $\left\langle {S}^{2}\right\rangle$ and the ideal s(s + 1) value, where s is the spin quantum number. Analysis of ⟨S²⟩ reveals that the spin contamination of the investigated molecules is small. For the systems studied here, moreover, we found no clear connection between the spin contamination at the starting point solutions and the error-rate of the calculated G₀ W₀–IEs. More details, explaining the equations to compute $\left\langle {S}^{2}\right\rangle$ along with the degree of spin contamination, are given in the supporting information [36].

As discussed in section 2.2, the solution to the quasiparticle equation might not be unique. We confirm for the majority of our test-set systems the graphical method on a fine frequency grid provides a single solution for the HOMO with a pronounced spectral weight of around 0.9. Among all, only for a few systems and only at G₀ W₀@B3LYP and @PBE0, we found multiple solutions. In such systems, however, the spectral weights are rather different, favoring in general one solution against all others. Here, we only reported the solution with the largest weight in the spectral function and discussed cases with multiple solutions in the supporting information [36].

To further examine the quality of the G₀ W₀ approximation, we calculate the EA⁺. Detailed results are listed in the supporting information [36], where figure C-1 presents a summary of the deviations between IEs and EA⁺s for all the starting points. Here, we provided the reader with the statistical deviations collected in table 2 as well as a color map representation of the absolute error with respect to the reference in figure 1. Clearly, the MAE of the EA⁺s computed at the G₀ W₀@UHF is improved as compared with the corresponding deviation in table 1. Moreover, EA⁺s are underestimated with an MSE of −0.23 eV, which contracts with the systematic overestimation of the IEs for the UHF starting point (MSE = 0.55). This is in line with the well-known tendency of the HF, inherited by the G₀ W₀ .

Similarly to the G₀ W₀@UHF case, we observe a sign reversal of the MSE of the EA⁺s calculated from B3LYP and PBE0 starting points as compared to the MSE of the IEs (see tables 1 and 2). However, in contrast to the case of UHF, for these two hybrid functionals the MAE increases with respect to that found for the IEs of the neutral molecules. This is consistent with a larger error in the energy position of the LUMO in the initial DFT calculations, as compared to that of the HOMO, and points to the importance of an improved description of exchange in correcting the energies of the occupied states. G₀ W₀@CAM–B3LYP, on the other hand, provides a similar MAE for EA⁺s (0.17 eV) and IEs (0.14 eV) in which the MSE is positive in both cases. This seems to indicate that range-separated CAM–B3LYP functional as a starting point for the G₀ W₀ provides a more balanced description in terms of localization and delocalization errors than the other hybrid functionals considered here [30, 33, 34].

To gain an insight into the role of the exact exchange, we benchmarked the IEs and EA⁺s of our test-set molecules within the tuned PBE0 functional, where the fraction of exchange α varies from 0.35 to 0.75. The corresponding data are given in the supporting information [36], table D-1. As compared with the results obtained by G₀ W₀ starting from standard PBE0 (α = 0.25), we find all statistical deviations have slightly improved when the content of exchange is in the range of 0.35 to 0.4. For example, a fraction of 0.35 results in an MSE of 0.02 and 0.09 eV for the IEs of neutral molecules and the EA⁺s of the cations, respectively. When the exchange fraction is more than 50%, on the contrary, all deviations increase and provide a positive (negative) MSE for the IE (EA⁺). This indeed reflects the general rule of thumb that a larger content of exact exchange usually leads to an overestimation (underestimation) energies of the occupied (unoccupied) states.

To sum up, this benchmark reveals that the G₀ W₀ scheme can treat open-shell molecules with a reasonable accuracy comparable to that of CCSD(T), while its performance in terms of the computational cost is favorable. G₀ W₀ on top of the optimal starting points, i.e. hybrid functionals, yields an MAE of 0.1–0.2 eV for IEs, while EA⁺s show a slightly larger MAE of 0.2–0.3 eV. For G₀ W₀ @UHF, on the other hand, the MAE was reduced from 0.55 eV to 0.4 eV when EA⁺ of the cations are computed instead of the IEs of the neutral molecules. We ascribed this apparent improvement to the overestimation of the role of exchange in UHF, that gives rise to too large IEs. Among the hybrid functionals studied here, long-ranged CAM-B3LYP provides the smallest discrepancy of 0.03 eV between the MAE of IE and MAE of EA⁺.

The interpretation of photo-electron spectroscopy is frequently performed in terms of effective independent electron theories where the concept of MO appears naturally. Thus, the energies and symmetries of Frontier MOs become central to understanding the interaction of molecules with light. The connection of such single-electron description with a many-body picture has been discussed by many authors [37], pointing out that the quasiparticle concept can fail in the presence of strong electron–electron interactions or high excitation energies for which many different excitations in the system (including those with the intrinsic collective character) can contribute [22, 37, 39, 83, 84]. In the present case, however, we focus on the lowest photo-ionization energies and explore the accuracy of our mean-field starting points and G₀ W₀ to describe them, in particular the relative energy positions of MOs with different character. Our connection to a full many-body description is provided by the Dyson orbitals obtained from CCSD calculations, which are expected to include a fair account of many-body effects both in their excitation energies and spatial distributions [84]. We find that, while the character of the MOs obtained by HF and KS is in general in good agreement with the Dyson orbitals, their energy ordering is different in several cases. Interestingly, as shown below, G₀ W₀ is able to recover the correct energy order of the Frontier MOs. In the following, we focus our attention on three specific cases, namely, CN, NH₂, and O₂ molecules and discuss the sequence of spin-up MOs.

4.2.1. First and second IEs of molecular radicals

The cyano radical (CN) is a paramagnetic system exhibiting one unpaired electron in the σ_2p HOMO with the fully occupied (degenerate) π_2p orbitals at lower energy. Figure 3(a) illustrates the energy eigenvalues of these orbitals calculated by different methods. Dyson orbital energies are also given as reference (solid lines). The reference energies show a separation of 0.88 eV between the σ_2p (HOMO) and π_2p (HOMO-1), which agrees with previous results [85–89]. This is clearly a challenging system for mean-field methods, and we see both HF and KS with different xc functionals present incorrect energy ordering of MOs.

Zoom In Zoom Out Reset image size

Within the UHF, σ_2p eigenvalue is ∼1.5 eV larger than the corresponding reference, while the energy of the degenerate π_2p level is surprisingly underestimated by about 0.6 eV. Such high errors lead to an inverted energy ordering of orbitals with a large σ_2p –π_2p separation of 1.2 eV. DFT using three different hybrid functionals, on the other hand, systematically underestimates the energies of both orbital types. This underestimation is significantly higher for the degenerate π_2p orbitals than for σ_2p , resulting in the wrong ordering, but with a smaller gap (∼0.3 eV) than that of UHF.

Regardless of the starting-point, G₀ W₀ correction (shaded area) not only leads to the correct MO ordering in CN but also estimated quasiparticle orbital energies are much closer to the reference. Particularly, energy separation between σ_2p –π_2p obtained from the G₀ W₀ starting from mean-field approximations with high contents of Fock exchange, namely, UHF and CAM-B3LYP, are very small (∼0.2 eV), since both methods tend to overestimate σ_2p and underestimate π_2p energies. In contrast, G₀ W₀ calculations starting from either PBE0 or B3LYP provide improved results, with an energy separation around 0.5 eV.

Next, we turn our attention to the study of the amino radical NH₂. The singly occupied HOMO of NH₂ corresponds to the 2p nitrogen orbital perpendicular to the molecular plane belonging to the B₁ irreducible representation of the C_2v molecular symmetry (1b₁ orbital). The lower orbital is obtained as the bonding interaction between the 2p_z nitrogen orbital (aligned with the twofold rotational axis) and the in-phase combination of the 1s orbitals on the two hydrogen atoms, which belongs to the totally symmetric representation (3a₁ orbital).

As depicted in figure 3(b), UHF predicts the inverted order for the 1b₁ and 3a₁ orbital energies. Like in CN, G₀ W₀@UHF corrects the UHF orbital energies with a final one-electron energy ordering in agreement with the reference Dyson orbital energies. DFT calculations, on the other hand, correctly predict the ordering of MOs. However, attributed eigen-energies are significantly underestimated (more than 3.7 eV for B3LYP). Here, the G₀ W₀ correction greatly improves the accuracy of orbital energies with an MAE of 0.22 and 0.18 eV for the HOMO and HOMO-1, respectively.

4.2.2. Swapping of IEs in spin-triplet molecules

In the triplet O₂ molecule, with a ³Σ_g symmetry ground-state, because of a relatively large gap between 2s and 2p states of the oxygen atom, no hybridization occurs between s and p atomic orbitals. Therefore, σ_2p orbital is energetically situated below the two-fold degenerate π_2p and ${\pi }_{2p}^{{\ast}}$ orbitals [90–93]. The same ordering is obtained by the reference, where Dyson orbital energies are determined to be −12.55, −17.73, and −19.46 eV for ${\pi }_{2p}^{{\ast}}$ (HOMO), π_2p (HOMO-1), and σ_2p (HOMO-2), respectively. As shown in figure 4, neither UHF nor DFT using hybrid functionals can capture the correct energy ordering of the HOMO-1 and HOMO-2 orbitals. UHF overshoots π-type orbital energies leading to a large separation between σ_2p and π_2p (2.1 eV) and wrong ordering that is not corrected even at the G₀ W₀@UHF level. Similarly, hybrid functionals used within the DFT considerably underestimate σ_2p energy and provide a vanishingly small energy spacing between σ_2p and π_2p orbitals.

Zoom In Zoom Out Reset image size

Quasiparticle energies obtained from G₀ W₀ started from hybrid functionals interchange the energy order of σ_2p and π_2p orbitals, providing the same as that given by the reference. From a quantitative point of view, G₀ W₀@CAM–B3LYP gives the quasiparticle σ_2p energy very close to the reference, while estimated IEs of π_2p and ${\pi }_{2p}^{{\ast}}$ are overestimated by 1 and 0.36 eV, respectively. Conversely, G₀ W₀@PBE0 yields the closest energies of both π-type orbitals to the reference values, whereas the energy of the σ_2p orbital is underestimated by 0.62 eV. Eventually, G₀ W₀@B3LYP and @HSE06 accurately estimate the HOMO energy while IEs of π_2p and σ_2p orbitals are underestimated by around 0.5 and 0.9 eV, respectively. More details are given in the supporting information [36].

It is clear that diverse MO types respond differently to the GW correction. It has been recently shown that the GW self-energy often provides a poor estimation of σ_p orbital energies [94]. The self-screening problem of the GWA can also affect the energy level of orbitals with different bonding types [17, 95]. Therefore, to quantify the origin of such errors in the computed orbital energies of the oxygen molecule, we analyze the behavior of the two GW self-energy components Σ(E) = Σ_x + Σ_c(E), the exchange and correlation terms, for the different starting points.

Figures 5(a)–(c) show the diagonal expectation values of the exchange contribution Σ_x to the self-energy for the three topmost majority-spin occupied orbitals, i.e. ${\pi }_{2p}^{{\ast}}$, π_2p , and σ_2p in the O₂ molecule in its triplet ground-state configuration. For each of the mean-field methods, we evaluate the expectation values of Σ_x as $\left\langle {\psi }^{m}\vert {{\Sigma}}_{\text{x}}\left[{\rho }^{m}\right]\vert {\psi }^{m}\right\rangle$, where ψ^m and ρ^m are the corresponding mean-field's wave function and one-particle density matrix, respectively. For each orbital, the reference is given by the expectation value of the exact-exchange with the CCSD density ρ^CCSD and the corresponding Dyson orbitals ψ^d in the bra-ket, $\left\langle {\psi }^{d}\vert {{\Sigma}}_{\text{x}}\left[{\rho }^{\text{CCSD}}\right]\vert {\psi }^{d}\right\rangle$. Notice that the density matrix obtained from CCSD calculations, taking into account the contributions from single and double excitations, contains high quality information about the electron correlations in the system. Therefore, the difference between the mean-fields' exchange energy and the reference for each orbital, as shown in figure 5, reveals the effect of electron correlations beyond the mean-field theory. A detailed comparison of the contribution of density matrices can be found in the supporting information [36], figure F-1, where $\left\langle {\psi }^{m}\vert {{\Sigma}}_{\text{x}}\left[{\rho }^{\text{CCSD}}\right]\vert {\psi }^{m}\right\rangle$ values are also presented.

Zoom In Zoom Out Reset image size

Results show that all exchange energies are overestimated with respect to the reference. While the expectation values of the exchange energy for σ_2p , π_2p , and ${\pi }_{2p}^{{\ast}}$ orbitals within the DFT solutions are overestimated with a mean error of 0.54, 1.22, and 0.66 eV, the UHF solution leads to an error of 0.53, 1.86, and 1.2 eV, respectively. Such large deviations, particularly for π-type orbitals, directly reflect the origin of the incorrect spectra at the UHF level, shown in figure 4.

Figure 5(d) features the diagonal expectation values of the dynamic correlation part of the G₀ W₀ self-energy $\left\langle {\psi }^{m}\vert {{\Sigma}}_{\text{c}}(E)\vert {\psi }^{m}\right\rangle$, where E is the corresponding MO energy at the G₀ W₀ level. The reference values are obtained by subtracting the Dyson orbital energies in equation (8) from Fock energies $\left\langle {\psi }^{d}\vert F\left[{\rho }^{\text{CCSD}}\right]\vert {\psi }^{d}\right\rangle$ for each orbital. Clearly, the correlation energies attributed to π_2p -type orbitals are higher than σ_2p . This can be ascribed to the larger extension and polarizability of the π_2p orbitals. Notice that the G₀ W₀ provides rather higher correlation energy for π_2p (HOMO-1) as compared to the ${\pi }_{2p}^{{\ast}}$ (HOMO). Adding these positive correlation energies pulls up π-type orbital energies in the spectra, leading to more accurate energies associated with the correct MO picture for hybrid orbitals, as shown in figure 4. In comparison with reference values, however, one sees that the correlation energy of the π_2p orbital is underestimated in direct relation to the fraction of the exact exchange. Therefore, hybrid functionals with a higher fraction of long-range exchange incorrectly provide smaller separation between σ_2p and π_2p . This also explains the dramatic failure of the G₀ W₀@UHF, where estimated correlation energy for UHF-π_2p orbital is drastically insufficient to provide the correct spectra.

It is worth mentioning here the multiple solutions obtained for HOMO-2 orbital of the triplet oxygen molecule. While removal one electron from ${\pi }_{2p}^{{\ast}}$ (HOMO) features a single pronounced peak in the spectral function, ionization of up-spin σ_2p orbital results in two peaks as described in the supporting information [36]. This is consistent with a previous study [14], where the two solutions are attributed to the splitting between two possible ionic states, namely, doublet (²Σ_g ) and quartet (⁴Σ_g ), which is experimentally determined to be 2.3 eV [96]. We find values of 2.03 and 2.56 eV for this splitting at the G₀ W₀@B3LYP and @PBE0, respectively. Notice that this does not qualitatively affect the spectra shown in figure 4 since the second solution occurs at lower energy (∼21 eV).

Summarizing the discussion, we emphasize that MOs with different symmetry respond with varying sensitivity to different approximations to describe the exchange and correlation effects. In the case of CN and NH₂ molecules, we showed that mean-field solutions might predict an incorrect sequence of the Frontier orbitals due to the significant over- or underestimation of the IEs, while the GW correction can provide reasonably accurate orbital energies associated with the correct ordering. Likewise, in the example of O₂, we found that the G₀ W₀ can correct the overestimation of the exchange energy by adding a well-set correlation, thanks to its energy-dependent screening. As a result, the dynamic contribution of the GW self-energy can result in an electronic structure that is qualitatively different from that of mean-field methods.