Intruder-free cumulant-truncated driven similarity renormalization group second-order multireference perturbation theory

Affordable and accurate multireference (MR) electronic structure theories are indispensable for constructing potential energy surfaces. Because of their low computational cost, multireference methods that treat dynamical correlation effects via low-order perturbation theory are the most frequently applied.^1–9 A major limitation is extending these approaches to treat many active orbitals. For example, the ideal active space for large conjugated hydrocarbons (e.g., acenes) would include all π orbitals, which becomes a bottleneck for systems with more than 10–12 carbon atoms.^10,11 Additionally, computations on polynuclear transition metal compounds, which require at least one d shell in the active space for each metal atom, can quickly become prohibitively expensive for more than three metal centers.^12,13 A common strategy to reduce the cost of multireference perturbation theories (MRPTs) is to express the perturber states using partially or fully internally contracted configurations. This approach reduces the scaling of a method with respect to the number of active orbitals (N_A) from exponential to polynomial; however, it still requires calculating and storing the three- and four-body reduced density matrices (RDMs) of the reference wavefunction, whose storage scales as $O(NA6)$ and $O(NA8)$⁠, respectively. This reduction in computational and storage costs has enabled MRPT computations with as many as 20 active orbitals. Extending multireference computations beyond this limit requires avoiding the costs introduced by the 3- and 4-RDMs.

Strategies for reducing the cost associated with high-order RDMs fall into two main categories: (1) modifying the existing algorithms to eliminate the bottlenecks and (2) approximating the RDMs. The first category encompasses methods such as second-order complete-active-space perturbation theory (CASPT2)^1,2 in combination with the density matrix renormalization group (DMRG-CASPT2) as proposed by Kurashige and Yanai¹⁴ and later optimized by Wouters and co-workers.¹⁵ These works take advantage of the matrix product structure of DMRG to reduce the cost of terms involving the 4-RDM down to that of the 3-RDM terms. Anderson, Shiozaki, and Booth¹⁶ reported CASPT2 and strongly contracted (SC) n-electron valence state perturbation theory (NEVPT2)^5,7,8 calculations that avoid storing the 4-RDM by combining these methods with full configuration interaction quantum Monte Carlo.¹⁷ An alternative formulation of NEVPT2 was proposed by Sharma and Chan, in which the Hylleraas functional is minimized in the space of wavefunctions described by matrix product states.¹⁸ While this approach does not require higher-body RDMs, its efficiency is contingent on the representation of the first-order correction of the wavefunction as a matrix product state with a low bond dimension. Chatterjee and Sokolov proposed an intermediate approach to partially contracted NEVPT2.¹⁹ This approach avoids the evaluation of the 4-RDM without introducing any approximation and reduces the scaling to that of computing the 3-RDM. Some NEVPT2-based approaches entirely avoid computing the 3-RDM, such as stochastic SC-NEVPT2^20,21 and the imaginary-time uncontracted NEVPT2 method.^22,23 Other approaches have been reported to reduce the cost of CASPT2 and NEVPT2 computations with large active spaces.^24–26 Nevertheless, it is important to note that not all multireference approaches require high-order RDMs. For instance, the second-order driven similarity renormalization group MRPT (DSRG-MRPT2) requires only up to the three-body RDM.²⁷ More recent developments like the multireference adiabatic connection formalism^28–30 and multiconfigurational pair density functional theory^31,32 require only the 2-RDM.

The second category of methods that can reduce the scaling with respect to the number of active orbitals uses approximate RDMs. RDMs can be expressed in terms of connected quantities (reduced density cumulants) of equal and lower rank.^33–38 If RDMs are replaced by their cumulant expansion, it is possible to truncate the rank of the cumulants while preserving the correct size extensivity properties of the full theory. Cumulant approximation has been used in the contracted Schrödinger equation method^39–42 and was later extended to the spin-free case.^43–45 Zgid et al. investigated the performance of the spin-free cumulant approximation within SC-NEVPT2.⁴⁶ They explored the possibility of approximating the 3- and 4-RDMs in two ways: (1) cumulant approximated SC-NEVPT2 (cu-SC-NEVPT2), and (2) cumulant with diagonals approximated SC-NEVPT2 (cud-SC-NEVPT2). Unfortunately, these two approximations introduce “false intruders.”⁴⁶ Recently, Guo et al.^47,48 proposed a fully internally contracted NEVPT2 and revisited cumulant approximations therein. They observed that the cumulant approximation could introduce negative eigenvalues in the Koopmans matrix. Guo et al. further developed a more stable variant named full rank NEVPT2 (FR-NEVPT2),⁴⁸ in which rank-reduction is avoided. FR-NEVPT2 was found to be more stable than NEVPT2 for incomplete active space reference wavefunctions, but new contractions involving the 5-RDM appear [$O(NA10)$ elements]. Unfortunately, cumulant approximations still result in inaccurate Koopmans matrix eigenvalues and introduce “false intruders” within FR-NEVPT2.⁴⁸ Note that the “false intruders” encountered in NEVPT2 arise due to the neglect of high-order RDMs in the energy denominators⁴⁶ and unlike classical intruders^49–51 do not reflect real near-degeneracies of the perturber and reference states. Indeed, the use of Dyall’s modified zeroth-order Hamiltonian⁵² in NEVPT2 is one of the reasons for the practical absence of intruder states in this method. One solution to the “false intruder” problem is to shift the energy denominators; Zgid et al. found that a large imaginary level shift (1.2 a.u.) is often required to obtain smooth curves,⁴⁶ leading to less accurate and, at times, even qualitatively wrong results, as observed in other contexts.^53,54 Various post-DMRG methods with approximated cumulants have been implemented. Saitow et al. proposed a DMRG-based multireference configuration interaction scheme with a cumulant approximated 4-RDM.^55,56 DMRG-CASPT2 combined with cumulant approximations has also been developed and extensively studied.^57–59 

In this work, we explore cumulant truncation in the DSRG-MRPT2 formalism to reduce the original scaling with respect to the number of active orbitals $[O(NA6)]$ and enable computations on large active spaces. Since truncation of the cumulants affects only the numerators entering the DSRG-MRPT2 amplitude and energy expressions, we expect that cumulant approximations would not introduce “false intruders” in DSRG-MRPT2. We test the two cumulant approximations considered by Zgid et al.⁴⁶ and propose two new schemes that retain a larger portion of the cumulants and scale as $O(NA4)$ and $O(NA5)$⁠. We study the energy errors these cumulant truncation schemes introduce on molecular properties, dissociation energies, and spin splittings. This paper proceeds as follows: In Sec. II A, we give an overview of the DSRG-MRPT2 method. In Sec. II B, we define four different levels of cumulant truncations and evaluate the potential for the appearance of “false intruders.” In Sec. II C, we examine the dependence of the energy of the cumulant truncation with respect to unitary transformations of the orbital basis. In Secs. III and IV, three different systems are used to examine the accuracy of cumulant-truncated DSRG-MRPT2. Finally, in Sec. V we summarize our results and discuss possible future research directions.

In this section, we summarize the MR-DSRG method and its second-order approximation. Starting with an MCSCF reference state (Ψ), the molecular spin orbitals {ϕ^p} are partitioned into core (C), active (A), and virtual (V) subsets of sizes N_C, N_A, and N_V, respectively. We use the indices m, n to denote core spin orbitals, u, v, w, x, y, z to denote active spin orbitals, and e, f to denote virtual spin orbitals. For convenience, we also define hole (H) and particle (P) spin orbital subsets as H = C ∪ A and P = A ∪ V of dimensions N_H = N_C + N_A and N_P = N_A + N_V, respectively. The labels i, j, k, l indicate hole spin orbitals and a, b, c, d indicate particle spin orbitals. For general spin orbitals (hole or particle), we use the labels p, q, r, s.

Since the DSRG-MRPT2 energy is always computed in a semicanonical basis, which uniquely defines orbitals up to an energy degeneracy, we focus on the degenerate invariance property of MR-DSRG methods.⁶⁷ This form of invariance implies that the energy remains unaffected under unitary rotations that exclusively mix degenerate orbitals within one of the three orbital subspaces. This is a departure from the traditional concept of invariance, which encompasses any unitary transformation compliant with Eq. (24). The source operator utilized in all MR-DSRG methods assumes a semicanonical orbital basis and is not invariant under general unitary transformations. However, since the source operator does not change form when two degenerate orbitals are mixed, it is possible to formulate a DSRG-MRPT2 procedure that, given any orbital basis, can reproduce the results in the semicanonical basis.^67,68 Degenerate invariance proves critical as degenerate semicanonical orbitals lack a unique definition, ensuring the DSRG-MRPT2 method consistently provides a unique solution.

When considering the degenerate invariance of the cumulant-truncated DSRG-MRPT2 methods, it is necessary to examine the consequence of cumulant truncation. In general, a unitary transformation of the form given in (24) does not preserve the structure of a cumulant in which certain blocks are set to zero. Consequently, the $λ3[d]$-, $λ3[1]$-, and $λ3[2]$-DSRG-MRPT2 methods inherently lack degenerate invariance. The exception is represented by the λ₂-DSRG-MRPT2 method, which retains degenerate invariance because it omits all three-body cumulant terms. We have tested the degenerate invariance of the cumulant-truncated DSRG-MRPT2 on two systems. The first one is the lowest singlet state of N₂ in the dissociation limit, where we use a full valence active space (ten electrons in eight orbitals) and test invariance by localizing the 2s and 2p orbitals. For N₂, we find that λ₃ is zero in the dissociation limit, and consequently, all truncated schemes yield the same energy as the DSRG-MRPT2 and are invariant with respect to unitary transformations of degenerate orbitals (tested by localizing the symmetry-adapted orbitals via the Pipek–Mezey procedure⁶⁹). The behavior observed for N₂ is in line with a previous report by Hanauer and Köhn,⁷⁰ which noted that the odd-rank cumulants of N₂ vanish in the dissociation limit. They explained the phenomenon with their statistical interpretation of cumulants: since all active orbitals are degenerate, the electron distribution is symmetric, resulting in all odd cumulants being equal to zero.⁷⁰ The second system we examine is a point along the Be + H₂ insertion reaction⁷¹ where the HOMO and LUMO become degenerate. Computations on this system used the cc-pVDZ basis set^72,73 and correlated all electrons. For Be + H₂, at the point of orbital degeneracy (ϵ_HOMO = ϵ_LUMO) obtained using a minimal active space (two electrons in two orbitals) also gives λ₃ = 0, and consequently, the energy is degenerate invariant for all cumulant-truncated DSRG-MRPT2 schemes. If we expand the active space to include the LUMO + 1 and LUMO + 2 orbitals, the three-body cumulant is no longer zero at the point of HOMO–LUMO orbital degeneracy. The impact of the various truncation schemes evaluated by computing the energy difference for the symmetry-adapted and Pipek–Mezey localized orbitals is shown in Table II. Both the DSRG-MRPT2 and λ₂-DSRG-MRPT2 are degenerate invariant (in agreement with our analysis), while the other truncated schemes show energy differences of various degrees. The largest impact of the lack of degenerate invariance is found for the $λ3[d]$-approximation (55 μE_h), a large portion of the three-body cumulant contribution to the correlation energy (75 μE_h). Including more elements of the three-body cumulant further reduces the energy difference between the symmetry-adapted and localized results by 8.2 (⁠$λ3[1]$-) and −0.3 (⁠$λ3[2]$-) μE_h.

As part of this work, we have developed pilot implementations of the λ₂-, $λ3[d]$-, $λ3[1]$-, and $λ3[2]$-DSRG-MRPT2 truncation schemes using Forte,⁷⁴ an open-source plugin for the Psi4 ab initio quantum chemistry package.⁷⁵ To thoroughly study the accuracy of the cumulant-truncated DSRG-MRPT2 methods, we applied them to compute properties of several systems: the singlet–triplet gap of CH₂ $(A11−B13)$⁠, the potential energy curves of N₂ (both the $XΣg+1$ state and $AΣu+3$ state), and the adiabatic singlet–triplet splitting of n-acenes (n = 2–4). The CH₂ and N₂ calculations used the flow parameter value $s=1Eh−2$⁠, following past work.²⁷ To be consistent with our previous paper,⁷⁶ we used the value $s=0.5Eh−2$ for the n-acenes calculations.

For CH₂, state-specific CAS self-consistent field (CASSCF) computations were performed to obtain the optimized molecular orbitals and reference state used in DSRG-MRPT2. The geometry and the basis set (double zeta with Cartesian H 2p polarization function and C 3d polarization function) for the CH₂ calculations were taken from Ref. 77. The singlet and triplet calculations used different C 3d polarization functions.⁷⁷ The CASSCF active space was taken from Ref. 46 and consisted of the 2a₁, 3a₁, 4a₁, 1b₁, 2b₁, and 1b₂ orbitals in C_2v symmetry. The 1a₁ core orbital was treated as a restricted doubly occupied orbital in the CASSCF procedure (optimized but with fixed occupation) and frozen (i.e., not included) in the correlated computations.

For N₂, we used the cc-pVQZ basis.^72,78 Optimized molecular orbitals were obtained via state-specific CASSCF, including all 2s and 2p orbitals in the active space (ten electrons and eight orbitals). The 1s core orbitals were treated as a restricted doubly occupied orbitals and, for consistency with previous work, were included in the correlated computations.⁴⁶ 

We first look at the accuracy of various methods. The predicted adiabatic singlet–triplet gaps of CH₂ are shown in Table III. As can be seen in Table III, the FCI error in DSRG-MRPT2 is larger (ΔE_FCI = 4.05 kcal mol⁻¹) than SC-NEVPT2 (ΔE_FCI = 1.72 kcal mol⁻¹) and CASPT2 (ΔE_FCI = 3.66 kcal mol⁻¹). CASSCF gives the smallest error (ΔE_FCI = 0.87 kcal mol⁻¹) among all the computations, likely due to error cancellation since it gives the largest absolute energy error for both singlet and triplet states. Next, we examine the error caused by cumulant truncation. Interestingly, the performance of SC-NEVPT2 and DSRG-MRPT2 differs significantly when applying the crudest approximation: The singlet–triplet gap cumulant error for cu-SC-NEVPT2 (ΔE_cu = 10.81 kcal mol⁻¹) is considerably larger than that for λ₂-DSRG-MRPT2 (ΔE_cu = 1.63 kcal mol⁻¹), which is smaller than the intrinsic DSRG-MRPT2 error (ΔE_cu = 4.05 kcal mol⁻¹). Higher-level cumulant truncation methods do not show this significant difference, e.g., the cumulant error for both cud-SC-NEVPT2 (ΔE_cu = 0.63 kcal mol⁻¹) and $λ3[d]$-DSRG-MRPT2 (ΔE_cu = 1.20 kcal mol⁻¹). The smaller sensitivity with respect to cumulant truncation of the DSRG-MRPT2, when compared to SC-NEVPT2, is consistent with the different ways the density cumulants enter into these methods, with the former depending only linearly with respect to the cumulants.

To examine in more detail the connection between cumulant structure and the accuracy of the approximate methods, in Fig. 1, we plot the nonredundant tensor elements of the three-body reduced density cumulant $(λuvwxyz)$ retained in the parent and truncated DSRG-MRPT2 schemes. This plot represents $λuvwxyz$ as a matrix indexed by the upper and lower tensor indices (compounded) for the ααα and ααβ spin cases for the singlet and triplet states. The βββ spin case is the same as the ααα case, and the αββ case is equivalent to the ααβ case, with elements ordered differently. See Figs. S1 and S2 in the supplementary material for the βββ and the αββ cases. λ₂-DSRG-MRPT2 is not included in the figure since all tensor elements are zero. These plots show that the largest elements are included in $λ3[2]$-DSRG-MRPT2. Still, many are not included in the simpler truncation schemes, explaining the small cumulant error introduced by $λ3[2]$-DSRG-MRPT2 for the CH₂ singlet–triplet gap. A quantitative argument can be made by analyzing the total number of significant elements of the three-body reduced density cumulant $(|λuvwxyz|>10−4)$⁠. In Table IV, we report this number for the full and truncated DSRG-MRPT2. In the full DSRG-MRPT2, only 3068 elements of $λuvwxyz$ are significant (less than 2% of the total number). The $λ3[2]$-DSRG-MRPT2 can retain a large fraction of these important elements (2452) even though it requires storing fewer elements [69% of the total, with $O(NA5)$ storage cost]. This truncation can therefore lead to significant storage and computational cost savings for large active spaces, bringing a potential reduction in cost compared to conventional computations of $O(NA)$⁠.

To evaluate the accuracy of the truncated DSRG-MRPT2 methods in describing potential energy curves, we computed the lowest singlet $(XΣg+1)$ and the lowest triplet $(AΣu+3)$ states of N₂ as a function of bond length. In cumulant-truncated SC-NEVPT2,⁴⁶ “false intruders” arise at many geometries along the potential energy curve for both bare (i.e., without level shift) cu-SC-NEVPT2 and bare cud-SC-NEVPT2. A large imaginary level shift (1.2 a.u.) is required to correct those singularities arising from the cumulant approximation. Therefore, this stricter test case is useful to assess if cumulant truncations introduce the intruder state problem in DSRG-MRPT2 and how they affect the potential energy curve.

Potential energy curves and the corresponding cumulant errors for the $XΣg+1$ and $AΣu+3$ states computed with the parent and cumulant-truncated DSRG-MRPT2 schemes are shown in Fig. 2. For all cumulant-truncated DSRG methods, the potential energy curve is everywhere continuous and smooth. Cumulant truncations are expected to impact the absolute energy the most when the neglected parts of the cumulants are the largest, which always happens in the recoupling region. The largest cumulant errors for the singlet and triplet states observed for λ₂-DSRG-MRPT2 (ΔE_cu = 11.20 and 3.05 kcal mol⁻¹) gradually decrease going from the $λ3[d]$ (ΔE_cu = 6.13 and 3.01 kcal mol⁻¹) to the $λ3[2]$ (ΔE_cu = 0.19 and 0.16 kcal mol⁻¹) approximations. It is encouraging to see that the accuracy of the absolute energy can be systematically improved by using more accurate cumulant truncation schemes and that the cumulant error for $λ3[2]$-DSRG-MRPT2 is always less than 0.2 kcal mol⁻¹. Cumulant errors in the absolute energy are smallest at the stretched geometry. We attribute this to the small magnitude of the three-body reduced density cumulant elements, as measured by the sum of the tensor Frobenius norm (F-norm) of the four different spin cases (ααα, βββ, ααβ, αββ). The F-norm is negligible at the stretched geometry for both the $XΣg+1$ (0.023) and the $AΣu+3$ (0.013) states but not for the other geometries (e.g., 0.437 for the $XΣg+1$ state in the recoupling region and 1.277 for the $AΣu+3$ state at the equilibrium geometry). As a consequence, at the stretched geometry, even the crudest cumulant truncation scheme, λ₂-DSRG-MRPT2, introduces relatively small absolute errors (0.25 and 0.08 kcal mol⁻¹ for the $XΣg+1$ and $AΣu+3$ states, respectively).

To evaluate how many significant tensor elements are used by each truncation method, we plot the nonredundant elements in the three-body reduced density cumulant as a matrix for both the $XΣg+1$ and the $AΣu+3$ states of N₂ using the method described in the CH₂ section. Figures S3–S10 (see the supplementary material) show the nonredundant tensor elements for the $XΣg+1$ and the $AΣu+3$ states of N₂ at the equilibrium geometry (1.10 and 1.29 Å for the $XΣg+1$ and the $AΣu+3$ states, respectively), in the recoupling region (1.75 Å), and at the stretched geometry (3.00 Å). The results are similar to those found for CH₂: The largest elements are included in $λ3[2]$-DSRG-MRPT2, but many are not included in the simpler truncation schemes. Table V shows the total number of significant elements $(|λuvwxyz|>10−4)$ and total number of elements of λ₃ entering in the $XΣg+1$ and the $AΣu+3$ states for DSRG-MRPT2 computations at different geometries. Unsurprisingly, most significant elements are missing in $λ3[d]$-DSRG-MRPT2 and $λ3[1]$-DSRG-MRPT2, while $λ3[2]$-DSRG-MRPT2 (storing 59% of total elements) can include more than 70% of them in most cases, except for the $AΣu+3$ state stretched geometry calculation (66%).

We now turn to the impact of our cumulant truncation schemes on the dissociation energy and spectroscopic constants. The results are shown in Table VI. Let us first consider the effect of cumulant truncation on the ground state D_e as predicted by SC-NEVPT2. The error with respect to experiment without truncation is 0.14 eV, and the cu approximation scheme with a 1.2 a.u. imaginary level shift reduces the error to only 0.03 eV. As shown in Ref. 46, this improvement mainly comes from error cancellation resulting from a significant vertical shift of the potential energy curve in the dissociation region. The untruncated DSRG-MRPT2 displays a larger error than SC-NEVPT2 and CASPT2 in the dissociation energy (0.62 eV). However, as shown in Fig. 2, results are consistent across the various cumulant truncation schemes, smoothly converging to the untruncated values as the cumulant truncation is improved.

Spectroscopic constants for the parent and cumulant-truncated DSRG-MRPT2 were obtained by using DSRG-MRPT2 analytic energy gradients⁹⁷ for the $XΣg+1$ and $AΣu+3$ states of N₂. The smoothness of potential energy curves of the parent and cumulant-truncated DSRG-MRPT2 theories is critical for obtaining meaningful properties. The maximum absolute differences for the bond length (r_e) and the harmonic vibrational frequency (ω_e) for all DSRG-MRPT2 methods examined are 0.0020 Å and 14 cm⁻¹ for the $XΣg+1$ state and 0.0027 Å and 5 cm⁻¹ for the $AΣu+3$ state. Our results suggest that the crudest cumulant-truncated DSRG-MPRT2 scheme, λ₂-DSRG-MRPT2, could already be an efficient method for geometry optimization and frequency calculations. In Table VII, we further compared the predicted adiabatic singlet–triplet gap of N₂ using the parent and cumulant-truncated DSRG-MRPT2 methods with the experimental value. All energies reported do not include zero-point energy corrections. While the gap predicted by the DSRG-MRPT2 deviates from experiment by ∼0.3 eV, all the cumulant truncation schemes that include part of the three-body density cumulant are within 0.01 eV from the full DSRG-MRPT2.

In the last test case, we study the performance of the cumulant truncations for the case of large active spaces (up to 18 electrons and 18 orbitals). Oligoacenes, or n-acenes (where n is the number of linearly fused benzene rings), have been studied extensively by chemists because of their semiconducting and optical properties.^99–109 In particular, their singlet–triplet splittings have been of interest to theorists.^{79,99,110–122} Table VIII provides adiabatic singlet–triplet splittings for the first three members of the acene series using both unrelaxed and relaxed references. Our results show that all cumulant-truncated DSRG-MRPT2 schemes, including the crudest λ₂-DSRG-MRPT2 method, can preserve the accuracy of the parent method across different levels of basis sets, with a maximum cumulant error of 0.2 kcal mol⁻¹ for the energy gap. Contrary to expectations, we observed no specific trend for the performance of different cumulant truncations in the singlet–triplet splittings. To better understand the origin of this phenomenon, we compared the absolute energies of all computations. Additional data, including absolute energies of the singlet and triplet states of the acenes, are provided in Table S1 in the supplementary material. Figure 3 shows cumulant errors in the absolute energy of the acenes for computations with the cc-pVQZ basis set. Consistent with our previous test cases, these data show that the absolute energy can be systematically improved using higher-level cumulant truncations and that near-perfect cancellation of the cumulant errors happens across the truncation schemes. As the most promising cumulant-truncated DSRG-MRPT2 method, $λ3[2]$-DSRG-MRPT2 shows significant computational savings by storing 50.8%, 39.7%, and 32.5% of total elements for 2-acene, 3-acene, and 4-acene, respectively. Table S2 compares the computed adiabatic singlet–triplet splittings of the acene series to experimental values using the cc-pVQZ basis set and relaxed reference. These data include ZPVE corrections taken from our previous paper⁷⁶ (−3.4, −2.3, and −1.8 kcal mol⁻¹ for 2-, 3-, and 4-acene, respectively). In this case, the DSRG-MRPT2 underestimates the experimental gaps, a phenomenon already observed in our previous study.⁷⁶ 

In this work, we have explored the possibility of truncating the three-body reduced density cumulant in the second-order multireference perturbation theory based on the driven similarity renormalization group (DSRG-MRPT2). We tested the two cumulant truncations examined in the work of Zgid et al.⁴⁶ (λ₂-DSRG-MRPT2 and $λ3[d]$-DSRG-MRPT2) and also proposed two new cumulant truncations (⁠$λ3[1]$-DSRG-MRPT2 and $λ3[2]$-DSRG-MRPT2). In contrast to SC-NEVPT2, which requires a large imaginary level shift to deal with “false intruders” along potential energy curves, we demonstrated that cumulant-truncated versions of the DSRG-MRPT2 preserve the smoothness of the potential energy curves. Our results also show that the accuracy of cumulant-truncated DSRG-MRPT2 can be systematically improved with higher-level truncations, and $λ3[2]$-DSRG-MRPT2 is particularly promising since the error introduced from this truncation scheme in energy differences is found to generally be less than 0.5 kcal mol⁻¹. At the same time, the memory cost decreases from $O(NA6)$ to $O(NA5)$⁠. Simpler approximation methods could be applied to reduce the memory and computational cost of energy gap calculations due to near-perfect error cancellations observed when the two states have similar molecular geometries. A careful analysis is necessary when applying cumulant truncation methods to evaluate potential energy surfaces as less expensive schemes may introduce non-negligible errors. The cumulant truncation schemes discussed in this work are particularly useful to reduce the storage and computational cost in applications with large active spaces (N_A > 10–20). Further extension of the DSRG-MRPT2 that will enable applications to large systems will also need to address the scaling of the DSRG-MRPT2 with respect to the number of core and virtual orbitals via methods that can take advantage of the locality of dynamical electron correlation. Combining these two strategies is expected to open the door to future applications of multireference methods and Hamiltonian downfolding for quantum computing¹¹⁰ to large systems with large active spaces.

See the supplementary material for (1) plots of the nonredundant tensor elements of the three-body cumulants of CH2 and N2 and (2) adiabatic singlet–triplet splittings of the oligoacenes corrected for zero-point vibrational energy effects.

This work was supported by the U.S. National Science Foundation under Award No. CHE-2312105 and by a Camille Dreyfus Teacher-Scholar Award (No. TC-18-045).

The authors have no conflicts to disclose.

Shuhang Li: Data curation (lead); Investigation (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Jonathon P. Misiewicz: Supervision (equal); Writing – review & editing (equal). Francesco A. Evangelista: Conceptualization (lead); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.