Restoring Size Consistency of Approximate Functionals Constructed from the Adiabatic Connection.

Approximate exchange-correlation functionals built by modeling in a nonlinear way the adiabatic connection (AC) integrand of density functional theory have many attractive features, being virtually parameter-free and satisfying different exact properties, but they also have a fundamental flaw: they violate the size-consistency condition, crucial to evaluate interaction energies of molecular systems. We show that size consistency in the AC-based functionals can be restored in a very simple way at no extra computational cost. Results on a large set of benchmark molecular interaction energies show that functionals based on the interaction strength interpolation approximations are significantly more accurate than second-order perturbation theory.

With applications that stretch from solid-state physics to biochemistry, Kohn–Sham density functional theory (KS-DFT)[1] is presently the most employed electronic structure method. Although the theory is in principle exact, any practical implementation of KS-DFT must rely on approximations for the exchange–correlation (XC) functional, which should capture all of the many-body effects beyond the simple Hartree theory. Despite the existence of hundreds of different XC approximations[2,3] and their widespread success in various disciplines,[2] KS-DFT still encounters open issues, which hamper its overall predictive power[2,4−6] and make the quest for better approximations a crucial research field for computational chemistry, solid-state physics, and materials science.[2,4−6]

The density-fixed adiabatic connection (AC) formalism[7,8] provides an exact expression for the XC energy functional Exc[ρ]where Wλ[ρ] is the AC integrandΨλ[ρ] is the Fermionic wave function with density ρ(r) that minimizes the sum of the kinetic energy T̂ and the electron–electron repulsion V̂ scaled by the coupling constant λ, and U[ρ] is the Hartree energy. For small systems, Wλ[ρ] has also been computed exactly through eq . However, this requires the solution of the many-body Schrödinger equation.[9−11] Thus, for all practical purposes, Wλ[ρ] must be approximated.

Equation has been a fundamental milestone in guiding the construction of approximations. Early AC-based XC functionals used forms that depend linearly on some chosen input ingredients, such as the exchange energy from Hartree–Fock (HF) theory as the value to be recovered at λ = 0, and semilocal approximations at some λ = λp between 0 and 1. These forms are commonly used for the construction of hybrid[12−14] and double-hybrid[15−17] density functionals, resulting in mixing a fixed fraction of HF exchange and second-order perturbation theory with semilocal functionals. They often work well for main-group chemistry, but they show important limitations for various other problems as, for example, the chemistry of transition metals[18] (where they even worsen the results with respect to simpler semilocal functionals), metal–molecule interfaces,[19] and even noncovalent bonding (unless an ad hoc van der Waals correction is used).[20,21] Their main disadvantage is that the mixing fractions are fixed and cannot adapt to different systems or to different parts of a system.

To address this problem, several models in which the input ingredients enter in a nonlinear way have been proposed.[11,22−26] These latter forms do not need to rely on empiricism and can adapt automatically to the peculiarities of the system under study. Along these lines, Ernzerhof had proposed Padé forms for the λ dependence of the AC integrand,[22] which later were used for construction of the MCY family of functionals that are constrained to be free of one-electron self-interaction error.[25,27] Another example of models that use input ingredients in a nonlinear way is provided by the interaction strength interpolation (ISI) functionals, which depend explicitly on the weak- and strong-coupling ingredients,[23,24,28−30] essentially extending to nonuniform densities Wigner’s[31,32] idea for approximating the energy of the homogeneous electron gas. Despite the advantages of the nonlinear forms over the linear ones, the former encounter a fundamental flaw: the XC functionals that are constructed from them are not size-consistent for systems composed of different species of fragments.[27,33] This depends on the fact that these methods employ as input ingredients global quantities (integrated over all space). A route that is currently being explored addresses this issue by modeling the AC at each given spatial position r, using energy densities wλ(r)[26,34−36] instead of quantities integrated over all space. This strategy is very promising but does not allow using in a straightforward way semilocal ingredients[26,33−39] because of the inherent ambiguity in the definition of energy densities, a problem shared with the construction of local hybrid functionals.[37,40,41]

In this Letter, we show that size consistency of the global (integrated over all space) AC forms in which the ingredients enter in a nonlinear way can be restored in a remarkably simple way, making it possible to obtain meaningful interaction energies at no additional computational cost.

Consider a system M (e.g., a molecular complex) composed of a set of fragments A, with i = 1, ..., N. The interaction energy is a key quantity in chemistry, and it is defined aswhere E(M) is the energy of the bound system M and E(A) are the energies of the individual fragments. If we now compute the energy of a system M* made of the same fragments A placed at a very large (infinite) distance from each other, any size-consistent method should give Eint(M*) = 0, or equivalently

We should stress at this point that size consistency in DFT is in general a very subtle issue, particularly when dealing with fragments with a degenerate ground state (e.g., open-shell atoms) as the (spin) density is no longer an intensive quantity.[42,43] To disentangle this more general DFT problem from the one of size consistency of the nonlinear AC models, here we focus on the cases where the fragments have a nondegenerate ground state,[44] considering noncovalent interactions.

The idea behind AC-based functionals is to use a certain number of input ingredients W[ρ], constructing a λ-dependent function that interpolates between them. For example, many standard hybrid functionals model Wλ[ρ] with functions of the kindwhere WλDFA[ρ] is a given density functional approximation (usually a semilocal functional), with its exchange component Wλ=0DFA[ρ] = ExDFA[ρ], and ExHF[ρ] is the HF exchange energy.

These kinds of expressions, when inserted into eq , yield a fixed fraction 1/n of HF exchange energy mixed with a semilocal density functional approximation. Because the input ingredients, in this case, WλDFA[ρ], ExHF[ρ], and ExDFA[ρ], enter linearly in the model of eq , the resulting XC functional automatically satisfies the size-consistency condition of eq if the individual ingredients do.

As examples of approximations in which the ingredients enter in a nonlinear way, consider first the Padé([1/1]) form introduced by Ernzerhof[22]with a[ρ] = W0[ρ] = ExHF[ρ], b[ρ] = W0′[ρ] (which can be obtained from second-order perturbation theory), and c[ρ] = λp–1 – W0′[ρ](W0[ρ] – WλDFA[ρ])−1, where WλDFA[ρ] could be a semilocal functional at a chosen value λp. We see immediately that in this case, even if the input quantities W0, W0′, and Wλ satisfy the size-consistency condition of eq , the resulting XC functional from eq does not, because it is given by a nonlinear function fPad of these ingredients, ExcPad = fPad(W0, W0′, Wλ).

Another example is the idea of Seidl and co-workers[23,24] to build approximate Wλ[ρ] by interpolating between its weak (λ → 0) and strong (λ → ∞) coupling expansionswhich allows one to avoid bias toward the weakly correlated regime and to include more pieces of exact information. The λ → 0 limit of eq is provided by the exact exchange and the second-order perturbation theory, while the functionals W∞[ρ] and W∞′[ρ] describe a floating Wigner crystal with a metric dictated by the density.[28,45]

Different formulas that interpolate between the limits of eqs and 8 are available in the literature.[23,24,26,28,46] As in the Padé example of eq , when these interpolation formulas are inserted into eq they give an XC energy that is a nonlinear function of the four ingredients (or a subset thereof), W0[ρ], W0′[ρ], W∞′[ρ], and W∞[ρ], appearing in eqs and 8.

It is clear from these examples that we can write a general XC functional obtained by modeling the AC aswhere fACM is a nonlinear function that results from the integration (via eq ) of the given adiabatic connection model (ACM) and W[ρ] = {W1[ρ], ..., W[ρ]} is a compact notation for the k input ingredients that have been used. Then we haveandThis equation is one of the main points in this work. Although conceptually simple, it shows that the energy of a set of infinitely distant fragments (M*) can be expressed as a function of the quantities of the isolated fragments. Notice that this holds in this special case because fACM is a function of global size-consistent[44] quantities. For size-inconsistent wave function methods, such as CISD, this is usually not true, and the energy of M* needs to be computed by performing an extra calculation with the fragments at a large distance, which might be tricky to do in practice.

Essentially, all of the models that have been proposed in the literature[27] satisfy the conditionmeaning that they are size-consistent when a system dissociates into equal fragments (size-extensivity). This is also a key difference with the size-consistency problem of wave function methods, which also arises in the case of equal fragments. However, when the A are of different species, eqs and 11 give in general different results, and attempts to make them equal for a nonlinear model have failed so far.[27]

As said, evaluating eq or 11 has exactly the same computational cost as both equations only need the input ingredients for the individual fragments. The idea behind the size-consistency correction (SCC) is thus extremely simple, and it is related to discussions reported in refs (47−49); it consists of using the difference between eqs and 10 to cancel the size-consistency error that is made when evaluating interaction energies from eq Note that this correction is fundamentally different from a direct calculation of ∑ E(A) – E(M*) because, due to the use of eq , only knowledge of the isolated fragments is required here, while there is no need to deal with the (possibly tricky) calculation of the supramolecular energy M*.

Adding ΔSCC to an interaction energy computed via eq is equivalent to always evaluating interaction energies with respect to eq instead of eq , i.e.

As an example of the performance of the SCC, we examine here ACMs that link the two limits of eqs and 8. As said, we focus on noncovalent interactions because the fragments A have a nondegenerate ground state, which should guarantee size consistency of the input ingredients.[42,43] Moreover, in this case, the interaction energy is small, and therefore, the correction can be relevant; for covalent interactions, in fact, the correction is of the same order of magnitude as that for noncovalent ones, but the interaction energy is at least 2 orders of magnitude larger. All calculations have been performed using a development version of the TURBOMOLE package,[51,52] with computational details similar to those of refs (48) and (53), in which the ISI-like functionals are evaluated on HF orbitals (see section IV.D of ref (53) for some discussion of this choice). Thus, in eq , W0[ρ] = ExHF[ρ] and W0′[ρ] is twice the second-order Möller–Plesset (MP2) correlation energy, W0′[ρ] = 2EcMP2. The strong-coupling functionals W∞[ρ] and W∞′[ρ] of eq are approximated with the point charge plus continuum (PC) semilocal model,[54] which is reasonably accurate.[28,45] We test different interpolation formulas that have been proposed in the literature, namely, SPL,[23] rev-ISI,[28] and LB.[46] Additionally, we also tested the Padé[1,1] formula of eq by using λp = ∞. The interpolation formulas and additional computational details are reported in the Supporting Information.

As a first example, in Figure , we show the absolute errors in the interaction energy for a set of dispersion complexes made of fragments of different species obtained from the rev-ISI and the Padé interpolation formulas, computed with and without the SCC. From this figure, we can see that in both cases the error is reduced by an order of magnitude when the correction is applied, i.e, when eq is used.

Figure 1

Absolute errors in the interaction energy (kcal/mol, log-scale) for a set of dispersion heterodimers containing noble gas atoms obtained with the rev-ISI and Padé([1/1]) functionals with and without inclusion of the SCC of eq (1: He–Ne; 2: He–Ar; 3: Ne–Ar; 4: Ar–Kr; 5: CH4–Ne; 6: C6H6–Ne).

In Figure , we also report the interaction energy curve for He–Ne obtained from the rev-ISI functional. We see that the rev-ISI curve has a very reasonable shape, but because of the size-consistency error when computed with eq , it goes to a positive value with respect to the sum of the fragment energies. Instead, when the SCC is applied, the correct asymptotic value of the dissociation curve (given by eq ) is used to compute interaction energies. Very similar figures are obtained when we consider other interpolation formulas and other systems, with the overall shift that is sometimes positive and sometimes negative.

Figure 2

Interaction energy curve for the HeNe heterodimer obtained with the rev-ISI functional with and without the SCC of eq . The MP2 curve is shown for comparison, and the accurate curve has been taken from ref (50).

Finally, we use the SCC to assess the accuracy of AC-based functionals for more noncovalent complexes relevant for chemistry and biology. For this purpose, we employ the well-established quantum-chemical data set for noncovalent interactions S66.[55]

In Figure , we report the values ofwhere the index i = 1, ..., 66 labels the various complexes, Eref is the reference interaction energy of the ith complex and Erev-ISI-SCC and Erev-ISI are the corresponding interaction energies calculated with rev-ISI with and without the SCC correction, respectively. Thus, a negative (positive) δSCC means that the SCC reduces (increases) the relative error of the interaction energies. Notice that in Figure the S66 complexes are sorted in ascending order according to the computed |ΔSCC| value (see inset of Figure ). Thus, one can see that for systems where |ΔSCC| is non-negligible (i.e., i ≳ 30) the inclusion of SCC brings an improvement to the results (δSCC < 0) and that the improvement can be as large as 10%. On the other hand, there are some systems (i.e., i ≲ 30) with a negligible SCC. This is not surprising as the S66 data set contains 17 homodimers for which the AC-based functionals are already size-consistent. Moreover, there is another case in which the size-consistency error becomes negligible: when the ratio q = W(A)/W(B) between the ith input ingredient of fragment A and of fragment B is roughly the same for all i, q ≈ q, a case that becomes mathematically equivalent to eq . In summary, Figure shows that the larger the ΔSCC value, the larger the reduction of the errors. This indicates that the inclusion of ΔSCC is significant and works correctly for most noncovalent complexes having different constituent units.

Figure 3

Difference δSCC of the absolute relative errors of interaction energies calculated with and without SCC for the complexes of the S66 test set (sorted with increasing |ΔSCC|; see the inset).

More generally, for all of the SCC-ISI-like functionals that we examined, the performance for noncovalent interactions is quite good, being comparable or better than state-of-the-art computational approaches (see Table ). Especially for dispersion and mixed complexes, all of the ACMs perform 7 or 8 times better than the B2PLYP double hybrid and twice better than MP2 (note that both of these methods have the same computational cost as the SCC-ISI-like functionals). This is quite relevant, considering that ISI-like functionals have not been explicitly constructed to model interaction energies and do not employ any empirical parameter (in contrast, e.g., to the approaches in the last three lines of Table ). Notice also that for AC-based functionals not only is the mean absolute error low but also the variance (last column of Table ). Therefore, these functionals can describe different types of interactions with similar accuracy.

Table 1

Mean Absolute Error (And Variance, Last Column), in kcal/mol, for the S66 Data Set and Some of Its Subsets for Different AC-Based Functionals Including SCC and Evaluated on HF Density and Orbitalsa

method	H-bonds	dispersion	mixed	total	variance
rev-ISI	0.35	0.44	0.20	0.33	0.08
ISI	0.37	0.42	0.19	0.33	0.09
SPL	0.42	0.42	0.19	0.35	0.11
LB	0.36	0.41	0.19	0.31	0.14
MP2	0.11	0.81	0.45	0.45	0.29
SCS-MI-MP2	0.19	0.45	0.20	0.19	0.10
SCS-CCSD	0.30	0.08	0.08	0.27	0.05
B2PLYP	0.72	2.79	1.63	1.71	1.26

For all results, see Table S2 in the Supporting Information. The last four lines report, for comparison, results from literature.[55,56]

We have shown that XC functionals built by approximating the AC integrand with functions in which the input ingredients enter in a nonlinear way can be made size-consistent at no extra computational cost. The starting idea is that size consistency is restored once we consider fragments that are infinitely far apart, whose energy, by virtue of eq , we compute from the sum of quantities of individual fragments. We focused here on only the case of noncovalent interactions, but the method is generally applicable also to covalent systems. We also remark that, even though in this work we only consider a few ACMs functionals, the SCC based on eq has a more general applicability to any functional built using the AC framework as well as to any functional depending nonlinearly on size-consistent global quantities.

We have shown that our SCC provides in many cases an important correction to the interaction energy and leads to considerable improvement of the accuracy of various ACMs. Thus, it is a simple and efficient way to correct one of the main drawbacks of actual ACMs, which can now be reliably used for different applications. This opens the quest for the development of improved ACMs. A promising route in this direction is the construction of approximations by interpolating energy densities along the AC, which requires nonlocal functionals for the strong-interaction limit[35,57] and/or for the λ = 1 case.[58]