Efficient Calculation of the Dispersion Energy for Multireference Systems with Cholesky Decomposition: Application to Excited-State Interactions

Accurate and efficient prediction of dispersion interactions in excited-state complexes poses a challenge due to the complex nature of electron correlation effects that need to be simultaneously considered. We propose an algorithm for computing the dispersion energy in nondegenerate ground- or excited-state complexes with arbitrary spin. The algorithm scales with the fifth power of the system size due to employing Cholesky decomposition of Coulomb integrals and a recently developed recursive formula for density response functions of the monomers. As a numerical illustration, we apply the new algorithm in the framework of multiconfigurational symmetry adapted perturbation theory, SAPT(MC), to study interactions in dimers with localized excitons. The SAPT(MC) analysis reveals that the dispersion energy may be the main force stabilizing excited-state dimers.

M odeling of dispersion forces is crucial for an accurate representation of noncovalent interactions in molecular systems 1,2 and materials. 3−7 Unfortunately, approaches to calculating the dispersion energy in excited-state complexes have been scarce. 8−13 Semiempirical dispersion energy corrections for density functionals for ground-state complexes generally fail for dimers in excited states. 13,14 So far, two pairwise dispersion approaches have been extended to excited states. First, the local response dispersion (LRD) model of Nakai and co-workers 8,15 was applied to exciton-localized complexes from the S66 data set. 16,17 Second, Feng et al. 18 used the exchange-hole dipole moment (XDM) method of Becke and Johnson 19,20 to obtain van der Waals C 6 coefficients in systems involving inter-and intramolecular charge transfer excitations. The proposed generalizations of both LRD and XDM rely on the excited-state electron density extracted from TD-DFT.
When ground-state interactions are concerned, accurate values of the dispersion energy can be obtained from single reference symmetry-adapted perturbation theory (SAPT) 21, 22 based on either coupled-cluster 23−25 or DFT description of the monomers. 26−29 These methods are not applicable to excited states. Recently, we have developed a wave function-based approach to the dispersion energy in ground and excited states, 30,31 which employs the extended random phase approximation (ERPA) for density response. 32 The dispersion energy can then be predicted for any molecular system with a local exciton. 12 A complete description of noncovalent interactions is accessible when combining our model with multiconfigurational SAPT, 33 SAPT(MC), or a supermolecular approach based on the multiconfigurational self-consistent field (MCSCF), in particular, complete active space (CASSCF), description of the dimer. 34 In the latter method, named CAS +DISP, the supermolecular CASSCF energy is corrected for the missing part of the dispersion energy. Both SAPT(MC) and CAS+DISP have already proven useful in studying excitedstate organic dimers. 12 Currently, the bottleneck in both SAPT(MC) and CAS +DISP is the calculation of coupled dispersion and exchangedispersion energy contributions, as the computational cost of both components grows formally with the sixth power of the system size.
The goal of this work is to extend the applicability of the SAPT(MC) and CAS+DISP methods to larger systems by reducing the scaling of the coupled dispersion energy from m 6 to m 5 . For this purpose, a novel algorithm is proposed. It employs a Cholesky decomposition technique and recently introduced recursive formula for the computation of density response functions. 35 The m 5 scaling is achievable if the interacting monomers are described with multiconfigurational (MC) wave functions, e.g., CASSCF, and the number of active orbitals is much smaller than that of virtual orbitals, which is typically the case. The new developments are applied to study molecular interactions in excited-state organic complexes of larger size than those affordable until now for multiconfigurational dispersion methods.
The approach for multireference functions parallels that in previous works focused on coupled dispersion energy computations for single-reference wave functions. In particular, SAPT based on the Kohn−Sham description of the monomers, SAPT(DFT), 28,29 may employ either the density-fitting (DF) 36 or Cholesky decomposition 37 techniques, both leading to the m 5 scaling of the dispersion energy computation. The algorithm of Bukowski et al., 36 recently improved by Xie et al., 38 is most general and applicable to computing the density response of the monomers from both local and hybrid functionals. In the case of the exchange-dispersion energy, the computational cost remains as large as m 5 in singlereference SAPT (m 6 in the multireference case) 31 even in the DF/Cholesky formulation. 38−42 The spin-summed second-order dispersion formula written in terms of monomer response properties obtained within the extended random phase approximation (ERPA) reads 30 where pqrs represents natural orbitals (NOs) of the monomers. Modified two-electron integrals in the NO representation are defined as The same partitioning also applies for the p′q′ multi-index. Thus, a straightforward implementation of ERPA requires steps that scale as n OCC is the number of generalized secondary orbitals (M s i denotes cardinality of the set s i ). An evaluation of the dispersion energy formula, eq 1, shares the same scaling behavior (indices μ and ν run over all n OCC n SEC eigenvectors). Below we propose an algorithm leading to a lowered m 5 scaling.
Using the integral identity and introducing the frequency-dependent matrix C X (ω) leads to another representation of E disp The C X (ω) matrix is equivalent to the real part of the density linear response function taken with the imaginary argument iω (see eq 33 in ref 44) and is obtained by solving the following equation 35 [ The modified two-electron integrals g pqrs of eq 2 can be represented via the decomposition where vectors D L are obtained by Cholesky decomposition of the AO Coulomb matrix followed by transformation to the natural-orbital representation and scaling by n p 1/2 + n q 1/2 factors. The matrix product of C X (ω) with D yields a reduced-dimension intermediate which allows one to write Notice that the obtained formula applied to excited-state systems would miss the so-called non-Casimir−Polder terms arising from negative transitions ω ν X < 0, i.e., transitions from higher-to lower-energy states, present in the density response function of an excited-state monomer. Extended RPA does yield negative excitations, which are considered spurious and discarded in practice, i.e., excluded from eq 1. In ref 12, we have shown how to account for negative excitations explicitly, The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter in a physically meaningful manner. Upon inspection, we found that such non-Casimir−Polder terms are negligible for the studied systems, and they are not discussed any further.
The key step in the proposed reduced-scaling algorithm for computing dispersion energy with a multireference description of monomer wave functions assumes partitioning of a monomer Hamiltonian H X into a partially correlated effective . Then, the parametric representation of the Hamiltonian is introduced where the H X operator is multiplied by the coupling constant α ∈ [0, 1]. There are two underlying requirements in the Hamiltonian partitioning. The first one is that the wave function describing monomer X is of zeroth order in α for H X . The other condition is that scaling of the ERPA equations corresponding to H X is lowered from m 6 to m 5 at α = 0. It has been shown that a group-product-function Hamiltonian 45,46 H X (0) satisfies both conditions for the MC wave function based on an ansatz assuming the partitioning of orbitals into inactive, active, and virtual ones. Notice that for a single-reference wave can be chosen as a noninteracting Hamiltonian, as in the Møller−Plesset (MP) perturbation theory.
After employing the partitioned Hamiltonian H X (eq 12), the resulting ERPA Hessian matrices ± ( ) become linear functions of α. (Notice that from now on the index X in Hessian and response matrices is dropped for simplicity.) By contrast, the response matrix C(α, ω) (eq 8) depends on α in a nonlinear fashion. Let us represent the projected response matrix C̃(α, ω) (eq 10) as a power series expansion in the coupling constant α around α = 0. Truncating the expansion at the nth order and setting α = 1, we obtain where C ( ) k ( ) follows from an efficient recursive scheme derived in ref 35 The required matrices are given by the ERPA matrices ± (0) and ± (1) Notice that by setting n = 0 in eq 14 for each monomer, the dispersion energy obtained from eq 11 will be equivalent to the uncoupled approximation introduced in ref 30. For a sufficiently large n, one recovers full dispersion energy, i.e., the coupled dispersion energy, 30 numerically equal to that following from eqs 1−3.
Since the dimensions of the Hessian matrices and the matrix D are m 2 × m 2 and m 2 × N Chol , respectively, matrix multiplications involved in eqs 15−17 scale as m 4 N Chol ∼ m 5 22) is negligible if the number of active orbitals, M s 2 , is much smaller than that of virtual orbitals, which is usually the case in practical calculations. The combination of eq 11 with the recursive scheme for is the main contribution of this work.
A valid concern is whether the expansion of the linear response function at α = 0 (eq 14) leads to a convergent series. Although definite proof cannot be given, it is reasonable to expect that if a monomer wave function leads to stable ERPA equations around α = 0 then the series converges. Our numerical tests on two data sets of small, weakly correlated dimers 47,48 have shown no convergence problems (Supporting Information). In most cases, expansion up to n = 8 was sufficient to achieve μE h accuracy, amounting to the mean absolute percentage error below 0.1% in the dispersion energy. Convergence tests carried out on larger dimers, selected from the S66 test set, in both ground and excited states have led to the same conclusions; see Table S1 in the Supporting Information. An example of the convergence of E disp (2) computed according to the procedure given by eqs 11−22 with the truncation order n in the range from n = 1 to 10 is presented for the benzene-cyclopentane complex in Figure S1 in the Supporting Information.
The evaluation of the α-expanded response matrix requires access to ± (1) Hessians together with three projected C( ) k ( ) matrices needed to carry out the recursion. Due to their size (m 2 × m 2 and m 2 × N Chol , respectively), for systems approaching 100 atoms, these quantities can no longer fit into memory and have to be stored on disk. In this regime, the disk storage and the number of I/O operations will become the main bottleneck of the proposed approach.
The scaling of second-order induction energy computations in SAPT(MC) can be reduced to m 4 using α-expansion of the response functions, accompanied by the Cholesky decomposition of two-electron integrals. However, it is possible to achieve such scaling without relying on the coupling constant expansion; see the Supporting Information for details. For the Hartree−Fock treatment of the monomers, such an alternative approach is identical to induction energy computations in the coupled Hartree−Fock scheme, as first proposed by Sadlej. 49 Numerical demonstration of the developed algorithm was carried out for both ground and excited states of noncovalent complexes selected from the S66 benchmark data set of Hobza and co-workers. 16,17 The dimers were divided into two sets according to their size. Smaller systems (up to eight heavy atoms in a dimer and ca. 500−600 basis functions with a basis The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter set of a triple-ζ quality) were analyzed in detail in our recent work. 12 Larger systems (up to 11 heavy atoms in a dimer and ca. 800−900 basis functions), which are beyond the capabilities of the m 6 -implementation of the MC dispersion energy, are analyzed for the first time. This group contains five complexes: benzene-cyclopentane, benzene-neopentane, AcOH-pentane, AcNH 2 -pentane, and peptide-pentane, where peptide refers to N-methylacetamide (see Figure 1).
Both ground-and excited-state calculations were performed using ground-state geometries taken from ref 16. All supermolecular calculations employed the Boys−Bernardi counterpoise correction. 50 The excitons were localized on benzene (π → π*), AcOH (n → π*), AcNH 2 (n → π*), and peptide (n → π*) molecules. The CCSD(T) results extrapolated to the complete basis set limit (CBS) 16 served as a benchmark for ground-state interaction energies. In the case of complexes involving excited states, reference results were taken from ref 8. They were obtained by combining CCSD(T)/CBS ground-state interaction energies with excitation energies calculated at the EOM-CCSD 51 / 6-31++G(d,p) 52−54 level of theory.
The Cholesky decomposition of the Coulomb integrals matrix, ⟨pr|qs⟩, was performed in the AO basis with a modified program developed for refs 55 and 56. The Cholesky vectors, R pq,L , were generated with the convergence criterion ∑ p≥q (⟨pp|qq⟩ − ∑ L R pq,L R pq,L ) < 10 −2 , which is the same as used in ref 35.
Second-order dispersion energies and SAPT(MC) 33 energy components based on CASSCF treatment of the monomers   were computed in the GammCor 57 program. From now on, SAPT(MC) based on CASSCF wave functions is denoted as SAPT(CAS). The frequency integration in eq 11 has been carried out using the eight-point Gauss−Legendre quadrature. The necessary integrals and reduced density matrices were obtained from the locally modified Molpro 58 package. Supermolecular CASSCF and DFT-SAPT calculations were performed in Molpro. All calculations employed the aug-cc-pVTZ basis set. 59,60 Excited-state wave functions were computed with two-stateaveraged CASSCF. We chose the same active spaces for both ground-and excited-state calculations. The active space for benzene included three π bonding and three π* antibonding MOs, which means six active electrons on six orbitals, labeled as CAS (6,6). 61 For AcOH, we chose the CAS(8,8) active space including two n, π, π*, two σ, and two σ* orbitals. 62 For AcNH 2 , the CAS(6,5) space was selected which involves σ, n, π, π*, and σ* orbitals. 63 The peptide (N-methylacetamide) active space, CAS (6,6), was composed of σ, π, π*, and σ* orbitals, and two lone-pair orbitals n located on the oxygen atom. 64 To improve the accuracy of the SAPT(CAS) interaction energy, especially for systems dominated by large polarization effects, we need to include higher than second-order induction terms. For ground-state systems which can be represented with a single Slater determinant, these terms can be approximated at the Hartree−Fock (HF) level of theory and represented as the δ HF correction 65 where E int HF corresponds to the supermolecular HF interaction energy and all terms are computed with Hartree−Fock wave functions. There is no straightforward way to account for higher-order polarization terms in excited-state computations.
To tackle this problem, we assume that the change of higherthan-second-order induction terms upon excitation is proportional to a corresponding shift in the second-order induction and define the δ CAS correction as (2) HF (24) where the labels GS/ES correspond to dimers in ground and excited states, respectively. Notice that in our previous work 12 a similar scaling expression involved sums of induction and exchange-induction (E exch−ind (2) ) terms. In this work, the latter is not computed directly but follows from an approximate scaling relation (see below). Such a treatment of E exch−ind (2) energy could introduce additional error into the δ CAS term, and we decided not to include it in eq 24.
Compared to single-reference SAPT schemes, in multiconfigurational SAPT it is not straightforward to apply the Cholesky decomposition to second-order exchange energy components, i.e., E exch−ind (2) and E exch−disp (2) . The difficulty follows from the necessity to obtain separately lower and upper triangles of transition density matrices (see the discussion in section 2 of ref 31); possible solutions will be addressed in our future work. To account for both exchange-induction and exchange-dispersion terms in this study, we propose a simple scaling scheme   where aVXZ = aug-cc-pVXZ, and it is assumed that the convergence of second-order polarization and exchange components with the basis set size is identical. Expressions for E exch−ind (2) and E exch−disp (2) terms are given in ref 33; the induction energy is computed with the m 5 -scaling algorithm presented in the Supporting Information.
All presented SAPT(CAS) results include δ HF and δ CAS corrections for ground-and excited-state complexes, respectively, as well as scaled second-order exchange components defined in eqs 25 and 26. CAS+DISP is a sum of the supermolecular CASSCF interaction energy and the dispersion energy, DISP = E disp (2) + E exch−disp (2) , computed in the same fashion as in SAPT(CAS), i.e., using the newly developed expression given in eq 11 and the scaling relation from eq 26.
In Table 1, we present SAPT interaction energy decomposition for ground-and excited-state complexes. Regardless of the electronic state of the dimer, all systems can be classified as dispersion-dominated with the E disp (2) /E elst (1) ratio ranging from 2.8 to 3.8.
As can be deduced from Table 1, the most significant change in dispersion interactions upon transition from the ground to the excited state occurs in complexes of benzene (benzene··· cyclopentane and benzene···neopentane). The effect amounts to ΔE disp (2) ≈ 0.2 kcal mol −1 , which corresponds to a decrease in the dispersion energy in the excited state. In both complexes, the redistribution of the electron density upon π → π* excitation on benzene is accompanied by a non-negligible drop in the electrostatic attraction. The latter energetic effect is, however, canceled by the simultaneous depletion of the exchange repulsion. Thus, a decline of the dispersion energy contributes in a major way to the weakened net attraction in the excited state. We observed the same trends in dimers of benzene with H 2 O, MeOH, and MeNH 2 studied in ref 12. Compared to benzene π → π* systems, complexes of npentane involve an n → π* exciton localized on the carbonyl group of the interacting partner (AcOH, AcNH 2 , and peptide). These systems exhibit an increase in the dispersion energy upon excitation, which ranges from −0.06 to −0.15 kcal mol −1 ( Table 1). In AcOH···pentane, the enhanced dispersion is comparable in magnitude with a concurrent decrease in the first-order Pauli repulsion, both of which contribute to the overall stabilization of the excited-state dimer. In contrast, in peptide···pentane and AcNH 2 ···pentane, the net repulsive components become stronger and outweigh the dispersion attraction so that both complexes are more strongly bound in the ground state. For peptide···pentane, the weakened interaction in the excited state can be attributed mainly to a significant increase in second-order exchange-induction (ΔE exch−ind (2) = 0.21 kcal mol −1 ). The other interaction energy components undergo relatively minor changes; the only stabilizing effect of −0.07 kcal mol −1 is due to dispersion. In the case of AcNH 2 ···pentane, increased static polarizability of acetamide in the excited state results in a stronger induction attraction. (The net change in induction and δ corrections amounts to −0.60 kcal mol −1 .) This, however, is countered by a steep rise in the repulsive components. In particular, exchange-induction and first-order exchange increase by 1.00 and 0.45 kcal mol −1 , respectively. Note that a similar pattern occurred in the methylamine···peptide (n − π*) interaction. 12 The use of eq 11 enables spatial visualization of the dispersion interactions. Following the work of Parrish et al., 67,68 Wuttke et al., 69 and our recent development, 70 we propose a spatially local dispersion density function , where Q X (r) are constructed from occupied natural orbital densities weighted by contributions to the dispersion energy. For details, see the Supporting Information. The Q AB (r) density function integrates to the dispersion energy The additional cost of the Q AB (r) computation is marginal, as all intermediate quantities are available from the calculation of E disp (2) . Changes in the E disp (2) components are visualized in Figure 2 by using the difference between the ground-and excited-state dispersion interaction density, Q AB (r). Both the sign and magnitude of the effect are correctly captured�one observes a notable depletion of the dispersion density in π − π* complexes compared to a weaker accumulation for n − π* dimers. In agreement with the character of the underlying excitons, in the π − π* case, the majority of the ΔE disp (2) term is delocalized over the benzene ring, while in n − π* dimers it is basically confined to the carbonyl group.
In Tables 2 and 3, we report total ground-and excited-state interaction energies, respectively, calculated at the CASSCF, CAS+DISP, and SAPT(CAS) levels of theory. Addition of the dispersion energy to supermolecular CASSCF changes the character of the interaction from repulsive to attractive, reducing mean errors by 2 orders of magnitude. As a consequence, CAS+DISP results closely match the coupledcluster (CC) reference 8 with mean absolute percentage errors (MA%Es) of 0.8% for the ground and 3.7% for excited states. SAPT(CAS) performs similarly to the CAS+DISP model (MA%E values of 2.5 and 3.3% for ground and excited states, respectively). The DFT-based LRD model of Nakai et al. 8 combined with the LC-BOP functional 71−73 is somewhat less accurate. The model underestimates interaction energies, which amounts to mean errors at the level of 10% (Tables 2  and 3). Note, however, that the DFT results were obtained with the 6-311++G(2d,2p) basis set.
Since the aug-cc-pVTZ basis set is not sufficient to saturate the dispersion energy with respect to the basis set size, the good agreement of both SAPT(CAS) and CAS+DISP with coupled cluster is partially due to error cancellation. (The CC values include CBS-extrapolated ground-state energies.) Indeed, individual SAPT(CAS) energy components for groundstate complexes are systematically underestimated with respect to their SAPT(DFT) counterparts (Tables S2−S4 in the  Supporting Information). This reflects the effective neglect of intramonomer electron correlation effects in the SAPT(CAS) approach. 30,31,33 To summarize, we have presented an algorithm for secondorder dispersion energy calculations with multiconfigurational wave functions that scales with the fifth power of the system size. Until now, m 5 scaling in coupled dispersion energy computations could be achieved only with a single-determinant description of the monomers. 29, 36,38,39,74 The prerequisite for m 5 scaling with a multiconfigurational reference is that the number of active orbitals in wave functions of the monomers is considerably smaller compared to the number of virtual orbitals. In practice, this condition is typically fulfilled in interaction energy calculations performed using augmented basis sets. The new m 5 dispersion energy algorithm was employed together with state-averaged CASSCF wave functions to examine interactions involving localized excitons of the π − π* and n − π* types. For the representation of the interaction energy, multiconfigurational dispersion energy was complemented either with SAPT(MC) 33 energy components or with supermolecular CASSCF interaction energy, the latter known as the CAS+DISP 34 method. In line with earlier investigations, 12 we have found out that in low-lying excited states the dispersion energy may be the driving force behind the stability of the complex. Hence, both accurate and efficient algorithms adequate for dispersion computations with multiconfigurational wave functions are mandatory. Spatial mapping of the dispersion energy density helps to identify the regions affected most by the exciton. Visualizing the remaining SAPT(MC) energy components could aid in the interpretation of energetic effects that occur upon electron density rearrangement in excited states. Work along this line is in progress.

■ ASSOCIATED CONTENT Data Availability Statement
The raw data are available in the Zenodo repository at 10.5281/zenodo.8131535.
Description of the m 4 algorithm for induction computations in the SAPT(MC) framework. Additional results including SAPT(DFT) and SAPT(CAS) energy decomposition, total CASSCF energies, and convergence of the dispersion energy with n for data sets in refs 16