Efficient Adiabatic Connection Approach for Strongly Correlated Systems: Application to Singlet-Triplet Gaps of Biradicals.

Strong electron correlation can be captured with multireference wave function methods, but an accurate description of the electronic structure requires accounting for the dynamic correlation, which they miss. In this work, a new approach for the correlation energy based on the adiabatic connection (AC) is proposed. The ACn method accounts for terms up to order n in the coupling constant, and it is size-consistent and free from instabilities. It employs the multireference random phase approximation and the Cholesky decomposition technique, leading to a computational cost growing with the fifth power of the system size. Because of the dependence on only one- and two-electron reduced density matrices, ACn is more efficient than existing ab initio multireference dynamic correlation methods. ACn affords excellent results for singlet-triplet gaps of challenging organic biradicals. The development presented in this work opens new perspectives for accurate calculations of systems with dozens of strongly correlated electrons.

Electron correlation energy is defined with respect to the energy of a model (a reference) used to describe a given system. In other words, given a Hamiltonian , if Ψref is the reference wave function and Eref is the corresponding energy, i.e.,then electron correlation comprises all electron interaction effects not accounted for by the chosen model, and the correlation energy pertains to the energy errorcomputed with respect to the exact energy Eexact (an eigenvalue of the Hamiltonian ). Strongly correlated molecular systems require model wave functions consisting of multiple configurations to capture static correlation effects. The complete active space (CAS) method assumes the selection of a number of (active) electrons and orbitals crucial to the static correlation followed by exact diagonalization in the active orbital subspace.[1,2] The CAS model is a base of the CASSCF wave function and is also frequently employed in density matrix renormalization group (DMRG) calculations. The DMRG method is one of the most promising tools for strongly correlated molecules[3−7] because of its favorable scaling, which enables the handling of much more extensive active spaces than CASSCF allows. The reference energy, Eref in eq , for all CAS-based methods does not include a substantial portion of the electron correlation, called dynamic correlation, Ecorr in eq . Even the inclusion of dozens of orbitals in the active space is not sufficient to achieve a reliable description, and the necessity to recover dynamic correlation remains the major challenge of DMRG.[6] Although there exist many post-CAS methods aimed at including dynamic correlation (e.g. ref (7)), none are satisfactory because of the limitations in both accuracy and efficiency. In particular, perturbation-theory-based approximations may suffer from the lack of size-consistency, intruder states, or the unbalanced treatment of closed- and open-shell systems, which must be cured by level-shifting.[8] The limitation of PT2 when combined with DMRG is the high scaling with the number of active orbitals resulting from the treatment of three- and four-electron reduced density matrices (RDMs). Efforts to reduce the cost of handling high-order RDMs in NEVPT2 are worth noticing. These include the stochastic strongly contracted scheme,[9,10] employing the cumulant expansion[11] or prescreening techniques.[12] However, the improved efficiency may come at the cost of introducing additional intruder states.[13] Alternative approaches for molecular systems with strongly correlated electrons, which might also be viewed as CAS-like methods, are represented by embedding schemes. These comprise the self-energy embedding theory,[14] active-space embedding,[15] or subsystem embedding subalgebras[16] leading to the active-space coupled-cluster downfolding techniques.[17]

The goal of this work is to address the challenge of recovering dynamic correlation and proposing an efficient and reliable computational method applicable to large active spaces. The presented approach builds upon the adiabatic connection formalism first introduced in the framework of Kohn–Sham DFT[18−21] and recently formulated for CAS models.[22,23]

Although the following discussion will pertain to a ground-state energy, the presented formalism is general and can be directly applied to higher states. Derivation of the formula for the correlation energy in the adiabatic connection (AC) formalism begins with assuming a model Hamiltonian (typically the electron–electron interaction is either reduced or removed from ) such that the reference function Ψref is its eigenfunctionThe AC Hamiltonian is introduced as a combination of and a scaled complementary operator The eigenequation for readswhere index ν pertains to the νth electronic state. The role of the coupling parameter α is to adiabatically turn on full electron correlation by varying α from 0 to 1. Namely, at α = 0 electron interaction is reduced according to the assumed model and the reference wave function is obtained as The α = 1 limit corresponds to electrons interacting at their full strength so that both the exact energy and wave function are obtainedExploiting the Hellmann–Feynman theorem , satisfied for α ∈ [0, 1], it is straightforward to show that the correlation energy, eq , is given exactly as

The choice for the Hamiltonian depends on the reference wave function. Our interest is in multireference CAS-based models which assume partitioning orbitals into sets of inactive (fully occupied), active (fractionally occupied), and virtual (unoccupied) orbitals and constructing Ψref as an antisymmetrized product of a single determinant comprising inactive orbitals and a multiconfigurational function utilizing active orbitals. Thus, we represent as a sum of group Hamiltonians (22,24)where I corresponds to an inactive, active, or virtual group and consists of one- and two-particle operatorsNotice that denotes a two-electron integral in the x1x2x1x2 convention and the effective one-electron Hamiltonian heff is a sum of kinetic and electron-nuclei operators and the self-consistent field interaction of orbitals in group I with the other groups (second term in eq ). Throughout the letter, it is assumed that indices p, q, r, and s denote natural spin orbitals of the reference (α = 0) model and are the corresponding natural occupation numbers. For this choice of , the α-dependent integrand in the correlation energy expression, eq , includes, among others, one-electron terms depending on the difference between 1-RDM at α = 1 and the reference term, γα – γα=0 = γα – γref. Such terms are set to 0 under the assumption that for the properly chosen multireference wave function for a strongly correlated system, the variation of γα with α can be ignored.

As has been shown in refs (22) and (23) and also in the Supporting Information (SI), choosing as a group Hamiltonian and assuming that 1-RDM stays constant with α turn eq into the following AC correlation energy expressionwhere γ0ν,α are one-electron transition reduced-density matrices (1-TRDM)It is important to notice a prime in the AC formula in eq , which indicates that terms pertaining to pqrs belonging to the same group are excluded. This implies that electron correlation already accounted for by the active-orbitals component of the reference wave function is not counted twice in .

We now briefly recapitulate developments presented in our earlier works[22,23,25,26] leading to approximate correlation energy methods called AC and AC0. To formulate a working expression for the AC correlation energy, we have used Rowe’s equation of motion[27,28] in the particle-hole RPA approximation, where the excitation operator generating a state ν, , is approximated by single excitation operators as . To distinguish this approximation from the conventional RPA,[27,29−31] which assumes a single determinant as a reference, we used the ph-RPA equationswhich are introduced for a general, multiconfigurational reference: the extended RPA (ERPA).[28,32] The ERPA equations have been written for the AC Hamiltonian (eq ) leading to defined asExplicit expressions of in terms of 1- and 2-RDMs are presented in the SI. Both and are symmetric and positive-definite at α = 0 and 1 for the Hellmann–Feynman reference wave function Ψref. Because the coupling constant dependence is passed to ERPA equations only via AC Hamiltonian , the matrices are linear in α, i.e.,In the ERPA model,[33] the α-dependent 1-TRDMs (eq ) are given by the eigenvectors as , which allows one to turn eq into a spin-free formula[25]where . Equations and 19 form the basis for practical correlation energy calculations. This, however, requires solving the ERPA problem which formally scales with the sixth power of the system size. In addition, using the reference wave function in which the choice of the active orbitals is not optimal could lead to developing instability in the ERPA problem for α ≫ 0.[25] To lower the computational cost and avoid potential instabilities, we introduced an AC0 variant, assuming linearization of the integrand in eq , namely, using , by keeping the linear terms in α and carrying out the α integration,[25]The low computational cost of AC0 stems from the fact that ERPA equations must be solved only at α = 0, and for this value of the coupling constant, the matrices are block diagonal. The largest block has dimensions of (Nact denotes the number of active orbitals), so the cost of its diagonalization is marginal even for dozens of active orbitals.

Despite the fact that encouraging results have been obtained with AC0 when combined with CASSCF[25,26,34] or DMRG,[35] α integration should in principle account for correlation more accurately than AC0. It is thus desirable to develop an AC method which on the one hand is exact at all orders of α and on the other avoids solving the expensive ERPA problem. Ideally, such a method would be free of potential instabilities that might occur when α approaches 1. A novel AC method satisfying all of the requirements is presented in this work.

Let us use the integral identity = 1 to express the AC correlation energy by means of the α-dependent dynamic density–density response matrix.[36] This can be attained by employing the relationsin eq , resulting in the formulawhereand the prime in eq indicates that when taking a product of matrices C and g, terms pqrs ∈ active are excluded. By using spectral representations of the matrices and in terms of the ERPA eigenvectors,[37] it is straightforward to show that the dynamic linear response matrix Cα(ω) follows from the linear equation given as (see the SI for details)To reduce the computational cost of solving eq , we introduce a decomposition of the modified two-electron integrals gwhere D are the natural-orbital-transformed Cholesky vectors of the Coulomb matrix multiplied by factors , cf. eq . We expand Cα(ω) at α = 0and solve eq iteratively in the reduced space by retrieving, in the nth iteration, the nth-order correction C( projected onto the space spanned by NChol transformed Cholesky vectors To account for the prime (exclusion of terms for all-active indices pqrs) in the AC correlation energy, eq , we define the auxiliary matrices of the transformed Cholesky vectors asAssuming an expansion of the response matrix Cα(ω), cf. eq , up to nth order in α and employing the Cholesky decomposition of integrals, eq , together with matrices D1 and D2 in eq leads to a new AC formula for the correlation energy readingThe matrices defined ashave dimensions of M2 × NChol, which are reduced comparing to the M2 × M2 dimensions of C(ω)( because by construction the number of Cholesky vectors is 1 order of magnitude smaller than M2, i.e., NChol ≈ M. Employing the linearity in α of the matrices , cf. eq , in eq one finds the following recursive formulas for the nth-order term where the required matrices are given by the ERPA matrices and (see the SI for their explicit forms in terms of 1- and 2-RDMs)

The correlation energy expression in eq together with the recursive relation in eqs –34 is the central achievement of this work. It allows one to compute the correlation energy for strongly correlated systems at the cost of scaling with only the fifth power of the system size. All matrix operations scale as M4NChol down from M6 scaling of the original ERPA problem in eq . Notice that the cost of computing the Λ(ω) matrix is marginal because the inverted matrix is block diagonal with the largest block having dimensions of .

By setting the maximum order of expansion of the response matrix C(ω) in eq to 1, the correlation energy AC becomes equivalent to the AC0 approximation, cf. eq . In the limit of n → ∞, the value approaches the AC energy given according to the formula in eq if the Taylor series is convergent. Numerically, this equality requires sufficient accuracy both in the frequency integration and in the Cholesky decomposition of two-electron integrals.

Going beyond the first-order terms in the coupling constant is potentially beneficial because higher orders gain importance as α approaches 1. Higher-order contributions are effectively maximized if the AC integrand Wα in eq , , is linearly extrapolated from Wα=1 to the exact limit Wα=0 = 0. Such an extrapolation method leading to the approximation has already been proposed in ref (22). If it is used together with the formula in eq , the expansion shown in eq , and the Cholesky decomposition of two-electron integrals, then one obtains the formulawhich will be denoted as AC1. Notice that the frequency-integrated kth-order term in eq contributes to the correlation energy by a factor of (k + 1)/2 greater than its counterpart in the expression given in eq .

The Cholesky decomposition of the Coulomb integrals matrix in the AO basis was carried out using a modified program originally used in refs (38) and (39). The implementation was carried out according to ref (40). The Cholesky vectors in the AO basis, R, were generated until the satisfaction of the trace condition . The convergence threshold was previously tested as a part of the default set of numerical thresholds in Table of ref (38).

Table 1

ST Gaps (ET – ES), Mean Errors (ME), Mean Unsigned Errors (MUE), and Standard Deviations (std dev) Computed with Respect to CC3 Reference Dataa

molecule	T state	CASSCFb	AC1n	ACn	AC0	NEVPT2c	CASPT2d	CC3d
ethene	13B1u	3.78	4.53	4.56	4.69	4.60	4.60	4.48
E-butadiene	13Bu	2.77	3.44	3.43	3.46	3.38	3.34	3.32
all-E-hexatriene	13Ag	2.66	2.83	2.81	2.80	2.73	2.71	2.69
all-E-octatetraene	13Bu	2.25	2.46	2.43	2.39	2.32	2.33	2.30
cyclopropene	13B2	3.78	4.42	4.44	4.56	4.56	4.35	4.34
cyclopentadiene	13B2	2.75	3.34	3.34	3.37	3.32	3.28	3.25
norbornadiene	13A2	3.07	3.92	3.89	3.86	3.79	3.75	3.72
benzene	13B1u	3.74	4.17	4.21	4.37	4.32	4.17	4.12
naphtalene	13B2u	2.93	3.19	3.21	3.29	3.26	3.20	3.11
furan	13B2	3.54	4.09	4.16	4.30	4.33	4.17	4.48
pyrrole	13B2	3.95	4.47	4.52	4.67	4.73	4.52	4.48
imidazole	13A′	4.42	4.70	4.74	4.85	4.77	4.65	4.69
pyridine	13A1	3.81	4.28	4.34	4.53	4.47	4.27	4.25
s-tetrazine	13B3u	2.43	2.27	2.05	1.51	1.64	1.56	1.89
formaldehyde	13A2	3.32	3.80	3.74	3.77	3.75	3.58	3.55
acetone	13A2	4.17	4.27	4.29	4.90	4.10	4.08	4.05
formamide	13A″	4.72	5.31	5.47	5.60	5.64	5.40	5.36
acetamide	13A″	4.77	5.46	5.57	5.73	5.52	5.53	5.42
propanamide	13A″	4.79	5.51	5.61	5.80	5.54	5.44	5.45
ME		– 0.38	0.08	0.10	0.18	0.10	0.00	-
MUE		0.45	0.13	0.13	0.24	0.14	0.07	-
std dev		0.35	0.15	0.11	0.23	0.13	0.12	-

All values are in eV.

Active spaces from ref (42).

Results from ref (47).

Results from ref (42).

For the ω integration in the AC correlation energy, we have used a modified Gauss–Legendre quadrature as described in ref (41). With the 18-point grid, the accuracy of the absolute value of energy achieves 10–2 mHa, which results in 10–2 eV accuracy in the singlet–triplet (ST) gaps.

To assess the accuracy of the AC approaches, we have applied them to two benchmark data sets of singlet–triplet energy gaps: the single-reference system set of Schreiber et al.[42] and the multireference organic biradicals studied by Stoneburner et al.[43] In the single-reference data set, we employed the TZVP[44] basis set and compared our data against the CC3[42] results. The aug-cc-pVTZ basis and the doubly electron-attached (DEA) equation-of-motion (EOM) coupled-cluster (CC) 4-particle–2-hole (4p–2h) reference[43] were used for the biradicals. All CASSCF calculations were performed in the Molpro[45] program. All AC methods were implemented in the GammCor program.[46]

Computing the correlation energy with the AC method requires either fixing the maximum order of expansion with respect to the coupling constant, n in eq , or continuing the expansion until a prescribed convergence threshold is met. The advantage of the former strategy is that size consistency is strictly preserved. For each system, we found that the AC correlation energy converges with n for the chosen active space. Typical convergence behavior for the singlet, triplet, and ST energies is presented in Figure . It can be seen that for n = 3 the AC ST gap deviates by only 0.01 eV from the AC value, computed using eq . For all other biradicals and single-reference systems, we found that setting n = 10 in eq is sufficient to converge ST gaps within 10–2 eV, thus n = 10 has been set for all systems.

Figure 1

Differences in AC and AC correlation energies for singlet (S) and triplet (T) states (left axis) and ST gaps (right axis) as a function of n for the C4H2-1,3-(CH2)2 biradical. Notice that black markers overlap with the red ones.

In Table , we present ST gaps for the subset of the ref (42) data set. The CASSCF method predicts ST gaps that are too narrow, with the mean error approaching −0.4 eV, which results from the unbalanced treatment of closed-shell singlet and open-shell triplet states. The addition of correlation energy using the adiabatic connection greatly reduces the errors. The mean unsigned error (MUE) of AC0 amounts to 0.24 eV. The performance is further improved by AC, which affords MUE of 0.13 eV. Maximizing the contribution from higher-order terms in α, attained in AC1, leads to ST gaps of the same unsigned error as that of AC. Noticeable, the signed error is reduced, which indicates that higher-order terms play a more important role in the open-shell states than in the closed-shell states. The accuracy of AC is on a par with NEVPT2 and only slightly worse than the best CASPT2 estimations from ref (42). The standard deviation of AC0, amounting to 0.23 eV, is reduced to 0.11 eV by AC, which parallels the standard deviation of the perturbation methods.

In ref (43), the systematic design of active spaces for biradicals based on the correlated participating orbital (CPO) scheme[48] is presented. Here, we take a different approach and identify the most appropriate CASs by means of single-orbital entropies and two-orbital mutual information.[49−51]

Figure shows the correlation measures for singlet and triplet states of prototypical biradicals, C4H4 and C5, obtained with, respectively, CAS(20,22) and CAS(14,16) active spaces (cf., the description in the SI). We observe that the π orbitals of C4H4 and C5 are well separated from the others in terms of their single-orbital entropies (s > 0.19; see the SI) and represent a natural choice of the active space selection. The largest values of s correspond to the singly occupied frontier orbitals in the singlet states. These orbital pairs also possess the largest values of mutual information, which stems from the strong correlation of the frontier orbitals due to the singlet-type coupling of these open shells. Notice that both single-orbital entropies and mutual information from the singly occupied orbitals are much lower in the case of the triplet states. This is due to the fact that the triplet states were calculated as high-spin projections and thus can be qualitatively described with a single determinant. However, when analyzing the triplet states, one can see that all of the π orbitals have similar values of their single-orbital entropies and that CAS(2,4) (the nCPO active space in ref (43)) is not a reasonable choice. In fact, for this imbalanced active space, we have experienced divergence of the AC series (last entry in Table in the SI).

Figure 2

DMRG mutual information (colored edges) and single-orbital entropies (colored vertices) of C4H4 and C5 for the lowest singlet and triplet states. Numbers in the graphs correspond to indices of the DMRG-SCF (C4H4) and CASSCF (C5) natural orbitals presented together with their occupation numbers in the SI. Blue circles represent the π orbitals with s > 0.19.

The analysis of mutual information and single-orbital entropies of prototypical biradicals has allowed us to define optimal active spaces: CAS(4,4) for C4H4, C4H3NH2, C4H3CHO, and C4H2NH2(CHO), CAS(4,5) for C5, and CAS(6,6) for C4H2-1,2-(CH2)2, and C4H2-1,3-(CH2)2. The choice of the orbitals in CAS is therefore such that all valence π orbitals on, or adjacent to, the carbon the carbon ring are included and only the mostly correlated orbitals, with occupancies in the range of (0.05, 1.95), enter the active space. The chosen active spaces are close to the πCPO scheme considered in ref (43), with the difference that nearly unoccupied orbitals in πCPO, shown to be uncorrelated according to our mutual information analysis, are excluded.

Similar to the single-reference case, the performance of the CASSCF method for the ST gaps in biradicals is seriously affected by the lack of dynamic correlation (Table ). Even though the CASSCF gaps of three systems (C5, 1,2- and 1,3-isomers) are in error by only 0.1 eV, the overall MUE is as large as 0.20 eV and the mean average unsigned percentage error (MU%E) exceeds 100%. The AC0 method overcompensates for the errors in CASSCF. For biradicals 1, 3, and 4, the excessive reduction of the ST gaps by AC0 results in a wrong ordering of states. Both AC and AC1 approaches capture correlation at high orders of α and greatly improve on AC0. The ordering of states is correct and the average error falls below 0.10 eV as compared to the 0.16 eV error in AC0. AC1 performs slightly better than AC in terms of MUEs, with the errors being 0.06 and 0.08 eV, respectively, and significantly better in terms of percentage errors. The improved MU%E of AC1 (14 vs 26%) is due to the good performance of this method for small gaps (systems 1, 3, and 4). These excellent results imply a crucial role of the higher-order terms in AC which should enter the correlation energy with high weights.

Table 2

ST Gaps (ET – ES) in eV and Errors with Respect to DEA-EOMCC[4p-2h][43] (ref) Predictionsa

system	CASSCF	AC1n	ACn	AC0	ftPBEb	RASPT2c	ref
1	0.44	0.18	0.13	0.00	0.11	0.19	0.18
2	–0.70	–0.71	–0.77	–0.86	–0.64	–0.65	–0.60
3	0.39	0.12	0.07	–0.08	0.03	0.11	0.12
4	0.41	0.17	0.12	–0.02	0.09	0.16	0.16
5	0.60	0.39	0.37	0.14	0.45	0.32	0.25
6	3.33	3.44	3.44	3.46	3.34	3.27	3.37
7	–0.92	–0.89	–0.90	–0.92	–0.66	–0.80	–0.80
ME	0.13	0.00	–0.03	–0.14	0.01	–0.01
MUE	0.20	0.06	0.08	0.16	0.09	0.04
MU%E	102.3	13.7	26.1	69.6	37.6	7.2
Std. Dev.	0.20	0.09	0.10	0.11	0.12	0.05

Labels: (1) C4H4, (2) C5H5+, (3) C4H3NH2, (4) C4H3CHO, (5) C4H2NH2(CHO), (6) C4H2-1,2-(CH2)2, and (7) C4H2-1,3-(CH2)2.

ftPBE results taken from ref (52).

RASPT2 (valence-π, πCPO active space) results taken from ref (43).

Table includes the ftPBE results from ref (52). The latter method performs better than other MC-PDFT[53] approaches for ST gaps of biradicals. Similarly to AC approximations, MC-PDFT is a post-CASSCF method relying on only 1- and 2-RDMs obtained from CAS. It employs density functional exchange-correlation functionals with modified arguments to describe electron correlation. As shown in Table , the accuracy of ST gap predictions by ftPBE does not match that of the AC1 method, with the percentage error nearly tripled and amounting to 38%. When comparing the computational efficiency of the adiabatic connection and MC-PDFT approximations, AC (or AC1) formally scale with the fifth power of the system size, which is one order more than scaling the MC-PDFT. (The timings of both methods are presented in the SI.) It should be noticed, however, that in the cases of both AC and MC-PDFT the major share of the total computational time is spent on the CASSCF calculation.

The accuracy achieved by AC1 comes close to that of the RASPT2 method. A comparison of RASPT2 (or CASPT2[43]) results with those of AC requires some care. These perturbation methods involve parameters to remove intruder states and to compensate for their tendency to underestimate gap energies between closed- and open-shell states.[54] The default value of the ionization potential-electron affinity shift[8] used in ref (43) improves ST gaps of biradicals predicted by CASPT2 and RASPT2 methods. In general, however, the shift may be problematic for strongly correlated systems, e.g., complexes with transition metals, and their tuning may be required.[55,56]

In summary, we have proposed a computational approach to the correlation energy in complete active space models. The novel AC formula for the correlation energy is based on a systematic expansion with respect to the adiabatic connection coupling constant α. Application to singlet–triplet gaps of single- and multireference systems revealed the need to account for higher-order terms in the α expansion. The AC/AC1 approaches showed a systematic improvement over the first-order AC0 method. The AC1 variant, which maximizes contributions from the higher-order terms, was identified as the best-performing AC approximation. Owing to the Cholesky decomposition technique, the AC methods achieve scaling of the computational time with the system size. Because they involve only 1- and 2-RDMs, they are well-suited to treat large active spaces. Importantly, the formalism used to derive AC is not limited to a particular form of the model Hamiltonian , thus further improvements in accuracy could be achieved with models other than that assumed in this work.

Compared to other correlation energy methods for strong correlation, AC emerges as having the most favorable accuracy to cost ratio. Advantages of AC over perturbation methods, such as CASPT2 or RASPT2, include not only the ability to treat dozens of active orbitals but also the lack of parameters and strict size consistency.[57] We believe that the presented development opens new perspectives for meeting the challenge of strong correlation, e.g., by DMRG[6] methods.