Introduction to the Fifth-rung Density Functional Approximations: Concept, Formulation, and Applications

DHAs are now the leading actor in quantum chemistry. The name of “doubly hybrid” was originally used by Truhlar and co-workers for their empirical multi-coefficient methods of MC3BB and MC3MPW [12], where the DFT total energy is linearly mixed with the total energy of the simplest wavefunction-based correlation method, i.e. the second-order Møller-Plesset perturbation theory (MP2). Later, Grimme introduced another kind of doubly hybrid strategy to connect MP2 and meta-GGAs [13]. B2PLYP is the first DHA of this family [13]. Popular B2PLYP-type DHAs include DSD-BLYP [14], $\omega$B97X-2 [15], and so forth [16].

In fact, DHAs can be derived in the generalized KS-DFT framework and viewed as an empirical extension of the Görling-Levy second-order perturbation theory (GL2) [17] along the adiabatic-connection (AC) approach [18]. It was first proposed by Zhang and co-workers in 2009 [19], which emphasizes on the use of accurate density and orbitals of KS non-interacting systems for evaluating the total energy of DHAs [20, 21]. XYG3 is the first DHA of this family [19]. Popular XYG3-type DHAs includes XYGJ-OS [22], xDH-PBE0 [23], $\omega$B97M(2) [24], xrevDSD-PBEB86-D4 [25], and so forth.

The rapid growth of computation capacity has been boosting the development of comprehensive databases for main-group chemistry with accurate benchmark reference. Two well-established databases are the MGCDB84 set of the Head-Gordon group with about 5000 test cases [24] and the GMTKN55 set of the Grimme group with 1505 test cases [26], respectively. A serial of recent benchmark works using these databases repeatedly demonstrated that the DHAs present significant and consistent improvement over the conventional DFAs on the first four rungs, and are competent for various kinds of chemical interactions of the main-group chemistry [25, 20]. The XYG3-type DHAs stand out in this round of benchmark according to benchmark results obtained independently by different research groups. The top-class performers all belong to the XYG3 family of DHAs, including the XYG7 of the Xu group, the $\omega$B97M(2) of the Head-Gordon group, and the xrevSDS-PBEB86-D4 of the Martin group. Their performance on main-group chemistry has been comparable to, or even better than the high-level wavefunction-based methods[20].

At present, the mainstream DHAs aforementioned use the second-order perturbation energy (PT2) to capture the electron correlation effect that involves the unoccupied orbitals. Recently, Santra and co-workers observed further notable improvement by introducing the third-order perturbation energy (PT3) in the context of XYG3-type doubly hybrid framework [27]. However, the computational cost increases significantly. Meanwhile, it is well-documented that the perturbation algorithms at any finite order (PTn) completely fail for describing the strong correlation effect, which is closely related to the so-called static correlation in quantum chemistry. Therefore, the DHAs including the PT2 and/or PT3 correlation contributions cannot be used to study the bond dissociation problem and the complex systems involving transition metals. A number of studies have begun to address this limitation of current DFAs by introducing some kinds of renormalization algorithms to the standard PT2 with little extra computational costs. As will be introduced in the next section, the random-phase approximation (RPA) is another successful fifth-rung DFA in materials science. The RPA-type correlation can be viewed as a renormalization of the direct term of PT2. RPA and its variants have also been taken as the candidate to replace the PT2 in the DHAs [28, 29, 30, 31, 32].

RPA is a fundamental concept that plays a central role in many-body physics. The concept was developed by Bohm and Pines [33, 34] in the early 1950’s in an endeavor to describe the cohesive properties of the so-called jellium – interacting electrons moving in a background of uniform positive charge. Using a Hamiltonian formulation of interacting many-electron system, Bohm and Pines were able to decouple the collective motion of electrons – the plasma oscillations – from their individual motions, a procedure named as RPA. It was soon recognized by Hubbard [35] that the original RPA formulation is equivalent to the infinite summation of ring diagrams from the viewpoint of diagrammatic many-body perturbation theory [36]. Since then, the RPA concept has gone beyond the realm of condensed matter physics and significantly influenced all branches of physics.

Although the RPA concept can be applied to any interacting many-particle systems (for its applications in nuclear physics, see e.g., Ref. [37]), the next key step towards applying RPA to real materials was the incorporation of RPA into the KS-DFT framework in the 1970’s [38, 39, 18]. This formulation turned RPA into a first-principles electronic-structure method, suitable for computing the ground-state energy of real materials. Within the KS-DFT framework, RPA can be viewed as a fifth-rung approximation to the exchange-correlation energy functional, according to Jacob’s ladder of DFT [3]. However, the application of RPA to realistic systems was impeded by its high computational cost and the lack of efficient algorithms at the time. The first application of RPA to small molecules only appeared in early 2000’s, carried out by Furche [40] and by Fuchs and Gonze [41]. Since then, RPA has been applied to a variety of systems including atoms [42, 43], molecules [40, 41, 44, 45, 46, 47, 48], solids [49, 50, 51, 52, 53, 54, 55, 56], surfaces [57, 58], interfaces [59, 60], layered materials [61], and defects [62, 63]. The consensus arising from these studies is that RPA is capable of describing the delicate energy differences in complex chemical environments [57, 58, 54, 56, 64, 65, 66], the correct asymptotic behavior of van der Waals (vdW) complexes [67, 68, 69] and layered materials [61, 70], and the correct dissociation limit of closed-shell molecules [40, 41]. Evidences show that RPA can provide unprecedented accuracy compared to lower-rung DFAs at tractable computational cost. As such, RPA is expected to play an increasingly more important role in computational materials science, with the rapid development of more efficient algorithms and the availability of more powerful computing resources.

In a review paper [71] published in 2012, Ren et al. discussed the history of the RPA concept, its formulation as a first-principles method , and its applications in quantum chemistry and computational materials science up to that time. These points were nicely summarized by David Pines in his recent review paper titled as “Emergent behavior in strongly correlated electron systems" [72]. In particular, Pines noted that,

Sixty-plus years later, the RPA continues to play a significant role in nuclear physics [66], bosonic field-theory [67], the quarkgluon plasma [68], many-fermion solvable models [69], and especially in computational chemistry and materials science. A recent review by Ren et al [70], to which the interested reader is referred, describes the impact of the RPA in the theoretical chemistry and materials science community, cites some thirty articles that indicate the renewed and widespread interest in the RPA during the period 2001-2011, discusses how it enables one to derive the $1/r^{6}$ interaction between spatially separated closed shell electron systems, and, shows, in some detail, how the RPA enables one to go beyond density functional theory in computing ground state energies.

This paragraph highlights the far-reaching impact of RPA in a variety of fields, in particular in computational chemistry and materials science. It is worthwhile to mention that several other, complementary RPA review papers [73, 74] also appeared around that time, where the theoretical foundation and applications of RPA are discussed from different perspective. More recent account of RPA of review character can be found in Ref. [75].

In this review, we will first give an account of the theoretical foundation of DHAs and RPA, highlighting their unique role in DFT. Computational schemes beyond PT2-based DHAs and RPA will be also briefly discussed in this review. This is followed by a description of the key algorithm of implementing DHAs and RPA using the resolution-of-identity (RI) technique, within an atomic-orbital basis-set framework. We then present some prototypical applications illustrating the usefulness of these fifth-rung functionals particularly in computational materials science. Note that, the general performance of DHAs have been extensively discussed for finite molecules, but only a limited number of their applications on solids were reported, which are introduced in this review. Most recent efforts devoted to further developing the theoretical and computational aspects of DHAs and RPA-based functionals as well as extending their capabilities will be briefly mentioned and commented, before we conclude the review.

Within the adiabatic connection (AC) approach to DFT [39, 18], one considers a continuous set of fictitious Hamiltonian’s,

$$\begin{split}\hat{H}(\lambda)&=\hat{T}+\hat{V}_{\lambda}+\lambda\hat{V}_{ee}=-\sum_{i}^{N}\frac{\nabla^{2}_{i}}{2}+\sum_{i}^{N}v^{\mathrm{aux}}_{\lambda}(\hat{{\bf r}}_{i})+\frac{\lambda}{2}\sum_{i\neq j}^{N}\frac{1}{|\hat{{\bf r}}_{i}-\hat{{\bf r}}_{j}|}\,,\end{split}$$

(17)

which connects $\hat{H}_{0}$ at $\lambda=0$ with $v^{\mathrm{aux}}_{\lambda=0}=v_{\mathrm{KS}}=v^{\mathrm{ext}}+v_{\mathrm{H}}+v_{\mathrm{xc}}$ (Eq. 12) and $\hat{H}$ at $\lambda=1$ with $v^{\mathrm{aux}}_{\lambda=1}=v^{\mathrm{ext}}$. The Hamiltonian $\hat{H}(\lambda)$ for $0<\lambda<1$ describes a collection of particles moving under the auxiliary external potential $v^{\mathrm{aux}}_{\lambda}({\bf r})$ and interacting with a scaled Coulomb interaction $\lambda/|{\bf r}-{\bf r^{\prime}}|$. The auxiliary potential $v^{\mathrm{aux}}_{\lambda}({\bf r})$ ($0<\lambda<1$) is chosen such that the density of $\lambda$-scaled systems is kept at the physical density, i.e., $n_{\lambda}({\bf r})=n({\bf r})$, along the AC path. By utilizing the coordinate scaling technique, it can be shown that

$$v^{\mathrm{aux}}_{\lambda}({\bf r})=v_{\mathrm{KS}}({\bf r})-\lambda\left[v_{\mathrm{H}}({\bf r})+v_{\mathrm{x}}({\bf r})+\lambda\frac{\delta E_{c}[n_{1/\lambda}]}{\delta n({\bf r})}\right].$$

(18)

Here $E_{c}[n_{1/\lambda}]$ is the correlation energy with respect to the scaled density (see Ref. [17, 92] for more detailed definition and derivation). $E_{c}[n_{1/\lambda}]=E_{c}[n]$ when $\lambda=1$; while in the high density limit with $\lambda=0$, $E_{c}[n_{1/\lambda}]=0$. For the fictitious system of $\hat{H}(\lambda)$ with fixed density $n({\bf r})$, the ground-state wavefunction can be denoted as $|\Psi_{\lambda}^{n}\rangle$,

$$\hat{H}(\lambda)|\Psi_{\lambda}^{n}\rangle=E_{\lambda}|\Psi_{\lambda}^{n}\rangle$$

(19)

with the normalization condition of $\langle\Psi_{\lambda}^{n}|\Psi_{\lambda}^{n}\rangle=1$. Here $E_{\lambda}$ is the ground-state energy along the AC path. $E_{0}=E_{\mathrm{KS}}$ is nothing but the ground-state energy of non-interacting KS system (Eq. 15) , while $E_{1}=E$ is the ground-state energy of the fully interacting system (Eq. 3).

The Hellman-Feynman theorem implies that

$$\begin{split}\frac{dE_{\lambda}}{d\lambda}&=\left\langle\Psi_{\lambda}^{n}\left|\frac{\partial\hat{H}(\lambda)}{\partial\lambda}\right|\Psi_{\lambda}^{n}\right\rangle=\left\langle\Psi_{\lambda}^{n}\left|\frac{\partial\hat{V}_{\lambda}}{\partial\lambda}\right|\Psi_{\lambda}^{n}\right\rangle+\left\langle\Psi_{\lambda}^{n}\left|\hat{V}_{ee}\right|\Psi_{\lambda}^{n}\right\rangle\\ &=\left\langle\Psi_{\lambda}^{n}\left|\sum_{i=1}^{N}\frac{\partial v^{\mathrm{aux}}_{\lambda}(\hat{{\bf r}}_{i})}{\partial\lambda}+\frac{1}{2}\sum_{i\neq j}^{N}\frac{1}{|\hat{{\bf r}}_{i}-\hat{{\bf r}}_{j}|}\right|\Psi_{\lambda}^{n}\right\rangle\\ &=\int d{\bf r}n({\bf r})\frac{\partial v^{\mathrm{aux}}_{\lambda}({\bf r})}{\partial\lambda}+\frac{1}{2}\iint d{\bf r}d{\bf r^{\prime}}\frac{\langle\Psi_{\lambda}^{n}\left|\hat{n}({\bf r})\hat{n}({\bf r^{\prime}})\right|\Psi_{\lambda}^{n}\rangle-n({\bf r})\delta({\bf r}-{\bf r^{\prime}})}{|{\bf r}-{\bf r^{\prime}}|}\end{split}$$

(20)

Here, we have used the expression for the density operator of the $N$-electron systems

$$\hat{n}({\bf r})=\sum_{i=1}^{N}\delta({\bf r}-\hat{{\bf r}}_{i})\,,$$

(21)

and the condition $\langle\Psi_{\lambda}|\hat{n}({\bf r})|\Psi_{\lambda}\rangle=n({\bf r})$. By integrating Eq. 20 over the whole adiabatic connection path, the ground-state energy of the interacting system can be obtained as

$$\begin{split}E=&E_{0}+\int_{0}^{1}\frac{dE_{\lambda}}{d\lambda}d\lambda\\ =&E_{0}+\int d{\bf r}~{}n({\bf r})\left[v^{\mathrm{aux}}_{\lambda=1}({\bf r})-v^{\mathrm{aux}}_{\lambda=0}({\bf r})\right]+\int_{0}^{1}\left\langle\Psi_{\lambda}^{n}\left|\hat{V}_{ee}\right|\Psi_{\lambda}^{n}\right\rangle d\lambda\\ =&T_{s}[\psi]+E_{\mathrm{ext}}[n]+\int_{0}^{1}\left\langle\Psi_{\lambda}^{n}\left|\hat{V}_{ee}\right|\Psi_{\lambda}^{n}\right\rangle d\lambda\\ =&T_{s}[\psi]+E_{\mathrm{ext}}[n]+E_{\mathrm{H}}[n]\\ &+\int_{0}^{1}d\lambda~{}\frac{1}{2}\iint d{\bf r}d{\bf r^{\prime}}~{}\frac{\langle\Psi_{\lambda}^{n}\left|\delta\hat{n}({\bf r})\delta\hat{n}({\bf r^{\prime}})\right|\Psi_{\lambda}^{n}\rangle-n({\bf r})\delta({\bf r}-{\bf r^{\prime}})}{|{\bf r}-{\bf r^{\prime}}|}\,,\end{split}$$

(22)

where we have introduced the fluctuation of the density operator $\delta\hat{n}({\bf r})=\hat{n}({\bf r})-n({\bf r})$, and used the fact that $\langle\Psi_{\lambda}|\delta\hat{n}({\bf r})|\Psi_{\lambda}\rangle=0$ and the definitions of $E_{\mathrm{H}}$, $E_{0}[n]$, and $V_{\mathrm{xc}}$ (Eqs. 6, 15, and 16). Comparing Eq. 3 to Eq. 22 yields the XC energy in terms of the coupling-constant integration

$$E_{\mathrm{xc}}[n]=\int_{0}^{1}W_{\lambda}[n]d\lambda,$$

(23)

with the XC potential along the AC path $W_{\lambda}$ defined as

$$\begin{split}W_{\lambda}[n]=&\left\langle\Psi_{\lambda}^{n}\left|\hat{V}_{ee}\right|\Psi_{\lambda}^{n}\right\rangle-E_{\mathrm{H}}[n]\\ =&\frac{1}{2}\iint d{\bf r}d{\bf r^{\prime}}~{}\frac{\langle\Psi_{\lambda}^{n}\left|\delta\hat{n}({\bf r})\delta\hat{n}({\bf r^{\prime}})\right|\Psi_{\lambda}^{n}\rangle-n({\bf r})\delta({\bf r}-{\bf r^{\prime}})}{|{\bf r}-{\bf r^{\prime}}|}.\\ \end{split}$$

(24)

It is easy to see that the XC potential at the starting point with $\lambda=0$ and $\Psi_{0}^{n}=\Phi_{0}$ is nothing but the exact exchange energy (Eq. 8) in the KS framework,

$$\begin{split}W_{0}[n]=&\left\langle\Phi_{0}\left|\hat{V}_{ee}\right|\Psi_{0}\right\rangle-E_{\mathrm{H}}[n]=E_{\mathrm{x}}^{\mathrm{EX}}[n].\\ \end{split}$$

(25)

If $W_{\lambda}$ were known for the whole AC path, we would be able to approach the exact XC functional via the above coupling-constant integration (Eq. 23). From this perspective, better DFAs can be designed by making a better description of the integrands $W_{\lambda}$ along the AC path.

Let us start by reformulating the Hamiltonian $\hat{H}_{\lambda}$ (Eq. 17) as

$$\hat{H}_{\lambda}=\hat{H}_{0}+\lambda\hat{H}^{\prime},$$

(26)

where $\hat{H}_{0}$ is the Hamiltonian of the KS non-interacting system (Eq. 11) and the perturbation $\hat{H}^{\prime}$ is therefore

$$\begin{split}\hat{H}^{\prime}=&\hat{V}_{ee}+\frac{1}{\lambda}\sum_{i}^{N}[v_{\lambda}^{\mathrm{aux}}({\bf r}_{i})-v_{\mathrm{KS}}({\bf r}_{i})]\\ =&\hat{V}_{ee}-\sum_{i}^{N}\left[v_{\mathrm{H}}({\bf r}_{i})+v_{\mathrm{x}}({\bf r}_{i})+\lambda\frac{\delta E_{c}[n_{1/\lambda}]}{\delta n({\bf r}_{i})}\right].\end{split}$$

(27)

Eq. 27 suggests that the perturbation along the AC path is almost linear, except for the correlation part which evolves nonlinearly with respect to the coupling-constant parameter $\lambda$. As a consequence, the perturbative $\hat{H}^{\prime}$ for $\lambda\rightarrow 0$ is written down as

$$\Delta=\lim_{\lambda\rightarrow 0}\hat{H}^{\prime}=\hat{V}_{ee}-\sum_{i}^{N}[v_{\mathrm{H}}({\bf r}_{i})+v_{x}({\bf r}_{i})].$$

(28)

According to the standard perturbation theory, the ground-state energy of $\hat{H}_{\lambda}$ can be expanded perturbatively based on the ground-state energy of the KS non-interacting system $E_{0}$ (Eq. 15)

$$E_{\lambda}=E_{0}+\lambda E^{(1)}+\lambda^{2}E^{(2)}+\lambda^{3}E^{(3)}+O(\lambda^{4}),$$

(29)

where $E^{(k)}$ refers to the $k$th-order energy correction to $E_{0}$. In particular, the first-order energy correction for $\lambda=0$ is

$$\left.E^{(1)}\right|_{\lambda=0}=\left.\frac{\partial E_{\lambda}}{\partial\lambda}\right|_{\lambda=0}=\left\langle\Psi_{0}^{n}\left|\Delta\right|\Psi_{0}^{n}\right\rangle=E_{x}^{\mathrm{EX}}-E_{\mathrm{H}}-V_{x},\\$$

(30)

and the second-order energy correction for $\lambda=0$ is

$$\begin{split}\left.E^{(2)}\right|_{\lambda=0}=\frac{1}{2}\left.\frac{\partial^{2}E_{\lambda}}{\partial\lambda^{2}}\right|_{\lambda=0}=E_{c}^{\mathrm{GL2}}=E_{c}^{\mathrm{GL-SE}}+E_{c}^{\mathrm{PT2}}\\ \end{split}$$

(31)

which is widely recognized as the Görling-Levy (GL) theory of the coupling-constant perturbation expansion to the second order [17]. The first term of $E_{c}^{\mathrm{GL2}}$ is associated with the single excitation (SE) contribution

$$E_{c}^{\mathrm{GL-SE}}=\sum_{m}^{occ}\sum_{a}^{vir}\sum_{\sigma}\frac{\left|\left\langle\psi_{m\sigma}\left|v_{x}^{\mathrm{HF}}-v_{x}\right|\psi_{a\sigma}\right\rangle\right|^{2}}{\varepsilon_{m\sigma}-\varepsilon_{a\sigma}},$$

(32)

where $v_{x}^{\mathrm{HF}}$ is the Hartree-Fock non-local exchange potential, which is the derivative of $E_{x}^{\mathrm{EX}}$ against the KS orbitals, while $v_{x}$ is the KS local exchange potential defined in Eq. 14. The second part of $E_{c}^{\mathrm{GL2}}$ counts the contributions from double excitations, namely $E_{c}^{\mathrm{PT2}}$,

$$E_{c}^{\mathrm{PT2}}=\frac{1}{4}\sum_{m,n}^{occ}\sum_{a,b}^{vir}\sum_{\sigma,\sigma^{\prime}}\frac{\left|\left\langle\psi_{m\sigma}\psi_{n\sigma^{\prime}}||\psi_{a\sigma}\psi_{b\sigma^{\prime}}\right\rangle\right|^{2}}{\varepsilon_{m\sigma}+\varepsilon_{n\sigma^{\prime}}-\varepsilon_{a\sigma}-\varepsilon_{b\sigma^{\prime}}},$$

(33)

where $\langle\psi_{m\sigma}\psi_{n\sigma^{\prime}}||\psi_{a\sigma}\psi_{b\sigma^{\prime}}\rangle$ is defined in Eq. 10. $E_{c}^{\mathrm{PT2}}$ shares the same formula as the second-order Møller-Plesset perturbation theory (MP2).

As discussed above, the initial integrand along the AC path is given by Eq. 25, which is simply the Hartree-Fock like exact-exchange energy (Eq. 8). Moreover, the second-order derivative against $\lambda$ based on Eq. 20 together with Eq. 31 defines the initial slope of $W_{\lambda}$

$$\begin{split}W_{0}^{\prime}[n]=&\left.\frac{\partial W_{\lambda}[n]}{\partial\lambda}\right|_{\lambda=0}=2E^{\mathrm{GL2}}.\end{split}$$

(34)

If the XC potential $W_{\lambda}$ along the AC path is linear (see Figure 1A), it is precisely determined by $W_{0}$ and $W_{0}^{\prime}$

$$W_{\lambda}^{\mathrm{linearAC}}=W_{0}+W_{0}^{\prime}*\lambda=E_{x}^{\mathrm{EX}}+2E_{c}^{\mathrm{GL2}}*\lambda.$$

(35)

As a result, the corresponding XC functional $E_{xc}^{\mathrm{linearAC}}$ can be obtained by applying Eq. 23,

$$E_{xc}^{\mathrm{linearAC}}=E_{x}^{\mathrm{EX}}+E_{c}^{\mathrm{GL2}}.$$

(36)

Eq. 36 looks similar to the standard MP2 method, which is composed of the HF exchange and the MP2 correlation. However, it has to be emphasized that the MP2 method is the lowest-level correlated wavefunction method for many-electron systems, while Eq. 36 is the exact (very accurate) exchange-correlation functional for any systems with linear (quasi-linear) $\lambda$-dependence. The key difference here is that, the Hartree-Fock orbitals are used for the standard MP2 calculations, while the KS non-interacting orbitals, which deliver the exact density of the interacting system, are employed for the evaluation of HF-like exchange and GL2 correlation in Eq. 36. The linear AC functional in the context of coupling-constant expansion indicates that, for the development of fifth-rung DFAs, the KS unoccupied orbitals can be introduced in the form of second-order perturbation theory. Equally important is the quality of input density and orbitals, which ensures the double of $E_{c}^{\mathrm{GL2}}$ properly represents the initial slope of the XC potential along the AC path.

In this consideration, XYG3-type DHAs were proposed to go beyond the linear model by introducing several empirical parameters to extrapolate the second-order perturbation to infinite order (see Figure 1B). For example, the XYG3-type DHAs utilizing the B3LYP orbitals (xDH@B3LYP)[20] has the form of

$$\begin{split}E_{xc}^{\mathrm{xDH@B3LYP}}=&a_{1}E_{x}^{\mathrm{EX}}+a_{2}E_{x}^{\mathrm{S}}+a_{3}E_{x}^{\mathrm{B88}}+a_{4}E_{c}^{\mathrm{VWN}}\\ &+a_{5}E_{c}^{\mathrm{LYP}}+a_{6}E_{c}^{\mathrm{osPT2}}+a_{7}E_{c}^{\mathrm{ssPT2}}.\end{split}$$

(37)

where $E_{x}^{\mathrm{S}}$ and $E_{x}^{\mathrm{B88}}$ are the exchange contribution from the Slater-type LDA and the Becke88 GGA. $E_{c}^{\mathrm{VWN}}$ is the local Vosko-Wilk-Nusair correlation, and $E_{c}^{\mathrm{LYP}}$ is the Lee-Yang-Parr correlation approximation. Here, $E_{c}^{\mathrm{osPT2}}$ and $E_{c}^{\mathrm{ssPT2}}$ are the opposite-spin and same-spin part of PT2 (Eq. 33), respectively.

$$E_{c}^{\mathrm{osPT2}}+E_{c}^{\mathrm{ssPT2}}=E_{c}^{\mathrm{PT2}}.$$

(38)

It is easy to find that in xDH@B3LYP, the single-excitation contribution ($E_{c}^{\mathrm{SE}}$) is not taken into account. It was argued that $E_{c}^{\mathrm{SE}}$ can be absorbed into the fitting parameters of xDH@B3LYP against the experimental data. Note that, the importance of $E_{c}^{\mathrm{SE}}$ was revisited by Ren et. al. [47] for improving the RPA method in describing vdW-bonded molecules. This will be discussed in Section 2.4.

Equation 37 suggests that there are at most 7 empirical parameters in the general formula of xDH@B3LYP methods. In fact, any xDH@B3LYP methods can be written down using this formula. The first approximation in this family, i.e., XYG3 proposed in 2009 [19], contains only 3 parameters, indicating that four constraints are imposed during the parameter optimization. Similarly, XYGJ-OS proposed in 2011 contains 4 parameters with 3 constraints [22]. Recently, the accuracy limit of xDH@B3LYP methods was explored by fully optimizing all 7 parameters. The resulting XYG7 functional [20] is among the top-class performers in the comprehensive benchmark for main-group chemistry, which can provide a balanced description of both covalent and noncovalent interactions. Its accuracy is comparable to or even better than the very expensive composite methods in wave function theory.

A key deficiency of PT2-based DHAs is in the description of strong correlation systems, where the degeneration of low-lying states could result in extremely large correlation contributions based on the perturbation treatment at any finite orders. In the context of coupling-constant expansion along the AC path, it means that the initial slope $W_{0}^{\prime}[n]=2E_{c}^{\mathrm{GL2}}$ becomes extremely large, such that the initial slope itself is of little relevance to the exact XC energy, which is mainly determined by the XC potential at the fully interacting system $W_{\lambda=1}[n]$ (see Figure 1C).

To address this problem, we should go beyond the finite-order perturbation treatment of the correlation effects. Among different kinds of solutions already proposed, RPA is a promising candidate, as will be introduced in the following section. It offers a systematic way to renormalize different contributions in PT2, and thus address the divergence problem of finite-order perturbation theory as well as DHAs for strong-correlation systems.

Let us return to the coupling-constant integration formula of the XC functional. From the viewpoint of zero-temperature fluctuation-dissipation theorem [93], the density-density correlation (fluctuations) in Eqs. 23 and 24 can be related to the imaginary part of the density response function (dissipation) of the system

$$\langle\Psi_{\lambda}|\delta\hat{n}({\bf r})\delta\hat{n}({\bf r^{\prime}})|\Psi_{\lambda}\rangle=-\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\text{Im}\chi^{\lambda}({\bf r},{\bf r^{\prime}},i\omega)\,.$$

(39)

Here the density response function $\chi^{\lambda}({\bf r},{\bf r^{\prime}},\omega)=\delta n_{\lambda}({\bf r},\omega)/\delta v^{ext}({\bf r^{\prime}},\omega)$ describes the variation of the density of the partially interacting system at the spatial point ${\bf r}$, up to the linear order, due to a change of the local external potential at ${\bf r^{\prime}}$. From Eqs. 23, 24 and 39, we arrive at the renowned ACFDT expression for the XC energy in DFT,

$$\begin{split}E_{\text{xc}}=&\frac{1}{2}\int_{0}^{1}d\lambda\iint d{\bf r}d{\bf r^{\prime}}\frac{1}{|{\bf r}-{\bf r^{\prime}}|}\left[-\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\text{Im}\chi^{\lambda}({\bf r},{\bf r^{\prime}},\omega)-\delta({\bf r}-{\bf r^{\prime}})n({\bf r})\right]\\ =&\frac{1}{2}\int_{0}^{1}d\lambda\iint d{\bf r}d{\bf r^{\prime}}\frac{1}{|{\bf r}-{\bf r^{\prime}}|}\left[-\frac{1}{\pi}\int_{0}^{\infty}d\omega\chi^{\lambda}({\bf r},{\bf r^{\prime}},i\omega)-\delta({\bf r}-{\bf r^{\prime}})n({\bf r})\right]\,.\end{split}$$

(40)

The fact that the above frequency integration can be performed along the imaginary frequency axis originates from the pole structure of $\chi^{\lambda}({\bf r},{\bf r^{\prime}},\omega)$ and the fact that it becomes purely real on the imaginary axis. Such a property simplifies the ground-state energy calculation within the ACFDT framework considerably. The ACFDT expression in Eq. 40 translates the problem of computing the XC energy to the computation of the response functions of a continuous set of fictitious systems along the AC path, which in practice have to be approximated as well.

Obviously, the exact density response function $\chi^{\lambda}({\bf r},{\bf r^{\prime}},i\omega)$ is not known either. However, according to linear-response time-dependent DFT (TDDFT), the interacting response function $\chi^{\lambda}$ for $\lambda>0$ is linked to the noninteracting response function $\chi^{0}$ via the Dyson-like equation

$$\begin{split}\chi^{\lambda}({\bf r},{\bf r^{\prime}},i\omega)=&\chi^{0}({\bf r},{\bf r^{\prime}},i\omega)+\int d{\bf r}_{1}\int d{\bf r}_{2}\chi^{0}({\bf r},{\bf r}_{1},i\omega)\\ &\times\left[\frac{\lambda}{|{\bf r}_{1}-{\bf r}_{2}|}+f_{\mathrm{xc}}^{\lambda}({\bf r}_{1},{\bf r}_{2},i\omega)\right]\chi^{\lambda}({\bf r}_{2},{\bf r^{\prime}},\omega)\,,\end{split}$$

(41)

where $f_{\mathrm{xc}}^{\lambda}$ is the so-called exchange-correlation kernel, which in time domain is given by

$$f_{\mathrm{xc}}^{\lambda}({\bf r}_{1},{\bf r}_{2},t-t^{\prime})=\frac{\delta v_{\mathrm{xc}}^{\lambda}({\bf r}_{1},t)}{\delta n({\bf r}_{2},t^{\prime})}=\frac{\delta^{2}A^{\lambda}_{\mathrm{xc}}[n]}{\delta n({\bf r}_{1},t)\delta n({\bf r}_{2},t^{\prime})}\,.$$

(42)

Here $A^{\lambda}_{\mathrm{xc}}$ is the exchange-correlation part of the action functional [94].

Obviously, the $f_{\mathrm{xc}}$ kernel is a very complex quantity to deal with, and there are considerable ongoing efforts to find better approximations for it [95]. In this context, the random phase approximation amounts to simply neglecting the $f_{\mathrm{xc}}$ kernel, and the resultant interacting response function is termed as the RPA response function,

$$\chi^{\lambda}_{\text{RPA}}({\bf r},{\bf r^{\prime}},i\omega)=\chi^{0}({\bf r},{\bf r^{\prime}},i\omega)+\int d{\bf r}_{1}d{\bf r}_{2}\chi^{0}({\bf r},{\bf r}_{1},i\omega)\frac{\lambda}{|{\bf r}-{\bf r^{\prime}}|}\chi^{\lambda}_{\text{RPA}}({\bf r}_{2},{\bf r^{\prime}},\omega)\,.$$

(43)

In physical terms, the RPA here corresponds to linearized time-dependent Hartree approximation, by which the variation of the exchange-correlation potential due to an external perturbation is neglected.

In Eqs. 41 and 43, $\chi^{0}({\bf r},{\bf r}_{1},i\omega)$ is the independent-particle response function of the KS reference system ($\lambda=0$) and is known explicitly in terms of the single-particle KS orbitals $\psi_{p\sigma}({\bf r})$, orbital energies $\epsilon_{p\sigma}$, and occupation factors $f_{p\sigma}$,

$$\chi^{0}({\bf r},{\bf r^{\prime}},i\omega)=\sum_{p,q,\sigma}\frac{(f_{p\sigma}-f_{q\sigma})\psi_{p\sigma}^{\ast}({\bf r})\psi_{q\sigma}({\bf r})\psi_{q\sigma}^{\ast}({\bf r^{\prime}})\psi_{p\sigma}({\bf r^{\prime}})}{\epsilon_{p\sigma}-\epsilon_{q\sigma}-i\omega}\,.$$

(44)

Based on Eqs. 40 and 43, the XC energy in RPA can be further decomposed into an exact exchange (EX) term and a RPA correlation term,

$$E_{\mathrm{xc}}^{\mathrm{RPA}}[n]=E_{\mathrm{x}}^{\mathrm{EX}}[\psi]+E_{\mathrm{c}}^{\mathrm{RPA}}[\epsilon,\psi],$$

(45)

where

$$\begin{split}E_{\text{x}}^{\text{EX}}[\psi]=&-\frac{1}{2}\sum_{p,q,\sigma}f_{p\sigma}f_{q\sigma}\iint d{\bf r}d{\bf r^{\prime}}\frac{\psi_{p\sigma}^{\ast}({\bf r})\psi_{q\sigma}({\bf r})\psi_{q\sigma}^{\ast}({\bf r^{\prime}})\psi_{p\sigma}({\bf r^{\prime}})}{|{\bf r}-{\bf r^{\prime}}|}\,\end{split}$$

(46)

and

$$\begin{split}E_{\text{c}}^{\text{RPA}}[\epsilon,\psi]=&-\frac{1}{2\pi}\iint d{\bf r}d{\bf r^{\prime}}\frac{1}{|{\bf r}-{\bf r^{\prime}}|}\int_{0}^{\infty}d\omega\left[\int_{0}^{1}d\lambda\chi^{\lambda}_{\text{RPA}}({\bf r},{\bf r^{\prime}},i\omega)-\chi_{0}({\bf r},{\bf r^{\prime}},i\omega)\right]\\ =&\frac{1}{2\pi}\int_{0}^{\infty}d\omega\text{Tr}\left[\text{ln}(1-\chi^{0}(i\omega)v)+\chi^{0}(i\omega)v\right],\end{split}$$

(47)

with $v({\bf r},{\bf r^{\prime}})=1/|{\bf r}-{\bf r^{\prime}}|$. Equation 46 defines the exact-exchange energy (Eq. 8) within the ACFDT context, and is extended here to fractional occupation numbers and to be evaluated with KS orbitals. Furthermore, in the second line of Eq. 47, for brevity the following convention

$$\text{Tr}\left[fg\right]=\iint d{\bf r}d{\bf r^{\prime}}f({\bf r},{\bf r^{\prime}})g({\bf r^{\prime}},{\bf r})$$

(48)

has been adopted.

So far, we have described how the RPA method is defined as an approximate XC energy functional in the KS-DFT context, within the ACFDT framework. As mentioned above, RPA can also be derived from the perspective of the coupled cluster theory and the Green-function based many-body perturbation theory. For instance, to link with the ring coupled cluster doubles (rCCD) theory, the RPA correlation energy can be obtained as,

$$E_{\text{c}}^{\text{RPA}}=\frac{1}{2}\text{Tr}\left(BT^{\text{rCCD}}\right)=\frac{1}{2}\sum_{mn}^{occ}\sum_{ab}^{vir}\sum_{\sigma,\sigma^{\prime}}B_{ma\sigma,nb\sigma^{\prime}}T_{nb\sigma^{\prime},ma\sigma}^{\text{rCCD}}\,.$$

(49)

where $T_{nb\sigma^{\prime},ma\sigma}^{\text{rCCD}}$ is the rCCD amplitude, to be determined by solving the so-called the Riccati equation [96], and $B_{ma\sigma,nb\sigma^{\prime}}=\langle\psi_{m\sigma}\psi_{n\sigma^{\prime}}|\psi_{a\sigma}\psi_{b\sigma^{\prime}}\rangle$. Due to the limited space, here we will not elaborate on the alternative formulations of RPA any further, and interested readers are referred to Ref. [71] and the original references [96, 97] for more details.

Finally, combining Eqs. 3 and 45, one obtains the total ground-state energy of the RPA method,

	$\displaystyle E[n,\psi]$	$\displaystyle=T_{\mathrm{s}}[\psi]+E_{\mathrm{ext}}[n]+E_{\mathrm{H}}[n]+E_{\mathrm{x}}^{\mathrm{EX}}[\psi]+E_{\mathrm{c}}^{\mathrm{RPA}}[\epsilon,\psi]$
		$\displaystyle=\langle\Psi_{\lambda=0}^{n}|\hat{H}|\psi_{\lambda=0}^{n}\rangle+E_{\mathrm{c}}^{\mathrm{RPA}}[\epsilon,\psi]$		(50)

where the first four terms can be grouped together, yielding an energy expression that corresponds to the Hartree-Fock energy evaluated with KS orbitals. This is simply given by the expectation value of the full interacting Hamiltonian within the KS Slater determinant.

In practice, RPA calculations are usually done perturbatively on top of a preceding calculation based on lower-rung DFAs. That is, the single-particle orbitals $\{\psi_{p\sigma}\}$ and orbital energies $\{\epsilon_{p\sigma}\}$ employed in Eq. 50 are generated using lower-rung approximations such as LDA, GGAs, or meta-GGAs. Thus, practical RPA calculations are often denoted as “RPA@DFA”, where “DFA” refers to the actual functional used in the preceding reference-state calculation.

Despite its appealing features, RPA does not go without shortcomings. The conventional wisdom is that RPA describes the long-range correlation very well, whereas it is not adequate for short-range correlations. This issue has been well known for homogeneous electron gas. For real materials, RPA was found to underestimate the cohesive energies of both molecules [40, 98] and solids [51]. This stimulated much research interest and several beyond-RPA schemes have been developed to fix this deficiency. Among these, three approaches are particularly noteworthy: 1) Hybridizing RPA with the lower-rung DFAs in the empirical double hybrid framework; 2) improving RPA by restoring the contribution of the $f_{\mathrm{xc}}$ kernel within the ACFDT formalism; and 3) improving RPA from the perspective of diagrammatic many-body perturbation theory. The first empirical approach gains increasing attention in quantum chemistry. The proposed RPA-based DHAs includes (a) dRPA75 [99] and PWRB95 [29] where the direct RPA correlation is used to replace the PT2 term in the double hybrid formula; (b) SCS-dRPA75 [100] and scsRPA [101] that utilize different kinds of spin-component scaled RPA; and (c) SCS-dRPA75rs [102] with range-separated RPA. Successful beyond-RPA schemes of the second type include the truncated adiabatic LDA kernel correction [103], and the more systematic construction of the XC kernel with respect to a series expansion of the coupling constant $\lambda$ within the ACFDT framework [104, 95].

The diagrammatic many-body expansion approach to correct RPA, on which we shall concentrate here, exploits the fact that RPA, in contrast with other density functional approximations (DFAs), has a clear diagrammatic representation – an infinite summation of the ring diagrams. The question arises if there is a simple and systematic way to incorporate the missing diagrams to arrive at an improved theory. The fermionic nature of electrons requires the many-electron wavefunction to be antisymmetric, and the diagrammatic representation of the correlation energy contains graphs describing both direct processes and exchange processes. For example, in the Møller-Plesset perturbation theory, both types of processes are included at each order, which ensures the theory to be one-electron “self-correlation free” – the correlation energy for an one-electron system being zero. The standard PT2 correlation (Eq. 33) contains one direct term and one exchange term, represented respectively by the leading diagram in the first two rows of Fig. 2. These two terms cancel each other for one-electron systems. RPA is represented by the ring diagrams summed up to infinite order, as represented by the first-row diagrams in Fig. 2, where only “direct” processes are accounted for ³³3The RPA discussed here is often referred to as direct RPA in quantum chemistry literature. In some literature, RPA in fact corresponds to the time-dependent Hartree-Fock theory, where the exchange terms are also included. Nowadays, RPA with exchange included is usually referred to as RPAX or full RPA. A simple way to include the exchange processes and eliminate the self-correlation error is to antisymmetrize the two-electron Coulomb integrals within the rCCD formulation of RPA [cf. Eq. 49]. By doing so, a second-order screened exchange (SOSEX) term [105, 106, 107] is added, and the resultant RPA+SOSEX correlation energy is given by

$$E_{\text{c}}^{\text{RPA+SOSEX}}=\frac{1}{2}\sum_{mn}^{occ}\sum_{ab}^{vir}\sum_{\sigma,\sigma^{\prime}}\left[B_{ma\sigma,nb\sigma^{\prime}}-\delta_{\sigma\sigma^{\prime}}B_{mb\sigma,na\sigma}\right]T^{\text{rCCD}}_{nb\sigma^{\prime},ma\sigma}\,,$$

(51)

with the two terms in Eq. 51 corresponding to RPA and SOSEX correlation energies respectively. The SOSEX contribution can be represented by the diagrams shown in the second row of Fig. 2, where the leading term corresponds to the exchange contribution in PT2. The RPA+SOSEX scheme has the interesting feature that it is one-electron self-correlation free, and improves substantially the total energy [106, 107]. The underbinding problem of RPA for chemically bonded molecules and solids is on average alleviated by the SOSEX correction.

As illustrated in Fig. 2, both RPA and SOSEX correlations can be interpreted as infinite-order summations of selected types of diagrams, with the PT2 terms in the leading order. This perspective is instructive for identifying important contributions that are still missing in the RPA+SOSEX scheme. In fact, at the second order, in addition to the direct and exchange terms, there is yet another type of contribution, arising from the SE, as already discussed above in the context of Görling-Levy perturbation theory (Eqs. 31 and 32). For beyond-RPA corrections, the SE-type correction was initially not derived according to the Görling-Levy perturbation theory in Sec. 2.2, but rather the Rayleigh-Schrödinger perturbation theory (RSPT) [108]. As shown in Ref. [47], within RSPT, the 2nd-order SE contribution is given by

$$\begin{split}E_{\text{c}}^{\text{SE}}=&\sum_{m}^{occ}\sum_{a}^{vir}\sum_{\sigma}\frac{|\langle\Phi_{0}|\hat{H}^{\prime}|\Phi_{m\sigma}^{a\sigma}\rangle|^{2}}{E_{0}-E^{(0)}_{m\sigma,a\sigma}}\\ =&\sum_{m}^{occ}\sum_{a}^{vir}\sum_{\sigma}\frac{|\langle\psi_{m\sigma}|\hat{v}^{\text{HF}}-\hat{v}^{\text{MF}}|\psi_{a\sigma}\rangle|^{2}}{\epsilon_{m\sigma}-\epsilon_{a\sigma}}\\ =&\sum_{m}^{occ}\sum_{a}^{vir}\sum_{\sigma}\frac{|\langle\psi_{m\sigma}|\hat{f}|\psi_{a\sigma}\rangle|^{2}}{\epsilon_{m\sigma}-\epsilon_{a\sigma}}\end{split}$$

where $\hat{v}^{\text{MF}}$ is the KS mean-field (MF) potential, according to which the reference state is generated, and $\hat{v}^{\text{HF}}$ is the non-local Hartree-Fock potential, evaluated using orbitals $\{\psi_{p\sigma}\}$ which themselves are generated according to the potential $\hat{v}^{\text{MF}}$. Furthermore, $\hat{f}=-\nabla^{2}/2+\hat{v}^{\text{ext}}+\hat{v}^{\text{HF}}$ is the single-particle Hartree-Fock Hamiltonian (also known as the Fock matrix in the quantum chemistry literature), and $\Phi_{0}$ and $\Phi_{m\sigma}^{a\sigma}$ are the ground-state Slater determinant formed by the lowest-occupied $\psi_{p\sigma}$’s and the singly-excited configuration generated by exciting one electron from an occupied state $m$ (with spin $\sigma$ and energy $\epsilon_{m\sigma}$) to a virtual state $a$ (with the same spin and energy $\epsilon_{a\sigma}$), respectively. A detailed derivation of Eq. (2.4) can be found in the supplementary material of Ref. [47]. Obviously, for a Hartree-Fock reference where $\hat{v}^{\text{MF}}=\hat{v}^{\text{HF}}$, Eq. 2.4 becomes zero, a fact known as Brillouin theorem [108]. Therefore, this term is not present in standard MP2 theory based on the Hartree-Fock reference.

One may notice that the SE expression derived within RSPT is different from the GL-SE given by Eq. 32, with the exact-exchange optimized effective potential (OEP) in the GL-SE is replaced in Eq. 2.4 by the mean-field potential, on which the reference state is based. In fact, RSPT can also be interpreted from the ACFDT perspective along a linear adiabatic connection path that linearly interpolates the non-interacting reference Hamiltonian and the full interacting Hamiltonian. However, in this case, the electron density is not (and cannot be) fixed any longer along the adiabatic connection path. Such a linear connection path was considered long ago by Harris and Jones in Ref. [109], and further discussed in Refs. [110, 111, 112]. As pointed by Klimeš et al. in Ref. [112], the SE term given by Eq. 2.4 accounts for the change of the mean-field Hartree and exchange energy as $\lambda$ goes from 0 to 1, stemming from the change of both the electron density and density matrix. In contrast, the GL-SE only accounts for the change of exchange energy, since the electron density is fixed, and so does the Hartree energy.

From the practical point of view, Eq. 2.4 is easier to implement, since one can use the readily available LDA- or GGA-type KS potential as $v^{\mathrm{MF}}$, without solving a rather cumbersome OEP equation at the exact change level. Yet, as discussed above, more physical effects missing in the standard RPA are accounted for by Eq. 2.4. Due to these advantages, the SE expression derived within RSPT is de facto the SE correction term for beyond-RPA schemes.

In Ref [47], it was shown that adding the SE term of Eq. 2.4 to RPA significantly improves the accuracy of vdW-bonded molecules, which the standard RPA scheme generally underbinds, and the SOSEX correction does not appreciably improve the situation. Similar to the RPA and SOSEX cases, one can identify a sequence of single-excitation processes up to infinite order, as illustrated in the third row in Fig. (2). Summing these single-excitation diagrams to infinite order represents a renormalization of the 2nd-order SE, and is therefore termed as renormalized single excitations (rSE) [71, 68]. Remarkably, this infinite summation can be translated into a closed expression similar to Eq. (2.4),

$$E_{\text{c}}^{\text{rSE}}=\sum_{m}^{occ}\sum_{a}^{vir}\sum_{\sigma}\frac{|\tilde{f}_{ma}^{\sigma}|^{2}}{\tilde{\epsilon}_{m\sigma}-\tilde{\epsilon}_{a\sigma}}\,,$$

(52)

where $\tilde{\epsilon}_{m\sigma}$ and $\tilde{\epsilon}_{a\sigma}$ are the eigenvalues obtained by diagonalizing separately the occupied-occupied and virtual-virtual subblocks of the Fock matrix,

	$\displaystyle\sum_{l}f_{ml}^{\sigma}{\cal O}_{ln}^{\sigma}$	$\displaystyle=\tilde{\epsilon}_{m\sigma}{\cal O}_{mn}^{\sigma}$		(53)
	$\displaystyle\sum_{c}f_{ac}^{\sigma}{\cal U}_{cb}^{\sigma}$	$\displaystyle=\tilde{\epsilon}_{a\sigma}{\cal U}_{ab}^{\sigma}\,,$

with $f_{pq}^{\sigma}=\langle\psi_{p\sigma}|\hat{f}|\psi_{q\sigma}\rangle$. Note that $f^{\sigma}_{pq}$ matrix is not diagonal since $\psi_{p\sigma}$’s are the reference KS orbitals and not the eigenfunctions of $\hat{f}$ – the single-particle Hartree-Fock Hamiltonian. Now, $\tilde{f}_{ma}^{\sigma}$ in the numerator of Eq. (52) are the “transformed" off-diagonal block of the Fock matrix

$$\tilde{f}_{ma}^{\sigma}=\sum_{n}^{occ}\sum_{b}^{vir}{\cal O^{\ast}}_{mn}^{\sigma}f_{nb}^{\sigma}{\cal U}_{ba}^{\sigma}\,,$$

(54)

where the eigenvectors obtained in (53) are used here as the transformation coefficients. The physical origin of the SE corrections is that the commonly used KS references for RPA and beyond-RPA calculations are not the optimal starting point. The SE corrections accounts for the “orbital relaxation" effect, which leads to a lowering of the ground-state total energy. Furthermore, it was found that SE and rSE corrections generally increases the binding energies, in particularly for weakly bonded systems. When the SE contributions are added, the electron density effectively contracts compared to that yielded by local and semilocal DFAs, and this reduces the Pauli repulsion between closed-shell atoms and molecules. Consequently, the underbinding behavior of the standard RPA is corrected for both weakly bonded molecules [47, 68] and solids [112]. Recently, the concept of rSE has been extended to Green’s function by Yang and coworkers [113, 114]. These authors showed that, by utilizing the so-called renormalized singles Green’s function, the starting point dependence in the usual perturbative $G_{0}W_{0}$-type calculations are significantly reduced [113, 114].