How Much Dispersion Energy Is Included in the Multiconfigurational Interaction Energy?

We demonstrate how to quantify the amount of dispersion interaction recovered by supermolecular calculations with the multiconfigurational self-consistent field (MCSCF) wave functions. For this purpose, we present a rigorous derivation which connects the portion of dispersion interaction captured by the assumed wave function model—the residual dispersion interaction—with the size of the active space. Based on the obtained expression for the residual dispersion contribution, we propose a dispersion correction for the MCSCF that avoids correlation double counting. Numerical demonstration for model four-electron dimers in both ground and excited states described with the complete active space self-consistent field (CASSCF) reference serves as a proof-of-concept for the method. Accurate results, largely independent of the size of the active space, are obtained. For many-electron systems, routine CASSCF interaction energy calculations recover a tiny fraction of the full second-order dispersion energy. We found that the residual dispersion is non-negligible only for purely dispersion-bound complexes.


INTRODUCTION
Calculations of the intermolecular interaction energy typically involve the supermolecular approach. The perturbation theory is also a common choice and is often used to expose the nature of the interaction. In hybrid approaches, the supermolecular results are supplemented with energy contributions obtained in a perturbative manner, the most popular example being density functional theory (DFT) corrected for the dispersion energy. 1, 2 A given wave function theory (WFT) approximation is useful in describing interactions either if it accounts for highorder components in the intermolecular interaction operator or when it is reasonably accurate in the lowest orders and it is known which higher-order components are missing. A representative example from the first group are the coupled cluster (CC) methods 3,4 which predict interaction energies with high accuracy when effects of single, double, and triple excitations are accounted for. The Hartree−Fock (HF) method belongs to the second group. The HF energy contains electrostatic, exchange, and induction interaction energy components and has been rigorously shown to miss the second-order dispersion component. 5 In spite of this deficiency, the supermolecular HF result can be corrected for the dispersion energy which leads to a low-cost hybrid approach. 6−8 The accuracy of the hybrid approach can be significantly improved by including intramonomer correlation effects in a perturbative manner which is known as the HFDc (1) method. 8,9 HF calculations assume a single-determinantal description of the interacting species. Although high-order single-reference CC methods may be applied to multireference systems, 10 the cost of going to full triple or quadruple excitations is too high in practice. The multireference CC extensions 11,12 are expensive, and their applications to intermolecular interactions remain scarce (e.g., see refs 11, 13, 14). Going beyond singledeterminantal methods in the intermolecular interaction theory is a challenge. The multiconfigurational (MC) wave function methods are a good starting point and often the required one for open-shell systems, systems dominated by strong correlation, excited states, and transition-metal complexes. Although including multiple configurations guarantees proper representation of static correlation and partial recovery of the intramonomer correlation effects, 15,16 it is not a viable way to restore the remaining part of the dynamic correlation which in the case of intermolecular interactions manifests as dispersion. Notwithstanding significant progress in this field, none of the existing MC methods have yet achieved the levels of both accuracy and efficiency comparable to what is now within reach of black-box interaction-energy calculations for single-reference systems. Therefore, the search for new and improved MC approaches remains active.
The most popular strategies to capture dispersion forces in the MC framework are the configuration interaction (CI) expansion and the perturbation theory. Each of these approaches involves certain limitations both in terms of computational cost and in terms of accuracy. A significant problem of CI is the lack of size consistency. Although it can be ameliorated with approximate corrections, 17,18 the residual size inconsistency typically still has to be addressed. This is done by obtaining the asymptotic interaction energy limit in a dimer calculation with the monomers separated at a large distance and shifting the results with respect to this value. 19 Note that in this scheme, the counterpoise correction (CP) 20 cannot be applied and the basis set superposition error (BSSE) is not removed. Additionally, because of the lack of triple excitations, multireference CI may underestimate the dispersion energy in van der Waals complexes by as much as 30%. 21,22 A common choice to account for dynamic correlation via the multireference perturbation theory are the complete active space (CAS) perturbation theory (CASPT2) 23,24 and nelectron valence state perturbation theory (NEVPT2) 25−27 methods. Truncated at the second order, both require access to four-electron density matrices, which quickly becomes a computational bottleneck with respect to the size of the active space. Moreover, CASPT2 and multireference formulations of the Møller−Plesset theory 28,29 are plagued with intruder states 30 which have to be removed to guarantee a smooth potential energy surface. Various techniques that eliminate intruders exist, yet one should be cautious when choosing the optimal value of the shift parameter. 31 The intruder-state problem is by design avoided in the NEVPT2 method. The accuracy of single-reference MP2 calculations is often insufficient for quantitative predictions. A major weakness is that MP2 restores the dispersion energy at the uncoupled level of theory 32−34 which leads, for example, to the well-known overestimation of interaction energy in π−π complexes 35,36 or underestimation in the case of alkali-metal dimers. 37 This implies that the accuracy of the multireference second-order perturbation theory is seriously limited.
Several variants of the multiconfiguration DFT (MC DFT) 38,39 capable of treating dispersion forces have also been developed. In MC DFT, the connection between the DFT and wave function descriptions of the system is realized by partitioning of the Coulomb interaction operator into the short-and long-range components. The short-range (sr) electron−electron correlation effects are recovered by a density functional approximation, whereas the multiconfiguration wave function model acts in the long range. The usual strategy is to restrict the wave function to a few determinants, which is sufficient to recover the static correlation at a relatively low cost but not enough to capture the long-range dynamic correlation. This translates into the absence of the dispersion contribution in the interaction energy. As a remedy, Fromager and coauthors 40 proposed a long-range adaptation of the NEVPT2 approach. Recently, Hapka et al. 41 introduced a less expensive alternative based on the multireference formulation of the adiabatic connection formalism 42 which relies only on one-and two-electron reduced density matrices of the system. The cost-efficient route in which the complete active space self-consistent field (CASSCF) srPBE functional 43,44 is simply corrected with the semi-classical D3 dispersion correction of Grimme 45 has been advocated by Stein and Reiher. 46 The direct addition of the dispersion energy to the supermolecular CASSCF results, without correcting for the double counting, was attempted several times in the past. In the study of the Hg···Hg dimer, Kunz et al. 15 used this approach to partially account for the dynamic correlation within the monomersan effect that cannot be recovered in the simple "HF plus dispersion" model. Rajchel et al. 47 verified the performance of the "CAS plus dispersion" method against the RCCSD(T) results for the potential energy curves of the Sc···Cr complex in the 8 Σ + , 8 Π, and 8 Δ states. The authors found that "CAS plus dispersion" consequently underestimated the magnitude of the interaction and predicted the location of the minima at ca. 0.2−0.3 bohr longer distances compared to the CC reference. The CAS-based approach was also applied to describe the unusually strong binding in the He···BeO( 1 Σ + ) dimer and showed excellent agreement with the MRCISD + Q potential. 48 In all of the mentioned examples, it was argued that the compact representation of the active space selected in CASSCF calculations is too small to recover the dispersion energy. Nevertheless, to the best of our knowledge, no systematic study of this claim has been presented to date and expressions for the residual dispersion energy recovered in the supermolecular CASSCF have not been derived.
In this work, we present how to rigorously extract the dispersion term in the interaction energy calculated with MC wave functions so that the double counting of the intermonomer correlation is avoided. As will be discussed, any wave function dissociating to MC wave functions for both monomers leads to recovering a small portion of the dispersion energy, denoted as E DISP MC . Being able to explicitly calculate the E DISP MC contribution, we propose the "CAS plus dispersion" approach which strictly avoids the double counting of the dispersion contribution and recovers the total interaction energy as where E int MC is the supermolecular interaction energy obtained within the assumed WFT and E DISP (2) is the full second-order dispersion energy. We present the results for wave functions of the CAS type and use the symmetry-adapted perturbation theory (SAPT) 49,50 to obtain the E DISP (2) energy contributions. It is worth noticing that Stein and Reiher 46 raise the concern of dispersion double counting in the context of large active spaces available with the density-matrix renormalization group wave functions. They propose a procedure that could approximately remove the residual dispersion inherently captured by an MC function. This approach requires orbital localization, followed by the removal of configurations (determinants) corresponding to excitations of a dispersionlike character. 51−53 In contrast, our protocol for E disp MC evaluation avoids orbital localization and stems directly from the exact definition of the dispersion energy.
In the next section, we derive an exact expression for the dispersion energy recovered by supermolecular MC WFT calculation and discuss a method for its approximation. In Section 3, we present the details of the implementation and calculations. Section 4 provides numerical demonstration of the proposed "CAS plus dispersion" method on the example of model He···He and He···H 2 dimers as well as representative By definition, the interaction energy vanishes when subsystems are infinitely separated where R AB measures the distance between A and B. Let ĤA B , ĤA, and ĤB denote the dimer Hamiltonian and the Hamiltonians of the isolated monomers, respectively. In the assumed WFT model, computation of the interaction energy in the supermolecular way requires finding wave functions for a dimer Ψ AB and isolated monomers Ψ A and Ψ B , and making use of the pertinent Hamiltonians as follows A primary condition which must be satisfied by a physically relevant WFT is size consistency, that is, vanishing of the supermolecular interaction energy in the infinite-separation limit, which can be written as which in turn imposes the separability condition for the (antisymmetric) wave function Ψ AB reading where the antisymmetrization operatorÂ is applied on the product Ψ A Ψ B of monomer functions. In the Rayleigh−Schrodinger (RS) perturbation theory applied to molecular interactions, a dimer N AB -electron Hamiltonian is written as a sum of N A -and N B -electron Hamiltonians of the monomers (N AB = N A + N B ) and a remainder, which is the interaction operator V̂i nt A sum of ĤA + ĤB is taken as the zero-order Hamiltonian, V̂i nt is taken as a perturbation, and the Ψ A Ψ B product of monomer functions is the zero-order dimer eigenfunction. It is important to notice that the antisymmetrized wave function given in eq 6 is an eigenfunction of ĤA B in the limit R AB → ∞, but it is not an eigenfunction of the zero-order Hamiltonian ĤA + ĤB. As is well known, the interaction energy, eq 4, follows from perturbation theory as a sum of the first-, second-, and higher-order terms in the perturbation 54 where E elst (1) denotes the electrostatic energy, while E ind (2) and E disp (2) denote the second-order induction and dispersion energy contributions, respectively. Each term is well defined and has a clear physical interpretation. In particular, the second-order dispersion energy is given by the expression where i (j) refers to states Ψ A,i (Ψ B,j ) of the energy E A,i (E B,j ) corresponding to ĤA (ĤB) and a prime indicates that summation indices run through all states with the exception of the states of interest Ψ A and Ψ B for which the dispersion is computed. The leading term of the multipole expansion of the dispersion energy of molecular systems vanishes as R AB −6 and dominates in the asymptotic regime of the total interaction energy in van der Waals complexes.
The problem that we address in the next section is whether one can rigorously show to what extent a given wave function model accounts for the dispersion interaction in the supermolecular approach. The positive answer will be provided by singling out from the interaction energy expression terms of a dispersion character, which decay asymptotically in the same (R AB −6 ) fashion as the full dispersion energy, cf. eq 9, and which are equivalent to the latter if WFT corresponds to the full CI (FCI) description of a dimer.
2.2. Dispersion Interaction in the MC Wave Function Description of the Interaction Energy. In this section, we present the expression for the dispersion energy, E disp MC , which is included in the supermolecular interaction energy, obtained with the MC wave function description of a dimer. We begin by considering a general N-electron wave function which can be expressed as the antisymmetrized product of wave functions ∏ Ψ =̂Ψ A P P (10) whereÂ is the antisymmetrization operator (including a normalization constant). Ψ P is an N P -electron antisymmetric wave function, and ∑ P N P = N. In addition, it is assumed that the whole set of spin orbitals is partitioned into disjoint subsets and a given function Ψ P is expanded in a set of Slater determinants built from spin orbitals belonging to a subset P. The assumed ansatz for Ψ is therefore a group product function. 55 We refer to Ψ P as a group function, where P denotes a group, indexes a group function, and refers to a pertinent subset of spin orbitals. Ψ P can be a single determinant or a combination of many determinants, for example, an MC wave function. MC self-consistent field (MCSCF) wave functions, which are in the focus of this work, conform to the ansatz given in eq 10. They can be written as an antisymmetrized product of the inactive part Ψ inact and the active part Ψ active where Ψ inact is a single Slater determinant built of one-electron groups, that is, spin orbitals φ p (12) while Ψ active is a combination of Slater determinants.
If Ψ of the form as in eq 10 is found variationally by minimizing the energy ⟨Ψ|Ĥ|Ψ⟩, where Ĥis the electronic Hamiltonian, then, for each group wave function Ψ P , one can form a group electronic Hamiltonian ĤP readinĝ Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article The intragroup Hamiltonian ĤP ,intra includes the kinetic energy operator, the interaction of N P electrons, assigned to a group P, with all nuclei in the system (denoted as M nuc ), and the appropriate electron−electron interaction (14) The other part of the group Hamiltonian given in eq 13 describes the intergroup mean field interaction of electrons assigned to a group P with electrons in the other groups. It arises as a result of self-consistent field (SCF) optimization and depends on the self-consistently obtained group spin orbitals. The mean field intergroup interaction operator, ĤP ,SCF , is a sum of one-electron Coulomb and exchange operators pertaining to all groups except for P, namely where r 1i = |r 1 − r i |, x = (r,σ) stands for a combined spatial-spin coordinate, and the Pî j operator interchanges respectively, the density and one-electron reduced density matrix functions corresponding to Ψ Q . A wave function Ψ P is one of the eigenfunctions of ĤP with the corresponding energy ε P , that is, if we write the eigenproblem for the group Hamiltonian ĤP then, ∃iΨ P = Ψ P,i . Within the assumed wave function approximate model, Ψ P,i is an MC wave function given as a combination of M P Slater determinants Φ P,K ) and the Hamiltonian ĤP is the diagonal in the M P -dimensional space spanned by Φ P,K functions. It is straightforward to notice that the total electronic energy is given as a sum of group energies, ε P = ⟨Ψ P |ĤP|Ψ P ⟩, with the intergroup Coulomb-exchange interaction energy subtracted to avoid its double counting Now, consider a perturbed electronic Hamiltonian, linear in the perturbation parameter λ A corresponding perturbed group Hamiltonian, cf. eqs 13 and 15, will include a direct perturbing term (linear in λ) and its SCF part will acquire a λ dependence because of the λ dependence of both the density and density matrices, namely The total electronic energy, eq 18, expanded in terms of λ, is a sum of terms in different orders of the perturbation. Because our interest lies in the dispersion energy, which appears in the second order of the interaction operator, we only focus on the second-order term reading where using the second Hellmann−Feynman theorem, 56 the second-order group energy is given as Notice that we have assumed hermiticity of the perturbing operators and that wave functions are real-valued.
Let us now turn to the molecular interaction theory, where the interaction operator V̂i nt , introduced in eq 7, is treated as a perturbation. The zero-order part is given as a sum of monomer Hamiltonianŝ 25) and the explicit form of the perturbation, Ĥ( 1) = V̂i nt , reads where the subscripts a and b label the nuclei and i and j denote the electrons of monomers A and B, respectively. A monomer X is composed of N X electrons and M nuc,X nuclei. Throughout the text, X will stand for either A or B. When monomers are infinitely separated, the size-consistency condition given in eq 5 imposes that the total dimer wave function acquires the form of a product of monomer functions, as seen in eq 6. This can only be satisfied if each N Pelectron wave function Ψ P dissociates into a product of monomer wave functions Ψ P A and Ψ P B , which are, respectively, N P A -and N P B -electron and N P = N P A + N P B . Neglecting the electron exchange, the zero-order group function reads Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article In the R AB → ∞ limit, the group density and the density matrix corresponding to Ψ P separate into monomer group functions, namely which is a consequence of vanishing overlap between Ψ P A and Ψ P B in the dissociation limit. Partitioning of the total electronic Hamiltonian assumed in eqs 25 and 26 results in the following decomposition of the group Hamiltonian ĤP, cf. eq 21 where ĤP A ,intra collects intragroup operators for a group P in monomer A analogously for ĤP B ,intra , and the interaction operator for a group P includes the interaction of electrons in A assigned to a group P with nuclei in monomer B, the interaction of electrons in monomer B assigned to a group P with nuclei in monomer A, and the interaction among electrons in monomers A and B within a group P. At infinite intermonomer separation (the λ = 0 limit), the SCF part of the group Hamiltonian, eq 22, is partitioned aŝ where ρ Q X and γ Q X are the monomer group density and the density matrix, as seen in eqs 28 and 29, respectively. This allows one to write the group Hamiltonian, eq 30, as a sum of monomer group Hamiltonianŝ Our goal is to identify the dispersion terms in the secondorder energy given in eq 23 corresponding to the intermolecular perturbation operator given in eq 26. First, one should notice that the last term in eq 23, which includes second-order perturbations of the Coulomb-exchange inter-group interactions, (E Coul-Exch PQ ) (2) , cannot contribute to the dispersion energy. These terms describe the interactions of the perturbed one-electron density and density matrices. Next, let us consider the second-order group energy ε P (2) , presented in eq 24, where the group interaction operator ĤP (1) = V̂P AB is given in eq 32. The derivatives of ĤP ,SCF (λ) with respect to λ, involved in ε P (2) , are given by the pertinent derivatives of densities and density matrices Thus, they are one-electron operators and as such cannot contribute to the dispersion interaction either. Therefore, the only term in eq 23 which may give rise to the dispersion interaction is the one involving the operator V̂P AB . Retaining solely this term in E (2) and denoting the resulting energy as E d (2) , we write To see under which condition dispersion terms emerge from E d (2) , we need to find the first-order correction for group wave functions, Ψ P (1) . For this purpose, consider zero-order group functions, Ψ P,i (0) , that are eigenfunctions of the ĤP (0) Hamiltonian, eq 35 describing a group P in the noninteracting, R AB → ∞, limit. Let Ψ P X ,i denote monomer functions corresponding to a group monomer Hamiltonian ĤP X given in eq 36 Then, any (N P -electron) ith eigenfunction Ψ P,i (0) is a product of (N P A -and N P B -electron) monomer-group functions with the corresponding energy , , The index i, labeling ε P (0) and Ψ P , is understood as a combined index (kl), where k and l label group states of A and B, respectively.
In the assumed multiconfiguration approximation, a monomer wave function can be written as , , X X X X (43) where {Φ P X ,K } is a set of single determinants built from spin orbitals belonging to a group P X . In the multireference approximations considered in this work, orbital and configurational spaces are finite, and the wave functions are obtained with the optimal orbitals and configuration expansion coefficients. Having defined a group zero-order Hamiltonian, eq 35, with the corresponding eigenfunctions shown in eq 41, applying the standard perturbation theory to eq 16 using the as the perturbation (see eqs 30, 32 and 37) leads to obtaining first-order perturbation to an ith state of the group P, Ψ P,i (1) , which can be written in a general form as , , , The first term results from admitting state mixing and the sum runs over all unperturbed states of the Hamiltonian ĤP (0) orthogonal to Ψ P,i (0) . The second and third terms in eq 44 are due to orbital optimization, that is, allowing for mixing of orbitals between groups. If, for example, φ p (x) is an unperturbed spin orbital belonging to a group P and assigned to monomer A, p ∈ P A , then, its first-order perturbation φ p (1) is expanded in the space of all spin orbitals except those belonging to P A . Such first-order perturbations contribute to Ψ P,i (1) as combinations of singly excited determinants Φ P g , A denotes a determinant built of N P A − 1 orbitals from a group P A and one spin orbital from another group Q ≠ P A (the same for monomer B). Assume from now on that the considered dimer is in its ground state, thus the indices i, l, and k in eq 44 will be set to 0. Notice, however, that the presented analysis is general, not limited to ground states. The explicit expression for C P,0j coefficients in eq 44 reads ,SCF (1) ,0 ,0 Coefficients {c P X ,0g } arising from orbital perturbations cannot be given explicitly. They follow as solutions of linear equations which are generalizations of the coupled-perturbed HF equations. 57,58 As in the case of the latter equations, c P X ,0g coefficients are linearly dependent on the elements ,0 ,SCF , which allows us to write ,0 ,SCF (1) ,0 ,0 , similar for c P B ,0g , where the functions F P X ,gf depend on the energies ε P X ,0 and the matrix elements ofĤ P Let us return to the energy expression in eq 38 obtained after eliminating these terms from the second-order energy of the dimer, which cannot give rise to the dispersion energy. The C P,0j coefficients (eq 45) which enter in the first term of eq 44 depend on the interaction operator V̂P AB and a one-electron operator ĤP ,SCF (1) . Only the former is retained because, as we discussed, the latter does not contribute to the dispersion interaction. The second and third terms in eq 44 (expression for Ψ P,i (1) ), which arise from the orbital perturbation in SCF methods, cannot recover the dispersion interaction either. Their insertion into eq 38 leads to products of terms of the form ⟨Ψ Ψ |̂|Φ Ψ ⟩ V P P P P g P ,0 ,0 , , in which the interaction operator V̂P AB does not connect the ground state with bothmonomer-excited states. Such terms, corresponding to one monomer remaining in the ground state, contribute to the induction components of the interaction energy, not to dispersion.
Thus, out of all terms in E d (2) , we retain only the ones involving C P,0j and the operator V̂P AB , cf. eq 45. After expressing both the zero-order energy and wave function in terms of its monomer counterparts, eqs 41 and 42, we arrive at the energy expression which is a part of the dimer energy and gives rise to the dispersion interaction. It reads where ω P X ,i denotes a transition energy for a group P and a monomer X and only the terms for which i ≠ 0 and j ≠ 0 have been kept (naturally, the terms for which i = 0 and j ≠ 0 and vice versa will be of the induction type; therefore, they have been excluded from eq 48). A more convenient form of E disp MC follows after introducing transition density matrices. Notice that in the assumed wave function, eq 10, a set of occupied spin orbitals is divided into disjoint subsets P corresponding to group functions Ψ P . In the dimer dissociation limit, each subset decomposes into sets P A and P B of the orbitals assigned to the pertinent monomers. A transition density matrix for a monomer X in the representation of the orbitals from P X reads where ⟨pa|qb⟩ are two-electron integrals in the 1212 convention. This is a central equation of this work. E disp MC yields the amount of dispersion interaction energy recovered by the employed MC wave function model conforming to the ansatz in eq 10. Let us list the most relevant conclusions that the obtained result leads to (a) A given group P has nonzero contributions to the interaction energy of the dispersion type only if in the R AB → ∞ limit, the wave function Ψ P dissociates into a product Ψ P A Ψ P B , where both Ψ P A and Ψ P B are MC functions. The reason is that only the states excited in the configurational space of the P A and P B groups enter E disp MC , which obviously implies that configurational spaces, cf. eq 43, must be at least two-dimensional.  In this work, we focus our attention on the CAS model. A dimer wave function takes a product form, as shown in eq 11, where the active function Ψ active is an n-electron full configuration function obtained from m orbitals called active. Such a function is denoted as CAS(n,m). The dispersion energy recovered by a dimer CAS function dissociating into CAS functions for monomers can be written as where N active,X is the number of active electrons in monomer X, M nuc,X denotes the number of nuclei in X, and ρ inact,X and γ inact,X are the density and the density matrix corresponding to inactive orbitals, respectively. We have shown that the supermolecular CAS interaction energy recovers a portion of the second-order dispersion energy, quantified by eq 52, which directly relates to the size of the active space. In the extreme case, when all orbitals are active and the CAS turns into the FCI, the E disp CAS term becomes identical to the full second-order dispersion energy shown in eq 9. In other cases, the CAS supermolecular interaction energy can be corrected for the missing portion of the dispersion interaction by subtracting E disp CAS and adding the full dispersion energy recovered by a method of choice, which leads to the "CAS corrected for the dispersion" interaction energy method When defining monomer CAS wave functions, which are dissociation products of the assumed CAS(n,m) dimer function, one should be cautious. It is clear that the number of active electrons in the monomers, n A and n B , should sum up to n to satisfy the size-consistency condition. However, if m = m A + m B active orbitals are considered, m A of which are localized on A and m B of which are localized on B, then in the limit R AB → ∞, a product of CAS(n A ,m) and CAS(n B ,m) functions yields the same energy as a sum of monomer energies computed with the pertinent CAS(n X ,m X ) functions. Thus, the aforementioned product is a valid CAS zero-order wave function of a dimer, satisfying the size-consistency condition. In practical terms, this implies that in the supermolecular computation of the CAS interaction energy, the energy of monomers should be obtained with CAS(n X ,m) functions, where m includes m X active orbitals on X and m − m X orbitals on the other monomer. This problem is similar to, but not the same as, the counterpoise correction for the BSSE. We call it a counterpoise CAS-superposition error (CASSE) correction. It vanishes when the CAS function turns into the FCI and a complete basis set is used.
It should be noticed that the dispersion correction, E disp CAS , discussed in this section follows from the standard RS perturbation theory, which misses exchange contributions. If perturbation theory with symmetry forcing 59,60 were applied, it would reveal that the CAS interaction energy includes only a fraction of the second-order exchange-dispersion term, whose magnitude depends on the size of the active space. One would therefore arrive at the conclusion that the CAS must be corrected for the exchange-dispersion energy in a way similar to that for the dispersion energy, that is, in analogy to eq 54. In the next section, an approximate exchange-dispersion correction is proposed.

COMPUTATIONAL DETAILS
The second-order dispersion energies calculated in this work include both the polarization and exchange components. We use the E DISP (2) notation for the energy terms calculated directly within the SAPT 49,50 framework Unlike E disp CAS , the second-order exchange-dispersion energy recovered in supermolecular CASSCF calculations, E exch-disp CAS , is not obtained directly but approximated by scaling of the full term In analogy to eq 55, we use the (57) notation to represent the complete second-order dispersion energy recovered by the CAS (the E disp CAS term is given in eq 52). In other words, E DISP CAS denotes the dispersion energy which is accounted for in the supermolecular interaction energy, eq 2.
The "CAS corrected for dispersion" energy takes the form For comparison, we also present the results obtained by taking the sum of the CASSCF interaction energy and the E DISP (2) term (the dispersion energy included in the CASSCF interaction energy is not removed), the approach which we denote as E int CAS+DISP . An analogous notation, E int HF+DISP , is used for the HF method supplemented with full dispersion.
Evaluation of the E disp CAS formula, eq 52, requires access to both transition energies and transition density matrices of the monomers corresponding to the active group Hamiltonian, eq 53. These may be obtained either by direct diagonalization of the Hamiltonian or from full linear response equations. In this work, we followed the latter approach for the model He···He and He···H 2 dimers using time-dependent linear response equations of ref 61. For these systems, therefore, E disp CAS is computed exactly, within a given basis set. For many-electron systems, we approximate the required response properties Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article using solutions of extended random phase approximation 62,63 (ERPA) equations. The ERPA eigenproblem for each monomer is formulated for the Ĥa ct X Hamiltonian, as seen in, for example, eqs 6−9 in ref 16, which leads to a Hessian matrix of a simple form 64 (for details, see the Supporting Information). Both construction and diagonalization of the Hessian for a CAS wave function require steps which scale as M act 6 , where M act denotes the number of the active orbitals. Thus, the cost of computing E disp CAS in the ERPA approximation is small considering that regular CAS calculations employ an active space that is considerably smaller compared to the virtual space.
The CASSE introduced in the previous section was accounted for in E disp CAS energy calculations according to a three-step protocol. First, monomer CAS(n A ,m A )SCF and CAS(n B ,m B )SCF calculations were converged in the dimercentered basis set. Next, a set containing inactive and active natural orbitals of a given monomer was extended with active orbitals from its interacting partner and orthonormalized according to the Schmidt procedure. The new sets containing m = m A + m B active orbitals for each monomer were completed by adding CASSCF secondary orbitals from monomer calculations. Finally, these orbital sets entered CAS-CI calculations.
The CAS wave function captures dynamic correlation effects only within the active space. This affects the quality of both the density and density matrices 16,65 and translates into a rather poor representation of the intramonomer correlation effects in the supermolecular CAS energy. To correct this deficiency, one may describe the first-order interaction energy terms, that is, the electrostatic and exchange energies, at the fully correlated level of theory instead of the CASSCF one. This may be realized in an approximate manner by subtracting the firstorder terms calculated with CASSCF density matrices and replacing them with more accurate ones. When applied together with the discussed dispersion correction, this leads to  (59) where the first term is given by eq 58. The first and second terms in square brackets represent corrections of electrostatic and exchange interaction components, respectively, where the E elst (1) and E exch (1) terms account for intramonomer correlation effects. The E elst (1) (CAS) term is computed according to the known expression for electrostatic electron interaction, as seen in, for example, ref 49, with the electron density obtained from the pertinent CAS wave function. It is assumed that first-order exchange terms are calculated in the S 2 approximation (e.g., see eq 9 in ref 66) which requires access only to one-and twoelectron reduced density matrices of the monomers. The CASc + DispCAS approach can be interpreted as a multireference counterpart of the HFDc (1) method. 67 For the sake of convenience, acronyms and notations used in figures and tables throughout the article have been listed in Table 1.
Two model systems studied in this work are the He···He dimer and the He···H 2 dimer. The reference total interaction energy curve for He···He was taken from the work of Przybytek et al. 68 The benchmark interaction energies for excited states of the He···H 2 dimer were obtained at the FCI level of theory in the aug-cc-pVTZ basis set. 69 For manyelectron systems from the A24 data set, the coupled-cluster The E disp CAS energy calculations were performed in the inhouse code. All necessary integrals and CASSCF one-and twoelectron reduced density matrices were obtained from the Molpro 71 program. The same software was used for supermolecular CASSCF calculations. The convergence threshold for SCF orbital optimization was set to 10 −10 in the gradient norm. All reported interaction energies were corrected for the BSSE using the CP scheme 20 of Boys and Bernardi. In contrast to E disp CAS , the supermolecular CASSCF results were not corrected for the CASSE.
CASSCF wave functions for systems from the A24 data set were obtained with MP2 natural orbitals employed as an initial guess for the CAS optimization. The supermolecular CAS interaction energy may suffer from the lack of size consistency if dimer active orbitals, obtained for the assumed CAS(n,m) model and yielding the lowest energy of a dimer in the dissociation limit, do not correspond to the optimal active orbitals of the monomers [it is assumed that monomers are described with CAS(n X ,m X ) functions such that n A + n B = n and m A + m B = m]. To eliminate the risk of size-consistency violation, for each dimer, the CASSCF energy was computed at an intermonomer distance of 200 bohr and it was checked that the dimer energy agrees with the sum of monomer energies within at least 10 −4 mHa.
For both the He···He and He···H 2 dimers, the reference values of SAPT energy components were calculated at the FCI level of theory in a house-developed code. 65,72 In the case of He···He, both dispersion and exchange-dispersion components were extrapolated to the CBS limit according to the aug-cc-pVQZ → aug-cc-pV5Z two-point scheme of Bak et al. 73 For He···H 2 , all calculations employed the aug-cc-pVTZ basis and CBS extrapolation was not performed. For systems from the A24 data set, the SAPT energy components entering eqs 58 and 59 were obtained at the DFT-SAPT 74−76 level of theory based on the asymptotically corrected 77 PBE0 78,79 functional using the dimer-centered basis for monomer calculations. The  (2) second-order dispersion energy, eq 9 E exch-disp (2) second-order exchange-dispersion energy E DISP (2) sum of second-order dispersion and exchange-dispersion energies, eq   In this section, we examine the performance of the introduced "CAS plus dispersion" models for the interaction between two-electron monomers. As mentioned, in this case, we use accurate CAS transition density matrices which enter the E disp CAS expression. In calculations of the latter, the CASSE is removed according to the protocol described in the previous section. The second-order dispersion energy is obtained from SAPT calculations at the FCI level of theory.
As we have rigorously shown in Section 2.2, CASSCF recovers a portion of the dispersion energy, as seen in eq 52. Therefore, a simple addition of the dispersion component, cf. eq 55, to the supermolecular CASSCF interaction energy leads to double counting of intermonomer correlation effects, the magnitude of which depends on the size of the active space. This is illustrated in Figure 1 (see "CAS + DISP" curves) on the example of the He···He dimer in the ground state. While a small CAS(4,4) wave function of the dimer [dissociating to a product of CAS(2,2) monomer functions] supplemented with dispersion deviates from the benchmark only by 6.6% in the van der Waals minimum (R eq = 5.6 bohr), going to the CAS(4,28) model leads to 72% overestimation of the interaction energy. Additionally, dispersion double counting shifts the minimum from the correct distance of 5.6 bohr for CAS(4,4) to 5.3 bohr for CAS (4,28). In the proposed CAS + DispCAS approach, we improve the CAS + DISP results by removing the partial dispersion contained in the supermolecular CASSCF energy, cf. eq 58. Consequently, the extension of the active space does not lead to overcorrelation and the errors in the He···He minimum region do not exceed 2% regardless of the active space size (Figure 1). The only exception is the CAS(4,4) + DispCAS curve where the error reaches 7%, which is mainly due to the poor representation of the first-order energy terms (see Figure S1 in the Supporting Information).
In the asymptotic regime the He···He interaction is dominated by the dispersion energy which decays with the inverse sixth power of the interatomic distance. Therefore, at large distances, the CASSCF interaction energy and the E disp CAS term calculated directly according to eq 52 should become identical. Indeed, in Figure 2 we demonstrate a perfect agreement on the example of CAS(4,10) and CAS(4,18) wave functions. Even in the larger active space, CASSCF recovers only approximately 50% of the C 6 coefficient, that is, 0.735 a.u. compared to the reference value of 1.461 a.u. 80 A closer comparison with the long-range behavior of the FCI curve reveals that neither CAS(4,10) nor CAS(4,18) recovers higher C n coefficients. This is not surprising because only p orbitals are included in the active space. A more challenging model system is the He···H 2 dimer in which the hydrogen molecule is in one of the two lowest excited states, that is, either the 1 Σ u + or the 1 Π u state. Note that the dispersion corrections employed in "CAS plus dispersion" approaches may be applied to excited states provided that the zero-order wave function is nondegenerate. We present the results for the T-shaped geometry where we vary the distance R between the center of mass of H 2 and the He atom. The equilibrium bond lengths in the H 2 molecule are 2.443 bohr for the 1 Σ u + state and 1.951 bohr the 1 Π u state. We begin by  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article comparing the performance of the "CAS plus dispersion" models. In Figure 3, we show the results for the He···H 2 ( 1 Π u ) complex obtained by choosing 2 active orbitals on the He atom and either 4 or 10 active orbitals on the H 2 molecule. Clearly, increasing the number of active orbitals leads to a severe overestimation of the interaction energy if the E DISP (2) term (eq 55) is directly added to the CASSCF interaction energy (cf. "CAS + DISP" curves in Figure 3). Because of the correlation double counting in the larger CAS space, the minimum is predicted at 2.8 bohr, while the smaller CAS + DISP model correctly locates it at 3.2 bohr. When the amount of the added dispersion is tailored to the employed multireference model, as proposed in the CAS + DispCAS approach, interaction energy curves greatly improve. For CAS (4,6), the error in the van der Waals minimum amounts to 2.3%. In CAS (4,12), it is larger and reaches 11%. Considering that the dispersion energy is calculated exactly, this should be attributed to the erroneous representation of the remaining energy components in the supermolecular CASSCF approach.
Because applications of multireference methods to interacting excited systems are scarce, it is instructive to compare the performance of CAS + DispCAS with that of other approaches.
In Figure 4, the results from CASPT2 and SAPT(FCI) calculations up to the second order in V̂i nt (the latter denoted as E (1) + E (2) ) are shown next to the CASSCF and "CAS plus dispersion" models. The He atom and H 2 ( 1 Σ u + ) molecule were described with the CAS(2,5) and CAS (2,4) wave functions, respectively. In the He···H 2 ( 1 Π u ) complex, we chose the CAS(2,2) active space for He and CAS(2,10) for H 2 . The corresponding active spaces assumed for a dimer CASSCF calculation were (4,9) and (4,12) for the considered states. By confronting Figure 4, one observes a strikingly poor performance of the second-order SAPT energy (see "E (1) + E (2) " curves). Although the corresponding curves account for the second-order dispersion energy and in consequence show a correct asymptotic behavior, at intermediate-and short intermonomer distances second-order SAPT becomes even less accurate than CASSCF. This indicates that higher than second order induction terms, accounted for in CASSCF but missed in SAPT, are comparable in magnitude to the considered E (1) + E (2) energy contributions. In the CASPT2 method the perturbation correction is added to the CASSCF energy, so the aforementioned higher-order terms are included. In spite of this, CASPT2 significantly underestimates the interaction energy for the 1 Σ u + state. In particular, the minimum is predicted at R = 4.1 bohr and 1.0 mE h deep, which is too shallow compared to the FCI benchmark that is 1.7 mE h deep and located at 3.7 bohr. For the 1 Π u state, CASPT2 performs better and the curve has a correct shape. The interaction energy at the minimum (R = 3.2 bohr) is underestimated by 12.3%, which is similar in magnitude to the −11.0% overestimation from CAS + DispCAS calculations.
As already discussed, direct addition of the dispersion energy in the CAS + DISP approach leads to severe overcorrelation in both the He···H 2 ( 1 Π u ) and He···H 2 ( 1 Σ u + ) states (Figure 4). Contrary to this, CAS + DispCAS avoids dispersion double counting and performs well in both cases. The positions of van der Waals minima match the FCI benchmark and are predicted with −9.6 and −11.0% relative percent errors for the Σ u + and Π u states, respectively. The improvement of the first-order interaction energy components according to the procedure proposed in eq 59, cf. the CASc + DispCAS curves, has a marginal effect on the CAS + DispCAS energies in the 1 Π u state. For example, the error in the minimum is reduced from −11% to −8.6%. This points toward the inaccurate treatment of higher than second order terms. In the 1 Σ u + state, the  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article difference is more pronouncedat the minimum geometry (R = 3.7 bohr) CASc + DispCAS deviates from the FCI by 2.2%, while CAS + DispCAS is off by − 9.6%.

Correcting Dispersion and First-Order Interaction Energy Components in the CAS Model for Many-Electron Systems.
For model, few-electron systems discussed in Section 4.1, it is possible to employ large active spaces and restore a significant portion of the dispersion energy in sizeconsistent supermolecular CASSCF calculations. To verify what amount of the dispersion energy is captured by CASSCF for many-electron systems and examine the effect of correcting first-order interaction energy according to eq 59, in this section, we focus on a set of small, noncovalently bound dimers selected from the A24 database. 70 (2), and C 2 H 2 ···C 2 H 2 (2)]. In Table 2, we report the interaction energy values following from "CAS plus dispersion" and "HF plus dispersion" calculations [CCSD-(T)/CBS is taken as the reference]. Table 2 also contains error statistics in terms of the mean error Δ̅ , the standard deviation σ, the mean absolute error Δ̅ abs , and the maximum absolute error Δ max . We found that for the selected many-electron systems, the effect of the CASSE on the E disp CAS energy was negligible and the results presented in this section do not account for it (see Table S6 in the Supporting Information).
As expected for single-reference systems, both the CASSCF and HF methods not corrected for dispersion show similar, poor performance with mean errors that amount to 128% and 141%, respectively ( Table 2). Addition of the E DISP (2) dispersion term from SAPT(DFT) calculations drastically improves the results (Δ̅ abs = 19% for CAS + DISP and Δ̅ abs = 25% for HF + DISP). Note that the CAS + DISP model clearly outperforms its HF counterpart for hydrogen-bonded dimers. In contrast, HF + DISP performs better for systems dominated by dispersion. The CAS-based approach may be corrected further by removing from CAS + DISP the dispersion energy recovered by CAS, that is, the E DISP CAS term, as done in the CAS + DispCAS model (see eq 58). This has a weak effect the values of |E DISP CAS | are smaller than 0.05 kcal/mol for the investigated systems. The only notable error reduction in CAS + DispCAS with respect to CAS + DISP occurs for the van der Waals dimers and does not exceed 5% (see also Table S4 in the Supporting Information). We conclude that in many-electron systems described with moderate active spaces, the double counting of the dispersion energy does not play a significant role, with the only exception being van der Waals complexes where it introduces discernible errors.
Both HF and CASSCF clearly benefit when the first-order electrostatic and exchange interaction energy components are treated at the SAPT(DFT) level of theory. The resulting HFc and CASc models provide a nearly 2-fold reduction in terms of both mean errors and the standard deviation when compared with their uncorrected HF and CAS equivalents ( Table 2). In particular, the CASc + DispCAS approach, as formulated in eq 59, performs the best in the group (Δ̅ = −5.7% and σ = 12.2%). Still, relatively large errors between 13 and 27% persist for dispersion-bounded dimers for which the interaction energy has a repulsive character.

CONCLUSIONS
We presented a way to quantify the amount of the dispersion energy recovered in the supermolecular interaction energy, the residual dispersion energy, calculated with MC wave functions. Applying the RS perturbation theory to a group product wave function, we arrived at a general expression that rigorously yields a portion of the dispersion energy in the MCSCF interaction energy (eq 51). It is different from 0 if a dimer wave function is such that it dissociates into monomer wave functions both of which are at least two-configurational. The dispersion-in-MCSCF energy increases upon active space  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article expansion and reaches the full second-order dispersion limit when all orbitals are active. By deriving the expression for the residual dispersion, we were able to rigorously address the problem of dispersion double counting which occurs when the dispersion energy is added directly to the supermolecular MCSCF interaction energy. We developed a dispersion correction for wave functions of the CAS type. The correction is tailored to the employed active space so that its addition to the CAS interaction energy, as shown in eq 54, accounts for exactly the amount of the second-order dispersion energy which is missed by the CAS. The resulting "CAS plus dispersion" method requires computing the E disp CAS energy, eq 52, either in an exact manner by diagonalization of the active Hamiltonians for the monomers or approximately, by solving ERPA equations formulated with the active Hamiltonian for each monomer. One should notice that the proposed route of avoiding dispersion double counting in supermolecular CAS calculations is rigorous and does not rely on orbital localization, contrary to the physically motivated scheme outlined in ref 46. In the latter approach, the degree to which double counting of dispersion is avoided will depend on how well orbitals are localized on fragments. The most expensive step of computing E disp CAS scales only with the sixth power of the number of active orbitals if ERPA equations are used, so it is modest comparing to the cost of CASSCF calculation. To assure the proper physical damping of the dispersion energy, we include exchange-dispersion terms approximated with second-order SAPT energy expressions, eqs 55−58. Therefore, in the CAS + DispCAS scheme proposed in this work, the CASSCF interaction energy is corrected for both dispersion and exchange-dispersion components.
On the example of the model He···He and He···H 2 systems in the ground and excited states, we demonstrated that the direct addition of the full second-order dispersion to the CAS interaction energy leads to overcorrelation which increases with the number of active orbitals. In contrast, the tailored dispersion correction applied in the CAS + DispCAS approach yields an accurate interaction energy regardless of the size of the active orbital set which validates the method. For the He··· H 2 ( 1 Σ u + ) and He···H 2 ( 1 Π u ) complexes, CAS + DispCAS proved to be superior to CASPT2 and second-order SAPT. This shows that the method may be successfully applied to study excited-state interactions provided that an accurate representation of the dispersion energy is available.
To examine the extent of the dispersion double counting in CASSCF calculations for many-electron systems, we studied several representative dimers selected from the A24 data set. We found out that in small active spaces feasible for sizeconsistent calculations [up to CAS (12,12) for the dimer], CASSCF recovers only 2−3% of the true dispersion. The overcorrelation becomes relevant in van der Waals (dispersiondominated) systems where it translates into interaction energy errors at the level of 4−5%. This effect may be larger in the case of systems in excited states, when a few active diffuse orbitals contribute significantly to the dispersion interaction. Taking into account the low cost of computing the residual dispersion (as shown in this work, it depends only on the number of active orbitals, which is always much lower than the total number of orbitals), it is advisable to remove it from the supermolecular energy, even when the employed active space is relatively small.
Applications of the proposed CAS + DispCAS method for truly multireference systems require an adequate representation of the dispersion energy. Currently, it is possible to calculate both the polarization and exchange components of the second-order dispersion energy for many-electron CAS wave functions within the ERPA formalism. 16,65 A substantial improvement of the CAS + DispCAS accuracy is achieved if both the electrostatic and exchange first-order components recovered by the CAS are corrected for intramonomer correlation effects, as shown in eq 59. This approach is analogous to the HF-based HFDc (1) method of ref 8. Improving the quality of 1-RDMs from CASSCF calculations to compensate for the missing intramonomer correlation is less straightforward compared to the single-reference HFDc (1) approach. The multiconfiguration DFT theory could be explored in this context. Works on further verification and improvement of the accuracy provided by the CAS + DispCAS method, particularly in complexes involving excited states, are being carried out in our group.
ERPA equations for the active group Hamiltonian; second-order E disp CAS(n,m) energy values for the He···He interaction energy curve; CASSCF, CASPT2, SAPT, and "CAS plus dispersion" interaction energy curves for the He···H 2 ( 1 Σ u + , 1 Π u ) dimers; and absolute CASSCF energies and individual relative percent errors of "CAS plus dispersion" models for the systems from the A24 data set (PDF)