We extend our recently developed heat-bath configuration interaction (HCI) algorithm, and our semistochastic algorithm for performing multireference perturbation theory, to calculate excited-state wavefunctions and energies. We employ time-reversal symmetry, which reduces the memory requirements by more than a factor of two. An extrapolation technique is introduced to reliably extrapolate HCI energies to the full CI limit. The resulting algorithm is used to compute fourteen low-lying potential energy surfaces of the carbon dimer using the cc-pV5Z basis set, with an estimated error in energy of 30-50 *μ*Ha compared to full CI. The excitation energies obtained using our algorithm have a mean absolute deviation of 0.02 eV compared to experimental values.

## I. INTRODUCTION

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.^{25} To overcome this limitation, other methods such as restricted^{26,27} and generalized^{28} 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,^{36} 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^{47} (HCI) algorithm, in conjunction with semistochastic perturbation theory.^{48} 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\u221eh$ 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 $\u03f51$ or assumed that our calculations were converged, here we extrapolate the energies with respect to the perturbation energy Δ*E*_{2}. We find that the extrapolation with respect to Δ*E*_{2} 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.

## II. HEAT-BATH CONFIGURATION INTERACTION (HCI)

### A. Overview

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

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.

Perturbative stage

- (a)
Identify the most important perturbative corrections to the variational energy.

- (b)
Sum the contributions found in step 2a.

- (a)

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.

### B. Variational stage

In HCI, the variational wavefunction at any iteration is given by $\psi \u2009=\u2009\u2211iVciDi$, 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 $\u03f51$ is a user-defined parameter that controls the size of the final variational space.^{80}

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*_{0}, for sufficiently large variational expansions.

### C. Perturbative stage

In SCI+PT algorithms, the perturbative correction Δ*E*_{2} 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^{50} and its efficient approximation using diagrammatic perturbation theory.^{57} 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 $\u2211i(\u03f52)Haici$ includes only terms for which $Haici>\u03f52$. 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.

### D. Heat-bath algorithm for acceleration of both stages

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 $\u03f51/|ci|$ or $\u03f52/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.

## III. SEMISTOCHASTIC PERTURBATION THEORY

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.^{48} Here, an initial deterministic perturbative correction $\Delta E2D[\u03f52d]$ is calculated using a relatively loose threshold $\u03f52d$. Then, the stochastic calculation is used to evaluate the bias in the deterministic calculation (due to using an insufficiently small $\u03f52d$) by calculating the two stochastic energies $\Delta E2[\u03f52]$ and $\Delta E2[\u03f52d]$ (the second-order perturbative energy calculated with $\u03f52$ and $\u03f52d$, respectively) with $\u03f52\u2009<\u2009\u03f52d$. The total second-order energy is given by the expression

Both $\Delta E2[\u03f52]$ and $\Delta E2[\u03f52d]$ 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 $\u03f52d\u2009=\u2009\u03f52$, 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.^{83} An unbiased estimate of the second-order perturbative correction to the energy^{48} 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_{2} dimer in the cc-pV5Z basis set in Table I. The variational plus the deterministic parts of the algorithm yield energies for the $\Sigma g+1$ ground and excited states accurate to 0.3 mHa in just 20 min on a single computer node.

Component . | E_{0} (Ha)
. | E_{1} (Ha)
. | Time (min) . |
---|---|---|---|

Variational energy | −75.805 98 | −75.715 73 | 9 |

Deterministic component of PT correction | −0.002 20 | −0.002 32 | 11 |

Stochastic component of PT correction | 0.000 03(1) | 0.000 03(1) | 18 |

Total HCI energy | −75.808 15(1) | −75.718 02(1) | 38 |

Extrapolated total HCI energy | −75.807 90(3) | −75.717 73(3) |

Component . | E_{0} (Ha)
. | E_{1} (Ha)
. | Time (min) . |
---|---|---|---|

Variational energy | −75.805 98 | −75.715 73 | 9 |

Deterministic component of PT correction | −0.002 20 | −0.002 32 | 11 |

Stochastic component of PT correction | 0.000 03(1) | 0.000 03(1) | 18 |

Total HCI energy | −75.808 15(1) | −75.718 02(1) | 38 |

Extrapolated total HCI energy | −75.807 90(3) | −75.717 73(3) |

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

## IV. EXCITED STATES IN HCI

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\u2009\u2208\u2009V$. 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).

## V. EXTRAPOLATION OF HCI ENERGIES

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: $\u03f51$, which controls the variational stage, and $\u03f52$, which accelerates the perturbative energy calculation by screening out the many tiny contributions. In the limit in which $\u03f51$ goes to zero, the HCI energy equals the FCI energy, and in the limit in which $\u03f52$ 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 $\u03f52$ = 10^{−8} Ha, which is sufficiently small to give near exact PT energies, and perform runs at several different values of $\u03f51$.

In the original HCI paper,^{47} we extrapolated to the FCI limit by extrapolating the HCI energy with respect to $\u03f51$. However, this is often nonlinear with a curvature that increases as $\u03f51$ = 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.

## VI. FURTHER IMPROVEMENTS TO HCI

Since our most recent HCI paper,^{48} 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\u221eh$ 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.

### A. Angular momentum symmetry

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.

### B. Time-reversal symmetry

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\alpha J\beta $, then its time-reversed partner is $(\u22121)n\alpha (n\beta +1)J\alpha I\beta $, where *n*_{α} and *n*_{β} are the number of alpha and beta electrons. Note $(\u22121)n\alpha (n\beta +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\u2260L$,

whereas if $I\u2260J$ and $K\u2260L$,

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.

## VII. RESULTS

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^{108} 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.^{109}

### A. Carbon dimer in cc-pVQZ basis

In order to compare to DMRG and FCIQMC energies in the literature, we first computed the potential energy surfaces of the three lowest $\Sigma g+1$ and two lowest $\Sigma 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^{106} 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.

R/Å
. | DMRG variational energy (Ref. 106) . | FCIQMC energy (Ref. 103) . | HCI variational energy (this work) . | HCI total energy (this work) . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

1.0 | – | – | – | −0.655 70 | −0.486 65 | −0.376 54 | −0.655 98 | −0.486 88 | −0.376 92 | −0.656 20 | −0.487 16 | −0.377 25 |

1.1 | −0.761 24 | −0.621 83 | −0.502 28 | −0.761 14 | −0.621 70 | −0.502 12 | −0.761 03 | −0.621 57 | −0.501 96 | −0.761 28 | −0.621 86 | −0.502 33 |

1.2 | −0.799 20 | −0.694 59 | −0.544 90 | −0.799 13 | −0.694 50 | −0.544 79 | −0.799 01 | −0.694 35 | −0.544 61 | −0.799 27 | −0.694 65 | −0.544 98 |

1.242 53 | −0.802 64 | −0.712 08 | −0.549 53 | −0.802 58 | −0.712 00 | −0.549 42 | −0.802 44 | −0.711 82 | −0.549 24 | −0.802 71 | −0.712 13 | −0.549 61 |

1.3 | −0.799 33 | −0.726 33 | −0.548 71 | −0.799 27 | −0.726 26 | −0.548 61 | −0.799 13 | −0.726 07 | −0.548 42 | −0.799 39 | −0.726 39 | −0.548 81 |

1.4 | −0.779 65 | −0.732 67 | −0.537 76 | −0.779 61 | −0.732 61 | −0.537 66 | −0.779 45 | −0.732 40 | −0.537 46 | −0.779 73 | −0.732 74 | −0.537 89 |

1.6 | −0.724 01 | −0.704 87 | −0.510 54 | −0.723 95 | −0.704 80 | −0.510 47 | −0.723 74 | −0.704 57 | −0.510 24 | −0.724 10 | −0.704 95 | −0.510 72 |

1.8 | − | − | − | −0.680 56 | −0.654 07 | −0.496 39 | −0.680 29 | −0.653 89 | −0.496 12 | −0.680 71 | −0.654 24 | −0.496 61 |

2.0 | −0.645 52 | −0.614 69 | −0.492 90 | −0.645 48 | −0.614 70 | −0.492 97 | −0.645 22 | −0.614 53 | −0.492 69 | −0.645 65 | −0.614 86 | −0.493 16 |

R/Å
. | DMRG variational energy (Ref. 106) . | FCIQMC energy (Ref. 103) . | HCI variational energy (this work) . | HCI total energy (this work) . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

1.0 | – | – | – | −0.655 70 | −0.486 65 | −0.376 54 | −0.655 98 | −0.486 88 | −0.376 92 | −0.656 20 | −0.487 16 | −0.377 25 |

1.1 | −0.761 24 | −0.621 83 | −0.502 28 | −0.761 14 | −0.621 70 | −0.502 12 | −0.761 03 | −0.621 57 | −0.501 96 | −0.761 28 | −0.621 86 | −0.502 33 |

1.2 | −0.799 20 | −0.694 59 | −0.544 90 | −0.799 13 | −0.694 50 | −0.544 79 | −0.799 01 | −0.694 35 | −0.544 61 | −0.799 27 | −0.694 65 | −0.544 98 |

1.242 53 | −0.802 64 | −0.712 08 | −0.549 53 | −0.802 58 | −0.712 00 | −0.549 42 | −0.802 44 | −0.711 82 | −0.549 24 | −0.802 71 | −0.712 13 | −0.549 61 |

1.3 | −0.799 33 | −0.726 33 | −0.548 71 | −0.799 27 | −0.726 26 | −0.548 61 | −0.799 13 | −0.726 07 | −0.548 42 | −0.799 39 | −0.726 39 | −0.548 81 |

1.4 | −0.779 65 | −0.732 67 | −0.537 76 | −0.779 61 | −0.732 61 | −0.537 66 | −0.779 45 | −0.732 40 | −0.537 46 | −0.779 73 | −0.732 74 | −0.537 89 |

1.6 | −0.724 01 | −0.704 87 | −0.510 54 | −0.723 95 | −0.704 80 | −0.510 47 | −0.723 74 | −0.704 57 | −0.510 24 | −0.724 10 | −0.704 95 | −0.510 72 |

1.8 | − | − | − | −0.680 56 | −0.654 07 | −0.496 39 | −0.680 29 | −0.653 89 | −0.496 12 | −0.680 71 | −0.654 24 | −0.496 61 |

2.0 | −0.645 52 | −0.614 69 | −0.492 90 | −0.645 48 | −0.614 70 | −0.492 97 | −0.645 22 | −0.614 53 | −0.492 69 | −0.645 65 | −0.614 86 | −0.493 16 |

### B. Carbon dimer in cc-pV5Z basis

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 $\u03f51$, 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 $\u03f51\u2009=\u20095\xd710\u22125$ Ha) at several points along the curves and found the maximum deviations from the correct total spin to be $S2\u2009\u2272\u20095\xd710\u22124$ for the singlets, $S2\u22122\u22721\xd710\u22124$ for the triplets, and $S2\u22126\u22722\xd710\u22124$ 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.

R/Å
. | $X\Sigma g+1$ . | a^{3}Π_{u}
. | $b3\Sigma g\u2212$ . | A^{1}Π_{u}
. | $c3\Sigma u+$ . | B^{1}Δ_{g}
. | $B\u20321\Sigma g+$ . | d^{3}Π_{g}
. | C^{1}Π_{g}/^{5}Π_{g}
. | $1\Sigma u+/\u22121$ . | 1^{3}Δ_{u}
. | $23\Sigma u+/\u2212$ . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

1.0 | −0.662 52 | −0.589 52 | −0.504 14 | −0.550 13 | −0.645 46 | −0.465 81 | −0.494 05 | −0.542 66 | −0.485 60 | −0.471 69 | −0.286 27 | −0.313 76 |

1.1 | −0.767 01 | −0.725 04 | −0.662 60 | −0.686 97 | −0.738 83 | −0.628 00 | −0.627 99 | −0.661 69 | −0.601 29 | −0.571 69 | −0.447 93 | −0.475 45 |

1.2 | −0.804 61 | −0.786 57 | −0.741 65 | −0.749 86 | −0.765 05 | −0.710 40 | −0.700 37 | −0.708 62 | −0.644 96 | −0.606 67 | −0.535 47 | −0.562 09 |

1.3 | −0.804 44 | −0.804 88 | −0.773 73 | −0.769 51 | −0.754 08 | −0.745 53 | −0.731 84 | −0.714 85 | −0.648 48 | −0.605 60 | −0.578 48 | −0.602 77 |

1.4 | −0.784 60 | −0.798 79 | −0.778 64 | −0.764 81 | −0.724 83 | −0.753 22 | −0.737 99 | −0.700 09 | −0.656 80 | −0.586 06 | −0.596 44 | −0.615 43 |

1.5 | −0.756 63 | −0.779 88 | −0.768 45 | −0.747 36 | −0.690 09 | −0.745 66 | −0.728 93 | −0.677 59 | −0.669 95 | −0.558 51 | −0.608 48 | −0.612 40 |

1.6 | −0.728 95 | −0.755 24 | −0.750 62 | −0.724 16 | −0.658 98 | −0.730 27 | −0.709 75 | −0.657 10 | −0.671 53 | −0.623 66 | −0.619 00 | −0.614 93 |

1.7 | −0.705 82 | −0.728 97 | −0.729 53 | −0.699 51 | −0.636 23 | −0.711 47 | −0.684 17 | −0.642 16 | −0.666 52 | −0.627 77 | −0.622 95 | −0.617 93 |

1.8 | −0.685 50 | −0.703 34 | −0.707 80 | −0.675 75 | −0.623 33 | −0.691 88 | −0.658 60 | −0.630 00 | −0.657 97 | −0.626 92 | −0.622 14 | −0.613 52 |

1.95 | −0.658 37 | −0.668 74 | −0.676 77 | −0.644 56 | −0.614 98 | −0.664 07 | −0.627 20 | −0.613 73 | −0.642 70 | −0.620 83 | −0.616 09 | −0.599 05 |

2.1 | −0.635 61 | −0.640 12 | −0.649 63 | −0.620 25 | −0.606 65 | −0.640 09 | −0.605 72 | −0.600 43 | −0.627 60 | −0.612 57 | −0.607 67 | −0.58904 |

2.4 | −0.604 33 | −0.601 82 | −0.609 40 | −0.592 36 | −0.591 48 | −0.606 43 | −0.585 35 | −0.585 25 | −0.604 04 | −0.597 23 | −0.592 05 | −0.579 99 |

R/Å
. | $X\Sigma g+1$ . | a^{3}Π_{u}
. | $b3\Sigma g\u2212$ . | A^{1}Π_{u}
. | $c3\Sigma u+$ . | B^{1}Δ_{g}
. | $B\u20321\Sigma g+$ . | d^{3}Π_{g}
. | C^{1}Π_{g}/^{5}Π_{g}
. | $1\Sigma u+/\u22121$ . | 1^{3}Δ_{u}
. | $23\Sigma u+/\u2212$ . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

1.0 | −0.662 52 | −0.589 52 | −0.504 14 | −0.550 13 | −0.645 46 | −0.465 81 | −0.494 05 | −0.542 66 | −0.485 60 | −0.471 69 | −0.286 27 | −0.313 76 |

1.1 | −0.767 01 | −0.725 04 | −0.662 60 | −0.686 97 | −0.738 83 | −0.628 00 | −0.627 99 | −0.661 69 | −0.601 29 | −0.571 69 | −0.447 93 | −0.475 45 |

1.2 | −0.804 61 | −0.786 57 | −0.741 65 | −0.749 86 | −0.765 05 | −0.710 40 | −0.700 37 | −0.708 62 | −0.644 96 | −0.606 67 | −0.535 47 | −0.562 09 |

1.3 | −0.804 44 | −0.804 88 | −0.773 73 | −0.769 51 | −0.754 08 | −0.745 53 | −0.731 84 | −0.714 85 | −0.648 48 | −0.605 60 | −0.578 48 | −0.602 77 |

1.4 | −0.784 60 | −0.798 79 | −0.778 64 | −0.764 81 | −0.724 83 | −0.753 22 | −0.737 99 | −0.700 09 | −0.656 80 | −0.586 06 | −0.596 44 | −0.615 43 |

1.5 | −0.756 63 | −0.779 88 | −0.768 45 | −0.747 36 | −0.690 09 | −0.745 66 | −0.728 93 | −0.677 59 | −0.669 95 | −0.558 51 | −0.608 48 | −0.612 40 |

1.6 | −0.728 95 | −0.755 24 | −0.750 62 | −0.724 16 | −0.658 98 | −0.730 27 | −0.709 75 | −0.657 10 | −0.671 53 | −0.623 66 | −0.619 00 | −0.614 93 |

1.7 | −0.705 82 | −0.728 97 | −0.729 53 | −0.699 51 | −0.636 23 | −0.711 47 | −0.684 17 | −0.642 16 | −0.666 52 | −0.627 77 | −0.622 95 | −0.617 93 |

1.8 | −0.685 50 | −0.703 34 | −0.707 80 | −0.675 75 | −0.623 33 | −0.691 88 | −0.658 60 | −0.630 00 | −0.657 97 | −0.626 92 | −0.622 14 | −0.613 52 |

1.95 | −0.658 37 | −0.668 74 | −0.676 77 | −0.644 56 | −0.614 98 | −0.664 07 | −0.627 20 | −0.613 73 | −0.642 70 | −0.620 83 | −0.616 09 | −0.599 05 |

2.1 | −0.635 61 | −0.640 12 | −0.649 63 | −0.620 25 | −0.606 65 | −0.640 09 | −0.605 72 | −0.600 43 | −0.627 60 | −0.612 57 | −0.607 67 | −0.58904 |

2.4 | −0.604 33 | −0.601 82 | −0.609 40 | −0.592 36 | −0.591 48 | −0.606 43 | −0.585 35 | −0.585 25 | −0.604 04 | −0.597 23 | −0.592 05 | −0.579 99 |

. | . | Excitation energy (eV) . | |
---|---|---|---|

State . | R_{eq} (Å)
. | Calculated . | Experimental . |

$X\Sigma g+1$ | 1.242 53 | 0 | 0 |

a^{3}Π_{u} | 1.312 | 0.07 | 0.09 |

$b3\Sigma g\u2212$ | 1.369 | 0.78 | 0.80 |

A^{1}Π_{u} | 1.318 | 1.03 | 1.04 |

$c3\Sigma u+$ | 1.208 | 1.16 | 1.13 |

B^{1}Δ_{g} | 1.385 | 1.49 | 1.50 |

$B\u2032\Sigma g+1$ | 1.377 | 1.90 | 1.91 |

d^{3}Π_{g} | 1.266 | 2.50 | 2.48 |

C^{1}Π_{g} | 1.255 | 4.29 | 4.25 |

. | . | Excitation energy (eV) . | |
---|---|---|---|

State . | R_{eq} (Å)
. | Calculated . | Experimental . |

$X\Sigma g+1$ | 1.242 53 | 0 | 0 |

a^{3}Π_{u} | 1.312 | 0.07 | 0.09 |

$b3\Sigma g\u2212$ | 1.369 | 0.78 | 0.80 |

A^{1}Π_{u} | 1.318 | 1.03 | 1.04 |

$c3\Sigma u+$ | 1.208 | 1.16 | 1.13 |

B^{1}Δ_{g} | 1.385 | 1.49 | 1.50 |

$B\u2032\Sigma g+1$ | 1.377 | 1.90 | 1.91 |

d^{3}Π_{g} | 1.266 | 2.50 | 2.48 |

C^{1}Π_{g} | 1.255 | 4.29 | 4.25 |

## VIII. CONCLUSIONS AND OUTLOOK

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.

## ACKNOWLEDGMENTS

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.

## REFERENCES

In practice, the selection of determinants is improved by using a larger $\u03f51$ during the early HCI iterations.