Communication: An adaptive configuration interaction approach for strongly correlated electrons with tunable accuracy

Most popular methods in electronic structure theory by some means attempt to exploit the sparsity of full configuration interaction (FCI) wave functions.¹ The exponential scaling of the number of determinants with respect to the number of orbitals required for FCI calculations prevents its use for all but trivially small systems, or for active space calculations no larger than 18 electrons in 18 orbitals. Recently, wave function factorization techniques such as the density matrix renormalization group (DMRG),^2–6 and stochastic CI approaches such as Monte Carlo CI (MCCI)^7–10 and FCI Quantum Monte Carlo (FCIQMC)^11–14 have risen as promising alternatives to FCI and complete active space CI (CASCI), allowing for the description of chemically interesting systems.^15,16

In this study, we propose a new adaptive configuration interaction (ACI) method that produces compact wave functions with tunable accuracy. The ACI is based on the framework of selected CI,^17–20 which recently has received renewed attention.^21–30 It uses two parameters to control the treatment of electron correlation. As will be shown, a remarkable property of the ACI is its ability to compute electronic energies with almost perfect control over the energy error. Additionally, we demonstrate that the ACI is a viable alternative to traditional complete active space (CAS) methods by performing ACI computations on active spaces that are outside the reach of CASCI.

Given a set of orthonormalized one-electron molecular orbitals, {ϕ_p}, the ACI requires the user to specify the number of electrons, the spin multiplicity, and two orbital subsets: doubly occupied orbitals and active orbitals. The latter are partially occupied in all determinants generated by the ACI. The ACI procedure is illustrated in Fig. 1 and consists of the following steps:

• At each iteration k, we define the space of reference determinants [P^(k)],

  $P ( k ) = { Φ μ : μ = 1 , 2 , … , d k } ,$

  (1)

  where d_k is the dimension of the P^(k) space. To this space, we associate the configuration interaction wave function $Ψ P ( k )$⁠, defined as

  $| Ψ P ( k ) 〉 = ∑ μ = 1 d k C μ | Φ μ 〉 ,$

  (2)

  where the coefficients C_μ are determined by diagonalizing the Hamiltonian in the space P^(k). In most cases, we begin the ACI process with an initial reference space, P⁽⁰⁾, that contains a single determinant, though a set of determinants can be used to speed convergence.
• From the reference space P^(k), all singly and doubly excited determinants are generated. For a given P^(k) space, we define the usual first-order interacting space (FOIS), F^(k), as the union of all unique singly [S^(k)] and doubly [D^(k)] excited determinants out of the reference space,

  $F ( k ) = ( S ( k ) ∪ D ( k ) ) ∖ P ( k ) .$

  (3)

  Denoting the occupied (virtual) orbitals of determinant Φ_μ ∈ P^(k) as i, j, … (a, b, …), then S^(k) and D^(k) may be written compactly as $S ( k ) = { a ˆ a † a ˆ i Φ μ : ∀ Φ μ ∈ P ( k ) }$ and $D ( k ) = { a ˆ a † a ˆ b † a ˆ j a ˆ i Φ μ : ∀ Φ μ ∈ P ( k ) }$⁠.
• To each determinant in F^(k), we associate an estimate of the energy contribution. Following degeneracy-corrected perturbation theory,³¹ we consider the 2 × 2 Hamiltonian for a determinant Φ_I ∈ F^(k) interacting with the P-space wave function $Ψ P ( k )$⁠,

  $H = 〈 Ψ P ( k ) | H ˆ | Ψ P ( k ) 〉 〈 Ψ P ( k ) | H ˆ | Φ I 〉 〈 Φ I | H ˆ | Ψ P ( k ) 〉 〈 Φ I | H ˆ | Φ I 〉 = E P V V ∗ E I .$

  (4)

  Diagonalization of H yields two real eigenvalues (λ₁, λ₂, λ₁ ≤ λ₂), the lower of which differs from the energy of $Ψ P ( k )$ (E_P) by

  $ϵ ( Φ I ) = λ 1 − E P = Δ 2 − Δ 2 4 + | V | 2 ,$

  (5)

  where Δ = E_I − E_P. Eq. (5) defines the energy importance criterion used in ACI to screen the first-order interacting space.
• Using the energy importance criterion, we define the secondary space Q^(k), the set of the most important determinants in F^(k). To build Q^(k), we use an aimed selection scheme.³² First, we sort the set F^(k) in decreasing order according to |ϵ(Φ_I)|, the absolute value of the energy importance criterion. Second, starting from the determinant with the lowest |ϵ(Φ_I)|, we exclude all those elements of F^(k) such that the cumulative energy error is less than a user-specified threshold σ expressed in units of mE_h,

  $∑ Φ I ∈ F ( k ) ∖ Q ( k ) | ϵ ( Φ I ) | ≤ σ .$

  (6)

  The determinants that are not discarded from F^(k) form the set Q^(k).
• With the Q^(k) space built, we can define the total model space at iteration k [M^(k)] as the union between the reference space and the secondary space,

  $M ( k ) = P ( k ) ∪ Q ( k ) ,$

  (7)

  and diagonalize the Hamiltonian in the space M^(k) to obtain the model space wave function,

  $| Ψ M ( k ) 〉 = ∑ Φ I ∈ M ( k ) C I | Φ I 〉 ,$

  (8)

  and the associated energy $E M ( k )$⁠. The model space energy may be corrected for the contributions of the determinants excluded from the secondary space [Φ_I ∈ F^(k)∖Q^(k)] using the second-order perturbative estimate,

  $E F ( k ) ≈ E M ( k ) + ∑ Φ I ∈ F ( k ) ∖ Q ( k ) ϵ ( Φ I ) .$

  (9)

  In the limit of V that goes to zero, the contribution of a single determinant to the second-order perturbative correction is approximately ϵ(Φ_I) ≈ − |V|²/Δ. Therefore, Eq. (9) may be viewed as a state-specific second-order correction based on the Epstein–Nesbet partitioning of the Hamiltonian.¹⁸ 
• Rather than directly augmenting the total model space as the iterations proceed, as is traditionally done in selected CI methods, we coarse grain the space M^(k) to form an updated reference space P^(k+1). Specifically, the M^(k)-space determinants are sorted according to the square of the CI coefficients (|C_I|²) in decreasing order. Determinants are progressively included in P^(k+1) until the sum of the squared coefficients is less than 1 − γσ, where γ is a constant that has units of (energy)⁻¹,

  $∑ Φ μ ∈ P ( k + 1 ) | C μ | 2 < 1 − γ σ .$

  (10)

• Steps i–vi are repeated until the energy of the M^(k) space is converged. This convergence of the energy coincides with the convergence of P^(k) and M^(k) with respect to the determinants included.

ACI improves upon previous selected CI methods like CIPSI¹⁸ and CI+PT²⁰ in a number of important ways. First, the aimed selection procedure gives the user a priori control over the absolute error in a computation. Additionally, the coarse-graining step (vi) increases the efficiency of the selection process (analogous to the initiator approximation of FCIQMC)¹² and decreases the dependence on the starting wave function guess.

For all ACI calculations, the parameters σ and γ are directly related to the desired energy accuracy. We found it convenient to assume a constant value of γ, and in this work all results were obtained using $γ=1m E h − 1$⁠. Accordingly, ACI results will be denoted as ACI(σ), while the ACI energy corrected for the determinant excluded from the secondary space [Eq. (9)] will instead be indicated as ACI(σ) + PT2.

To illustrate the ability of ACI to determine molecular energies with nearly constant accuracy along a potential energy surface, we examine the dissociation of N₂.³³ Figure 2 shows the error with respect to FCI for the ground-state potential energy curve of N₂ computed with ACI using canonical restricted Hartree–Fock (RHF) orbitals. In addition, we plot results for the internally contracted multireference CISD (MR-CISD),³⁴ and MR-CISD with Davidson’s correction (MR-CISD+Q)³⁵ based on a CAS self-consistent-field reference with six electrons in six orbitals [CASSCF(6,6)]. MR-CISD and MR-CISD+Q data from Ref. 36 were used. Figure 2(a) illustrates a distinguishing factor of the ACI: the absolute error at each point along the curve is reliably estimated by the energy threshold σ. Moreover, while the ACI(10) curve displays noticeable microscopic discontinuity, the ACI(5) and ACI(1) curves are progressively smoother. The inclusion of the second-order perturbative correction (see Fig. 2(b)) leads to curves that are approximately within 1 mE_h from the FCI energy. In comparison, MR-CISD shows fairly constant error throughout the dissociation, but with a noticeable increase near 1.6 Å. With the +Q correction, the error is fairly constant across the potential, though with a slight decrease in accuracy towards dissociation. Additionally, these energies are not variational.

Table I compares the energy error with respect to FCI (ΔE) and the size of the ACI determinant space for N₂ at the bond distances 1.1 and 3 Å. In both cases, ACI energy errors with respect to FCI show very good correlation with the value of σ. For a given value of σ, the energy difference |ΔE(r = 3) − ΔE(r = 1.1)| is only a fraction of the absolute error, showing the ability of the ACI method to describe both static and dynamic correlation in a balanced way. With the perturbative correction, the absolute energy errors are further reduced but the non-parallelism errors (NPEs) remain virtually unchanged. When we use natural orbitals from second-order Møller–Plesset perturbation theory or CASSCF (see Table I), the ACI(1) gives a more compact model space, with improved energy error with respect to RHF orbitals. This result suggests that the parameter σ effectively controls the ACI error regardless of the molecular orbital basis.

To demonstrate the ability of ACI in accurately exploiting the sparsity of CASCI spaces, we turn to the polyacene series.^37–43 The polyacene geometries from Ref. 37 were used in our calculations, and only the π bonding and antibonding pairs included in the STO-3G basis were correlated. This corresponds to a CAS(4n + 2,4n + 2) wave function, where n is the number of fused benzene rings. Such a problem is intractable with CASCI for n ≥ 4, but it is well suited for DMRG³⁷ and the two-electron reduced density matrix (v-2RDM) method.^43–45 To directly compare with previous results,^37,43 all ACI computations use canonical RHF orbitals.

Table II shows the vertical singlet-triplet splittings (ΔE_ST = E^S=1 − E^S=0) and the required number of determinants for various values of σ, in addition to comparison with DMRG³⁷ and v-2RDM⁴³ results. Note that to guarantee sub-kcal mol⁻¹ accuracy, a σ value less than 1.6 mE_h is in principle required. However, already for σ = 10 mE_h, the ACI error with respect to DMRG is consistently less than 1 kcal mol⁻¹ through hexacene. For σ = 50 mE_h, we see the maximum error at anthracene, and in general, the errors in the singlet-triplet splitting decrease with increasing n. A similar trend is seen with the v-2RDM data, where the maximum in error is at n = 5. The singlet-triplet splittings were also computed with ACI+PT2, and they on average deviated from the corresponding ACI values on the order of 0.01 kcal mol⁻¹ with a maximum of 0.4 kcal mol⁻¹ for the σ = 100 calculation of naphthalene. Our current pilot ACI code can be used to perform computations with up to about 5 × 10⁶ determinants, which currently limits the ACI(10) to hexacene and the ACI(100) to decacene. However, we anticipate that a production-level implementation of the ACI method that can take advantage of distributed memory architectures will be able to routinely target 10⁷–10⁸ determinants.⁴⁶ Furthermore, we anticipate that like in the case of the DMRG,^3,47 in the context of the ACI a localized molecular orbital basis will be crucial to significantly compress the number of variational parameters, and in turn, expand the applicability of this method to larger active spaces.

In summary, the major benefits of the ACI method are that: (i) electron correlation can be treated in a balanced way without a priori knowledge of a system’s electronic structure and (ii) that the energy error is precisely controlled by one user-specified parameter. The most promising application of ACI is as a replacement for prohibitively large CASCI reference wave functions used in perturbative and non-perturbative treatments of electron correlation. However, to overcome the deficiencies of CASCI, ACI can be coupled with orbital optimization schemes to produce CASSCF-like references for active spaces much larger than the CAS(18,18) limit. The straightforward computation of ACI reduced density matrices enables these extensions. Given that the ACI: (i) is a variational deterministic procedure, (ii) does not require orbital ordering, (iii) can be straightforwardly extended to excited states, and (iv) can be easily implemented on distributed memory architectures, it is an interesting alternative to DMRG, MCCI, FCIQMC, and v-2RDM methods.

All ACI results were obtained using our pilot code (Forte),⁴⁸ which is a suite of multireference methods written as a plugin to the open-source quantum chemistry package Psi4.⁴⁹ This work was supported by start-up funds provided by Emory University.