Excited states using semistochastic heat-bath configuration interaction

The accurate ab initio calculation of low-energy excited states is of great importance in many fields, including spectroscopy and solar energy. Unfortunately, excited-state calculations are complicated by the fact that they often exhibit strong multireference character, and a simple and accurate variational ansatz does not exist for them. As a result, the commonly used quantum chemical techniques such as density functional theory^1–3 and coupled cluster theory^4–7 become unreliable. Even when the excited states are qualitatively well described by a single determinant, single reference method, such equation of motion coupled cluster singles and doubles^8–10 (EOM-CCSD) are unable to describe states that have double-excitation character.

Complete active space self-consistent field^11–14 (CASSCF) and its extensions, including multireference configuration interaction^15,16 (MRCI), complete active space perturbation theory^17,18 (CASPT2), n-electron valence perturbation theory^19–21 (NEVPT2), and variants of multi-reference coupled cluster^22–24 (MRCC), are presently the most reliable methods for dealing with such problems. The shortcoming of these methods lies in their inability to treat problems that require more than 18 electrons and 18 orbitals in their active space because of the need for performing exact diagonalization, though recently it has become possible to treat systems with 20 electrons in 20 orbitals on massively parallel machines.²⁵ To overcome this limitation, other methods such as restricted^26,27 and generalized²⁸ active space (RASSCF/GASSCF) methods further subdivide the active space orbitals and put restrictions on their occupation pattern. However, these methods quickly become unaffordable as one relaxes these restrictions to calculate exact active-space energies. Methods such as the density matrix renormalization group^29–33 (DMRG) algorithm, Full Configuration Interaction Quantum Monte Carlo^34,35 (FCIQMC), and its semistochastic improvement,³⁶ when used as active-space solvers, are able to go well beyond the restriction imposed by CASSCF and represent a significant advance. They are nevertheless exponentially scaling, with the exception of DMRG for systems with a linear topology. Although other approaches such as variational Monte Carlo,^37–39 various flavors of projector Monte Carlo,^40–43 and reduced density matrix based methods^44–46 have shown promise, they have not yet become widely used in quantum chemistry.

In this paper, we present a new, efficient excited-state algorithm using our recently developed Heat-bath Configuration Interaction⁴⁷ (HCI) algorithm, in conjunction with semistochastic perturbation theory.⁴⁸ The HCI method is in the category of Selected Configuration Interaction followed by Perturbation Theory (SCI+PT) algorithms. The first SCI+PT algorithm, called Configuration Interaction by Perturbatively Selecting Iteratively (CIPSI), was developed over four decades ago by Malrieu and co-workers^49,50 and extended earlier selected CI algorithms that did not use perturbation theory.^51,52 Many variations of CIPSI have been developed over the years.^53–77 These methods perform a pruned breadth-first search to explore Slater determinant space aiming to identify those determinants that have the largest overlap with the targeted ground and excited states. This step is followed by Epstein-Nesbet perturbation theory^78,79 to include the energy contributions from the first-order interacting space.

HCI distinguishes itself from other SCI+PT methods in two respects. First, it changes the selection criterion, thereby allowing it to explore only the most important determinants without ever having to consider unimportant determinants (see Sec. II). Second, it performs the perturbation theory using a semistochastic algorithm that eliminates the severe memory bottleneck of having to store the large number of determinants in the first-order interacting space (although memory-efficient deterministic variants are also possible, we have so far found them to be computationally much more expensive than the semistochastic algorithm). These two ingredients make it more efficient than other SCI+PT approaches. It is worth mentioning that, similar to FCIQMC and DMRG, the cost of the method scales exponentially with the number of electrons. However, in our experience, HCI is much faster than FCIQMC and is also faster than DMRG, at least for molecules that do not have a large number of spin-orbitals with occupancies far from 0 and 1.

Here we show that HCI can be made more efficient by utilizing time-reversal symmetry and angular momentum symmetry, which is the largest Abelian subgroup of the full $D∞h$ point group. We also present a new method for extrapolation of HCI energies to the full configuration interaction (FCI) limit. In contrast to the first two publications^47,48 describing HCI, in which we either extrapolated with respect to the variational parameter $ϵ1$ or assumed that our calculations were converged, here we extrapolate the energies with respect to the perturbation energy ΔE₂. We find that the extrapolation with respect to ΔE₂ is often nearly linear and more reliable than the previous extrapolation procedure.

The outline of the paper is as follows. In Secs. II and III, we set the stage by describing the salient features of the HCI algorithm and semistochastic perturbation theory. Next, in Sec. IV, we show how the HCI algorithm can be extended to calculate excited states. In Sec. V, we present an improved method for extrapolating unconverged HCI energies to the FCI energies. In Sec. VI, we describe further improvements to the algorithm, including the incorporation of angular momentum symmetries and time-reversal symmetry. We finish the paper with some concluding remarks in Sec. VIII.

HCI is an efficient SCI+PT algorithm, which can be broken down into the following steps:

1. Variational stage

   • Identify the most important Slater determinants.

   • Find a variational wavefunction and energy by computing the ground state within the space spanned by determinants found in step 1a.

2. Perturbative stage

   • (a)

     Identify the most important perturbative corrections to the variational energy.

   • (b)

     Sum the contributions found in step 2a.

Step 1 is performed as an iterative process, which alternates between adding new determinants to the selected space and finding the lowest-energy wavefunction within the current selected space. HCI improves over other SCI+PT algorithms by using a very fast algorithm for selecting the important determinants in steps 1a and 2a and also by using a semistochastic algorithm for performing the summation in step 2b, which eliminates the need for storing all the determinants that contribute to the perturbative correction.

In HCI, the variational wavefunction at any iteration is given by $ψ=∑iVciDi$⁠, and the new determinants $Da$ that are added to the variational space are those for which

Here H_ai are matrix elements connecting states $Di$ within the current variational space $V$ to states $Da$ outside $V$⁠, and $ϵ1$ is a user-defined parameter that controls the size of the final variational space.⁸⁰ 

This criterion is used because it can be implemented efficiently without checking the vast majority of the determinants that do not meet the criterion (see Sec. II D). The determinants chosen using this scheme are approximately those chosen by the CIPSI algorithm, which chooses determinants $Da$ for which

is sufficiently large.

The determinants chosen by the two criteria are nearly the same because the terms in the numerator of Eq. (2) span many orders of magnitude, so the sum is highly correlated with the largest-magnitude term in the sum [Eq. (1)]. The denominators of Eq. (2) also vary with D_a, but to a much lesser extent, since the determinant energies, E_a, are much larger than the current variational energy, E₀, for sufficiently large variational expansions.

In SCI+PT algorithms, the perturbative correction ΔE₂ is typically computed using Epstein-Nesbet perturbation theory,

In the original CIPSI algorithm, this expression is computed by evaluating and summing all of the terms in the double sum. However, the vast majority of the terms in the sum are negligible, so this approach is not very efficient. Various schemes for improving the efficiency have been implemented, including only exciting from a rediagonalized list of the largest-weight determinants⁵⁰ and its efficient approximation using diagrammatic perturbation theory.⁵⁷ However, this is both more complicated than necessary (requiring a double extrapolation with respect to the two variational spaces to reach the full CI limit) and is more computationally expensive than necessary since even the largest weight determinants have many connections that make small contributions to the energy.

Instead, HCI uses a “screened sum,”

where $∑i(ϵ2)Haici$ includes only terms for which $Haici>ϵ2$⁠. Note that the vast number of terms that do not meet this criterion are never evaluated. Even with this screening, the number of connected determinants can be sufficiently large to exceed computer memory. This is addressed in Sec. III.

The key to the efficiency of the heat-bath scheme is as follows. The vast majority of the Hamiltonian matrix elements correspond to double excitations, and their values do not depend on the determinants themselves but only on the four orbitals whose occupancies change during the double excitation. Therefore, before performing an HCI run, for each pair of orbitals, the set of all double-excitation matrix elements obtained by exciting from that pair of orbitals is computed and stored in decreasing order by magnitude, along with the corresponding pairs of orbitals the electrons would excite to. Once this is done, at any point in the HCI algorithm, from a given reference determinant, all double excitations whose Hamiltonian matrix elements exceed a cutoff (either $ϵ1/|ci|$ or $ϵ2/ci$ for the variational and perturbative stages, respectively) can be generated efficiently, without having to loop over all double excitations. This is achieved by looping over all pairs of occupied orbitals in the reference determinant and traversing the list of sorted double-excitation matrix elements for each pair until the cutoff is reached.

This screening algorithm is utilized in both steps 1a and 2a of the algorithm and is a significant reason why the HCI algorithm is faster than other selected CI algorithms that do not truncate the search for double excitations or skip over the large number of determinants making negligible contributions to the energy.

The evaluation of Eq. (4) with a low computational cost requires the simultaneous storage of all included terms, indexed by a, and can easily exceed memory limitations for challenging problems.

To overcome this memory bottleneck, we introduced a semistochastic evaluation of the PT sum, in which the most important contributions (found in the same way as in the original HCI algorithm) are evaluated deterministically and the rest are sampled stochastically.⁴⁸ Here, an initial deterministic perturbative correction $ΔE2D[ϵ2d]$ is calculated using a relatively loose threshold $ϵ2d$⁠. Then, the stochastic calculation is used to evaluate the bias in the deterministic calculation (due to using an insufficiently small $ϵ2d$⁠) by calculating the two stochastic energies $ΔE2[ϵ2]$ and $ΔE2[ϵ2d]$ (the second-order perturbative energy calculated with $ϵ2$ and $ϵ2d$⁠, respectively) with $ϵ2<ϵ2d$⁠. The total second-order energy is given by the expression

Both $ΔE2[ϵ2]$ and $ΔE2[ϵ2d]$ are calculated using the same set of samples, and thus there is significant cancellation of stochastic error. Furthermore, because these two energies are calculated simultaneously, the incremental cost of performing this calculation, compared with a fully stochastic summation, is very small. Clearly, in the limit in which $ϵ2d=ϵ2$⁠, the entire perturbative calculation becomes deterministic.

The stochastic piece of the semistochastic summation algorithm is evaluated by sampling only N_d variational determinants at a time. Each variational determinant D_i is sampled, with replacement, with probability

using the Alias method,^81,82 which has previously been used by us to efficiently sample double excitations in determinant-space quantum Monte Carlo algorithms.⁸³ An unbiased estimate of the second-order perturbative correction to the energy⁴⁸ is given by the expression

where $w$_i denotes the number of times determinant D_i is sampled, and the summation is only over the $Nduniq$ unique determinants in the sampled N_d determinants. A minimum of just two determinants is sufficient to perform this calculation, thus completely eliminating the memory bottleneck; however, the statistical noise is greatly diminished if larger samples are chosen.

For large systems, the stochastic part of the algorithm is essential, but for small systems, it is possible to get within 1 mHa of the FCI energy using just the deterministic part of the algorithm. This is demonstrated for the C₂ dimer in the cc-pV5Z basis set in Table I. The variational plus the deterministic parts of the algorithm yield energies for the $Σg+1$ ground and excited states accurate to 0.3 mHa in just 20 min on a single computer node.

We note that very recently a different semistochastic perturbation theory has been proposed⁷⁵ wherein the statistical error decreases much faster than the inverse square root of the number of Monte Carlo samples.

When computing ground and excited states, all states are expanded in the same set of variational determinants,

where s denotes the state. This is akin to performing state-average variational calculations, rather than state-specific ones where each state would have its own set of determinants. In the excited-state algorithm, at each iteration of the variational stage, we add to the variational space $V$ the union of the new determinants that are important for each of the states. Thus, at each iteration, new determinants D_j are added to $V$ if

for at least one $Di∈V$⁠. Eq. (9) ensures that when several states are targeted, the variational space will include more determinants than when only the ground state is targeted since there will be determinants relevant to all targeted states. Note that this formula is different from the one used in state-average CASSCF theory where the density matrix is averaged, which is closer in spirit to taking the square root of the sum of the squares of the coefficients. This distinction becomes important in the event of degeneracies among the targeted variational states. In such a situation, rotations within the degenerate subspaces are arbitrary, and the value of the maximum magnitude of the coefficients is not invariant to such rotations, but the square root of the sum of the squares of coefficients is invariant. For such a situation, we recommend using the invariant criterion, for example, if the ground state is nondegenerate but the first two excited states are degenerate, one should use $maxci(0),ci(1)2+ci(2)2$ in place of $maxci(0),ci(1),ci(2)$⁠. In the applications considered in this paper, there are no exact degeneracies among the targeted states, so we use the simple formula in Eq. (9).

Apart from the generalization to excited states, the most important modification to HCI in this paper is a new procedure for extrapolation of the HCI total energy to the FCI limit. The HCI energy is a function of two parameters: $ϵ1$⁠, which controls the variational stage, and $ϵ2$⁠, which accelerates the perturbative energy calculation by screening out the many tiny contributions. In the limit in which $ϵ1$ goes to zero, the HCI energy equals the FCI energy, and in the limit in which $ϵ2$ goes to zero, the perturbative correction is exactly equal to the Epstein-Nesbet perturbation correction. In the calculations in this paper, we use a fixed $ϵ2$ = 10⁻⁸ Ha, which is sufficiently small to give near exact PT energies, and perform runs at several different values of $ϵ1$⁠.

In the original HCI paper,⁴⁷ we extrapolated to the FCI limit by extrapolating the HCI energy with respect to $ϵ1$⁠. However, this is often nonlinear with a curvature that increases as $ϵ1$ = 0 is approached. Consequently, it is difficult to choose a function that provides a good fit to the computed energies. Instead, in this paper, we extrapolate with respect to the perturbative correction to the energy. In the limit in which this perturbative correction is zero, both the variational and the total HCI energies equal the FCI energy. If the fitting function contains a linear term, the variational and the total HCI energies extrapolate to precisely the same value. For many systems, this extrapolation is empirically found to be close to linear, as shown in Fig. 1, but more generally the points lie on a sufficiently smooth curve for which a reliable low-order polynomial extrapolation is possible.

Since our most recent HCI paper,⁴⁸ in which we introduced the semistochastic algorithm for evaluating the HCI perturbative correction to the energy, we have improved the algorithm in several ways. First, we have introduced the ability to employ angular momentum symmetry, which is the largest Abelian subgroup of the $D∞h$ point group (for all but the smallest basis sets, it has orbitals, as well as many-body states, of a larger number of irreducible representations than does the D_2h point group that is commonly used for linear molecules). Second, we have implemented time-reversal symmetry, which can be used to perform separate calculations of the singlet and triplet states, thereby reducing the Hilbert space of the problem by nearly a factor of two and reducing the memory requirement of the Hamiltonian in the variational space by a factor of between two and four.

For real orbitals, the 2-electron integrals have 8-fold permutational symmetry. Hence only slightly more than an eighth of the integrals need to be stored. For linear molecules, the orbitals can be chosen to be eigenstates of the z-component of angular momentum, $L^z$⁠, and the orbitals are complex. In that case, there is only 4-fold permutational symmetry. However, with the usual choice of phase, the integrals are real, and four of the eight are zero since they violate L_z conservation. Hence it is still possible to store only an eighth of the integrals, provided a check is performed to ensure L_z conservation. This enables us to use L_z symmetry to reduce the storage required for the Hamiltonian without increasing the storage required for the integrals.

The time-reversal operator exchanges the spin labels of the electrons. States with S_z = 0 are symmetric/antisymmetric under time reversal if S is even/odd. Consequently, the basis states can be chosen to be symmetric or antisymmetric linear combinations of time-reversed pairs of Slater determinants.

Consider two spatial orbital configurations, I and J. If a determinant is formed by assigning the α electrons to I and the β electrons to J, i.e., $IαJβ$⁠, then its time-reversed partner is $(−1)nα(nβ+1)JαIβ$⁠, where n_α and n_β are the number of alpha and beta electrons. Note $(−1)nα(nβ+1)$ is always equal to 1 for a system containing an equal number of alpha and beta electrons, so we will ignore this phase from now on. We choose to work in the basis of states $SIJ$⁠, where

where z is the eigenvalue of the time-reversal operator, which is either 1 for even S states or −1 for odd S states. Note that basis states for which I = J can occur only when z = 1.

The matrix elements between pairs of these time-reversal symmetrized states are straightforwardly evaluated. For example, if I = J and $K≠L$⁠,

whereas if $I≠J$ and $K≠L$⁠,

We use time-reversal symmetry only for the variational stage. Upon completion of the variational stage, we convert back to the determinant basis and perform Epstein-Nesbet perturbation theory in this basis.

Using time-reversal symmetrized states has two benefits. First, it shrinks the size of the Hilbert space, so that a larger variational manifold can be treated with a given memory size. Second, it allows one to target different irreducible representations separately. For example, if the ground state is a singlet and the first excited state is a triplet, then one can target the lowest triplet state as a ground state and avoid using the excited-state algorithm.

We consider the excited states of the strongly correlated carbon dimer. Despite its small size, the carbon dimer has strong multireference character even in its ground state and has been the focus of many experimental and theoretical studies.^41,85–107 Here we perform excited-state calculations in Dunning’s cc-pVQZ basis¹⁰⁸ to compare to calculations from other methods in the literature. Then, we compute fourteen low-lying potential energy surfaces in the larger cc-pV5Z basis, extrapolating to the full CI limit.

All integrals used in these calculations were obtained using the PySCF quantum chemistry package.¹⁰⁹ 

In order to compare to DMRG and FCIQMC energies in the literature, we first computed the potential energy surfaces of the three lowest $Σg+1$ and two lowest $Σg5$ states in the cc-pVQZ basis with a frozen core. These states were targeted by imposing L_z = 0 (Σ states) and using a basis of linear combinations of Slater determinants, which is symmetric under time-reversal symmetry (singlets and quintets).

The HCI variational and total energies are shown in Table II. Note that even in a relatively large cc-pVQZ basis, the variational energies are within 0.5 mHa of the converged total energies for both ground and excited states. The HCI total energies are lower than the (bond dimension 4000) DMRG and FCIQMC energies by 40-260 μHa and 120-710 μHa, respectively. For the ground state, DMRG energies that are in better agreement with HCI energies were obtained¹⁰⁶ by targeting just the ground state energy, e.g., the equilibrium energy is −75.802 69(1) Ha, consistent with the extrapolated HCI total energy. The discrepancy between the FCIQMC energies and the HCI total energies is likely due to the initiator bias in FCIQMC.

We next computed the potential energy surfaces of fourteen low-lying states of the carbon dimer in the cc-pV5Z basis,

Besides spatial symmetry, time-reversal symmetry was used to further reduce the size of the Hilbert space by targeting even S and odd S states separately. Thus, for example, the three Σ_g states were computed in two runs: one that targeted the two lowest-energy singlets and one that targeted only the lowest energy triplet.

To accelerate convergence of the HCI total energy with respect to $ϵ1$⁠, HCI natural orbitals were used. Within each of the symmetry sectors (a symmetry sector is defined by spatial and time-reversal symmetry), the natural orbitals corresponding to the state-averaged 1-RDM of the lowest variational states of interest were computed. Thus, at each geometry, ten sets of natural orbitals were computed, one for each of the ten symmetry sectors. After the natural orbitals were obtained, for each of the ten symmetry sectors, at each geometry, at least two HCI runs were performed with those natural orbitals in order to enable extrapolation to the FCI limit.

Each HCI calculation was performed starting from a single basis state of the target irreducible representation, found automatically using the following algorithm. Start from a low-energy determinant found by filling the orbitals whose diagonal 1-body integrals are smallest. This is the current best guess for a good HCI starting state, but it need not have the desired symmetry. Next, repeat the following step until convergence: Replace the current HCI starting determinant with the determinant of the target irreducible representation with the lowest energy out of the set of determinants that are no more than a double excitation away from the current HCI starting determinant.

Although this algorithm does not necessarily result in either the lowest-energy basis state of the target symmetry sector or the one with maximum overlap with the full CI ground state, it resulted in a good enough starting point for the HCI runs in this paper.

It should be mentioned that truncated CI methods such as HCI can produce states with spin contamination. Time-reversal symmetry ameliorates but does not completely eliminate this problem. We computed the total spin of each variational wavefunction (with $ϵ1=5×10−5$ Ha) at several points along the curves and found the maximum deviations from the correct total spin to be $S2≲5×10−4$ for the singlets, $S2−2≲1×10−4$ for the triplets, and $S2−6≲2×10−4$ for the quintet.

The energies of these fourteen states are shown in Table III and Fig. 2. In addition, we have calculated the excitation energies of the eight lowest-lying excited states, as shown in Table IV. These excitation energies have a mean absolute deviation of 0.02 eV relative to the experimental values.

We have presented an efficient excited-state method using Heat-bath Configuration Interaction and a method for extrapolating the resulting energies to the FCI limit. We incorporated symmetries including time-reversal symmetry and angular momentum conservation, enabling us to target excited states in different symmetry sectors. We then used the method to calculate fourteen low-lying potential energy surfaces of the carbon dimer in the large cc-pV5Z basis.

We are exploring including one more symmetry: the analog of time-reversal symmetry for orbital angular momentum (to separately target Σ⁺ and Σ⁻ states). For challenging problems, for which extrapolation to the full CI limit is not possible, we are also extending our extrapolation procedure to the case where the variational and perturbative steps have different active space sizes, resulting in an extrapolation to the limit of an uncontracted multireference perturbation theory with a complete active space (CAS) reference, rather than to the full CI limit.

The calculations in this paper were performed using the University of Colorado’s Research Computing cluster. S.S. and A.A.H. were supported by the startup package from the University of Colorado. The research was also supported in part by NSF Grant No. ACI-1534965. C.J.U. thanks Garnet Chan for valuable discussions.