Generalizing Double-Hybrid Density Functionals: Impact of Higher-Order Perturbation Terms.

Connections between the Görling-Levy (GL) perturbation theory and the parameters of double-hybrid (DH) density functional are established via adiabatic connection formalism. Moreover, we present a more general DH density functional theory, where the higher-order perturbation terms beyond the second-order GL2 one, such as GL3 and GL4, also contribute. It is shown that a class of DH functionals including previously proposed ones can be formed using the present construction. Based on the proposed formalism, we assess the performance of higher-order DH and long-range corrected DH formed on the Perdew-Burke-Ernzerhof (PBE) semilocal functional and second-order GL2 correlation energy. The underlying construction of DH functionals based on the generalized many-body perturbation approaches is physically appealing in terms of the development of the non-local forms using more accurate and sophisticated semilocal functionals.

Introduction

Nowadays, density functional theory (DFT)[1,2] becomes a standard framework of performing the electronic structure calculations of atoms, molecules, and solids, being an indispensable tool for quantum chemists and solid-state physicists. Despite its theoretical exactness, in practical calculations, the many-body interactions are included in DFT through density functional approximations (DFA) of the exchange-correlation (XC)[3] energy term. Therefore, the accuracy of the DFT depends on the accuracy of the XC approximations constructed from the exact quantum mechanical constraints satisfaction[4−9] and/or the empirical way.[10−13] The quantum mechanical constraints that have been used to construct the semilocal XC approximations are also equally important for developing higher-order accurate hybrids and double-hybrid (DH) functionals. Within the known exact constraints, we find uniform or non-uniform coordinate transformations,[5,14−16] density gradient expansions,[17−22] high and low-density limit of the correlation energy functionals,[23−25] correct asymptotic behavior of the XC energy density,[26−33] dimensional crossover or the quasi-2D behavior of the XC energy,[34−37] bounds of exchange and XC energies,[38−44] and exact properties of the XC hole.[8,45−47] All of these have been used regularly to construct the different rungs of non-empirical functionals[48] including higher-order accurate wave function methods.

The first three rungs of the XC DFAs,[48] which are recognized as the local density approximation (LDA),[1] generalized gradient approximation (GGA),[49] and meta-generalized gradient approximation (meta-GGA),[50] are constructed using the local or semilocal quantities (i.e., density, the gradient of density, Laplacian of density, and Kohn–Sham (KS) kinetic energy density), being very accurate for the diverse nature of the molecular[51,52] and solid-state properties.[9,53−61] However, there are important limitations in semilocal functionals performance.[62] Several resolutions are adopted to improve the functional performance, such as the inclusion of Hartree–Fock (HF) exchange within the semilocal approximations. Functionals constructed in this way are known as hybrid functionals, being extensively used for the chemical and solid-state properties.[51,63−82] Despite the attempt of constructing the best hybrid density functional, an important dynamic correlation, which is the key of the ab initio wave functional theory (WFT), are missing in density functional correlation energy functionals. Inclusion of non-local virtual orbital-dependent terms is shown to overcome many problems of DFAs including the best hybrid density functionals. The resulting functionals developed by mixing a part of second-order Møller–Plesset (MP2)[83] correlation energy expression are known as the DH DFA, which connect the density functional world with that of the ab initio WFT. Beginning from the Grimme’s B2-PLYP[62] DFA, several DH functionals have been proposed (see, e.g., refs (84) and (85) and references therein).

In general, the construction of a DH functional starts from the global hybrid DFA (also known as DFA0) having the form

Also, the range-separated (RS) scheme is applied in the construction of the hybrid density functionals using the RS screening approach of the exchange hole. Such XC functionals are known as the screened hybrids (SH), and they can achieve the correct asymptotic behavior of exchange potentials, important for the barrier heights, charge transfer excitation, and many-electron self-interaction related problems.[86−122] The screening of the exchange energy part is achieved through the decoupling of the electron–electron interaction aswhere ω is the range-separation parameter. Despite the plethora of considerable successes, there are limitations of both the DFA0 and SH functionals, which have been further improved by mixing the semilocal correlation functional with the non-local, second-order Görling–Levy (GL2), (i.e., MP2-like expression neglecting the single excited term) correlation energy expression

Therefore, the DH DFA is considered as the fifth rung functional where the two parameters {a, a} are fixed either empirically (with respect to the benchmark data) or in the non-empirical way. In particular, there exist two known relations between {a, a} parameters, namely, a = (a)2 proposed in ref (123), and used by B2-PLYP[62] or DS1DH[123] XC functionals, and the cubic scaling relation[124]which is adopted in the hereto study. The DH functionals achieved their accuracy over a large range of molecular properties beyond the limitations of the DFA0 and SH functionals.[85] In particular, the improvement of many-electron self-interaction problems and non-covalent interactions are noticeable.

The formal construction of the DH was first attempted in ref (125) using the adiabatic connection (AC)-based coupling constant integral formalism, which is the basis of the density functional many-body perturbation theory. Physically, the DH functionals are constructed from the perturbation theory of the XC functionals in the weakly interacting limit where the exchange part is treated by HF, and correlation is approximated by the GL2 expression. The higher-order terms of the perturbation series, i.e., GL3, GL4, etc., can also be used through the perturbation theory to construct the general mth-order hybrids (mH); however, due to substantial numerical cost of evaluation of these terms, they never have been considered. Nevertheless, in general, a form of the mH functional can be proposed, which allows one to include all the higher-order perturbation terms. This is the main motivation of the present study. In this work, we “generalized the double-hybrid density functional” derived from the many-body perturbation theory, which connects the higher-order GL perturbation terms via the adiabatic connection formalism. The novelty of the present method is also attributed via the assessment of the series of double hybrids based on mth-order integrand and RS scheme.

In the following, we will present the key equations and features of the prime DH functionals proposed so far. Following this, we will present the theory and results obtained from constructed functionals. We will conclude with future prospects.

Rungs of Double-Hybrid Density Functions

Adiabatic Connection View on DH Functionals

The starting point of the construction of DH functionals can be encapsulated through the AC formalism, where the coupling constant integral is used to construct the exact expression of the XC functional aswhere V̂, U[ρ], and Ψ[ρ] are the Coulomb operator, the Hartree energy, and the anti-symmetric wave function that yields the correct many-electron interacting density ρ(r) and minimizes the expectation value of ⟨T̂ + λV̂⟩, where T̂ is the kinetic energy operator and λ is the coupling constant. In the AC formalism, two limiting cases are important for the DH DFA construction, namely, the λ → 0 (weak-interaction limit) and λ → 1 (full-interaction limit). When λ → 0, the GL2[23−25,126] becomes the dominant correction to EXX energy. These two terms are most commonly employed to construct the XC energy integrand in this regime:[24]where W0[ρ] = E[ρ] and . We remark, however, that for any small value of λ (λ > 0), higher-order terms contribute as well. Therefore, one should also consider their inclusion in the W integrand. In the full interacting limit, in turn, the W[ρ] is usually approximated by non-local or semilocal formula.[14,128] Thus, the DH functional connects, by construction, the XC worlds of both the semilocal DFT and WFT approaches.[128] The utilization of the GL2 correlation energy term allows one to include non-local dynamic correlation effects into the XC DFA going beyond the semilocal correlation DFA description.

The AC formalism was applied in the construction of many semi-empirical or non-empirical DH DFAs among which we can recall the DFA0-DH,[129] DFA0-2,[130] DFA-CIDH,[127] or DFA-CIDH[127] DFA-QIDH[131] to be the most successful. In particular, the quadratic integrand DH (QIDH)[131] DFA was constructed considering W as the quadratic function[132] of λ, i.e., W[ρ] = a1[ρ] + a2[ρ]λ + a3[ρ]λ2. In the following, we adopt a general mth-order integrand function, which allows one to go beyond the QIDH formalism.

Beyond the DFA-QIDH Model

Let us consider the XC m-hybrid (mH) functional defined aswhere ξ, k = 1, ..., m + 2 are parameters. The E and E denote standard semilocal DFAs of exchange and correlation functionals, respectively, whereas the E is the HF (or exact) exchange energy:where are the two-electron integrals. The expression for the GL2 correlation energy reads[23]where ⟨ij∥ab⟩ = ⟨ij | ab⟩ – ⟨ij | ba⟩ are the antisymmetrized two-electron integrals, and ε is the energy of the KS orbital p, whereas v and v̂ denotes the local, multiplicative KS and the non-local HF exchange potentials, respectively. The higher-order GL energy terms[23,133] are denoted generally as EGL and the indices i, j, ... and a, b, ... are used to denote occupied and virtual KS orbitals, respectively.

The ξ parameters in eq can be, in principle, determined using adiabatic connection formalism, which defines, via eq , the exact XC functional. The coupling-constant integrand of eq has the well-known small-λ expansion[24,134]

Hence, truncating the above summation at n = m, it reads

Now, as in ref (131), we can define W integrand aswhere a (k = 1, ..., m) are the density or orbital dependent functionals. By comparing eq and eq at λ → 0, we get

The a[ρ] can be set from similar reasoning as in ref (131). At the upper limit (λ = 1), the integrand can be approximated bywhere the parameter α (together with m) ranging between [ – m,1] controls the amount of EXX in eq . Comparing eq and eq at λ = 1, we obtain

Now, employing eq into eq leads, after some algebra, to

Direct comparison with eq allows to connect (in the quantize manner) the ξ parameters with orders of the GL perturbation theory as

Note that, when m → 1, we get

On the other hand, when m → ∞, eq recoverswhich, in general, exhibits a diverging behavior.[135,136] For an arbitrary (but small) value of m, the kth GL terms in eq are always scaled by the coefficient ξ < < 1. Thus, the diverging behavior of DFA-mIH DFA should be avoided.

Generalized mth Hybrid with Range-Separation Model

Having eq , the obvious viable next step is to extend it to the domain of the range-separated DH (RSDH) framework. The background of the perturbative RS functionals, it was initially adopted by Ángyán et al.(137) Later, the RSDH was adopted by Chai et al.(138) in ωB97X family of functionals (ωB97X-2) and by Brémond et al.(139,140) in QIDH-based functionals (RSX-QIDH). Both the ωB97X-2 and RSX-QIDH showed an improvement over the base ωB97X and QIDH functionals. Particularly, a significant improvement has been observed for the problems related to the self-interaction corrections, barrier heights, and ionization potential, which are important and viable for RSDH.

In this respect, we extend the generalized formalism presented in Section to the RS scheme. By using W[ρ] ≈ E(ω) + E(ω) + 2E, we obtain the following expression:

Similarly to the DFA-mIH family of functionals, the LRC-ωDFA-mH functionals also depend on the choice of two parameters, namely, m and ω. This is done by disregarding the dependence on the α parameter in eq . Note that, in this case, the present RS functional form is closely related to the LRC-ωPBEh decomposition proposed by Rohrdanz et al.,[94] which reduces to the standard PBEh hybrid for a certain choice of the two parameters, namely, mixing (a) and range-separated (ω) parameters.

We underline, however, that by choosing W[ρ] ≈ αE + (1 – α)[E(ω) + E(ω)] + 2E, we can define also the family of RSX-mIH functionals, where for m = 2 we recover the RSX-QIDH[139,140] DFA.

Choice of Functionals and α and ω Parameters

Based on the generalized consideration presented in Sections and 2.3, we consider two forms of the functionals to assess their performance for the quantum chemistry test sets. However, due to the substantial computational cost of the inclusion of higher-order GL terms, we restrict our assessment to DFAs whose form include only the second-order GL2 term (in eqs and ) and m = 2,3,4. The impact of higher-order terms is discussed in Section .

In eq , for the semilocal part of XC approximation, we consider the PBE and thus we denote this family of functionals as PBE-mIDH (m, α) as it is based on the mth-order integrand limited to second order and the m and α parameters. In the second case, we investigate LRC-ωDFA-mH functionals defined via eq where, as the SR part of DFA, we have chosen the ωPBE[92] functional and therefore the functional family is denoted as LRC-ωPBE-DH (m, ω) as it is based on the m and empirically fitted ω.

In order to make the PBE-mIDH and LRC-ωPBE-DH methods feasible, we need to set the values of α and ω parameters, respectively. In the case of α, we can directly utilize the cubic scaling relation (ξ2 = ξ13) from ref (124) usually adopted in the construction of novel DH DFAs. This allows one to express α as a function of m:

Thus, the PBE-mIDH family of functionals can be considered as “non-empirical” in this sense. The present PBE-mIDH model not only correctly recovers the PBE-QIDH[131] for m = 2, but also it recovers the PBE0-2[130] functional for m = 3. This allows one to view the construction of the PBE0-2 functional from an entirely new physical point of view. As m = 2 and m = 3 are the two previously proposed functionals, here, we additionally assess the functional performance for m = 4.

In the case of the LRC-ωPBE-DH (m, ω) family of functionals, we cannot establish a direct relation between m and ω. Thus, in order to determine the ω screened parameter, several schemes can be utilized, such as empirical fitting,[89,92,95,99,119,141] optimal tuning,[142] or imposing an additional constraint to the XC functional.[139,143]

In this study, in order to be consistent with the LRC-ωPBEh functional, we have tuned ω considering the mean absolute error (MAE) of the small representative test sets AE6 and BH6.[144] The AE6 is the atomization energies of six molecules simulating the error statistics of the full G2/148 molecular test set. Meanwhile, the BH6 is the hydrogen transfer barrier heights of six small molecules. In the optimization procedure, we have used the AE6 and BH6 geometries from Lynch and Truhlar,[144] and for the reference values, we consider the non-relativistic coupled-cluster singles-doubles with perturbative triples (CCSD(T)[145] values [FC-CCSD(T)/cc-pVQZ-F12] of ref (146).

In Figure , we report the MAE of AE6 and BH6 for few values of m. We observe that, for m = 2, the minimum of the MAE becomes at ω = 0.2 bohr–1. This choice of {m = 2, ω = 0.2} is certainly a good choice for BH6 also, where the MAE is obtained to be about 1 kcal/mol smaller than the LC-ωPBE and LRC-ωPBEh. Interestingly, in ref (94), the choice of ω is also 0.2 bohr–1, which is certainly matching with our present LRC-ωPBE-DH functional.

Figure 1

MAE (in kcal/mol) of AE6 and BH6 test sets for different values of m and ω.

In Table , we present the parameters {m, α} and {m, ω} of the studied functionals. Note that, for LRC-ωPBE-DH (m = 3, ω), we chose ω = 0.2 bohr–1, which gives a balanced description of both AE6 and BH6, while for LRC-ωPBE-DH (m = 4, ω), we chose ω = 0.5 bohr–1, considering solely the optimization of the AE6 test set.

Table 1

Choice of Different Sets of {m, α} and {m, ω} as Used in the Calculationsa

PBE-mIDH (m, α)		LRC-ωPBE-DH (m, ω)
m	α1		m	ω	MAE (kcal/mol)
2	0.080	AE6/BH6	2	0.2	5.5/1.08
3	0.175	AE6/BH6	3	0.2	5.3/1.382
4	0.217	AE6/BH6	4	0.5	4.2/2.84

The ω values are in bohr–1 unit.

Obtained via Eq. (18).

Minimum of AE6 is obtained at ω = 0.35 bohr–1 (3.6 kcal/mol) but choice of ω = 0.20 bohr –1 is the most preferable for both the AE6 and BH6.

Results

In this section, we discuss the general trends found in the performance of both families of functionals. The details about the computational setup can be found in Section .

Impact of the Higher-Order Terms on the Functional Performance

In Table , we report the MAE and mean absolute relative error (MARE) obtained for several ab initio WFT approaches together with the results obtained for both families of DH DFAs. The errors have been calculated with respect to CCSD(T) data obtained using the same computational setup. Additionally, to investigate the impact of higher-order terms (HOT) on the quality of predictions, we report the triple-hybrid (TH) (PBE-TH, LRC-ωPBE-TH) and quadruple-hybrid (QH) (PBE-QH, LRC-ωPBE-QH) results. To this end, the energies have been calculated via eqs and taking into account all terms up to mth order, on top of GKS orbitals obtained for m = 3 and m = 4 for TH and QH, respectively. As aforementioned, in order to evaluate the higher-order GLk energy contributions (GL3, GL4), we utilize the MP-like formulas (MP3 and MP4) (evaluated on the same GKS reference state), which give the dominant energy contribution. This is mostly dictated by the lack of quantum chemistry codes that are capable to calculate full higher-order GLk terms (see Section for more details).

Table 2

Error Statistic on Total Energies, Atomization Energies of AE6, and Binding Energies for Several Ab Initio Wave-Function Theory and Density Functional Theory Methodsa

					PBE-mIDH					LRC-ωPBE-DH
	CCSD	MP2	MP3	MP4	m = 2	m = 3	m = 4	PBE-TH	PBE-QH	m = 2	m = 3	m = 4	LRC-TH	LRC-QH
Total energies (TE)1
MAE [mHa]	8.93	16.81	11.82	2.42	5.57	6.14	7.11	5.08	3.11	5.25	5.86	6.84	4.81	4.92
MARE [%]	0.01	0.09	0.03	0.01	0.04	0.04	0.05	0.03	0.02	0.03	0.04	0.03	0.03	0.01
AE61
MAE3	10.3	12.4	6.7	5.9	4.0	4.7	6.2	3.8	2.7	3.5	4.5	2.3	3.4	4.3
MARE [%]	3.01	3.58	2.36	1.59	1.34	1.23	1.68	1.08	0.61	1.10	1.08	0.71	1.03	1.63
Non-covalent interactions (NCI)2
MAE3	0.17	0.13	0.05	0.04	0.08	0.09	0.09	0.07	0.06	0.09	0.06	0.05	0.07	0.07
MARE [%]	16.27	20.76	10.17	3.49	20.21	19.90	19.39	19.16	15.22	41.40	32.06	34.64	32.46	33.27

The full results can be found in ref (147). The LRC-TH and LRC-QH denotes LRC-ωPBE-TH and LRC-ωPBE-QH, respectively.

Calculation performed in cc-pVQZ basis set, all electrons are correlated.

Calculation performed in aug-cc-pVQZ basis set, all electrons are correlated.

In kcal/mol.

In the case of total energies (TE), both MAE and MARE indicate that all DH DFAs give the results at least twice better than MP2 total energies, which yield a MAE of 16.81 mHa (MARE = 0.09%). The addition of HOT leads to further improvement. In the case of TH, the inclusion of the GL3 term leads to a decrease in MAE/MARE of about 1 mHa/0.01%, whereas for QH, the decrease is slightly larger (2 ÷ 3 mHa/0.02 ÷ 0.03%) in comparison to their m = 3 and m = 4 DFA counterparts. Moreover, one can note that TH and QH DFAs give results that lie between MP3 and MP4 predictions. This indicates that the inclusion of HOT might lead to quite substantial improvement of the results than the standard DH expressions.[148]

In the case of the AE6 atomization energies, we see that all DH DFAs outperform all reported ab initio WFT methods. Here, the worst performance is observed for the MP2 method (MAE = 12.4 kcal/mol, MARE = 3.58%) and surprisingly also for the CCSD method (MAE = 10.3 kcal/mol, MARE = 3.01%). In the case of the PBE-mIDH family of DFA, we see the worsening, and for LRC-ωPBE-DH, the improvement of the predictions along with an increase in m. For the latter species, this can be easily linked with the optimization of the ω parameter, which was performed with respect to the AE6 set. One can note that inclusion of HOT in both families of DFAs lead, in general, to improvement in the prediction. The PBE-TH and PBE-QH functionals improve upon their DH counterparts, giving here a MAE of 3.8 kcal/mol (MARE = 1.08%) and a MAE of 2.7 kcal/mol (MARE = 0.61%), respectively. In the case of LRC-ωPBE-TH, we observe a similar trend. The only exception is seen in the case of LRC-ωPBE-DH (m = 4) and LRC-ωPBE-QH DFAs where the latter exhibits almost twice worse performance with respect to its DH variant. This can be most probably related to the ω value, which is not optimal for the LRC-ωPBE-QH energy expression.

In order to check how these trends change for AE6 with basis set size, in Fig. S5 of Ref (147)., we report the MARE calculated with respect to CCSD(T) data obtained using the cc-pVXZ[149] basis set family (where X = D, T, Q). As one can note, for most of DFAs we observe, a drastic improvement between DZ and TZ results. The differences between TZ and QZ results are several times smaller. One exception can only be noted for LRC-ωPBE-DH (m = 4) and LRC-ωPBE-QH functionals. Here, the DZ results are quite similar to TZ and QZ counterparts. As previously noted, this behavior can be related to the ω value, which was optimized solely with respect to the AE6 data set. Nevertheless, we conclude that, for the DZ basis set, all DFAs give the worst performance, whereas TZ and QZ results do not differ significantly from each other.

Finally, we have tested all DFAs with respect to a small set of NCI systems.[150] Similar to AE6 data, the worst performance among WFT approaches is given by CCSD (MAE = 0.17 kcal/mol, MARE = 16.27%) and MP2 (MAE = 0.13 kcal/mol, MARE = 20.76%) methods. The MP4, in turn, gives the best predictions here, yielding a MAE of 0.04 kcal/mol (MARE = 3.49%). All PBE-mIDH functionals show almost the same performance with MAEs in the range of 0.08–0.09 kcal/mol (MARE between 19.39% and 20.21%). The addition of HOT terms, also here, leads to an improvement mostly visible for PBE-QH DFA, which gives a MAE of 0.06 kcal/mol (MARE = 15.22%). In the case of the LRC-ωPBE-DH family of functionals, the inspection of MAEs shows almost identical performance (MAE between 0.5 and 0.09 kcal/mol) for all DFAs. However, MAREs exhibit a substantial increase in error (almost twice), which is caused by the very bad performance of LRC-ωPBE-DH (large MAE) for weakly interacting systems (see Tab. S2 in ref (147)). In ref (147), we also report MAE and MARE without considering WI systems. One can immediately note the substantial decrease in MARE values mostly visible for the LRC-ωPBE-DH family of functionals. This indicates that WI binding energies can be considered as very sensitive tools for testing the quality of quantum chemistry methods.

To conclude, the inclusion of HOT leads in general to the improvement of the result. However, the substantial cost of MP3-like and MP4-like correlation contribution makes it almost impossible to apply, especially for large molecular systems.

Performance without Considering HOT

Molecular Properties

To assess the PBE-mIDH and LRC-ωPBE-DH families (m = 2,3,4) of functionals, we report in Table the MAEs for several standard benchmark tests. For comparison, we also supplied the results obtained from PBE0-DH, PBE0, LC-ωPBE, and RSX-QIDH methods. Note that, except LC-ωPBE and PBE0, all other functionals belong to DH DFA species.

Table 3

Mean Absolute Errors (MAE) (in kcal/mol) for the Benchmark Tests Obtained Using Various Double-Hybrid Functionalsa

	MP2	PBE0-DH	PBE-QIDH (m = 2)	PBE0-2 (m = 3)	PBE-mIDH (m = 4)	LRC-ωPBE-DH (m = 2)	LRC-ωPBE-DH (m = 3)	LRC-ωPBE-DH (m = 4)	PBE0	LC-ωPBE (ω = 0.40)[1]	RSX-QIDH
Main group thermochemistry (MGT)
AE6	11.6	5.8	5.3	5.3	6.0	5.5	5.3	4.2	5.4	5.0	7.0
G2/148	10.7	5.7	5.4	5.6	6.2	5.3	5.4	4.9	4.5	4.4	7.8
G21EA	4.02	3.39	3.03	2.75	2.75	3.30	2.72	3.52	2.91	2.98	3.39
G21IP	3.37	3.31	2.64	2.23	2.11	2.85	2.21	2.62	3.78	4.78	3.12
PA26	1.52	2.55	2.04	1.60	1.38	2.20	1.60	1.80	2.52	2.60	1.80
BH76RC	2.50	1.44	1.52	1.68	1.82	1.51	1.69	1.78	1.59	3.62	2.12
SIE4×4	1.62	7.36	3.33	1.68	0.99	2.67	2.00	1.88	14.12	9.42	1.91
TMAE	5.04	4.22	3.32	2.97	3.13	3.33	2.98	2.95	4.97	4.68	3.87
Barrier heights (BH)
BH6	3.54	1.85	0.88	1.39	1.82	1.08	1.38	2.84	4.59	1.39	1.42
HTBH38	3.61	1.53	0.95	1.50	1.89	1.11	1.48	2.97	4.21	1.20	1.28
NHTBH38	5.49	1.70	1.78	2.55	3.04	2.20	2.63	4.82	3.76	2.22	2.57
TMAE	4.21	1.69	1.20	1.81	2.25	1.46	1.83	3.54	4.18	1.60	1.75
Non-covalent interactions (NCI)
HB6	0.24	0.31	0.29	0.31	0.32	0.26	0.24	0.25	0.32	0.71	0.91
DI6	0.57	0.30	0.24	0.25	0.32	0.28	0.22	0.20	0.36	0.81	0.30
CT7	0.56	0.44	0.40	0.42	0.44	0.33	0.31	0.27	0.84	1.12	0.57
PPS5	1.48	1.12	0.36	0.34	0.55	0.46	0.26	0.63	1.89	1.96	0.27
WI7	0.04	0.16	0.14	0.11	0.10	0.18	0.14	0.13	0.16	0.28	0.15
S22	1.39	1.60	0.74	0.26	0.37	0.92	0.34	0.28	2.37	2.75	0.71
TMAE	0.71	0.66	0.36	0.28	0.35	0.41	0.25	0.29	0.99	1.27	0.49
Overall
TMAE	3.32	2.19	1.63	1.69	1.91	1.73	1.69	2.26	3.38	2.52	2.04

The least MAE in each category is marked in bold. The total MAE (TMAE) is also mentioned at the bottom of each category.

ω = 0.40 bohr−1 as suggested in ref. (141)

One can immediately note that all DH functionals outperform the global hybrid PBE0 and RS hybrid LC-ωPBE in all cases. This indicates that the addition of the scaled MP2 term leads, in general, to the improvement of the results with respect to its global or RS hybrid predecessor.

The overall performance of all DH DFAs is quite similar. The total MAEs (TMAEs) can be found in the range of 1.63–2.26 kcal/mol. The best TMAE of 1.63 kcal/mol is given by PBE-QIDH (m = 2), which is closely followed by PBE0-2 (m = 3, TMAE = 1.69 kcal/mol) and LRC-ωPBE-DH (m = 3, TMAE = 1.69 kcal/mol) functionals. Slightly worse results are provided by LRC-ωPBE-DH (m = 2, TMAE = 1.73 kcal/mol) and PBE-mIDH (m = 4), which gives a TMAE of 2.26 kcal/mol. One can also identify general trends in the results depicted in Figure where we report the ratio of MAE (for m = 2,3,4) with respect to the MAE of the MP2 method (MAE), which represents the m = ∞ case for both families of functionals. In the case of PBE-mIDH DFAs (top of Figure ), increasing the amount of HF and MP2 terms leads, in general, to the worsening of the predictions. However, all functionals perform much better than the MP2 method. The change in overall performance between the m = 2 and m = 3 cases is very small. Nevertheless, an inspection of subset performance reveals some interesting trends. For the BH subset, we observe evident worsening of the results with m, whereas for MGT, the worse performance between DH cases is observed for m = 2. Here, along with m, we note worsening (AE6, G2/148, BH76RC) and improvement (G21EA, G21IP, PA26, SIE4×4) of the predictions. Almost the same is valid for the NCI subset, where the best performance is obtained for the m = 3 variant. This behavior is, however, dominated by the PPS5 and S22 results, which for m = 3 give the least MAE. For other NCI subsets, we observe the increasing (HB6, DI6, CT7) or decreasing (WI7) trend in MAE with an increase in m. Similar trends can be observed for the LRC-ωPBE-DH family of functionals (bottom of Figure ). The amount of HF and MP2 terms governs the general behavior of DFAs. With increasing m, we observe the improvement of the results for MGT and worsening for BH. In the case of NCI, likewise for the PBE-mIDH family of functionals, with increasing m, we observe an improvement in the description of most types of interactions especially visible for S22 subsets. This is in more detail shown in Table , where we report separately the MAEs for different subsets of complexes in the S22 test set. The PBE0-2 (m = 3) performs significantly better for vdW and mixed complexes. Meanwhile, LRC-ωPBE-DH (m = 3) is better for H-bonded complexes. This behavior can be probably justified by the fact that, for larger values of m, the SCF orbitals (optimized in the presence of 100% LR HF exchange term) give a better starting point (by reducing self-interaction error) for the description of long-range dispersion interaction included in the MP2 term. This is clearly seen in the case of the LRC-ωPBE-DH (m = 4) functional, which gives good performance for NCI and MGT benchmark sets, yielding here MAEs of 0.29 and 2.95 kcal/mol, respectively.

Figure 2

Overall and detailed (MGT, BH, and NCI subsets) trends in the results obtained for (top) PBE-mIDH and (bottom) LRC-ωPBE-DH families of functionals.

Table 4

Mean Absolute Errors (in kcal/mol) of the Subsets (H-Bond, Dispersion, and Mixed Complexes) of the S22 Test

interactions	MP2	PBE0-DH	PBE-QIDH (m = 2)	PBE0-2 (m = 3)	PBE-mIDH (m = 4)	LRC-ωPBE-DH (m = 2)	LRC-ωPBE-DH (m = 3)	LRC-ωPBE-DH (m = 4)	RSX-QIDH
H-bond	0.44	0.54	0.35	0.38	0.42	0.41	0.28	0.35	1.18
vdW	2.49	2.99	1.38	0.28	0.43	1.66	0.53	0.33	0.74
mixed	1.05	1.08	0.41	0.12	0.25	0.57	0.19	0.18	0.23

Looking at the MP2 results (representing the m = ∞ case), one can note that it gives the worst MAE for all benchmark sets (MAE = 3.32 kcal/mol). This finding suggest that the relative good performance of DH functionals rely on the mutual error cancellation effect between non-local HF and MP2 and semilocal DFA counterparts[148] (see also Section ).

As to the RSX-QIDH DFA performance, it gives an almost identical behavior as PBE0-DH. They both perform similarly for BH and NCI benchmark sets. One can note that RSX-QIDH gives much better performance only for MGT where it yields a MAE of 3.87 kcal/mol, what is still larger in error than the results given by LRC-ωPBE-DH (m = 2) DFA. It is also worth mentioning that RSX-QIDH improves over PBE-QIDH only for PA26 and SIE4×4 subsets while its performance remains unsatisfactory for all other benchmark sets. This is also the case for vertical ionization potentials (VIP) calculated as an energy difference (see G21IP results and discussion in Section ).

What is worth noting that, for the SIE4×4 test set, LRC-ωPBE-DH (m = 4) gives twice worse results (MAE = 1.88 kcal/mol) than the PBE-mIDH (m = 4) counterpart, which gives here MAE = 0.99 kcal/mol. The identical trend is observed for the m = 3 case, whereas for m = 2, we observe the improvement of the results (PBE-QIDH gives a MAE of 3.33 kcal/mol, whereas LRC-ωPBE-DH gives MAE = 2.67 kcal/mol). This behavior can be justified from the percentage of the (LR-)HF mixing with the semilocal DFAs. Note that the SIE4×4 test set consists of H2+, He+, NH3+, and H2O+ molecules (at several points along a dissociation curve of each of the complexes) important to study the problem related to (many-electron) self-interacting error ((ME)SIE).[151−156] The PBE-mIDH family of functionals mixes HF more upon increasing m values and shows less many-electron self-interaction error. In other sense, more HF mixing with DFA eliminates the delocalization error of DFAs as the HF is more localized. Similar trends are observed for the LRC-ωPBE-DH family of functionals. For MP2, we observe overestimation in the case of H2O+, indicating that the error cancellation from DFAs and HF mixing is important for H2O+. However, the error of RSX-QIDH is slightly larger than that of LRC-ωPBE-DH (m = 4).

Here, a bit of analysis of the screened functionals in terms of the dependence on both m and ω is required from the standpoint of performances of the LRC-ωPBE-DH family of functionals. This can be somehow explained by analyzing Figure and Table . From Figure , one can note that, with increasing ω values, we observe the significant worsening of the results especially visible for the BH6 set. The full MP2 method results can be recovered by going with m → ∞ and ω → ∞. From eq , in turn, we see that, when ω → 0, we should recover most PBE-mIDH family of functional results. From Table , it is also noticed that the main differences in the performance of SIE4 × 4 comes from the different mixtures of the HF percentage, where LRC-ωPBE-DH (m = 4) gives the least error within the LRC-ωPBE-DH (m) family of functionals. Overall, LRC-ωPBE-DH (m = 3) performs better within this family of methods, which may come from the balanced treatment of the HF exchange and MP2 correlation with semilocal parts of the functional.

Dissociation Curves of H2+, He2+, and ArKr+

In Figure , we show the dissociation curves of three representative test cases where the errors of functionals are dominated by the many-electron self-interaction problem.

Figure 3

Dissociation curves of H2+ (upper panel), He2+ (middle panel), and ArKr+ (lower panel).

The dissociation of the H2+ molecule (upper panel) can be considered as a paradigm in quantum chemistry,[157] describing the delocalization errors of XC functionals. This system is composed by one electron, such that E = 0, and the exchange energy compensates the Hartree energy.[158] We observe that PBE and PBE0 fail badly, showing huge delocalization errors at large internuclear distances, due to their XC hole densities.[159] The LC-ωPBE improves over PBE and PBE0 but is still considerably worse than PBE-mIDH (m = 4) and especially LRC-ωPBE-DH (m = 4), which is very close to the exact HF curve. We recall that the u-metaGGA exchange functional, which uses as an additional ingredient the Hartree potential, is exact for this system.[158]

Next, we consider the dissociation curves of rare gas cation dimers,[156] the He2+ (middle panel), and the asymmetric ArKr+ molecule (lower panel). In both cases, the PBE and PBE0 are showing delocalization errors. On the other hand, PBE-mIDH (m = 4) and LRC-ωPBE-DH (m = 4) are very accurate, being close to the CCSD(T) reference, at any internuclear distance R. Moreover, for the He2+ molecule, the HF and even MP2 show localization errors at large R,[157] while LRC-ωPBE-DH (m = 4) is almost exact. We also note that for H2+ and He2+ molecules LRC-ωPBE-DH (m = 4) improve significantly over the RSX-QIDH functional (reported in Fig. 1 in ref (140)).

These results, corroborating with the ones for the SIE4×4 test shown in Table , prove that PBE-mIDH (m = 4) and LRC-ωPBE-DH (m = 4) are good candidates for real and difficult applications of quantum chemistry.

Correlation Potentials and Densities

In Figure and Figure , we report the correlation potentials of PBE-mIDH and LRC-ωPBE-DH families of functionals, respectively, for two representative cases, namely, Ne atom and CO molecule.

Figure 4

Correlation potentials of the Ne (uncontracted ROOS-ATZP[160]) atom and CO molecule (uncontracted cc-pVTZ,[149]R0 = 1.128Å) for various PBE-mIDH methods (m = 2 (PBE-QIDX), m = 3 (PBE0-2), and m = 4) obtained using a one-step scheme via the optimized effective potential (OEP) method from refs (161−163). The reference KS[CCSD(T)] has been obtained using a method from ref (164) using the same computational setup.

Figure 5

Similar as in Figure , but for the LRC-ωPBE-DH (m = 2,3,4) family of functionals.

One can note that all DH DFAs reproduce quite reasonable the physical features of reference correlation potential obtained from the self-consistent ab initio OEP2-sc[148,165] method and the one reconstructed from CCSD(T) relaxed density (KS[CCSD(T)]) via a WY inverse approach.[164]

The differences mostly occur in the core (due to Laplacian of the density[166,167] term in GGA potentials[168]) and asymptotic regions (due to the differences in the asymptotic behavior between OEPx (−1/r) and semilocal DFAs[169]). The oscillation in both regions are significantly diminished for LRC-ωPBE-DH potentials, which are quite similar to those obtained with OEP2-SOS[170] or interaction-strength interpolation methods.[171] In the asymptotic region, this is due to the incorporation of a proper long-range asymptotic (−1/r) behavior in the XC potential via LR part of the EXX term. This feature has a visible impact on the quality of VIP potentials obtained directly from the HOMO energy (see Section for more details). In the core region, in turn, the LR EXX potential introduces small but non-vanishing contribution,[101,163] which, combined with mutual error cancellation between exchange and correlation part of DFA, leads to a reduction of nonphysical oscillations.

As to the accuracy of methods, one can note that with increasing m the DFAs seem to better reproduce the reference CCSD(T) correlation potentials. Note also that both the PBE0-DH and RSX-QIDH v are quite similar in shape with the m = 2 variants differing only in the asymptotic region. For the m = ∞ case, this is the MP2 method, we should obtain the potential similar in shape to the OEP2-sc method,[172] which slightly overestimates the reference one. The opposite trend is observed when the calculations are performed in the full KS framework[162,163] (see Figures S1 and S3 in the Supporting Information[147]). The fully local potential affect also the unoccupied KS orbitals what in general have a large impact on the functional performance (see the discussion in this context in Sec 4.1. of ref (163)). The increasing percentage of EXX and MP2 leads to increasing discrepancies between the OEP DH v and the one obtained from CCSD(T) density. In this case, for m = ∞, we will recover the OEP-GL2 method, which, as it is generally known,[165,170,172−176] significantly overestimates correlation effects. This indicates that the full KS realization of DH functionals might lead to a quite different behavior (trends) for systems from Table .

In Figures and , we report the correlation densities, which support the correlation potential analysis. Also here, PBE0-DH DFA gives almost the same behavior as PBE-QIDH (m = 2) except the core region where it exhibits larger overestimation. This can be simply explained by the larger percentage of exchange DFA mixed in the PBE0-DH functional. Almost the same behavior is observed for the RSX-QIDH density (in the case of RS hybrids), which agrees quite nicely with LRC-ωPBE-DH (m = 2) data. With increasing m, we observe an improvement of densities also in the core region where nonphysical oscillations (inherited from GGA terms[176]) start to be less pronounced. Thus, for m = 4, the agreement with the reference CCSD(T) densities is much better. This analysis somehow also confirms the findings from ref (177) where it was shown that the accuracy of the PBE-based DH densities increases with the amount of EXX and MP2 terms included in the functional (PBE0-DH > PBE-QIDH > PBE0-2 > PBE-mIDH (m = 4)). As previously, the full self-consisted KS densities (see Figs. S2 and S4 in ref (147)) lead to an opposite trend, and thus the better agreement with reference is obtained for m = 2. The larger the m, the bigger contribution comes from the GL2 term, which leads to an overestimation of correlation effects. This might indicate that DH functionals can benefit from the utilization of orbitals obtained within the GKS framework.

Figure 6

Correlation densities, Δρ = ρ – ρEXX, of the Ne atom (uncontracted ROOS-ATZP[160]) and CO molecule (uncontracted cc-pVTZ,[149]R0 = 1.128Å) for various PBE-mIDH methods (m = 2 (PBE-QIDX), m = 3 (PBE0-2), and m = 4). The reference CCSD(T) correlation density was calculated with respect to HF density as Δρ = ρCCSD(T) – ρHF using the same computational setup.

Figure 7

Similar as in Figure , but for the LRC-ωPBE-DH (m = 2,3,4) family of functionals.

Ionization Potentials

In Figure , we report the mean absolute relative errors (MARE) for VIP calculated in two manners: (i) using total energies of neutral and ionic species VIP = E(N) – E(N – 1); (ii) as VIP = – εHOMO what is directly related with the quality of the XC potential[148,169−171,176,178] (see Section for more details). We note that, for exact DFT calculations, the two computed quantities should be identical.

Figure 8

Mean absolute relative errors (MAREs) calculated with respect to CCSD(T), computed with different methods for the vertical ionization potentials of several atoms and molecules. The systems are listed in Table S1 of ref (148). The MP2 and OEP2-sc results are taken from ref (148).

Let us first turn our attention to the results obtained as energy differences. In this case, all the methods perform quite similarly, giving MAREs in the range of 1.19%–2.77%. The best results are obtained by the PBE-QIDH (m = 2) functional with a MARE of 1.19%, closely followed by the OEP2-sc method (MARE = 1.27%). Surprisingly, the RSX-QIDH functional gives here slightly worse performance with a MARE of 1.62% (see also the G21IP results in Table. ). The worst VIP is provided by PBE0-DH (MARE of 2.77%) DFA. With increasing of m parameter, we also observe an increase in MARE for both families of functionals, which for m = ∞ should reach the MP2 results, which display a MARE of 2.1%.

Now let us focus on the VIPs obtained from the opposite of the energy of HOMO. Here, the worst performance is provided by PBE0-DH (MARE = 13.31%) followed by PBE-QIDH, which yields a MARE of 9.13%. These two functionals actually exhibit the worse correlation (and thus XC) potentials reported in Section . As previously, the increase in EXX and MP2 terms leads to the improvement of VIP. For PBE-mIDH (m = 4), we get a MARE of 5.6%, which is about 0.9–1.0% off the MP2 (m = ∞) and OEP2-sc results. These predictions are still quite bad with respect to other second-order approaches such as the OEP2-SOSc[179] method, which for the same VIP set yields a MARE of 1.03% (similar to the state-of-art IP-EOM-CCSD[180,181] method, which gives here a MARE of 1.07%).

In the case of RS DH functionals, we observe a drastic improvement in the quality of HOMO energies. Already for m = 2, the inclusion of the LR EXX term reduced the MARE to 3.18%. Surprisingly, for m = 4, we get a MARE of 1.84%, which is identical with MARE calculated for VIP from energy differences (MAEs in this cases are 0.26 and 0.23 eV for −εHOMO and E(N) – E(N-1), respectively). This indicates the good quality of both the XC potentials and total energies obtained with these DFAs. Similar to LRC-ωPBE-DH (m = 4) DFA, the RSX-QIDH HOMO energies agree quite nicely with VIP calculated as energy differences, giving here a MARE of 1.86%. This can be justified by the optimization/proper selection of the ω parameter in RSX-QIDH DFA, which, as recently shown,[143] is almost equivalent to optical tuning.[142]

Conclusions

We have introduced a generalized theory of the DH density functional using the AC formalism and higher-order many-body perturbation theory. The flavors of the present construction are quite physical and allow one to construct mH functionals of any order. We have shown that the present generalization recovers already existing DFA, namely, the PBE-QIDH and PBE0-2 functionals for certain choices of the parameter m.

Based on the here developed theory, generalized PBE-mIDH and long-range corrected DH functionals have been investigated. These functionals are applied to well-known molecular test cases, and their performances showed an improvement in many cases compared to other DH functionals. In particular, it has been shown that the constructed PBE-mIDH and long-range DH (mixing a larger amount of HF and MP2 terms with increasing value of m parameter) may offer a guiding tool in other areas where mitigating self-interaction error is important. The constructed functionals may also be attractive from the standpoint of application in time-dependent density functional theory.

As the final reflection, the inclusion of HOT of the many-body perturbation theory such as GL3/MP3 and GL4/MP4 leads in general to the improvement of the results, seen in the here investigated total, atomization, and binding energies. However, the substantial computational cost does not allow us to construct a usable functional with the higher perturbation terms such as triple and quadruple hybrids. Those can be considered in the future development of functional construction.

Computational Details

For molecular properties we consider selection of well-defined test sets from the Minnesota 2.0 and GMTKN55.[182] The test sets are divided into three subsets, namely, main group thermochemistry (MGT), barrier heights (BH), and non-covalent interactions (NCI). For the thermochemistry test set, we consider the following: (1) AE6, atomization energies of 6 molecules;[144,146] (2) G2/148, atomization energies of 148 molecules;[183] (3) G21EA, 25 electron affinities;[182,184] (4) G21IP, 36 ionization potentials;[182,184] (5) PA26, 26 proton affinities;[185−187] (6) BH76RC, 30 reaction energies of the BH76 test;[182] and (7) SIE4×4, 16 single-point self-interaction correction problems.[182] For barrier heights, we consider the following: (1) BH6, 6 barrier heights;[144] (2) HTBH38, 38 hydrogen barrier heights;[144] and (3) NHTBH38, 38 non-hydrogen barrier heights.[188] For the non-covalent interaction test set, we consider the following: (1) HB6, 6 hydrogen bond test set;[189] (2) DI6, 6 dipole interactions;[189] (3) CT7, 7 charge transfer molecules;[189] (4) PPS5, 5 π–π stacking molecules; (5) WI7, 7 weak interaction complexes, bound by dispersion-like interactions;[189] and (6) S22, 22 non-covalent interaction molecules including H-bond, dispersion interactions, and mixed bonds.[190,191] We underline that the here considered test sets represent the overall important chemical properties. Further, we also consider the dissociation energies of H2+, He2+, and ArKr+ to assess the functional performance. The calculations of all the functionals are performed using the developed version of the Q-CHEM code[192] except the RSX-QIDH functional calculations, which are done in the NWChem package.[193] The def2-QZVP basis set is used for our present calculations, except for the G21EA test set, for which the def2-QZVPD basis set is used. All quantities in this part have been calculated without counterpoise (CP) corrections for the basis set superposition error (BSSE).

Additionally, to assess the accuracy of PBE-mIDH and LRC-ωPBE-DH (m = 2,3,4) family of functionals, we have considered several test cases:

in Section , we investigate the impact of higher-order terms (HOT) in eqs and onto the quality of results. In this respect, we define the triple-hybrid (TH) (m = 3, including GL2 and GL3 correlation contributions) and quadruple-hybrid (QH) (m = 4, including GL2, GL3, and GL4 correlation contributions) functionals. In the assessment, we have examined the total energies of several atoms and molecules,[150] interaction energies of a small set of non-covalently bonded molecules from[150] and atomization energies of AE6[144] benchmark set. This is mostly dictated by the substantial computational cost of the inclusion of GL3 and GL4 terms.

We remark that, because for the HF reference state, the Taylor series around λ = 0 has almost the same small-λ expansion[134] expressed throughout MPk terms; usually in the construction of DH DFAs, one utilized an MP2-like correlation energy expression (the same expression as for the HF theory, but with KS reference orbitals and orbital energies). This is further justified by the small contribution coming from the single excited term (i.e., the last term of eq ) in the GL2 energy expression.[170] Higher-order GLk energy contribution differ much more from MPk ones (evaluated on the same KS reference state), mostly due to non-vanishing off-diagonal contributions in H0 (see ref (181) and ref (133) for the diagrammatic representation of third- and fourth-order energy terms). Nevertheless, the dominant energy contribution comes from the MPk-like formula. This is especially valid in our generalized formalism where the HF EXX energy enters the XC expression with at least 66% (HF dominant) contribution. Taking into account that the DH is generally treated within the GKS framework, we expect that the difference between GLn and MPn terms should be rather small. Thus, in the current study, we approximate GL2 and higher-order GLk terms by their MPk counterparts, namely, GL2 ≈ MP2, GL3 ≈ MP3, GL4 ≈ MP4, etc.

The calculations in Section have been carried out within the GKS framework using the NWChem package.[193] All interaction energies have been calculated with counterpoise (CP) corrections to remove the basis set superposition error (BSSE).

using the formalism developed in refs (162). and (163) we have computed the correlation potentials for the Ne atom and CO molecule (see Section ). Due to the explicit orbital dependence of EXX and GL2 energy terms, the corresponding correlation potential (v(r) = δEc [ρ]/δρ(r)) was computed using the optimized effective potential (OEP) formalism in a post-SCF fashion for fixed reference orbitals, orbital energies, and densities coming from generalized KS (GKS)[194,195] calculations in a one-step procedure.[161] Because eqs and , mix semilocal DFAs with WFT energy expressions, in the following, we define the DH correlation functionals as

We note that a similar approach was already successfully used in several studies[161,171,178,196−198] to investigate the most relevant features of functional derivatives. As shown in ref (161), the behavior of correlation potential depends on the reference orbitals utilized in a one-step procedure. Therefore, the correlation potentials might be considered as a first indicator of the quality of GKS orbitals and densities. Additionally, for comparison, in ref (147), we report the full self-consisted OEP correlation potentials obtained using the formalism developed in ref (162) and taking into account all functional derivative terms. Likewise, in our previous studies,[148,168,170,171,179] in order to solve the OEP equations, we have employed the finite-basis set procedure from refs (199) and (200).

In the same section, we report correlation densities[170,172,176,201,202] defined aswhere ρxref is the density obtained for the exact exchange only (OEPx)[200,203] (xref = EXX) and HF (xref = HF) calculations for DFT and conventional WFT methods, respectively. We underline, however, that for all investigated cases, the difference between HF and EXX densities is very small (see, e.g., Fig. S4 in ref (169)); thus, the mutual comparison is consistent. In the case of DH DFAs, the ρ was obtained as follows. As one knows in the DH GKS scheme, the KS equations are solved with disregarded v contribution. This means that the self-consistent density does not take into account the orbital relaxation due to this missing correlation term. However, using the converged GKS orbitals, we can include the missing contribution describing this effect. The re-scaled MP2 part[204] can be obtained in a post-SCF fashion from the relaxed MP2 density matrices[205−207] constructed using the Lagrangian approach.[208−210] In principle, the corrected DH density (ρ) should allow to generate via one of the inverse methods (e.g., the WY method[164]) the same correlation potentials as the one obtained via the one-step procedure. We note that a similar strategy is adopted within the orbital-optimized DH framework[211,212] or EKT-DH[213] method to include the effect of missing second-order contribution. As previously, the full self-consistent OEP results are reported in ref (147).

To compare our results, we have utilized the CCSD(T) correlation densities obtained from the relaxed density matrices[205−207] constructed using the Lagrangian approach.[208−210] The reference correlation potential have been obtained using the method from ref (164) taking as a starting point the relaxed CCSD(T) densities. All calculations have been carried out with a locally modified version of the ACES II[214] program. The Ne and CO OEP calculations have been performed in fully uncontracted ROOS-ATZP[160] and cc-pVTZ[149] basis sets, respectively. For more technical details, we refer the reader to refs (170) and (172).

as in ref (148), we have obtained the vertical ionization potentials (VIP)[179] for several atoms and molecules computed in two manners: (i) as the energy difference between the neutral and the ionic species[215] (VIP = E(N) – E(N – 1)); (ii) the opposite of energy of the highest occupied molecular orbital (HOMO) (VIP = – εHOMO). In the case of DH, the HOMO energies are calculated as in refs[148], (162), and (163), where the second-order self-energy correction[133] is calculated in a post-SCF fashion on top of GKS DH orbitals (likewise for HF orbitals[215,216]). Like in ref (148), we did not take into account the additional orbital relaxation term due to the second-order correlation; however, as pointed out in ref (217), this contribution should be rather small. The computational setup, namely, basis sets and molecular geometries are identical as in ref (179) (for more details, see Section ). The calculations have been performed with a locally modified version of PSI4 quantum chemistry packages.[218]