# Coupling Natural Orbital Functional Theory and Many-Body Perturbation Theory by using non-dynamically correlated canonical orbitals

Mauricio Rodriguez, Ion Mitxelena, Fabien Bruneval, Mario Piris

► To cite this version:

Mauricio Rodriguez, Ion Mitxelena, Fabien Bruneval, Mario Piris. Coupling Natural Orbital Functional Theory and Many-Body Perturbation Theory by using non-dynamically correlated canonical orbitals. Journal of Chemical Theory and Computation, 2021, 17 (12), pp.7562. 10.1021/acs.jctc.1c00858 . cea-03516535

HAL Id: cea-03516535

https://cea.hal.science/cea-03516535

Submitted on 7 Jan 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Coupling Natural Orbital Functional Theory and
Many-Body Perturbation Theory by using
non-dynamically correlated canonical orbitals

M. Rodríguez-Mayorga,*,†,‡ I. Mitxelena,¶
F. Bruneval,† and M. Piris¶,§

†DEN, Service de Recherches de Métallurgie Physique, CEA, Université Paris-Saclay,
F-91191 Gif-sur-Yvette, France

‡Department of Theoretical Chemistry, VU University Amsterdam, De Boelelaan 1083,
1081 HV Amsterdam, The Netherlands.

¶Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU) and Donostia
International Physics Center (DIPC), 20018 Donostia, Euskadi, Spain.

§IKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Euskadi, Spain.

E-mail: marm3.14@gmail.com

**Abstract**

We develop a new family of electronic structure
methods for capturing at the same time the
dynamic and non-dynamic correlation effects.
We combine natural orbital functional theory
(NOFT) and many-body perturbation theory
(MBPT) through a canonicalization procedure
applied to the natural orbitals to gain access to
any MBPT approximation. We study three dif-
ferent scenarios: Corrections based on second-
order Møller-Plesset (MP2), Random-Phase
Approximation (RPA), and coupled-cluster sin-
gles doubles (CCSD). Several chemical prob-
lems involving different types of electron cor-
relation in singlet and multiplet spin states
have been considered. Our numerical tests re-
veal that RPA-based and CCSD-based correc-
tions provide similar relative errors in molecu-
lar dissociation energies (D<sub>e</sub>) to the results ob-
tained using a MP2 correction. With respect
to the MP2 case, the CCSD-based correction
improves the prediction, while the RPA-based
correction reduces the computational cost.

**1 Introduction**

Many body perturbation theory¹ (MBPT) is
one of the preferred methodologies among the
chemists and physicists to account for electron
correlation energies (E<sub>c</sub>). While physicists usu-
ally treat MBPT by building Green's functions
and applying the GW approximation as pro-
posed by Hedin,² chemists prefer the *n*-order
Møller-Plesset³ (MP*n*) perturbation theory or
the coupled-cluster (CC) approaches.<sup>4–6</sup> These
methods lead to different correlation energy ex-
pressions, but have in common that they are
applied as a one-shot correction on top of the
result obtained with a single determinant (SD)
wavefunction from Hartree-Fock (HF) or den-
sity functional calculations.

Correlation energies from the GW approx-
imation<sup>7,8</sup> are generally evaluated using the
Galitskii-Migdal<sup>9</sup> (GM) or the Klein func-
tional.<sup>10</sup> The latter produces the Random-
Phase Approximation (RPA),<sup>10–14</sup> whose func-
tional expression can be obtained from the adi-
abatic connection formalism, whereas the GM
is an explicit functional of the Green’s func-
tion where the so-called self-energy can be ap-
proximated by the GW product. In chemistry, it is often sufficient to obtain $E_c$ using MP2 or the CC approximation including only single and double excitations (CCSD).$^{15–18}$ Lately, the physical and chemical approaches of MBPT are spreading from one community to other;$^{19–25}$ improving the applicability of MBPT and building bridges between the two perspectives.$^{26,27}$

A shortcoming of MBPT-based methods is the poor description of the so-called non-dynamic electron correlation. The latter cannot always be overlooked and is present in many physical and chemical processes, such as bond breaking$^{28,29}$ or Mott-Hubbard insulators.$^{30–32}$ To account for the non-dynamic electron correlation, any method must be capable of describing the multi-configurational character of the corresponding physical state. Although multi-configurational methods, such as complete active space self-consistent field (CASSCF), are suitable for retrieving non-dynamic electron correlation energies,$^{33}$ they cannot be considered as black-box methods. Indeed, a proper definition of the active space is crucial to obtain a good description of any system.$^{34,35}$

An alternative approach is the one-particle reduced density matrix (1RDM) functional theory. In 1975, Gilbert laid the groundwork by extending the Hohenberg-Kohn theorem to non-local potentials.$^{36}$ A few years later, Valone proved$^{37}$ the existence of the functional for ensembles using the Levy’s functional$^{38}$ in the domain of all ensemble N-representable 1RDMs$^{39}$. Actually, we must only reconstruct the electron-electron potential energy $V_{ee}$ in terms of the 1RDM since the non-interacting part of the electronic Hamiltonian is a one-particle operator. $V_{ee}$ is defined independently of the external potential under consideration and is therefore a universal functional.

The 1RDM can be expressed in an orthonormal basis set {$\phi_p$} of dimension $M$ as

$$ \gamma(\mathbf{r}, \mathbf{r}') = \sum_{pq}^{M} D_p^q \phi_p^*(\mathbf{r}') \phi_q(\mathbf{r}), \quad (1) $$

whose diagonal part corresponds to the elec-

tronic density: $\rho(\mathbf{r}) = \gamma(\mathbf{r}, \mathbf{r})$. The diagonalization of the 1RDM $^1D$ leads to the so-called natural orbital (NO) representation

$$ \gamma(\mathbf{r}, \mathbf{r}') = \sum_{p}^{M} n_p \chi_p^*(\mathbf{r}') \chi_p(\mathbf{r}) \quad (2) $$

where {$\chi_p$} are known as the natural orbitals (NOs) and {$n_p$} are the corresponding occupation numbers (ONs). The 1RDM functional theory is referred to as Natural Orbital Functional (NOF) theory (NOFT)$^{40,41}$ when the 1RDM is expressed in its NO representation, hence $V_{ee}$ [$^1D$] becomes $V_{ee}$ [$\{\chi_p, n_p\}$].

At the moment, the explicit reconstruction $V_{ee}$ [$\{\chi_p, n_p\}$] has resulted in an unattainable goal, and we can only make approximations. It is well known that the ground-state energy of a system with a Hamiltonian involving not more than two-body interactions can be cast as an exact functional $E$ [$^2D$] of the two-particle reduced density matrix (2RDM) $^2D$.42 Thus, the typical approach is to employ a reconstruction functional $^2D$ [$\{n_p\}$]. Such NOF approximations$^{43,44}$ have proven to be efficient to account for the multi-configurational character of electronic systems. An important advantage is that they do not require the definition of an active space, which facilitates their applicability.

Nevertheless, the functional N-representability problem arises;45 thus, the reconstructed $^2D$ must be N-representable.$^{46}$ Otherwise, an approximate NOF can lead to non-physical energy values. So far, only PNOFs$^{47–50}$ are based on the reconstruction of $^2D$ subject to necessary N-representability conditions. These functionals are capable of achieving chemical accuracy in many cases,$^{51,52}$ however, they suffer from an important lack of dynamic correlation. To recover this correlation, orbital-invariant second-order perturbative corrections have been implemented$^{53,54}$ with significant results.$^{55,56}$ Note that an orbital-invariant formulation is used because the orbitals that serve as reference are the NOs resulting from an approximate NOF which are usually localized.

In this work, a new path is proposed based on a modification of the standard MBPT formulated for canonical HF orbitals. A canonicaliza-
tion procedure is applied to the NOs in order to produce suitable orbitals that can be used in any MBPT approximation to compute the dynamic energy correction. Our canonical orbitals should not be confused with the canonical representation obtained<sup>57</sup> by diagonalization of the matrix of Lagrange multipliers built from the optimized NOs. Our resulting orbitals are similar to HF ones. They diagonalize a generalized Fock matrix that includes electronic correlation but does not contain interaction terms in the Fockian between occupied and virtual orbitals, as occurs in the HF approximation. This is crucial because it allows us to connect our proposal with the NOF-oiMP2 method, just as the orbital-invariant formulation of MP2<sup>58</sup> connects with the standard MP2 theory.

We pursue two objectives: On the one hand, set the theoretical foundations for a new family of NOF-MBPT methods able to capture at the same time the non-dynamic and dynamic electron correlation and, on the other hand, test the proposed NOF-MBPT approximation by using different MBPT corrections.

This article is organized as follows. In section 2, we briefly review NOFT and explain the canonicalization procedure. Next, we develop energy corrections for the special cases of describing dynamical correlation effects with MP2, RPA or CCSD approximations. In section 3, the performance of these approaches is tested in detail with six different systems, namely H<sub>2</sub>, (CH<sub>4</sub>)<sub>2</sub>, CO, NH, F<sub>2</sub> and N<sub>2</sub>, as well as with a set of small molecules, in comparison with state-of-the-art electronic structure methods.

## 2 Theory

The success of NOF-MBPT methods is determined by the NOs used to generate the reference. The functional PNOF7s proved<sup>53</sup> to be the functional of choice for the NOF-oiMP2 method. The “s” emphasizes that PNOF7s takes into account only the static correlation between pairs, and therefore avoids double counting in the regions where the dynamic correlation predominates, already in the NOF optimization.

Along this work we employ PNOF7s to obtain the NOs and ONs.

### 2.1 PNOF7s

In NOFT, the total energy is approximated in terms of the ONs and NOs as

$$E = 2 \sum_{p}^{M} n_p \mathcal{H}_{pp} + \sum_{pqrs}^{M} {}^{2}D[n_p, n_q, n_r, n_s] \langle pq | rs \rangle \quad (3)$$

where

$$\mathcal{H}_{pp} = \int d\mathbf{r} \chi_p^*(\mathbf{r}) \left( -\frac{\nabla_r^2}{2} + v_{\text{ext}}(\mathbf{r}) \right) \chi_p(\mathbf{r}) \quad (4) \\ \langle pq | rs \rangle = \iint d\mathbf{r}_1 d\mathbf{r}_2 \frac{\chi_p^*(\mathbf{r}_1)\chi_q^*(\mathbf{r}_2)\chi_r(\mathbf{r}_1)\chi_s(\mathbf{r}_2)}{|\mathbf{r}_2 - \mathbf{r}_1|} \quad (5)$$

In Eq.(4), $v_{\text{ext}}(\mathbf{r})$ is the external potential defined by the geometry within the Born-Oppenheimer approximation and in the absence of any external field. The ONs are constrained to lie in the interval $0 \le n_p \le 1$ to have an ensemble N-representable 1RDM.<sup>39</sup>

For a given spin $S$, there are $(2S+1)$ energy degenerate pure states $|SM\rangle$, so a quantum mixed state formed with them is defined by the following $N$-particle density matrix statistical operator

$$\hat{\mathcal{D}} = (2S + 1)^{-1} \sum_{M=-S}^{S} |SM\rangle \langle SM| \quad (6)$$

According to the recent extension of the NOFT to multiplets,<sup>54</sup> and considering that the expected value of the projection of the total spin ($\hat{S}_z$) for the whole ensemble is zero,

$$\langle \hat{S}_z \rangle = (2S + 1)^{-1} \sum_{M=-S}^{S} M = 0 \quad (7)$$

we can adopt a spin-restricted formalism in which a single set of orbitals is used for $\alpha$ and $\beta$ spins. All spatial orbitals will be then double occupied, so that ONs for particles with $\alpha$ and $\beta$ spins are equal: $n_p^\alpha = n_p^\beta = n_p$.

Let us consider that $U$ single electrons determine the spin of the system, $S = U/2$, and the
rest of electrons, $L = N - U$, are spin-paired providing zero spin. Accordingly, we are going to divide the orbital space $\Omega$ into two subspaces: $\Omega = \Omega_U \oplus \Omega_L$. Both $\Omega_U$ and $\Omega_L$ are composed of $U$ and $L/2$ mutually disjoint subspaces $\Omega_g$, respectively. Each subspace $\Omega_g \in \Omega_U$ contains only one orbital $g$ with $n_g = 1/2$ which is individually occupied, but we do not know whether the electron has $\alpha$ or $\beta$ spin. In contrast, each $\Omega_g \in \Omega_L$ contains one orbital with $g \le L/2$, and $N_g$ orbitals with $p > N_\Omega = L/2 + U$. In Fig. 1, an illustrative example is shown.

Taking into account the spin, the total occupancy for a given subspace $\Omega_g \in \Omega_L$ is 2, which is reflected in additional sum rule, namely,

$$ \sum_{p \in \Omega_g} n_p = 1, \quad \Omega_g \in \Omega_L \qquad (8) $$

It follows that

$$ 2 \sum_{p \in \Omega_L} n_p = 2 \sum_{g=1}^{L/2} \sum_{p \in \Omega_g} n_p = L \qquad (9) $$

Consequently, the trace of the 1RDM is verified equal to the number of electrons:

$$ 2 \sum_{p \in \Omega} n_p = L + U = N \qquad (10) $$

The reconstruction of the 2RDM in terms of ONs implies intra- and inter-subspace contributions subject to some necessary $N$-representability conditions of the 2RDM, also known as (2,2)-positivity conditions.⁴⁶ The explicit expressions for nonzero 2RDM blocks can be found in Ref. 54. Substituting these expressions in Eq. (3) we obtain PNOF7s, namely

$$ E^{\text{PNOF7s}} = \sum_{g=1}^{L/2} E_g + \sum_{g=L/2+1}^{L/2+U} \mathcal{H}_{gg} + \sum_{f,g=1;f \neq g}^{L/2+U} E_{fg} \qquad (11) $$

$$ E_g = 2 \sum_{p \in \Omega_g} n_p \mathcal{H}_{pp} + \sum_{p,q \in \Omega_g} \Pi_{pq}^{\text{tra}} \mathcal{K}_{pq} \qquad (12) $$

$$ \Pi_{pq}^{\text{tra}} = \begin{cases} -\sqrt{n_p n_q}, & p \text{ or } q \le L/2 \\ +\sqrt{n_p n_q}, & p = q \text{ or } p, q > L/2 \end{cases} \qquad (13) $$

Figure 1: Splitting of the orbital space $\Omega$ into subspaces. In this example, $S = 3/2$ (quadruplet) and $U = 3$, so three orbitals make up the subspace $\Omega_U$, whereas ten electrons ($L = 10$) distributed in five subspaces $\{\Omega_1, \Omega_2, \dots, \Omega_5\}$ make up the subspace $\Omega_L = \Omega_L^a + \Omega_L^b$. Note that $N_g = 2$ and $N_\Omega = L/2 + U = 8$. The arrows depict the values of the ensemble occupation numbers, alpha (↓) or beta (↑), in each orbital.

<figure><img src="image_2.png" /></figure>

$$ E_{fg} = \sum_{p \in \Omega_f} \sum_{q \in \Omega_g} [n_p n_q (2\mathcal{J}_{pq} - \mathcal{K}_{pq}) - 4n_p h_p n_q h_q \mathcal{K}_{pq}] \qquad (14) $$

where $h_p = 1 - n_p$ is the hole in orbital $g$. $\mathcal{J}_{pq} = \langle pq | pq \rangle$ and $\mathcal{K}_{pq} = \langle pq | qp \rangle$ are the Coulomb and exchange integrals. To simplify the expressions, the spatial orbitals were assumed to be real. The first term in Eq. (11) draws the system as independent $L/2$ electron pairs. It is worth noting that this term recovers the Löwdin-Shull functional⁵⁹ with fixed phases that is almost exact for any two-electron system ($L = 2, U = 0$).⁶⁰

The solution is established by optimizing the energy (11) with respect to the ONs and NOs, separately. The conjugate gradient method is used for performing the optimization of the energy with respect to auxiliary variables that
enforce automatically the N-representability bounds of the 1RDM. A self-consistent procedure⁶¹ yields the NOs by an iterative diagonalization procedure in which orbitals are not constrained to remain fixed along the orbital optimization process.

## 2.2 NOF-MBPT

An important shortcoming of approach (11) - (14) is the absence of inter-subspace dynamic electron correlation. Indeed, the last term of Eq. (14) has significant values only when the ONs differ substantially from 1 and 0. Consequently, PNOF7s is able to recover the complete intra-subspace, but only the static inter-subspace correlation. To add the missing dynamic correlation, Piris proposed⁴⁹ the NOF-MP2 method, where the total energy is approximated as

$$E^{\text{NOF-MP2}} = E^{\text{SD}} + E^{\text{dyn}} + E^{\text{nd}} \quad (15)$$

The energy (15) is computed employing the NOs and ONs that minimize a given NOF. $E^{\text{SD}}$ is a HF-like energy obtained from the Slater determinant formed by the strongly-occupied NOs, the dynamic energy ($E^{\text{dyn}}$) is derived from the modified MP2 perturbation correction, while the non-dynamic energy ($E^{\text{nd}}$) comes from the static component of the functional. An important feature of the method is that double counting is avoided by taking the amount of dynamic and static correlation in each orbital as a function of its occupancy. Consequently, attenuation coefficients are included into $E^{\text{dyn}}$ and $E^{\text{nd}}$ to eliminate the double counting in the correlation energy contributions.

Despite the improvements achieved by NOF-MP2 using directly the NOs, especially for predicting dissociation energies, it was found that this approach fails to retrieve sufficient dynamic correlation energy in noble-gas dimmers or when the number of electrons increases in poly-atomic systems. One way to overcome this issue was to formulate the dynamic correction $E^{\text{dyn}}$ using the orbital-invariant MP2 energy⁵⁸ which led to the NOF-oiMP2 method.⁵³

In this work, we return to the original idea of using the standard dynamic correction with canonical orbitals, and generalize the Eq. (15) to define the dynamic energy by any MBPT correction. Consequently, the total energy of all $E^{\text{NOF-c-X}}$ methods developed below will also consist of three contributions, namely

$$E^{\text{NOF-c-X}} = E_{\text{can}}^{\text{SD}} + E_{\text{can}}^{\text{nd}} + E_{\text{can}}^{\text{dyn}}(X) \quad (16)$$

where $X = \text{MP2}$, RPA, CCSD, etc. The “can” label is introduced to highlight that the new reference orbitals are not the NOs but those that will result from a canonicalization process, unlike those used in the NOF-MP2 method.

### 2.2.1 Orbital canonicalization

We now proceed with the canonicalization procedure in order to obtain a diagonal effective Fock matrix that takes into account the effects of electron correlation and with which we can use the standard MBPT.

Following Ref. 49, we first define the intra-subspace and inter-subspace attenuation coefficients $C_p^{\text{tra}}$ and $C_p^{\text{ter}}$, respectively, as follows

$$C_p^{\text{tra}} = \begin{cases} 1 - 4h_p^2, & p \le N_{\Omega} \\ 1 - 4n_p^2, & p > N_{\Omega} \end{cases} \quad (17)$$

$$C_p^{\text{ter}} = \begin{cases} 1, & p \le N_{\Omega} \\ 1 - 4h_p n_p, & p > N_{\Omega} \end{cases} \quad (18)$$

where $N_{\Omega} = L/2 + U$. These coefficients represent the amount of dynamic electronic correlation in each $p$ orbital as a function of its occupancy. Indeed, $C_p^{\text{tra/ter}}$ goes from zero for half occupied orbitals to one if the orbital is empty or fully occupied. The exception is $C_p^{\text{ter}}$ for an orbital below $N_{\Omega}$ for which this dependence is not considered.

Next, we keep the diagonal elements of the Fock matrix unmodified, that is $\tilde{\mathcal{F}}_{pp} = \mathcal{F}_{pp}$, and use the attenuation coefficients as the case may be (intra-subspace and inter-subspace), to define the off-diagonal attenuated-Fock matrix elements and attenuated-two-electron repulsion integrals as
$$ \tilde{\mathcal{F}}_{pq} = \mathcal{F}_{pq} C_p^{\text{tra/ter}} C_q^{\text{tra/ter}} \quad (19) $$

$$ \langle pq|rs \rangle = \langle pq|rs \rangle C_p^{\text{tra/ter}} C_q^{\text{tra/ter}} C_r^{\text{tra/ter}} C_s^{\text{tra/ter}} \quad (20) $$

Recall that in the standard MBPT the canonical HF orbitals that diagonalize the Fock matrix are used to build the reference Slater determinant. On the other hand, the orbital invariant MBPT uses non-canonical orbitals that are related to the canonical ones through a unitary transformation without mixing the occupied and virtual spaces, which ensures that the reference energy remains fixed. Taking this into account, the canonicalization of NOs will therefore consist of the following steps:

1. Build $\tilde{\mathcal{F}}_{pq}$ in NO representation, and set all virtual-occupied ($\tilde{\mathcal{F}}_{ia}$) and occupied-virtual ($\tilde{\mathcal{F}}_{ai}$) elements equal to zero.

2. Diagonalize the resulting matrix to produce the non-dynamically correlated canonicalized orbitals $\tilde{\phi}_p$ with their respective energies $\tilde{e}_p$.

3. Arrange the canonicalized orbitals $\{\tilde{\phi}\}$ in ascending order of energies. Note that the $N_\Omega$ orbital is considered the last occupied orbital.

4. Transform $\langle pq|rs \rangle$ and $\langle pq|rs \rangle$ from NO basis $\{\chi\}$ to the canonicalized basis $\{\tilde{\phi}\}$.

Note that the $i, j$ indices refer to the strongly occupied NOs and $a, b$ to weakly occupied ones. Steps 1 to 4 apply only once, that is, the procedure is not self-consistent. It is called canonicalization since it produces a basis $\{\tilde{\phi}\}$ with a diagonal effective Fock matrix. Importantly, the new canonical orbitals are close to those of HF when dynamic correlation predominates, but differ substantially from these when there is non-dynamic correlation.

### 2.2.2 $E_{\text{can}}^{\text{dyn}}(X)$ contribution

To achieve the MBPT corrections, we consider the Slater determinant $\Psi^{\text{SD}}(\tilde{\phi})$ as the zeroth-order wavefunction $\Psi^{(0)}$, and the zeroth-order Hamiltonian $\hat{H}^{(0)} = \sum_i \tilde{e}_i |\tilde{\phi}_i\rangle\langle\tilde{\phi}_i|$ where $\tilde{e}_i$ is

the $i$th diagonal element of the effective Fock matrix. Accordingly, the first-order energy perturbation correction affords the reference energy in Eq. (16) as

$$ E_{\text{can}}^{\text{SD}} = \langle \Psi^{\text{SD}}(\tilde{\phi}) | \hat{H} | \Psi^{\text{SD}}(\tilde{\phi}) \rangle \quad (21) $$

It is worth noting that $E_{\text{can}}^{\text{SD}}$ differs from the HF energy since the $\tilde{\phi}$ orbitals used to construct $\Psi^{\text{SD}}$ are those obtained from the canonicalization rather than canonical HF orbitals.

Let us start with $X = \text{MP2}$. $E_{\text{can}}^{\text{dyn}}(\text{MP2})$ is obtained from the second-order correction $E^{(2)}$ of the MP2 method.⁶² The first-order wavefunction is a linear combination of all doubly excited configurations from $\Psi^{\text{SD}}(\tilde{\phi})$. The dynamic energy correction takes the form

$$ E_{\text{can}}^{\text{dyn}} = \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{virt}} \mathcal{A}_i \mathcal{A}_j \frac{2\langle ij|ab \rangle - \langle ij|ba \rangle}{\tilde{e}_i + \tilde{e}_j - \tilde{e}_a - \tilde{e}_b} \langle ij|ab \rangle \quad (22) $$

where the coefficients $\mathcal{A}_i$ were recently introduced⁵⁴ to consider one electron with $\alpha$ or $\beta$ spin in the subspaces $\Omega_U$ present in spin multiplets as

$$ \mathcal{A}_i = \begin{cases} 1, & i \le L/2 \\ 1/2, & i \in [L/2 + 1 : L/2 + U] \end{cases} \quad (23) $$

Notice that the NOF-c-MP2 represents an advantage over the NOF-oiMP2 since a single diagonalization of the attenuated Fock matrix is required instead of the solution of the residual equations.⁵³ The correlation energy recovered by NOF-c-MP2 is equivalent to the one obtained with NOF-oiMP2 because these two formulations are related by the unitary transformation mentioned above. In any case, both methods extend the MP2 method to systems that exhibit a multi-configurational character.

The attenuated two-electron integrals and orbital energies also provide access to other dynamic energy corrections given by the MBPT. Next, we explore the link between NOF-c and RPA, denoted as NOF-c-RPA, as well as with the CCSD, hereafter referred to as NOF-c-CCSD.

To introduce NOF-c-RPA, we must apply the
RPA approximation, i.e. the time-dependent Hartree approximation to build the screened Coulomb interaction (W). This quantity can be readily obtained by solving the Casida-like equations.⁶³,⁶⁴ However, since in our approach attenuated orbital energies and two-electron integrals are inserted into Casida-like equations, the RPA correlation energy is modified too. Indeed, an attenuated version of the Casida-like equations produces an attenuated-screened Coulomb interaction ($\tilde{W}$), which is switched off when non-dynamic electron correlation effects are dominant. More details concerning the role of the Ctra/ter coefficients in the Casida-like equations can be found in the supporting information.

Defining the attenuated-averaged-screened Coulomb interaction matrix elements in the particle-hole basis and using the imaginary frequencies we arrive to $\tilde{W}_{ia}^{jb}(i\omega)$.⁶⁵ Then, following Ref. 14, the resulting RPA correlation energy reads as

$$E_{\text{can}}^{\text{dyn}} = -\frac{2}{\pi} \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{virt}} \mathcal{A}_i \mathcal{A}_j \langle ij | ab \rangle \times \int_0^{\infty} d\omega \tilde{W}_{ia}^{jb}(i\omega) \tilde{P}_{ia}(i\omega) \tilde{P}_{jb}(i\omega) \quad (24)$$

where

$$\tilde{P}_{ia}(i\omega) = \frac{2(\tilde{e}_i - \tilde{e}_a)}{(\tilde{e}_i - \tilde{e}_a)^2 + \omega^2} \quad (25)$$

Replacing the MP2 correlation energy by the CCSD one leads to the NOF-c-CCSD approximation. In fact, we insert the attenuated two-electron integrals ($\langle pq|rs \rangle$) and attenuated orbital energies ($\tilde{e}_p$) into the CCSD amplitudes equations.¹⁵,¹⁸ The resulting dynamic contribution reads as

$$E_{\text{can}}^{\text{dyn}} = \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{virt}} \mathcal{A}_i \mathcal{A}_j \langle ij | ab \rangle, \\ \times (2[\tilde{t}_{ij}^{ab} + \tilde{t}_i^a \tilde{t}_j^b] - [\tilde{t}_{ji}^{ab} + \tilde{t}_j^a \tilde{t}_i^b]) \quad (26)$$

where $\tilde{t}_i^a$ and $\tilde{t}_{ij}^{ab}$ are the attenuated spinless $t_1$ and $t_2$ amplitudes.¹⁵ The meaning of the $\mathcal{A}_i$ coefficients is the same as in NOF-c-MP2 and NOF-c-RPA. Recall that a restricted formalism

is employed, therefore, the spin dependence has already been accounted.

It is worth to mention that we have also explored other MBPT flavors to define NOF-c-X methods, such as GW evaluated as a one-shot correction ($G_0W_0$) in the Galitskii-Migdal equation correlation energy.⁹ The dynamic correlation energy contribution of this method (NOF-c-$G_0W_0$@GM) can be written as

$$E_{\text{can}}^{\text{dyn}} = 2 \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{virt}} \sum_{s}^{N_{ab}} \mathcal{A}_i \mathcal{A}_j \frac{\tilde{w}_{ia} (\tilde{X}_{jb}^s + \tilde{Y}_{jb}^s)}{\tilde{e}_i - \tilde{e}_a - \Omega_s} \langle ij | ab \rangle \quad (27)$$

where the sum over $s$ runs over all possible excitations ($Nab = \text{occ} \times \text{virt}$). {$\tilde{w}_{ia}$} are the attenuated residues, {$\Omega_s$} refer to the excitation energies, and {$\tilde{X}_{jb}^s + \tilde{Y}_{jb}^s$} are the coefficients of the eigenvectors associated to excitations (see the Eqs. 3 and 6 in the supporting information for their definition).

We have also analyzed approximations including second-order exchange corrections.¹⁴,⁶⁶⁻⁶⁹ These are usually referred to as SOSEX corrections applied to RPA and $G_0W_0$, which can be introduced by antisymmetrizing the Coulomb integrals.⁶⁶,⁷⁰ The dynamic correlation energy of the NOF-c-RPA+SOSEX method reads as

$$E_{\text{can}}^{\text{dyn}} = -\frac{1}{\pi} \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{virt}} \mathcal{A}_i \mathcal{A}_j [2\langle ij | ab \rangle - \langle ij | ba \rangle] \\ \times \int_0^{\infty} d\omega \tilde{W}_{ia}^{jb}(i\omega) \tilde{P}_{ia}(i\omega) \tilde{P}_{jb}(i\omega) \quad (28)$$

while for the NOF-c-$G_0W_0$+SOSEX@GM method it can be written as

$$E_{\text{can}}^{\text{dyn}} = \sum_{i,j}^{\text{occ}} \sum_{a,b}^{\text{virt}} \sum_{s}^{N_{ab}} \mathcal{A}_i \mathcal{A}_j \frac{\tilde{w}_{ia} (\tilde{X}_{jb}^s + \tilde{Y}_{jb}^s)}{\tilde{e}_i - \tilde{e}_a - \Omega_s} \times [2\langle ij | ab \rangle - \langle ij | ba \rangle] \quad (29)$$

where the spin has been already accounted.

### 2.2.3 $E_{\text{can}}^{\text{nd}}$ contribution

In general, all X corrections introduced in the previous section lacks non-dynamic correlation, which is well recovered by PNOF7s. Neverthe-
less, we cannot simply add these contributions since double counting occurs.

Similar to the attenuation coefficients $C_p$, the intra-subspace attenuation coefficients $\Lambda_p$ are defined$^{49}$

$$ \Lambda_p = 1 - |1 - 2n_p| \quad (30) $$

$\Lambda_p$ represents the amount of non-dynamic electron correlation in each orbital as a function of its occupancy. Note that $\Lambda_p$ goes from zero for empty or fully occupied orbitals to one if the orbital is half occupied. Analogously, the inter-subspace attenuation coefficients are directly taken as $2n_p h_p$.

Using these functions, the $E_{\text{can}}^{\text{nd}}$ correction is obtained from the pure non-dynamic component of the PNOF7s as a sum of the static intra-subspace and inter-subspace correlation energies (see Ref. 53), namely

$$ E_{\text{can}}^{\text{nd}} = \sum_{g=1}^{L/2} \sum_{p,q \in \Omega_g, p \neq q} \sqrt{\Lambda_p \Lambda_q} \Pi_{pq}^{\text{tra}} K_{pq} - \sum_{g=\frac{L}{2}+1}^{N_{\Omega}} \frac{K_{gg}}{4} - 4 \sum_{f \neq g} \sum_{p \in \Omega_f} \sum_{q \in \Omega_g} n_p n_q h_p h_q K_{pq} \quad (31) $$

In Eq. (31), the second term eliminates the intra-pair $\alpha\beta$ contribution of the singly occupied orbitals to the energy (21) because in each pure state $|SM\rangle$ of the ensemble there is no such interaction.$^{54}$

# 3 Results and discussions

## 3.1 Computational details

All calculations have been performed with DoNOF code,$^{71}$ while we have followed the implementation of MOLGW package$^{72}$ for the evaluation of MBPT quantities and Refs. 14 and 68. For NOF-c-CCSD, we have implemented the CCSD amplitudes equations following Ref. 18, where the restricted $\{\tilde{e}_p\}$ and $\{\langle pq|rs \rangle\}$ integrals are employed. Let us remark that for multiplets we employ the standard singlet restricted MBPT expressions, because the $A_i$ factors are only used in the evaluation of $E_{\text{can}}^{\text{dyn}}$.

Thus, any other amplitudes optimization procedure valid for restricted singlets could have also been employed. 40-50 frequencies were employed for the integration of the frequency dependence in Eq. (24) using a Gauss-Hermite quadrature. Reference values obtained at CCSD(T) and UCCSD levels were computed with Gaussian 03.$^{73}$ Open MOLCAS$^{74}$ was employed for CASSCF calculations including second-order perturbative corrections (CASPT2). The CASPT2 dissociation energies ($D_e$) were validated using the UCCSD ones, which produces comparable values for the studied systems (see Table 1). Let us remark that UCCSD $D_e$ values presented in this work were obtained as the difference between the energy of the optimized geometry and the isolated atoms. In the case of $H_2$, total energies are used to represent potential energy curves (PECs). For the rest of systems PECs are presented in terms of $D_e$. All bonds are assumed to be broken at an inter-nuclear separation equal to 10.0 Å, except for LiF where 12.0 Å is considered.

## 3.2 $H_2$: A simple covalent single-bond case

Our first analysis focuses on $H_2$ bond cleavage process in its ground state ($X^1\Sigma_g^+$). To evaluate the performance and facilitate the analysis of NOF-c-X methods, we have also included the PECs obtained using the corresponding X methods (i.e. the standard MP2 and RPA results applied as a one-shot correction on top of a restricted HF calculation).

In Fig. 2, PECs are shown using solid lines for Full Configuration Interaction (FCI) and NOF-c-X methods, and dotted-lines for X methods. For this particular case, the CCSD and PNOF7s are exact, therefore they coincide with the FCI PEC. We notice that in the non-dynamic correlation limit (when the bond is fully broken) we obtain the correct asymptotic behaviour for the system regardless of the NOF-c-X method employed, in contrast to the case of MP2 and RPA. From the error curves, i.e. $E^X - E^{\text{FCI}}$ against $R_{H-H}$, we arrive to the same conclusion as errors tend to zero when the bond
<figure><img src="image_3.png" /><figcaption>Figure 2: Potential energy curves and absolute errors of energies of H₂ computed using cc-pVTZ basis set with FCI, NOF-c-X, and X methods (X=MP2, RPA, CCSD). PNOF7s and CCSD coincide with FCI.</figcaption></figure>

is stretched. This example shows how the attenuation procedure switches off the $E_{\text{can}}^{\text{dyn}}$ contribution and only the $E_{\text{can}}^{\text{nd}}$ contribution remains active when the bond is stretched. Dynamic correlation effects are dominant near the equilibrium geometry and neither NOF-c-MP2 nor NOF-c-RPA are able to describe this region exactly. Regarding X methods, it is clear that when the ONs are close to 0 or 1 all NOF-c-X methods tend to their corresponding X counterparts. This result suggests that the attenuation coefficients switch off the $E_{\text{can}}^{\text{nd}}$ contribution near the equilibrium geometry, and thereby only the $E_{\text{can}}^{\text{dyn}}$ term remains active. In the intermediate region the nature of electron correlation changes from dynamic to non-dynamic as the bond is stretched. Here there is a slight difference between NOF-c-X and FCI, so all NOF-c-X methods introduce artifacts. For H₂ the NOF-c-MP2 and the NOF-c-CCSD methods show a tiny maximum around R<sub>H-H</sub> = 2.2Å. Indeed, the maximum corresponding to NOF-c-CCSD is due to an underestimation of ~0.007 a.u. of the electronic correlation energy. Similarly, the NOF-c-MP2 method depicts a similar maximum.

Finally, we observe that NOF-c-RPA overestimates the $D_e$ of H₂ while NOF-c-MP2 underestimates it. Only NOF-c-CCSD provides almost the exact value, which slightly deviates from the FCI result.

## 3.3 (CH₄)₂: The weak-interactions problem

The $(CH_4)_2$ dimer has been widely used$^{75}$ to study the nature of weak-interactions.$^{76,77}$ For this systems, different orientations lead to different PECs depending, especially, on the Hydrogen bonds formed. In this work, we have chosen the eclipsed orientation where only two Hydrogen atoms approach each other, as shown in Fig. 3. Initially, we performed the geometry optimization employing the aug-cc-pVDZ basis set of an isolated CH₄ specie at the MP2 level. Then, PECs were constructed by keeping frozen all coordinates except the inter-molecular distance represented by the Carbon-Carbon distance.

For systems formed by weak-interactions, the PEC of PNOF7s shows a purely repulsive profile. This fact was already observed using the PNOF7s method on He₂, HeNe, and Ne₂ systems.$^{53}$ Here, we also confirm this observation and show that all NOF-c-X methods tend to their pure X counterparts.

The small difference between NOF-c-X and X is due to two major effects: a) the attenuation coefficients modifies the integrals, and b) the {$\tilde{\phi}$} basis does not exactly reproduce the HF one. As the computational time required to perform the optimization with PNOF7s functional is longer than the time required by doing HF, the X methods should be preferred over the
<figure><img src="image_4.png" /><figcaption>Figure 3: Potential energy curves ($\Delta E = E_{R_{C-C}} - E_{R_{C-C}=10\AA}$) of $(CH_4)_2$ and absolute errors computed with aug-cc-pVDZ basis set and CCSD(T), PNOF7s, and NOF-c-X (X=MP2, RPA, CCSD).</figcaption></figure>

NOF-c-X ones when systems are dominated by weak interactions.

The error curves, i.e. $\Delta E^X - \Delta E^{\text{Ref.}}$, confirm that NOF-c-X methods account for the missing dynamic correlation energy. For this particular example, the NOF-c-X methods always underestimate the CCSD(T) reference correlation energies.

## 3.4 CO: Breaking a multiple bond

When single bonds are broken and neutral species are formed, the multi-configurational character of a state is enhanced. Then, SD wavefunctions fail, and in some cases an unrestricted formalism can overcome this issue, but at the price of introducing spin contamination.⁷⁸ When multiple bonds are broken (e.g. on CO molecule in the $X^1\Sigma^+$ state) even the unrestricted formalism is unable to give a proper description.

NOFs are suitable for treating the multi-configurational character. In Fig. 4, we can observe that the PNOF7s PEC exhibits a good qualitative shape when the bonds are broken. However, the missing dynamic electron correlation in the equilibrium region makes PNOF7s underestimate $D_e$. NOF-c-X approximations produce an improvement w.r.t. PNOF7s on the predicted value of $D_e$ for any X approximation.

From Fig. 4 we observe that NOF-c-RPA underestimates $D_e$, a behavior that is also observed in most of the systems studied in this work. Nevertheless, the NOF-c-RPA approximation provides the lowest error in this case. Regarding NOF-c-MP2 and NOF-c-CCSD, we have observed that the NOF-c-CCSD total energies are lower than NOF-c-MP2 ones, but dissociation energies do not show any clear tendency.

Finally, as it happened in the $H_2$ case, the major deviations w.r.t. the reference PEC occur in the region where the character of the correlation energy changes (i.e. in the intermediate region). Indeed, the non-dynamic correlation energy character is enhanced around the $R_{C-O} = 1.8$ Å geometry as the bond is stretched, which leads to important changes on the attenuation coefficients as the ONs tend to 0.5 rapidly. Moreover, the error curves confirm that NOF-c-X methods fail to describe the intermediate region as large differences are shown between 1.5 Å and 2.5 Å, where NOF-c-X methods underestimates the electronic correlation energy. Despite producing the lowest errors in the PEC profile (~0.03 a.u.), NOF-c-CCSD shows slight deviations in the way it reaches the horizontal asymptote.
<figure><img src="image_5.png" /><figcaption>Figure 5: Potential energy curves ($ΔE = E_{R_{N-H}} - E_{R_{N-H}=10Å}$) of NH (on X$^3\Sigma^-$ state) and absolute errors computed using cc-pVDZ basis set with CASPT2 (including 4 active electrons on 4 active orbitals), PNOF7s, and NOF-c-X (X=MP2, RPA, and CCSD) methods.</figcaption></figure>

$D_e$ values are observed (especially for multiplet state cases).

We have computed relative errors for the predicted $D_e$. The corresponding average relative error is ~10% when NOF-c-CCSD is employed to compute singlet state systems. NOF-c-MP2 and NOF-c-RPA approximations provide average relative errors of ~12% and ~13%, respectively. The latter points out that NOF-c-CCSD overcomes the approximations based on MP2 and RPA for singlet state systems. In the case of spin multiplets, larger average relative errors were obtained (~19% for NOF-c-CCSD, ~16% for NOF-c-MP2, and ~21% for NOF-c-RPA), but still an important improvement w.r.t PNOF7s errors (~36%) is observed. NOF-c-X methods do not reach the accuracy of UCCSD (see Table 1). Nevertheless, let us remark that UCCSD is unable to produce reasonable PECs when multiple bonds are broken while NOF-c-X approximations do so.

We have analyzed the predicted equilibrium distances, and all NOF-c-X approximations provided accurate results. The maximum relative error (~5%) was obtained by using NOF-c-RPA on $He_2$. Certainly, the geometries predicted by PNOF7s are already accurate enough, and PNOF7s only fails to predict systems formed by weak interactions.

Finally, let us recall that the definition of the attenuation coefficients, Eqs. (18) and (17) was fixed to the original values presented in Refs.

53 and 54. Consequently, NOF-c-X approximations fail in general to reproduce the curvature of the PEC from the equilibrium distance to the horizontal asymptote. In other words, the interval on the PEC where both dynamic and non-dynamic correlation effects are important, is not properly described by any NOF-c-X approximation. Hence, NOF-c-X approximations also fail to describe a critical system such as $Be_2$, where both electron correlation flavors need to be well balanced.

## 3.7 More NOF-c-X approximations

Concerning other NOF-c-X methods, all these approximations provide accurate equilibrium geometries but have troubles to reproduce dissociation energies (the results obtained with these approximations can be found in the supplementary material). Below we sum up these results. We have observed that NOF-c-$G_0W_0$@GM overestimates the $D_e$ for singlet state systems because it provides too large $E_{dyn}$ values in the dynamic electron-correlation regime (near the equilibrium distance). This can be explained from the nature of the $G_0W_0$@GM approximation, which misses the kinetic correlation energy$^{80}$ (usually labeled as $T_c$) and also the second-order exchange correction.$^{13}$

NOF-c-$G_0W_0$@GM provides reasonable re-
<figure><img src="image_6.png" /><figcaption>Figure 6: Potential energy curves ($\Delta E = E_{R_{X-X}} - E_{R_{X-X}=6Å}$ for X=N or F) of N₂ and F₂; absolute errors w.r.t. FCI computed using cc-pVTZ basis set with PNOF7s and NOF-c-X (X=MP2, RPA, and CCSD).</figcaption></figure>
sults only for multiplet state systems due to er-
ror compensation. Indeed, in the non-dynamic
correlation limit the $E_{can}^{nd}$ contribution is max-
imal (see Eq. 31) while in the dynamic cor-
relation regime $G_0W_0$@GM slightly overesti-
mates the $E_{can}^{dyn}$ energy, thus, these two ef-
fects compensate making this method applica-
ble for multiplet state species. The inclusion of
the second-order exchange contribution, as it
is done in NOF-c-$G_0W_0$+SOSEX@GM, makes
this approximation perform better than NOF-
c-$G_0W_0$@GM for singlet state species.

In NOF-c-$G_0W_0$+SOSEX@GM the effect of
the second-order exchange reduces the amount
of $E_{can}^{dyn}$ retrieved in the dynamic electron-
correlation regime making the $D_e$ to lie close
to the NOF-c-RPA one. Lastly, the NOF-
c-RPA+SOSEX approximation underestimates
all $D_e$ because of the presence of the second-
order exchange. As it happened with NOF-
c-$G_0W_0$+SOSEX@GM, the second-order ex-
change reduces the amount of $E_{can}^{dyn}$ retrieved
and, hence, it is not possible to compensate the
large amount of non-dynamic electron correla-
tion energy. That is to say, the dynamic and
non-dynamic correlation energy contributions
are not balanced.

## 3.8 Computational aspects

In terms of computational gain, let us first focus
on the NOF-c-MP2 method, which retrieves the
NOF-oiMP2 method<sup>53</sup> by construction. Since
the NOF-c-MP2 requires a single diagonaliza-
tion of the attenuated Fock matrix, the canon-
icalization represents a computational advan-
tage over the NOF-oiMP2 formulation that re-
quires solving the residual equations.

The canonicalization also gives access to less
computationally demanding MBPT approxi-
mations; such as RPA.⁸¹ Indeed, NOF-c-RPA
produces competitive results compared to other
NOF-c-X approximations. In principle, the so-
lution of Casida-like equation to compute the
RPA approximation scales as $M^6$, which is usu-
ally the limiting step in using this method.
However, the inclusion of imaginary times in-
stead of frequencies permits to reach a cu-
bic scaling,²¹,⁸²⁻⁸⁴ which makes NOF-c-RPA a

competitive approach for retrieving dynamic
and non-dynamic electron correlation effects in
large systems (e.g. polymers, periodic systems,
among others.), where MP2 and CCSD calcula-
tions may become computationally prohibitive.
To attain this scaling the density fitting tech-
nique for computing two-electron repulsion in-
tegrals is required. This technique was recently
implemented in DO NOTF⁷¹ for the energy mini-
mization of PNOF approximations.⁸⁵

Lastly, NOF-c-CCSD calculations are more
expensive than NOF-c-MP2 and thereby than
NOF-oiMP2, since the amplitudes equations
need to be solved for the CC correction.
Notwithstanding, much effort has been put
to make CCSD computationally affordable⁸⁶⁻⁸⁸
and it is considered a routine method in many
computational chemistry packages.⁸⁹,⁹⁰ From
Table 1 we notice that the small gain on the
precision of the $D_e$ makes it worth to use NOF-
c-CCSD over any other NOF-c-X approxima-
tion (especially for singlets).

# 4 Closing remarks

In this work, we have introduced a canon-
icalization procedure to establish a link be-
tween NOFT and MBPT corrections (such as
MP2, RPA, CCSD, etc.). As a result, a fam-
ily of methods called as NOF-c-X approxima-
tions is proposed. Our numerical tests reveal
that RPA-based and CCSD-based corrections
provide similar relative errors and MAE (see
Table 1) to the results obtained using a MP2-
based correction (~12% for singlets and ~20%
for multiplets), which proves that NOF-c-X ap-
proximations yield competitive approximations
to NOF-oiMP2 and can reduce the computa-
tional cost, as well as improve the method ac-
curacy.

We have explored some MBPT corrections
based on the Galitiskii-Migdal functional;
also including second-order exchange correc-
tions. Our results suggest that only NOF-c-
$G_0W_0$+SOSEX@GM predicts reasonable $D_e$
values that are comparable with NOF-c-RPA
ones. Further studies could be performed with
more MBPT approximations, including for ex-
ample the recently developed Bethe-Salpeter
correlation energies.²⁷

Finally, let us highlight that PNOF7s pro-
duces a correlated density matrix while the ref-
erence system employed by NOF-c-X methods
is a single-determinant wave function for all
approximations, which facilitates the use of a
standard MBPT formulations. On the contrary,
the direct application of CCSD, for example, as
an on-top correction of PNOF7s may require a
modified CCSD version as proposed by Hollet
and Loos.⁹¹ In their case, a finite temperature
formulation of CCSD⁹² is required.

**Acknowledgement** The authors acknowl-
edge P.-F. Loos and J.W. Hollett for sharing
their FCI data. M.R.M. acknowledges Mr. S.
Sitkiewicz for his aid on CASPT2 calculations.
F.B. and M.R.M. acknowledge the financial
support of the Cross-Disciplinary Program on
Numerical Simulation of CEA, the French Al-
ternative Energies and Atomic Energy Commis-
sion. M.P. acknowledges the financial support
of MCIU/AEI/FEDER, UE (PGC2018-097529-
B-100) and Eusko Jaurlaritza (Ref. IT1254-
19). I.M. acknowledges a DOKBERRI 2020 II
grant from the University of the Basque Coun-
try (UPV/EHU).

## Supporting Information Available

We have collected in the supporting informa-
tion a deeper analysis of the attenuation coeffi-
cients, the expressions used to solve the Casida-
like equations, and some comments about the
RPA and G₀W₀@GM correlation energies. We
also provide a spreadsheet containing the re-
sults obtained in this work (i.e. the *De* and
equilibrium distances).
<table><caption>Table 1: Absolute errors of the $D_e$ in KJ/mol.</caption><thead><tr><th>Molecule</th><th>PNOF7s</th><th>NOF-c-MP2</th><th>NOF-c-RPA</th><th>NOF-c-CCSD</th><th>UCCSD</th><th>Basis set</th><th>Ref.</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>0</td><td>16</td><td>35</td><td>4</td><td>0</td><td>cc-pVTZ</td><td>FCI</td></tr><tr><td>He<sub>2</sub></td><td>0.08</td><td>0.03</td><td>0.05</td><td>0.01</td><td>0.01</td><td>aug-cc-pVTZ</td><td>CCSD(T)</td></tr><tr><td>(CH<sub>4</sub>)<sub>2</sub></td><td>2.99</td><td>0.50</td><td>0.24</td><td>0.63</td><td>0.47</td><td>aug-cc-pVDZ</td><td>CCSD(T)</td></tr><tr><td>LiH</td><td>5</td><td>8</td><td>33</td><td>14</td><td>5</td><td>cc-pVTZ</td><td>FCI</td></tr><tr><td>LiF</td><td>113</td><td>42</td><td>7</td><td>29</td><td>7</td><td>cc-pVTZ</td><td>FCI<sup>a</sup></td></tr><tr><td>Li<sub>2</sub></td><td>0</td><td>10</td><td>14</td><td>14</td><td>4</td><td>cc-pVDZ</td><td>CASPT2(2,2)</td></tr><tr><td>HF</td><td>117</td><td>31</td><td>51</td><td>26</td><td>10</td><td>cc-pVTZ</td><td>FCI<sup>a</sup></td></tr><tr><td>F<sub>2</sub></td><td>73</td><td>67</td><td>5</td><td>60</td><td>25</td><td>cc-pVTZ</td><td>FCI<sup>a</sup></td></tr><tr><td>CO</td><td>183</td><td>42</td><td>33</td><td>35</td><td>47</td><td>cc-pVDZ</td><td>CASPT2(6,6)</td></tr><tr><td>N<sub>2</sub></td><td>188</td><td>135</td><td>52</td><td>130</td><td>12</td><td>cc-pVTZ</td><td>FCI<sup>a</sup></td></tr><tr><td>H<sub>2</sub>O<sup>b</sup></td><td>126</td><td>3</td><td>25</td><td>15</td><td>2</td><td>cc-pVDZ</td><td>CASPT2(4,4)</td></tr><tr><td>O<sub>2</sub>(b<sup>1</sup>&Sigma;<sup>+</sup><sub>g</sub>)<sup>c</sup></td><td>215</td><td>9</td><td>79</td><td>5</td><td>32</td><td>cc-pVDZ</td><td>CASPT2(8,6)</td></tr><tr><td>NH(1&Sigma;<sup>+</sup>)<sup>c</sup></td><td>61</td><td>35</td><td>35</td><td>21</td><td>1</td><td>cc-pVDZ</td><td>CASPT2(4,4)</td></tr><tr><td>BeH(X<sup>2</sup>&Sigma;<sup>+</sup>)<sup>c</sup></td><td>42</td><td>4</td><td>11</td><td>25</td><td>3</td><td>cc-pVTZ</td><td>FCI</td></tr><tr><td>O<sub>2</sub>(X<sup>3</sup>&Sigma;<sup>+</sup><sub>g</sub>)<sup>c</sup></td><td>325</td><td>168</td><td>192</td><td>175</td><td>14</td><td>cc-pVTZ</td><td>CASPT2(8,6)</td></tr><tr><td>NH(X<sup>3</sup>&Sigma;<sup>-</sup>)<sup>c</sup></td><td>42</td><td>8</td><td>14</td><td>8</td><td>0</td><td>cc-pVDZ</td><td>CASPT2(4,4)</td></tr><tr><td>CN(X<sup>2</sup>&Sigma;<sup>+</sup>)<sup>c</sup></td><td>276</td><td>158</td><td>237</td><td>172</td><td>1</td><td>cc-pVDZ</td><td>CASPT2(9,8)</td></tr><tr><td>MAE (singlet)</td><td>83</td><td>31</td><td>28</td><td>27</td><td>15</td><td></td><td></td></tr><tr><td>MAE</td><td>104</td><td>43</td><td>48</td><td>43</td><td>13</td><td></td><td></td></tr></tbody></table>

<footnote><sup>a</sup> FCI results were taken from Ref. 91; <sup>b</sup> The symmetric stretching of the two Hydrogen atoms was scanned keeping a fixed angle (H-O-H of 103°; <sup>c</sup> In parenthesis we have specified the spatial and spin state computed (for the rest of systems we used the ground state); <sup>d</sup> The reference number of electrons and the size of the active space of the CASSCF calculation used as reference for CASPT2 is specified in parenthesis in the usual convention (number of electrons, number of orbitals).</footnote>
References

(1) Kelly, H. P. Many-body perturbation theory applied to atoms. *Phys. Rev.* **1964**, *136*, B896.

(2) Hedin, L. New Method for Calculating the One-Particle Green's Function with Application to the Electron-Gas Problem. *Phys. Rev.* **1965**, *139*, A796.

(3) Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. *Phys. Rev.* **1934**, *46*, 618.

(4) Čížek, J.; Paldus, J. Correlation problems in atomic and molecular systems III. Rederivation of the coupled-pair many-electron theory using the traditional quantum chemical methodst. *Int. J. Quant. Chem.* **1971**, *5*, 359.

(5) Bartlett, R. J.; Purvis, G. D. Many-body perturbation theory, coupled-pair many-electron theory, and the importance of quadruple excitations for the correlation problem. *Int. J. Quant. Chem.* **1978**, *14*, 561.

(6) Bartlett, R. J. Many-body perturbation theory and coupled cluster theory for electron correlation in molecules. *Annu. Rev. Phys. Chem.* **1981**, *32*, 359.

(7) Baym, G. Self-consistent approximations in many-body systems. *Phys. Rev.* **1962**, *127*, 1391.

(8) Holm, B. Total energies from GW calculations. *Phys. Rev. Lett.* **1999**, *83*, 788.

(9) Galitskii, V. M.; Migdal, A. B. Application of quantum field theory methods to the many body problem. *Sov. Phys. JETP* **1958**, *139*, 96.

(10) Klein, A. Perturbation Theory for an Infinite Medium of Fermions. II. *Phys. Rev.* **1961**, *121*, 950.

(11) Bohm, D.; Pines, D. A collective description of electron interactions: III. Coulomb

interactions in a degenerate electron gas. *Phys. Rev.* **1953**, *92*, 609.

(12) Eshuis, H.; Bates, J. E.; Furche, F. Electron correlation methods based on the random phase approximation. *Theor. Chem. Acc.* **2012**, *131*, 1084.

(13) Ren, X.; Rinke, P.; Joas, C.; Scheffler, M. Random-phase approximation and its applications in computational chemistry and materials science. *J. Mater. Sci.* **2012**, *47*, 7447.

(14) Ren, X.; Rinke, P.; Scuseria, G. E.; Scheffler, M. Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks. *Phys. Rev. B* **2013**, *88*, 035120.

(15) Scuseria, G. E.; Scheiner, A. C.; Lee, T. J.; Rice, J. E.; Schaefer III, H. F. The closed-shell coupled cluster single and double excitation (CCSD) model for the description of electron correlation. A comparison with configuration interaction (CISD) results. *J. Chem. Phys.* **1987**, *86*, 2881.

(16) Scuseria, G. E.; Janssen, C. L.; Schaefer III, H. F. An efficient reformulation of the closed-shell coupled cluster single and double excitation (CCSD) equations. *J. Chem. Phys.* **1988**, *89*, 7382.

(17) Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. *Rev. Mod. Phys.* **2007**, *79*, 291.

(18) Stanton, J. F.; Gauss, J.; Watts, J. D.; Bartlett, R. J. A direct product decomposition approach for symmetry exploitation in many-body methods. I. Energy calculations. *J. Chem. Phys.* **1991**, *94*, 4334.

(19) Bruneval, F.; Marques, M. A. Benchmarking the starting points of the GW approximation for molecules. *J. Chem. Theory Comput.* **2013**, *9*, 324.

(20) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.;
Scheffler, M. Ab initio molecular simulations with numeric atom-centered orbitals. *Comput. Phys. Commun.* **2009**, *180*, 2175.

(21) Förster, A.; Visscher, L. Low-Order Scaling $G_0W_0$ by Pair Atomic Density Fitting. *J. Chem. Theory Comput.* **2020**, *16*, 7381.

(22) Förster, A.; Visscher, L. GW100: A Slater-Type Orbital Perspective. *J. Chem. Theory Comput.* **2021**, *17*, 5080.

(23) McClain, J.; Sun, Q.; Chan, G. K.-L.; Berkelbach, T. C. Gaussian-based coupled-cluster theory for the ground-state and band structure of solids. *J. Chem. Theory Comput.* **2017**, *13*, 1209.

(24) Bruneval, F.; Rodriguez-Mayorga, M.; Rinke, P.; Dvorak, M. Improved One-Shot Total Energies from the Linearized GW Density Matrix. *J. Chem. Theory Comput.* **2021**, *17*, 2126.

(25) Förster, A.; Visscher, L. Low-Order Scaling Quasiparticle Self-Consistent GW for Molecules, doi: 10.3389/fchem.2021.736591. *Front. Chem.* 698.

(26) Scuseria, G. E.; Henderson, T. M.; Sorensen, D. The ground state correlation energy of the random phase approximation from a ring coupled cluster doubles approach. *J. Chem. Phys.* **2008**, *129*, 231101.

(27) Loos, P.-F.; Scemama, A.; Duchemin, I.; Jacquemin, D.; Blase, X. Pros and Cons of the Bethe-Salpeter Formalism for Ground-State Energies. *J. Phys. Chem. Lett.* **2020**, *11*, 3536.

(28) Cioslowski, J. Density-driven self-consistent-field method: Density-constrained correlation energies in the helium series. *Phys. Rev. A* **1991**, *43*, 1223.

(29) Valderrama, E.; Ludeña, E.; Hinze, J. Analysis of dynamical and nondynamical components of electron correlation energy by means of local-scaling density-functional theory. *J. Chem. Phys.* **1997**, *106*, 9227.

(30) Liechtenstein, A.; Anisimov, V. I.; Zanen, J. Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators. *Phys. Rev. B* **1995**, *52*, R5467.

(31) Sharma, S.; Dewhurst, J.; Shallcross, S.; Gross, E. Spectral density and metal-insulator phase transition in mott insulators within reduced density matrix functional theory. *Phys. Rev. Lett.* **2013**, *110*, 116403.

(32) Mitxelena, I.; Piris, M.; Rodríguez-Mayorga, M. On the performance of natural orbital functional approximations in the Hubbard model. *J. Phys. Condens. Matter* **2017**, *29*, 425602.

(33) Helgaker, T.; Jorgensen, P.; Olsen, J. *Molecular electronic-structure theory*; John Wiley & Sons, 2014.

(34) Keller, S.; Boguslawski, K.; Janowski, T.; Reiher, M.; Pulay, P. Selection of active spaces for multiconfigurational wavefunctions. *J. Chem. Phys.* **2015**, *142*, 244104.

(35) Jeong, W.; Stoneburner, S. J.; King, D.; Li, R.; Walker, A.; Lindh, R.; Gagliardi, L. Automation of Active Space Selection for Multireference Methods via Machine Learning on Chemical Bond Dissociation. *J. Chem. Theory Comput.* **2020**, *16*, 2389.

(36) Gilbert, T. L. Hohenberg-Kohn theorem for nonlocal external potentials. *Phys. Rev. B* **1975**, *12*, 2111.

(37) Valone, S. M. Consequences of extending 1 matrix energy functionals pure-state representable to all ensemble representable 1 matrices. *J. Chem. Phys.* **1980**, *73*, 1344–1349.
(38) Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. *Proc. Natl. Acad. Sci. USA* **1979**, *76*, 6062–6065.

(39) Coleman, A. J. Structure of Fermion Density Matrices. *Rev. Mod. Phys.* **1963**, *35*, 668–687.

(40) Piris, M. Natural Orbital Functional Theory. *Adv. Chem. Phys.* **2007**, *134*, 387–427.

(41) Piris, M.; Ugalde, J. M. Perspective on natural orbital functional theory. *Int. J. Quantum Chem.* **2014**, *114*, 1169–1175.

(42) Husimi, K. Some Formal Properties of the Density Matrix. *Proc. Phys. Math. Soc. Jpn.* **1940**, *22*, 264–314.

(43) Pernal, K.; Giesbertz, K. J. H. In *Density-Functional Methods for Excited States*; Ferré, N., Filatov, M., Huix-Rotllant, M., Eds.; Springer International Publishing, 2016; pp 125–183.

(44) Mitxelena, I.; Piris, M.; Ugalde, J. M. In *State of The Art of Molecular Electronic Structure Computations: Correlation Methods, Basis Sets and More*; Ancarani, L. U., Hoggan, P. E., Eds.; Advances in Quantum Chemistry; Academic Press, 2019; Vol. 79; pp 155–177.

(45) Piris, M. In *Many-body approaches at different scales: a tribute to N. H. March on the occasion of his 90th birthday*; Angilella, G. G. N., Amovilli, C., Eds.; Springer: New York, 2018; Chapter 22, pp 283–300.

(46) Mazziotti, D. A. Structure of Fermionic Density Matrices: Complete N-Representability Conditions. *Phys. Rev. Lett.* **2012**, *108*, 263002.

(47) Piris, M. A natural orbital functional based on an explicit approach of the two-electron cumulant. *Int. J. Quantum Chem.* **2013**, *113*, 620–630.

(48) Piris, M. Interacting pairs in natural orbital functional theory. *J. Chem. Phys.* **2014**, *141*, 044107.

(49) Piris, M. Global Method for Electron Correlation. *Phys. Rev. Lett.* **2017**, *119*, 063002–5.

(50) Mitxelena, I.; Rodríguez-Mayorga, M.; Piris, M. Phase Dilemma in Natural Orbital Functional Theory from the N-representability Perspective. *Eur. Phys. J. B* **2018**, *91*, 109.

(51) Mitxelena, I.; Piris, M. An efficient method for strongly correlated electrons in one dimension. *J. Phys. Condens. Matter* **2020**, *32*, 17LT01.

(52) Mitxelena, I.; Piris, M. An efficient method for strongly correlated electrons in two-dimensions. *J. Chem. Phys.* **2020**, *152*, 064108.

(53) Piris, M. Dynamic electron-correlation energy in the natural-orbital-functional second-order-Møller-Plesset method from the orbital-invariant perturbation theory. *Phys. Rev. A* **2018**, *98*, 022504.

(54) Piris, M. Natural orbital functional for multiplets. *Phys. Rev. A* **2019**, *100*, 32508.

(55) Lopez, X.; Piris, M. Performance of the NOF-MP2 method in hydrogen abstraction reactions. *Theor. Chem. Acc.* **2019**, *138*, 89.

(56) Mercero, J. M.; Ugalde, J. M.; Piris, M. Chemical reactivity studies by the natural orbital functional second-order Møller-Plesset (NOF-MP2) method: water dehydrogenation by the scandium cation. *Theor. Chem. Acc.* **2021**, *140*.

(57) Piris, M.; Matxain, J. M.; Lopez, X.; Ugalde, J. M. The one-electron picture in the Piris natural orbital functional 5 (PNOF5). *Theor. Chem. Acc.* **2013**, *132*, 1298.
(58) Pulay, P.; Saebø, S. Orbital-invariant formulation and second-order gradient evaluation in Møller-Plesset perturbation theory. *Theo. chim. acta* **1986**, *69*, 357.

(59) Löwdin, P.-O.; Shull, H. Natural orbitals in the quantum theory of two-electron systems. *Phys. Rev.* **1956**, *101*, 1730–1739.

(60) Rodríguez-Mayorga, M.; Ramos-Cordoba, E.; Via-Nadal, M.; Piris, M.; Matito, E. Comprehensive benchmarking of density matrix functional approximations. *Phys. Chem. Chem. Phys.* **2017**, *19*, 24029.

(61) Piris, M.; Ugalde, J. M. Iterative diagonalization for orbital optimization in natural orbital functional theory. *J. Comput. Chem.* **2009**, *30*, 2078.

(62) Moller, C.; Plesset, M. S. Note on an approximate treatment for many-electron systems. *Phys. Rev.* **1934**, *46*, 618–622.

(63) Casida, M. E. In *Recent Advances in Density Functional Methods, Part I*; Chong, D., Ed.; World Scientific: Singapore, 1995; p 155.

(64) Marques, M.; Gross, E. Time-dependent density functional theory. *Annu. Rev. Phys. Chem.* **2004**, *55*, 427.

(65) By imposing the trapezoidal rule for the integration along the coupling parameter $\lambda$ we may write $\overline{W}_{ia}^{jb}(i\omega) = \int_{0}^{1} d\lambda \overline{W}_{ia}^{jb}(i\omega, \lambda) \approx \overline{W}_{ia}^{jb}(i\omega, \lambda = 1)/2$.

(66) Grüneis, A.; Marsman, M.; Harl, J.; Schimka, L.; Kresse, G. Making the random phase approximation to electronic correlation accurate. *J. Chem. Phys.* **2009**, *131*, 154115.

(67) Paier, J.; Janesko, B. G.; Henderson, T. M.; Scuseria, G. E.; Grüneis, A.; Kresse, G. Hybrid functionals including random phase approximation correlation and second-order screened exchange. *J. Chem. Phys.* **2010**, *132*, 094103.

(68) Ren, X.; Tkatchenko, A.; Rinke, P.; Scheffler, M. Beyond the random-phase approximation for the electron correlation energy: The importance of single excitations. *Phys. Rev. Lett.* **2011**, *106*, 153003.

(69) Ren, X.; Marom, N.; Caruso, F.; Scheffler, M.; Rinke, P. Beyond the GW approximation: A second-order screened exchange correction. *Phys. Rev. B* **2015**, *92*, 081104.

(70) Eshuis, H.; Furche, F. A parameter-free density functional that works for noncovalent interactions. *J. Phys. Chem. Lett.* **2011**, *2*, 983–989.

(71) Piris, M.; Mitxelena, I. DoNOF: an open-source implementation of natural-orbital-functional-based methods for quantum chemistry. Version available in Github at the url: https://github.com/marm314/DoNOFsw/tree/mbpt. *Comput. Phys. Commun* **2020**, *259*, 107651.

(72) Bruneval, F.; Rangel, T.; Hamed, S. M.; Shao, M.; Yang, C.; Neaton, J. B. molgw 1: Many-body perturbation theory software for atoms, molecules, and clusters. *Comput. Phys. Commun* **2016**, *208*, 149.

(73) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.;
Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision C.02. Gaussian, Inc., Pittsburgh, PA, 2003.

(74) Aquilante, F.; Autschbach, J.; Baiardi, A.; Battaglia, S.; Borin, V. A.; Chibotaru, L. F.; Conti, I.; De Vico, L.; Delcey, M.; Fdez. Galván, I., et al. Modern quantum chemistry with [Open] Molcas. *J. Chem. Phys.* **2020**, *152*, 214117.

(75) Tsuzuki, S.; Uchimaru, T. Magnitude of attraction in CF4-CH4 interactions: Are CF4-CH4 interactions weaker than average of CF4-CF4 and CH4-CH4 interactions? *J. Fluor. Chem.* **2020**, *231*, 109468.

(76) Via-Nadal, M.; Rodríguez-Mayorga, M.; Matito, E. Salient signature of van der Waals interactions. *Phys. Rev. A* **2017**, *96*, 050501.

(77) Via-Nadal, M.; Rodríguez-Mayorga, M.; Ramos-Cordoba, E.; Matito, E. Singling out dynamic and nondynamic correlation. *J. Phys. Chem. Lett.* **2019**, *10*, 4032-4037.

(78) Let us remark that NOF-c-X methods employ a restricted formalism and, therefore, they do not introduce any contamination.

(79) Pd, R.; Chandra, P. Ground and excited states of HNC, NH, and NH2 transients: Ab initio geometries, electronic structures, and molecular properties. *J. Chem. Phys.* **2001**, *114*, 7450.

(80) Hellgren, M.; Caruso, F.; Rohr, D. R.; Ren, X.; Rubio, A.; Scheffler, M.;

Rinke, P. Static correlation and electron localization in molecular dimers from the self-consistent RPA and GW approximation. *Phys. Rev. B* **2015**, *91*, 165110.

(81) Let us comment that for RPA and GW approximations, the Casida equation can be solved in any basis representation,⁹³ although the best converge rates are obtained when the canonical basis is employed.

(82) Rojas, H. N.; Godby, R. W.; Needs, R. J. Space-Time Method for *Ab Initio* Calculations of Self-Energies and Dielectric Response Functions of Solids. *Phys. Rev. Lett.* **1995**, *74*, 1827.

(83) Foerster, D.; Koval, P.; Sánchez-Portal, D. An O(N3) implementation of Hedin's GW approximation for molecules. *J. Phys. Chem.* **2011**, *135*, 074105.

(84) Kaltak, M.; Klimeš, J.; Kresse, G. Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si. *Phys. Rev. B* **2014**, *90*, 054115.

(85) Lew-Yee, J. F. H.; Piris, M.; del Campo, J. M. Resolution of the identity approximation applied to PNOF correlation calculations. *J. Chem. Phys.* **2021**, *154*, 064102.

(86) Asadchev, A.; Gordon, M. S. Fast and flexible coupled cluster implementation. *J. Chem. Theory Comput.* **2013**, *9*, 3385.

(87) Kaliman, I. A.; Krylov, A. I. New algorithm for tensor contractions on multicore CPUs, GPUs, and accelerators enables CCSD and EOM-CCSD calculations with over 1000 basis functions on a single compute node. *J. Comput Chem.* **2017**, *38*, 842.

(88) Piecuch, P.; Kucharski, S. A.; Kowalski, K.; Musiał, M. Efficient computer implementation of the renormalized coupled-cluster methods: the r-ccsd [t], r-ccsd (t),
cr-ccsd [t], and cr-ccsd (t) approaches.
*Comput. Phys. Comm.* **2002**, *149*, 71.

(89) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L., et al. Psi4: an open-source ab initio electronic structure program. *Wiley Interdisciplinary Reviews: Computational Molecular Science* **2012**, *2*, 556.

(90) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S., et al. PySCF: the Python-based simulations of chemistry framework. *Wiley Interdisciplinary Reviews: Computational Molecular Science* **2018**, *8*, e1340.

(91) Hollett, J. W.; Loos, P.-F. Capturing static and dynamic correlation with $\Delta$ NO-MP2 and $\Delta$ NO-CCSD. *J. Chem. Phys.* **2020**, *152*, 014101.

(92) White, A. F.; Chan, G. K.-L. A time-dependent formulation of coupled-cluster theory for many-fermion systems at finite temperature. *J. Chem. Theory Comput.* **2018**, *14*, 5690.

(93) Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. *Chem. Rev.* **2005**, *105*, 4009.
Graphical TOC Entry

<figure><img src="image_7.png" /></figure>