We introduce vibrational heat-bath configuration interaction (VHCI) as an accurate and efficient method for calculating vibrational eigenstates of anharmonic systems. Inspired by its origin in electronic structure theory, VHCI is a selected CI approach that uses a simple criterion to identify important basis states with a pre-sorted list of anharmonic force constants. Screened second-order perturbation theory and simple extrapolation techniques provide significant improvements to variational energy estimates. We benchmark VHCI on four molecules with 12–48 degrees of freedom and use anharmonic potential energy surfaces truncated at fourth and sixth orders. When compared to other methods using the same truncated potentials, VHCI produces vibrational spectra of tens or hundreds of states with sub-wavenumber accuracy at low computational cost.

## I. INTRODUCTION

Precise computational predictions of the vibrational structure of molecules and solids require the inclusion of anharmonic effects, which correspond to interactions between harmonic normal modes. When treated quantum mechanically, this requires the accurate solution of the vibrational Schrödinger equation on a high-dimensional potential energy surface (PES). Like in electronic structure theory, a hierarchy of wavefunction-based methods are commonly employed to avoid the exponential cost associated with an exact quantum solution; such methods include vibrational self-consistent field (VSCF) theory,^{1–6} vibrational perturbation theory,^{7–9} vibrational coupled-cluster theory,^{10–12} and vibrational configuration interaction (VCI).^{13–15} These methods rely on the accuracy of the harmonic approximation, and alternative methods have recently been developed to achieve quantitative accuracy for strongly anharmonic vibrations. Such methods include the nonproduct quadrature approach;^{16} reduced rank block power method (RRBPM), which uses a tensor factorization of the vibrational wavefunction;^{17} and adaptive VCI (A-VCI), which accelerates the VCI process using nested basis functions.^{18}

For strongly correlated electronic systems, the past two decades have seen the development of powerful computational methods that could be brought to bear on strongly anharmonic vibrational systems. For example, a vibrational density matrix renormalization group approach was recently developed and applied to systems with up to 16 atoms.^{19} In this work, we are primarily interested in selected configuration interaction,^{20–22} which began in the 1970s with CI by perturbatively selecting iteratively (CIPSI)^{20} and has experienced renewed interest in the form of adaptive configuration interaction,^{23} adaptive sampling CI (ASCI),^{24,25} and heat-bath CI (HCI),^{26–29} among others. In CIPSI, many-body basis states are added to the variational CI space based on a first-order approximation of their wavefunction coefficients. Unlike CIPSI, which considers all states belonging to the first-order interacting space, ASCI only considers connections from the variational basis states that have sufficiently large wavefunction amplitudes. HCI uses a different selection criterion that allows it to exploit the fact that many Hamiltonian matrix elements are identical in magnitude; a presorting of these matrix elements (the two-electron repulsion integrals) enables fast and efficient identification of new basis states. In almost all selected CI calculations, once the variational space is suitably converged, a second-order perturbation theory (PT2) correction is added to the variational energy. Numerous recent studies have demonstrated their power: HCI was selected as the most accurate method among 20 for a study of transition metal atoms and their oxides,^{30} HCI was used to provide almost exact energies of the Gaussian-2 dataset,^{31} and both HCI and ASCI have been used in a recent comparative study of the ground-state energy of benzene.^{32}

A variety of selected CI approaches have been developed for vibrational problems,^{33–36} including a vibrational CIPSI^{37} and a recent vibrational ASCI (VASCI).^{38} Motivated by HCI’s strong performance and computational efficiency, in this work, we present a vibrational HCI (VHCI). The differences between electronic and vibrational problems, especially their Hilbert spaces and Hamiltonian forms, necessitate a number of new developments in order for VHCI to enjoy the same advantages as electronic HCI; we describe these developments in Sec. II. In Sec. III, we present VHCI plus perturbation theory results for molecules containing between 12 and 48 degrees of freedom, calculating tens to hundreds of excited state energies for fourth- and sixth-order truncated potentials. By comparing to other literature results for these potentials, we demonstrate that VHCI is a highly accurate and efficient approach for large molecular systems with strong anharmonicity, typically achieving sub-wavenumber accuracy with modest computational resources.

## II. THEORY

### A. Vibrational heat-bath configuration interaction

Expressed in terms of mass-weighted normal mode coordinates *Q*_{i} with frequencies *ω*_{i}, the nuclear Hamiltonian is given by

where *D* = 3*N* − 6 is the number of normal mode degrees of freedom and *V*_{an} is the anharmonic part of the electronic ground state potential energy surface (PES). For comparison with previous works, here we neglect the Coriolis rotational coupling, but such a term can be straightforwardly included.

In vibrational HCI, we choose to work in the basis of Hartree product states |** n**⟩ = |

*n*

_{1},

*n*

_{2}, …,

*n*

_{D}⟩ formed from the harmonic oscillator orbitals $\varphi ni(Qi)=\u2009Qi|ni\u2009$, leading to the wavefunction

In VCI, the variational space $V$ is commonly truncated by the excitation level, for example, by limiting the number of vibrational quanta per product state or per mode. However, this approach is not always efficient and may lead to a large Hilbert space with many product states that contribute minimally to the VCI solution. This issue of inefficient addition of states is addressed by selected VCI methods, which select states for inclusion in $V$ by criteria other than the excitation level as discussed in the Introduction, and here, we focus on the HCI criterion.

We briefly recall that in the variational stage of HCI, basis state |** m**⟩ is added to $V$ if |

*H*

_{mn}

*c*

_{n}| ≥

*ε*

_{1}, where

*ε*

_{1}is a user-defined convergence parameter. In electroinic structure theory, the matrix elements

*H*

_{mn}arising from double excitations depend only on the identity of the four orbitals involved. Therefore, many of the Hamiltonian matrix elements are identical in magnitude. This is not true in vibrational structure theory, and the explicit enumeration and testing of connected

*H*

_{mn}matrix elements are expensive. In order to devise an analogous VHCI screening criterion, we consider an approximate Hamiltonian $H\u0303$ obtained from a truncated normal mode expansion of the anharmonic part of the PES,

where the anharmonic force constants are partial derivatives of the PES, e.g., $Vijk=(\u22023Van/\u2202Qi\u2202Qj\u2202Qk)Q=0$. In the second line of Eq. (3), we define $Q\xafi=(ai\u2020+ai)$ and

where *p* is the order of the anharmonicity. We emphasize that while we limit ourselves to using $H\u0303$ throughout this study, in principle, this normal mode expansion is only necessary for the screening criterion used to identify important basis states; any representation of the anharmonic PES could be used for the subsequent evaluation of Hamiltonian matrix elements within this space.

To illustrate the diversity of nonzero Hamiltonian matrix elements, consider an anharmonic PES expanded to include cubic and quartic terms. The Hamiltonian matrix elements between product states |** m**⟩ and |

**⟩ that differ in their occupancy by one quantum in mode**

*n**i*, by two quanta in modes

*i*and

*j*(including

*i*=

*j*), and so on are given by

The factors $\u2009m|Q\xafi\u2026|n\u2009$ that are not zero are products of terms of the form $ni$, and these factors are responsible for the differentiation of most matrix elements in the Hamiltonian, precluding an efficient evaluation of the usual HCI criterion. For example, if the product states |** m**⟩ and |

**⟩ differ in the occupancy of three modes (**

*n**i*,

*j*,

*k*), then

and so on, where the notation indicates that the occupation of a given mode in |** m**⟩ is bigger or smaller than in |

**⟩. Therefore, in addition to the**

*n**identity*of the modes with different occupations, the VCI matrix elements of a truncated PES depend on (1) whether the occupations are bigger/smaller on the bra/ket side and (2) the overall excitation level. However, because these latter factors $\u2009m|Q\xafi\u2026|n\u2009$ are products of terms of the form $ni$, they are of order unity (at least when the excitation level is reasonably low) and can be reasonably neglected for the purpose of estimating the magnitude of the Hamiltonian matrix element. Doing so leads to a significant reduction in the number of unique numbers to be considered.

Specifically, let us define *W*_{i} = ∑_{j}(2 + *δ*_{ij})*W*_{ijj} and *W*_{ij} = ∑_{k}(2 + *δ*_{jk} + *δ*_{ij}*δ*_{jk})*W*_{ijkk}, which approximately account for single and double mode excitations. Then, given an approximate wavefunction of form (2) expanded over some space $V$, we propose to add state |** m**⟩ to the variational space if |

**⟩ and |**

*m***⟩ differ in their occupancy by one quantum in mode**

*n**i*, by two quanta in modes

*i*and

*j*(including

*i*=

*j*), and so on, and

The above criterion can be implemented very efficiently by pre-sorting the matrix elements *W*_{ijk}…, of which there are at most a polynomial number, e.g., *O*(*D*^{4}) for a quartic PES.

The ground-state VHCI algorithm is performed as follows:

Sort the list of

*W*_{ijk…}, largest to smallest.Initialize the space $V$ according to a small total number of quanta.

Build the Hamiltonian matrix in $V$ and calculate its ground-state eigenvalue and eigenvector.

Add states to $V$ via criterion (7) as follows. For each state |

⟩, traverse the sorted list of*n**W*_{ijk…}.If |

*W*_{i}*c*_{n}| ≥*ε*_{1}, add all states |⟩ that are included in*m**Q*_{i}|⟩.*n*If |

*W*_{ij}*c*_{n}| ≥*ε*_{1}, add all states |⟩ that are included in*m**Q*_{i}*Q*_{j}|⟩.*n*If |

*W*_{ijk}, …,*c*_{n}| ≥*ε*_{1}, add all states |⟩ that are included in*m**Q*_{i}*Q*_{j}*Q*_{k}, …, |⟩. At anharmonic order*n**p*, there are*O*(2^{p}) such states to add.If |

*W*_{ijk}, …,*c*_{n}| <*ε*_{1}, break and go on to the next |⟩.*n*

Return to step 3 until the calculation is converged.

This line is part of the algorithm. Specifically, it should appear in item 4, after item (d). Close space up above

A number of possible convergence criteria can be devised; here, following the original HCI paper,^{26} we consider the calculation converged when the total number of states added in step 4 is less than 1% of the current number of variational states. Importantly, most elements of $H\u0303mn$ are zero due to the properties of harmonic oscillators; the states |** m**⟩ and |

**⟩ can only differ by a few quanta, depending on the maximum anharmonic order of the PES. The same would not be true in the basis of product states obtained after a vibrational self-consistent field procedure,**

*n*^{1–6}which is why we intentionally work in the noninteracting normal mode basis. However, future work will be dedicated to the development of VHCI in alternative mean-field basis sets such as the one generated by the VSCF procedure.

Crucially, the presorting of scaled and/or summed anharmonic force constants (*W*_{i}, *W*_{ij}, *W*_{ijk}, …, ) combined with criterion (7) means that *most states* |** m**⟩

*are never even tested for addition*. Just like in electronic structure theory, this construction is the key to the efficiency of VHCI. Unlike in electronic structure theory, the ground state of the vibrational Schrödinger equation, referred to here as the zero-point energy (ZPE), is rarely of interest by itself. Following Ref. 27, we can straightforwardly modify the above VHCI algorithm to allow the simultaneous calculation of many excited states. In step 3, we find all eigenstates of interest (typically the

*N*

_{s}lowest in energy). Then, we perform the addition step 4 for each of those states, combining all of their added basis states |

**⟩ in the updated $V$. Clearly, this adds many more basis states at each iteration, but many of them are duplicates, and so the overall variational space is observed to grow sublinearly with**

*m**N*

_{s}.

The above procedure can be applied based on an arbitrary-order expansion of the PES. However, high-order anharmonic interactions with repeated mode indices will contribute to lower-order excitations beyond the single and double excitations described in Eqs. (5a) and (5b) for the case of a quartic PES. For example, a sixth-order PES will yield triple excitations due to cubic anharmonicity *and* fifth-order anharmonicity (when a mode index is repeated) and similarly quadruple excitations due to quartic and sixth-order anharmonicity. In this work, we only test one sixth-order PES in Sec. III D; although we could define new effective force constants for the screening procedure, e.g., $Wijk\u2032$ = *W*_{ijk} + ∑_{k}*W*_{ijkll}, instead we make the approximation ∑_{l}*W*_{ijkll} ≈ max_{l}*W*_{ijkll}. In other words, we allow fifth-order anharmonicities to propose triple excitations but only based on the magnitude of individual anharmonic force constants and not the sum of all such constants contributing to a given triple excitation. Our testing suggests that the error incurred is negligible.

### B. Epstein–Nesbet perturbation theory

To the variational energy of each state *E*_{var}, we add the second-order perturbation theory (PT2) correction

In exact PT2, the summation over ** m** includes all basis states that are connected to the variational space $V$. For large variational spaces, this perturbative space is enormous and the cost of the PT2 correction is prohibitive. To address this, HCI uses the same screening procedure as in the variational stage to efficiently include only the most important perturbative states.

^{26}For VHCI, we again use criterion (7), with a cutoff

*ε*

_{2}<

*ε*

_{1}, to determine whether basis state |

**⟩ should be included in the perturbative space. When**

*m**ε*

_{2}= 0, the PT2 calculation is exact within the first-order interacting space.

Throughout this work, we calculate the PT2 correction deterministically as described above, which ultimately limits the size of the systems that we can accurately study. In the future, we will pursue the stochastic or semistochastic evaluation of Eq. (8), as is now common practice in HCI,^{28} in order to study the vibrational structure of even larger systems.

## III. RESULTS

### A. Software and simulation details

All simulations besides those on naphthalene were performed on a 4-core (8-thread) Intel Core i7-6700 3.4 GHz desktop central processing unit (CPU) using up to 16 GB of RAM. Naphthalene calculations were performed on a cluster with up to two 12-core Intel Xeon Gold 6126 2.6 GHz CPUs and using up to 768 GB of RAM.

Our code is based on the Ladder Operator Vibrational Configuration Interaction package^{39} with extensive modifications. Our VHCI code is available on GitHub.^{40} The Hamiltonian matrix is stored in a sparse format using the Eigen linear algebra library.^{41} The eigenvalues and eigenvectors of interest are calculated with the Lanczos algorithm as implemented in the SPECTRA linear algebra library.^{42} We verified our code by comparing to the literature and to the results obtained with the PyVCI software package.^{36}

All VHCI calculations were initialized with a basis of all states containing up to two vibrational quanta in order to ensure that two-quantum overtones and combination bands, which account for many of the low-lying target states, are present from the first iteration. Our preliminary testing indicates that using a larger initial basis does not qualitatively improve the convergence of the results.

In order to fairly and directly compare to previous works, we use the exact same PES as used in those works. Moreover, these previous works used truncated forms of the PES, and thus, the true Hamiltonian *H* used in the final calculations is the same as the approximate Hamiltonian $H\u0303$ used in the VHCI screening criterion. The truncation of the PES, neglect of Coriolis terms, and various levels of electronic structure theory used to parameterize the PES preclude direct comparison to experiment.

### B. Acetonitrile: A standard benchmark

We first present results for acetonitrile, CH_{3}CN, a 6-atom, 12-dimensional system that has become one of the canonical benchmarks for new methods in vibrational structure theory.^{17,19,38,43,44} We use the quartic PES described in Refs. 17 and 45, for which normal mode frequencies were calculated using CCSD(T)/cc-pVTZ and higher-order force constants were calculated using B3LYP/cc-pVTZ. This PES was also used in highly accurate reference calculations.^{16,46} All PESs used in this study are available as input files on our GitHub repository.^{40}

We calculated the energies of the first 70 states of acetonitrile, of which the first 20 are reported in Table I. Data for all 70 states can be found in the supplementary material. We compare our results to those obtained using A-VCI,^{46} which we consider to be numerically exact (similar to those obtained with the nonproduct quadrature approach^{16}); for reference, the A-VCI results are obtained with ∼2.5 × 10^{6} variational states. Mode assignments here and throughout indicate the character of the product state with the largest weight in the VCI eigenvector; for acetonitrile, we use the mode-numbering convention of Refs. 16 and 38. Because acetonitrile contains several degeneracies, we also include labels for the total vibrational angular momentum and symmetry group. To quantitatively assess the accuracy of our results, we report the maximum absolute error and the root mean squared error of the lowest 70 states. In addition to the numerically exact A-VCI results, we also compare to results obtained recently using VASCI,^{38} which is a selected CI technique that is similar in spirit to VHCI.

Method . | . | ε_{1} = 1.0
. | ε_{1} = 0.1
. | VASCI^{38}
. | A-VCI^{18}
. | ||
---|---|---|---|---|---|---|---|

$NV$ . | . | 29 900 . | 125 038 . | 30 038 . | 2 488 511 . | ||

Assign. . | Sym. . | Var . | Full PT2 . | Var . | Full PT2 . | Full PT2 . | Var . |

ZPE | 9 837.43 | 9 837.41 | 9 837.41 | 9 837.41 | 9 837.41 | 9 837.41 | |

$\omega 11\xb11$ | E | 361.01 | 360.99 | 360.99 | 360.99 | 361.01 | 360.99 |

361.01 | 360.99 | 360.99 | 360.99 | 361.01 | 360.99 | ||

$2\omega 11\xb12$ | E | 723.22 | 723.18 | 723.18 | 723.18 | 723.22 | 723.18 |

723.25 | 723.19 | 723.18 | 723.18 | 723.23 | 723.18 | ||

$2\omega 110$ | A_{1} | 723.90 | 723.84 | 723.83 | 723.83 | 723.89 | 723.83 |

$\omega 4\xb11$ | A_{1} | 900.70 | 900.66 | 900.66 | 900.66 | 900.68 | 900.66 |

$\omega 9\xb11$ | E | 1 034.18 | 1 034.13 | 1 034.13 | 1 034.12 | 1 034.14 | 1 034.12 |

1 034.19 | 1 034.13 | 1 034.13 | 1 034.12 | 1 034.14 | 1 034.13 | ||

$3\omega 11\xb13$ | A_{1} | 1 086.65 | 1 086.56 | 1 086.56 | 1 086.55 | 1 086.64 | 1 086.55 |

$3\omega 11\xb13$ | A_{2} | 1 086.66 | 1 086.56 | 1 086.56 | 1 086.55 | 1 086.64 | 1 086.55 |

$3\omega 11\xb11$ | E | 1 087.88 | 1 087.79 | 1 087.78 | 1 087.78 | 1 087.88 | 1 087.77 |

1 087.88 | 1 087.79 | 1 087.78 | 1 087.78 | 1 087.88 | 1 087.77 | ||

$\omega 4+\omega 11\xb11$ | E | 1 259.89 | 1 259.82 | 1 259.81 | 1 259.81 | 1 259.86 | 1 259.81 |

1 259.94 | 1 259.83 | 1 259.82 | 1 259.81 | 1 259.86 | 1 259.81 | ||

ω_{3} | A_{1} | 1 389.10 | 1 388.99 | 1 388.98 | 1 388.97 | 1 388.99 | 1 388.97 |

$\omega 9\xb11+\omega 11\xb11$ | E | 1 394.86 | 1 394.71 | 1 394.69 | 1 394.68 | 1 394.76 | 1 394.68 |

1 394.89 | 1 394.71 | 1 394.70 | 1 394.68 | 1 394.79 | 1 394.68 | ||

$\omega 9\xb11+\omega 11\u22131$ | A_{2} | 1 395.08 | 1 394.93 | 1 394.91 | 1 394.90 | 1 395.00 | 1 394.90 |

$\omega 9\xb11+\omega 11\u22131$ | A_{1} | 1 397.83 | 1 397.70 | 1 397.69 | 1 397.68 | 1 397.77 | 1 397.68 |

Max. AE(70) | 0.61 | 0.14 | 0.05 | 0.01 | 0.50 | … | |

RMSE(70) | 0.34 | 0.05 | 0.03 | 0.00 | 0.20 | … |

Method . | . | ε_{1} = 1.0
. | ε_{1} = 0.1
. | VASCI^{38}
. | A-VCI^{18}
. | ||
---|---|---|---|---|---|---|---|

$NV$ . | . | 29 900 . | 125 038 . | 30 038 . | 2 488 511 . | ||

Assign. . | Sym. . | Var . | Full PT2 . | Var . | Full PT2 . | Full PT2 . | Var . |

ZPE | 9 837.43 | 9 837.41 | 9 837.41 | 9 837.41 | 9 837.41 | 9 837.41 | |

$\omega 11\xb11$ | E | 361.01 | 360.99 | 360.99 | 360.99 | 361.01 | 360.99 |

361.01 | 360.99 | 360.99 | 360.99 | 361.01 | 360.99 | ||

$2\omega 11\xb12$ | E | 723.22 | 723.18 | 723.18 | 723.18 | 723.22 | 723.18 |

723.25 | 723.19 | 723.18 | 723.18 | 723.23 | 723.18 | ||

$2\omega 110$ | A_{1} | 723.90 | 723.84 | 723.83 | 723.83 | 723.89 | 723.83 |

$\omega 4\xb11$ | A_{1} | 900.70 | 900.66 | 900.66 | 900.66 | 900.68 | 900.66 |

$\omega 9\xb11$ | E | 1 034.18 | 1 034.13 | 1 034.13 | 1 034.12 | 1 034.14 | 1 034.12 |

1 034.19 | 1 034.13 | 1 034.13 | 1 034.12 | 1 034.14 | 1 034.13 | ||

$3\omega 11\xb13$ | A_{1} | 1 086.65 | 1 086.56 | 1 086.56 | 1 086.55 | 1 086.64 | 1 086.55 |

$3\omega 11\xb13$ | A_{2} | 1 086.66 | 1 086.56 | 1 086.56 | 1 086.55 | 1 086.64 | 1 086.55 |

$3\omega 11\xb11$ | E | 1 087.88 | 1 087.79 | 1 087.78 | 1 087.78 | 1 087.88 | 1 087.77 |

1 087.88 | 1 087.79 | 1 087.78 | 1 087.78 | 1 087.88 | 1 087.77 | ||

$\omega 4+\omega 11\xb11$ | E | 1 259.89 | 1 259.82 | 1 259.81 | 1 259.81 | 1 259.86 | 1 259.81 |

1 259.94 | 1 259.83 | 1 259.82 | 1 259.81 | 1 259.86 | 1 259.81 | ||

ω_{3} | A_{1} | 1 389.10 | 1 388.99 | 1 388.98 | 1 388.97 | 1 388.99 | 1 388.97 |

$\omega 9\xb11+\omega 11\xb11$ | E | 1 394.86 | 1 394.71 | 1 394.69 | 1 394.68 | 1 394.76 | 1 394.68 |

1 394.89 | 1 394.71 | 1 394.70 | 1 394.68 | 1 394.79 | 1 394.68 | ||

$\omega 9\xb11+\omega 11\u22131$ | A_{2} | 1 395.08 | 1 394.93 | 1 394.91 | 1 394.90 | 1 395.00 | 1 394.90 |

$\omega 9\xb11+\omega 11\u22131$ | A_{1} | 1 397.83 | 1 397.70 | 1 397.69 | 1 397.68 | 1 397.77 | 1 397.68 |

Max. AE(70) | 0.61 | 0.14 | 0.05 | 0.01 | 0.50 | … | |

RMSE(70) | 0.34 | 0.05 | 0.03 | 0.00 | 0.20 | … |

We report results for variational VHCI (“Var”), as well as VHCI+PT2 without VHCI screening of the perturbative space (“Full PT2”); these results are given for two values of the variational energy cutoff *ε*_{1} that controls the number of variational states $NV$. Using *ε*_{1} = 1 cm^{−1} results in a variational space with $NV=29\u2009900$ basis states, which enables comparison to the largest reported VASCI calculation, with $NV=30\u2009038$. We find that VHCI+PT2 achieves a maximum absolute error and a root mean squared error that is less than half that of VASCI, which is somewhat surprising because the CIPSI-style selection criterion used in VASCI is typically thought to more rigorously identify important states. However, the precise variational space generated by both VHCI and VASCI can be tuned by their parameters (*ε*_{1} and the number of core/target states, respectively), so a direct comparison is not straightforward. In any event, the accuracy of VHCI is remarkable, given the modest computational cost; for example, the *ε*_{1} = 1 cm^{−1} calculation reported here took less than 3 min for the variational stage and less than 6 min for the full PT2 correction (on an 8-core desktop CPU).

In Table I, we also present results of a larger calculation with *ε*_{1} = 0.1 cm^{−1}, resulting in $NV=125\u2009038$. We obtain an extremely accurate spectrum even without the PT2 correction, with a maximum absolute error well below 0.1 cm^{−1}. This variational calculation took less than 30 min on the same 8-core desktop CPU. This larger calculation still benefits from minor corrections with PT2, which takes about 90 min and produces a spectrum in which no value differs from the exact result by more than 0.01 cm^{−1} across all 70 states. In summary, we have shown that VHCI can produce near-exact results for a 12-dimensional system with an extremely small computational effort.

### C. Ethylene: Sixth-order potential

Next, we study ethylene, C_{2}H_{4}, which is the same size as acetonitrile (6 atoms, 12 dimensions), but here we use a more realistic *ab initio* potential energy surface and consider anharmonic expansions up to sixth order. Specifically, we use the PES of Ref. 47, which was calculated entirely using CCSD(T)/cc-pVQZ, and we use the PyPES software package^{48} to convert from internal to Cartesian normal mode coordinates and generate an anharmonic expansion up to sixth order. In contrast to the quartic PES of acetonitrile that contains 299 nonzero anharmonic force constants, this sixth-order PES for ethylene contains 2651 force constants, 2375 of which are fifth or sixth order. Here, we compare results and performance when the potential is truncated at fourth order and at sixth order. We calculated the first 100 states; results for the first 16 states are shown in Table II and those for the additional 84 states can be found in the supplementary material.

PES . | Fourth-order . | Sixth-order . | Untruncated . | |||
---|---|---|---|---|---|---|

Method . | VHCI+Full PT2 . | . | VHCI+PT2 (ε_{2} = 0.01)
. | . | VCI(PyVCI)^{48}
. | Pruned VCI^{47}
. |

Assign. . | ε_{1} = 0.75
. | vDMRG^{19}
. | ε_{1} = 0.5
. | vDMRG^{19}
. | N_{tot} = 8
. | N_{tot} = 13
. |

ZPE | 11 006.11 | 11 006.19 | 11 011.61 | 11 016.15 | 11 011.63 | 11 014.91 |

ω_{10} | 808.61 | 809.03 | 819.99 | 831.17 | 820.11 | 822.42 |

ω_{8} | 914.87 | 915.29 | 926.33 | 933.47 | 926.45 | 934.29 |

ω_{7} | 927.87 | 928.31 | 941.65 | 948.26 | 941.78 | 949.51 |

ω_{4} | 1 006.74 | 1 007.13 | 1 017.45 | 1 018.26 | 1 017.56 | 1 024.94 |

ω_{6} | 1 216.94 | 1 217.17 | 1 222.16 | 1 227.05 | 1 222.23 | 1 224.96 |

ω_{3} | 1 338.46 | 1 338.87 | 1 341.95 | 1 343.46 | 1 342.01 | 1 342.96 |

ω_{12} | 1 429.93 | 1 430.47 | 1 438.31 | 1 441.52 | 1 438.39 | 1 441.11 |

ω_{2} | 1 606.41 | 1 622.11 | 1 622.78 | 1 629.04 | … | 1 624.43 |

2ω_{10} | 1 631.47 | 1 625.56 | 1 655.21 | 1 682.18 | … | 1 658.39 |

ω_{8} + ω_{10} | 1 718.27 | 1 722.77 | 1 748.05 | 1 770.02 | … | 1 757.70 |

ω_{7} + ω_{10} | 1 733.98 | 1 729.53 | 1 766.19 | 1 786.99 | … | 1 778.34 |

ω_{4} + ω_{10} | 1 809.39 | 1 810.47 | 1 837.68 | 1 852.60 | … | 1 848.61 |

2ω_{8} | 1 821.71 | 1 826.01 | 1 858.36 | 1 878.42 | … | 1 873.73 |

ω_{7} + ω_{8} | 1 821.96 | 1 827.96 | 1 871.42 | 1 886.16 | … | 1 885.12 |

2ω_{7} | 1 850.39 | 1 852.23 | 1 887.03 | 1 906.02 | … | 1 901.61 |

PES . | Fourth-order . | Sixth-order . | Untruncated . | |||
---|---|---|---|---|---|---|

Method . | VHCI+Full PT2 . | . | VHCI+PT2 (ε_{2} = 0.01)
. | . | VCI(PyVCI)^{48}
. | Pruned VCI^{47}
. |

Assign. . | ε_{1} = 0.75
. | vDMRG^{19}
. | ε_{1} = 0.5
. | vDMRG^{19}
. | N_{tot} = 8
. | N_{tot} = 13
. |

ZPE | 11 006.11 | 11 006.19 | 11 011.61 | 11 016.15 | 11 011.63 | 11 014.91 |

ω_{10} | 808.61 | 809.03 | 819.99 | 831.17 | 820.11 | 822.42 |

ω_{8} | 914.87 | 915.29 | 926.33 | 933.47 | 926.45 | 934.29 |

ω_{7} | 927.87 | 928.31 | 941.65 | 948.26 | 941.78 | 949.51 |

ω_{4} | 1 006.74 | 1 007.13 | 1 017.45 | 1 018.26 | 1 017.56 | 1 024.94 |

ω_{6} | 1 216.94 | 1 217.17 | 1 222.16 | 1 227.05 | 1 222.23 | 1 224.96 |

ω_{3} | 1 338.46 | 1 338.87 | 1 341.95 | 1 343.46 | 1 342.01 | 1 342.96 |

ω_{12} | 1 429.93 | 1 430.47 | 1 438.31 | 1 441.52 | 1 438.39 | 1 441.11 |

ω_{2} | 1 606.41 | 1 622.11 | 1 622.78 | 1 629.04 | … | 1 624.43 |

2ω_{10} | 1 631.47 | 1 625.56 | 1 655.21 | 1 682.18 | … | 1 658.39 |

ω_{8} + ω_{10} | 1 718.27 | 1 722.77 | 1 748.05 | 1 770.02 | … | 1 757.70 |

ω_{7} + ω_{10} | 1 733.98 | 1 729.53 | 1 766.19 | 1 786.99 | … | 1 778.34 |

ω_{4} + ω_{10} | 1 809.39 | 1 810.47 | 1 837.68 | 1 852.60 | … | 1 848.61 |

2ω_{8} | 1 821.71 | 1 826.01 | 1 858.36 | 1 878.42 | … | 1 873.73 |

ω_{7} + ω_{8} | 1 821.96 | 1 827.96 | 1 871.42 | 1 886.16 | … | 1 885.12 |

2ω_{7} | 1 850.39 | 1 852.23 | 1 887.03 | 1 906.02 | … | 1 901.61 |

Through testing, we confirmed that our VHCI+PT2 results are converged with *ε*_{1} = 0.75 cm^{−1} for the quartic case (with full PT2) and *ε*_{1} = 0.5 cm^{−1} for the sixth-order case (with *ε*_{2} = 0.01 cm^{−1}), resulting in $NV=$ 153 935 and 161 338 basis states, respectively. In fact, it was not possible to go to larger variational spaces for the quartic case because of unphysical divergences in the truncated PES. Similar behavior of truncated PESs studied at high excitation levels has been observed before.^{36,47,49}

In Table II, we also compare our results to those obtained with the vibrational density-matrix renormalization group (vDMRG)^{19} using the same fourth- and sixth-order truncated PESs. As discussed in the Introduction, DMRG and selected CI are similarly competitive methods in electronic structure theory. Surprisingly, we find that for both truncations, our energies are noticeably lower than the vDMRG results. For the quartic potential, our zero-point energy is 0.8 cm^{−1} lower than the vDMRG result, while for the sixth-order potential, it is 4.5 cm^{−1} lower. Figure 1 shows the convergence of the ZPE as a function of the size of the variational space, $NV$. We see that VHCI (obtained with *N*_{s} = 1) converges smoothly and quickly; for the fourth-order potential, it exceeds the accuracy of vDMRG with about 5000 basis states. We also show results obtained with our own VCI code, which includes basis states according to their total number of vibrational quanta. For both the fourth- and sixth-order potentials, VCI converges to the same ZPE as VHCI, although it requires significantly more basis states for comparable accuracy. The discrepancy with vDMRG may come from the latter using insufficient bond dimension, getting trapped in local minima during sweeps, or be simply due to slight differences in the PES. As a check on the latter, we have also included in Table II the VCI excitation energies previously reported by the authors of the PyPES software package^{36} from which the PES parameters were obtained. The results are in excellent agreement with our own VHCI+PT2 results, indicating that we are using the exact same PES as those authors.

In Fig. 1, we also plot the percent sparsity of the Hamiltonian matrix as a function of the number of variational basis states. We present the sparsity for VHCI (optimized for the ground state) and conventional VCI. For both the fourth-order and sixth-order potentials, VHCI produces a much sparser Hamiltonian matrix than VCI. These results indicate that VHCI is extremely effective at capturing the connectivity between important basis states (as demonstrated by the accurate ground-state energy) while also ignoring the connectivity to less important basis states (as demonstrated by the increased sparsity).

In the final column of Table II, we also show the results obtained in Ref. 47 using a pruned VCI basis with up to 13 vibrational quanta for the untruncated PES. As can be seen, the agreement improves significantly when going from the fourth-order PES to the sixth-order PES. For low-lying states, the disagreement at fourth order is on the order of 10 cm^{−1}–20 cm^{−1} and at sixth order is on the order of 5 cm^{−1}, indicating that the normal mode expansion is sensible and systematically improvable. Nonetheless, this error incurred due to the truncated potential is comparable to the difference in methods presented here, and thus, future work will focus on incorporating more accurate representations of the potential into the VHCI method.

### D. Ethylene oxide: Convergence and extrapolation

Next, we study ethylene oxide, C_{2}H_{4}O, a molecule with 7 atoms and 15 normal mode degrees of freedom. Compared to the two previous molecules, the three additional degrees of freedom make numerically exact VCI calculations much more difficult. We use the quartic PES of Bégué *et al.*,^{50} with normal mode frequencies calculated at the CCSD(T)/cc-pVTZ level and anharmonic force constants calculated using B3LYP/6-31+G(d,p). We use the same version of the PES as Refs. 18, 38, 43, and 44, which is available along with our source code on our GitHub repository.^{40} Like acetonitrile, ethylene oxide represents a well-studied system that is suitable for benchmarking new approaches to solving the vibrational structure problem. In Table III, we present results for the ten lowest and ten highest states from a VHCI calculation targeting *N*_{s} = 200 states. We performed the calculations with two different variational cutoff values *ε*_{1} = 2 cm^{−1} and 1 cm^{−1}, producing variational spaces with 132 163 and 259 070 basis states, respectively. We present results with just the variational stage (“Var”) and with heat-bath screened PT2 correction with *ε*_{2} = 0.01 cm^{−1}.

Method . | ε_{1} = 2.0
. | ε_{1} = 1.0
. | Linear ΔE_{2}
. | VASCI^{38}
. | A-VCI^{18}
. | ||
---|---|---|---|---|---|---|---|

$NV$ . | 132 163 . | 259 070 . | Extrap. . | 150 000 . | 7 118 214 . | ||

Assign. . | Var . | ε_{2} = 0.01
. | Var . | ε_{2} = 0.01
. | ε_{2} = 0.01
. | Full PT2 . | Var . |

ZPE | 12 461.63 | 12 461.50 | 12 461.55 | 12 461.48 | 12 461.47 | 12 461.6 | 12 461.47 |

ω_{1} | 793.11 | 792.73 | 792.88 | 792.69 | 792.65 | 792.8 | 792.63 |

ω_{2} | 822.30 | 821.98 | 822.07 | 821.94 | 821.91 | 822.2 | 821.91 |

ω_{3} | 878.62 | 878.33 | 878.41 | 878.30 | 878.28 | 878.4 | 878.27 |

ω_{4} | 1 017.70 | 1 017.24 | 1 017.42 | 1 017.20 | 1 017.15 | 1 017.4 | 1 017.14 |

ω_{6} | 1 121.87 | 1 121.30 | 1 121.53 | 1 121.24 | 1 121.19 | 1 120.8 | 1 121.17 |

ω_{5} | 1 124.37 | 1 123.75 | 1 124.00 | 1 123.70 | 1 123.65 | 1 123.8 | 1 123.62 |

ω_{7} | 1 146.42 | 1 145.84 | 1 146.08 | 1 145.78 | 1 145.73 | 1 146.0 | 1 145.72 |

ω_{8} | 1 148.68 | 1 148.07 | 1 148.31 | 1 148.02 | 1 147.97 | 1 148.3 | 1 147.96 |

ω_{9} | 1 271.43 | 1 270.89 | 1 271.06 | 1 270.83 | 1 270.78 | 1 271.3 | 1 270.78 |

ω_{1} + ω_{5} + ω_{9} | 3 175.50 | 3 170.13 | 3 172.33 | 3 169.19 | … | 3 170.9 | 3 167.97 |

ω_{1} + ω_{6} + ω_{9} | 3 178.91 | 3 173.60 | 3 175.76 | 3 172.68 | … | 3 173.9 | 3 171.53 |

2ω_{4} + ω_{8} | 3 181.83 | 3 177.51 | 3 179.28 | 3 176.82 | … | 3 178.0 | 3 175.93 |

ω_{2} + ω_{3} + ω_{11} | 3 187.19 | 3 182.12 | 3 183.76 | 3 181.06 | … | 3 182.1 | 3 180.16 |

3ω_{1} + ω_{2} | 3 197.27 | 3 191.75 | 3 192.18 | 3 187.44 | … | 3 190.7 | 3 184.85 |

ω_{1} + ω_{8} + ω_{9} | 3 198.00 | 3 189.61 | 3 193.89 | 3 190.78 | … | 3 192.7 | 3 189.65 |

ω_{2} + ω_{5} + ω_{9} | 3 205.08 | 3 200.30 | 3 202.16 | 3 199.48 | … | … | 3 198.56 |

ω_{2} + ω_{6} + ω_{9} | 3 206.11 | 3 201.03 | 3 202.96 | 3 200.15 | … | 3 200.9 | 3 199.21 |

ω_{1} + ω_{7} + ω_{9} | 3 210.44 | 3 205.96 | 3 207.74 | 3 205.19 | … | 3 206.7 | 3 204.35 |

2ω_{1} + 2ω_{2} | … | … | 3 218.48 | 3 214.59 | … | … | 3 212.77 |

Max AE(50) | 4.00 | 0.98 | 2.09 | 0.50 | 0.06 | 2.2 | … |

RMSE(50) | 2.40 | 0.51 | 1.22 | 0.26 | 0.03 | 0.8 | … |

Max. AE(200) | … | … | 7.33 | 3.23 | … | … | … |

RMSE(200) | … | … | 2.96 | 0.84 | … | … | … |

Core hours | 5.9 | 18.9 | 23.3 | 54.4 | … | 67.1 | 1 756.4 |

Cores | 8 | 8 | 8 | 8 | … | 2 | 24 |

Method . | ε_{1} = 2.0
. | ε_{1} = 1.0
. | Linear ΔE_{2}
. | VASCI^{38}
. | A-VCI^{18}
. | ||
---|---|---|---|---|---|---|---|

$NV$ . | 132 163 . | 259 070 . | Extrap. . | 150 000 . | 7 118 214 . | ||

Assign. . | Var . | ε_{2} = 0.01
. | Var . | ε_{2} = 0.01
. | ε_{2} = 0.01
. | Full PT2 . | Var . |

ZPE | 12 461.63 | 12 461.50 | 12 461.55 | 12 461.48 | 12 461.47 | 12 461.6 | 12 461.47 |

ω_{1} | 793.11 | 792.73 | 792.88 | 792.69 | 792.65 | 792.8 | 792.63 |

ω_{2} | 822.30 | 821.98 | 822.07 | 821.94 | 821.91 | 822.2 | 821.91 |

ω_{3} | 878.62 | 878.33 | 878.41 | 878.30 | 878.28 | 878.4 | 878.27 |

ω_{4} | 1 017.70 | 1 017.24 | 1 017.42 | 1 017.20 | 1 017.15 | 1 017.4 | 1 017.14 |

ω_{6} | 1 121.87 | 1 121.30 | 1 121.53 | 1 121.24 | 1 121.19 | 1 120.8 | 1 121.17 |

ω_{5} | 1 124.37 | 1 123.75 | 1 124.00 | 1 123.70 | 1 123.65 | 1 123.8 | 1 123.62 |

ω_{7} | 1 146.42 | 1 145.84 | 1 146.08 | 1 145.78 | 1 145.73 | 1 146.0 | 1 145.72 |

ω_{8} | 1 148.68 | 1 148.07 | 1 148.31 | 1 148.02 | 1 147.97 | 1 148.3 | 1 147.96 |

ω_{9} | 1 271.43 | 1 270.89 | 1 271.06 | 1 270.83 | 1 270.78 | 1 271.3 | 1 270.78 |

ω_{1} + ω_{5} + ω_{9} | 3 175.50 | 3 170.13 | 3 172.33 | 3 169.19 | … | 3 170.9 | 3 167.97 |

ω_{1} + ω_{6} + ω_{9} | 3 178.91 | 3 173.60 | 3 175.76 | 3 172.68 | … | 3 173.9 | 3 171.53 |

2ω_{4} + ω_{8} | 3 181.83 | 3 177.51 | 3 179.28 | 3 176.82 | … | 3 178.0 | 3 175.93 |

ω_{2} + ω_{3} + ω_{11} | 3 187.19 | 3 182.12 | 3 183.76 | 3 181.06 | … | 3 182.1 | 3 180.16 |

3ω_{1} + ω_{2} | 3 197.27 | 3 191.75 | 3 192.18 | 3 187.44 | … | 3 190.7 | 3 184.85 |

ω_{1} + ω_{8} + ω_{9} | 3 198.00 | 3 189.61 | 3 193.89 | 3 190.78 | … | 3 192.7 | 3 189.65 |

ω_{2} + ω_{5} + ω_{9} | 3 205.08 | 3 200.30 | 3 202.16 | 3 199.48 | … | … | 3 198.56 |

ω_{2} + ω_{6} + ω_{9} | 3 206.11 | 3 201.03 | 3 202.96 | 3 200.15 | … | 3 200.9 | 3 199.21 |

ω_{1} + ω_{7} + ω_{9} | 3 210.44 | 3 205.96 | 3 207.74 | 3 205.19 | … | 3 206.7 | 3 204.35 |

2ω_{1} + 2ω_{2} | … | … | 3 218.48 | 3 214.59 | … | … | 3 212.77 |

Max AE(50) | 4.00 | 0.98 | 2.09 | 0.50 | 0.06 | 2.2 | … |

RMSE(50) | 2.40 | 0.51 | 1.22 | 0.26 | 0.03 | 0.8 | … |

Max. AE(200) | … | … | 7.33 | 3.23 | … | … | … |

RMSE(200) | … | … | 2.96 | 0.84 | … | … | … |

Core hours | 5.9 | 18.9 | 23.3 | 54.4 | … | 67.1 | 1 756.4 |

Cores | 8 | 8 | 8 | 8 | … | 2 | 24 |

On the left-hand side of Fig. 2, we plot the energy of the first, 50th, and 200th state with respect to the variational cutoff *ε*_{1} for a variety of values of the perturbative screening parameter *ε*_{2}. We see that the curves with *ε*_{2} = 1 cm^{−1}, 0.1 cm^{−1}, and 0.01 cm^{−1} agree very closely for all *ε*_{1}, with *ε*_{2} = 0.1 cm^{−1} and 0.01 cm^{−1} being indistinguishable, confirming convergence with respect to *ε*_{2}. In Table III and Fig. 2, we compare our results to A-VCI calculations^{18} that were obtained from an optimized variational space containing over 7 × 10^{6} basis states, which we take to be numerically exact. In Fig. 2, we see that all variational calculations tend monotonically toward the exact energies as *ε*_{1} becomes small. Addition of the PT2 correction leads to more rapid convergence with respect to *ε*_{1}. For example, for the ground state, we see that the results obtained with *ε*_{1} = 20 cm^{−1} and converged PT2 give an answer that is closer to exact than that with *ε*_{1} = 5 cm^{−1} and no PT2 correction. The benefits of perturbation theory become less pronounced in high-lying excited states due to their large multiconfigurational character. In fact, for the 200th state, our calculation only produces the correct band assignment for *ε*_{1} ≤ 1 cm^{−1}. Table III shows the band assignment (obtained from the variational calculation) following the mode labeling convention of Ref. 18, and we leave the energy blank if we do not have an assignment that corresponds to the exact result. The 200th state is omitted for all methods except VHCI at *ε*_{1} = 1 cm^{−1} and the A-VCI reference; several of the other ten highest excitations are also omitted from the VASCI results.^{38} Because the incorrect band assignments of high-lying states prevent us from comparing level-by-level with exact results, we include statistics on just the first 50 states. Variational VHCI produces good agreement for the first 50 states, with a RMSE on the order of 1 cm^{−1}–2 cm^{−1} and maximum AE below 5 cm^{−1}; the addition of PT2 yields accuracy better than 1 cm^{−1} for both *ε*_{1} = 1 cm^{−1} and 2 cm^{−1}. In comparison, VASCI+Full PT2 produces a maximum AE of 2.2 cm^{−1} and RMSE of 0.8 cm^{−1} for the first 50 levels compared to the exact results. All of the correct band assignments are present in the first 50 states of both VHCI and VASCI at $NV=150\u2009000$, but the closely spaced states 47 and 48 have the wrong energetic ordering for VASCI and for VHCI with *ε*_{1} = 2 cm^{−1}. For *ε*_{1} = 1 cm^{−1}, we obtain all 200 states with the correct assignments, enabling direct comparison to A-VCI. We see very good agreement over all states and VHCI+PT2 has a maximum AE around 3 cm^{−1} and RMSE below 1 cm^{−1}. Remarkably, this entire calculation took less than 7 h on an 8-core desktop CPU. As a rough comparison to other methods, we also present the timings for all calculations in Table III.

Finally, following standard HCI practice in the electronic structure, we attempt to approximate the exact energies via extrapolation. On the right-hand side of Fig. 2, we plot the VHCI+PT2 energy of each state as a function of the PT2 correction Δ*E*_{2} seeking to extrapolate to Δ*E*_{2} = 0, following Ref. 27. We see that extrapolation is reasonable and successful for low-lying states, but not for high-lying states. In general, extrapolation of high-lying states is difficult because the level spacing becomes very small and the states are highly multiconfigurational such that tracking a single level as *ε*_{1} changes is challenging. We performed extrapolation for the first 50 states using a linear fit to the results obtained with *ε*_{1} = 1 cm^{−1} and 2 cm^{−1}, as shown graphically on the right-hand sides of Fig. 2 for the first and 50th states; we tested other polynomial fits and found linear extrapolation to be the simplest and best performing, although other schemes could be considered in the future. In Table II, we present the extrapolated energies of the first ten states as well as the statistics of the first 50 states; all energies can be found in the supplementary material. The results are nearly exact, with a maximum AE below 0.1 cm^{−1} for all 50 states, and obviously require no additional computing effort. We conclude that linear extrapolation of high-quality VHCI+PT2 energies is a powerful way to achieve nearly exact energies for the ground state and many low-lying excited states.

### E. Naphthalene: A 48-dimensional system

As a final test of VHCI, we consider naphthalene, C_{10}H_{8}, with 18 atoms and 48 normal modes, making it more than three times larger than any of the previous test systems. We use the quartic PES of Cané *et al.*,^{51} which was calculated at the B97-1/TZ2P level of theory and includes 4125 nonzero anharmonic force constants. Large variational calculations were previously performed with this PES using the Hierarchical Intertwined Reduced-Rank Block Power Method (HI-RRBPM).^{44} We compare to the affordable HI-RRBPM (A) results and the most expensive HI-RRBPM (G) results, which we consider to be the most accurate. Finally, we compare to variational VASCI calculations from Ref. 38 with 1 × 10^{6} and 1.5 × 10^{6} basis states, obtained using 25 000 and 15 000 core states, respectively. Following Refs. 38 and 44, we calculate the 128 lowest states of naphthalene, of which the ten lowest and ten highest are presented in Table IV; results for all states can be found in the supplementary material. We show results with VHCI variational cutoff values *ε*_{1} = 1.5 cm^{−1} and 1 cm^{−1}, producing ∼1.3 × 10^{6} and 2.3 × 10^{6} basis states, respectively. Although a fully converged PT2 correction is intractable, we calculate an approximate PT2 correction with heat-bath screening; we used *ε*_{2} = 0.2 cm^{−1}, producing a perturbative space containing ∼15 × 10^{6} states.

Method . | ε_{1} = 1.5
. | ε_{1} = 1.0
. | VASCI^{38}
. | HI-RRBPM^{44}
. | ||||
---|---|---|---|---|---|---|---|---|

$NV$ . | 1 322 334 . | 2 270 672 . | 1 000 000 . | 1 500 000 . | A . | G . | ||

Assign. . | Var . | ε_{2} = 0.2
. | Var . | ε_{2} = 0.2
. | Var . | Var . | Var . | Var . |

ZPE | 31 772.71 | 31 764.77 | 31 769.90 | 31 764.34 | 31 774.4 | 31 773.6 | 31 782.20 | 31 766.03 |

ω_{48} | 168.17 | 165.20 | 166.89 | 164.63 | 166.4 | 166.3 | 165.84 | 164.60 |

ω_{13} | 182.57 | 179.28 | 181.28 | 178.79 | 179.9 | 179.9 | 184.90 | 178.18 |

2ω_{48} | 335.22 | 329.48 | 333.17 | 328.69 | 336.4 | 336.1 | 338.21 | 329.41 |

ω_{13} + ω_{48} | 349.25 | 342.84 | 346.48 | 341.76 | 349.5 | 349.4 | 365.68 | 342.02 |

ω_{24} | 358.59 | 356.26 | 357.83 | 355.97 | 361.4 | 363.8 | 372.86 | 354.44 |

2ω_{13} | 362.91 | 357.23 | 360.75 | 356.49 | 363.4 | 366.6 | 397.32 | 357.66 |

ω_{16} | 392.38 | 389.05 | 391.19 | 388.64 | 394.1 | 396.3 | 405.35 | 387.71 |

ω_{28} | 468.62 | 464.50 | 467.08 | 463.93 | 470.9 | 475.3 | 468.20 | 463.47 |

ω_{47} | 477.41 | 472.97 | 475.83 | 472.45 | 479.7 | 483.1 | 477.10 | 472.41 |

3ω_{48} | 503.57 | 493.98 | 499.45 | 492.38 | 502.5 | 501.8 | 506.60 | 495.50 |

ω_{13} + ω_{23} | 978.99 | 972.19 | 976.88 | 971.60 | 991.6 | 1 007.1 | 1 091.20 | 974.99 |

ω_{24} + ω_{36} | 982.41 | 977.49 | 981.37 | 977.64 | 1 006.6 | 1 015.4 | 1 099.87 | 982.87 |

ω_{12} + ω_{24} | 984.01 | 978.50 | 982.39 | 977.99 | 1 002.2 | 1 015.7 | 1 100.56 | 985.06 |

ω_{9} + ω_{28} | 984.89 | 978.67 | 982.83 | 978.03 | 997.7 | 1 010.7 | 1 098.74 | 981.31 |

ω_{44} + ω_{47} | 986.66 | 979.81 | 984.23 | 979.12 | 1 001.8 | 1 013.3 | 1 106.07 | 987.92 |

ω_{9} + ω_{47} | 993.93 | 987.08 | 991.49 | 986.44 | … | … | 1 121.57 | 994.66 |

ω_{35} | 1 013.18 | 1 008.63 | 1 011.79 | 1 008.21 | … | … | 1 134.08 | 1 011.99 |

ω_{16} + ω_{36} | 1 017.62 | 1 011.96 | 1 015.91 | 1 011.51 | … | … | 1 138.09 | 1 013.55 |

2ω_{44} | 1 017.51 | 1 012.86 | 1 016.19 | 1 012.70 | … | … | 1 138.46 | 1 015.09 |

ω_{9} + ω_{44} | 1 024.79 | 1 020.17 | 1 023.43 | 1 019.82 | … | … | 1 148.79 | 1 018.96 |

Max AE(25) | 11.78 | 2.85 | 6.32 | 3.69 | 16.1 | … | 54.64 | … |

RMSE(25) | 6.62 | 1.35 | 4.18 | 1.68 | 9.0 | … | 30.44 | … |

Core hours | 1 234.4 | 2 088.0 | 2 620.8 | 3 874.8 | 1 584.7 | 1 834.9 | 1 167.4 | 63 590.4 |

Cores | 20 | 20 | 20 | 20 | 40 | 40 | 128 | 64 |

Method . | ε_{1} = 1.5
. | ε_{1} = 1.0
. | VASCI^{38}
. | HI-RRBPM^{44}
. | ||||
---|---|---|---|---|---|---|---|---|

$NV$ . | 1 322 334 . | 2 270 672 . | 1 000 000 . | 1 500 000 . | A . | G . | ||

Assign. . | Var . | ε_{2} = 0.2
. | Var . | ε_{2} = 0.2
. | Var . | Var . | Var . | Var . |

ZPE | 31 772.71 | 31 764.77 | 31 769.90 | 31 764.34 | 31 774.4 | 31 773.6 | 31 782.20 | 31 766.03 |

ω_{48} | 168.17 | 165.20 | 166.89 | 164.63 | 166.4 | 166.3 | 165.84 | 164.60 |

ω_{13} | 182.57 | 179.28 | 181.28 | 178.79 | 179.9 | 179.9 | 184.90 | 178.18 |

2ω_{48} | 335.22 | 329.48 | 333.17 | 328.69 | 336.4 | 336.1 | 338.21 | 329.41 |

ω_{13} + ω_{48} | 349.25 | 342.84 | 346.48 | 341.76 | 349.5 | 349.4 | 365.68 | 342.02 |

ω_{24} | 358.59 | 356.26 | 357.83 | 355.97 | 361.4 | 363.8 | 372.86 | 354.44 |

2ω_{13} | 362.91 | 357.23 | 360.75 | 356.49 | 363.4 | 366.6 | 397.32 | 357.66 |

ω_{16} | 392.38 | 389.05 | 391.19 | 388.64 | 394.1 | 396.3 | 405.35 | 387.71 |

ω_{28} | 468.62 | 464.50 | 467.08 | 463.93 | 470.9 | 475.3 | 468.20 | 463.47 |

ω_{47} | 477.41 | 472.97 | 475.83 | 472.45 | 479.7 | 483.1 | 477.10 | 472.41 |

3ω_{48} | 503.57 | 493.98 | 499.45 | 492.38 | 502.5 | 501.8 | 506.60 | 495.50 |

ω_{13} + ω_{23} | 978.99 | 972.19 | 976.88 | 971.60 | 991.6 | 1 007.1 | 1 091.20 | 974.99 |

ω_{24} + ω_{36} | 982.41 | 977.49 | 981.37 | 977.64 | 1 006.6 | 1 015.4 | 1 099.87 | 982.87 |

ω_{12} + ω_{24} | 984.01 | 978.50 | 982.39 | 977.99 | 1 002.2 | 1 015.7 | 1 100.56 | 985.06 |

ω_{9} + ω_{28} | 984.89 | 978.67 | 982.83 | 978.03 | 997.7 | 1 010.7 | 1 098.74 | 981.31 |

ω_{44} + ω_{47} | 986.66 | 979.81 | 984.23 | 979.12 | 1 001.8 | 1 013.3 | 1 106.07 | 987.92 |

ω_{9} + ω_{47} | 993.93 | 987.08 | 991.49 | 986.44 | … | … | 1 121.57 | 994.66 |

ω_{35} | 1 013.18 | 1 008.63 | 1 011.79 | 1 008.21 | … | … | 1 134.08 | 1 011.99 |

ω_{16} + ω_{36} | 1 017.62 | 1 011.96 | 1 015.91 | 1 011.51 | … | … | 1 138.09 | 1 013.55 |

2ω_{44} | 1 017.51 | 1 012.86 | 1 016.19 | 1 012.70 | … | … | 1 138.46 | 1 015.09 |

ω_{9} + ω_{44} | 1 024.79 | 1 020.17 | 1 023.43 | 1 019.82 | … | … | 1 148.79 | 1 018.96 |

Max AE(25) | 11.78 | 2.85 | 6.32 | 3.69 | 16.1 | … | 54.64 | … |

RMSE(25) | 6.62 | 1.35 | 4.18 | 1.68 | 9.0 | … | 30.44 | … |

Core hours | 1 234.4 | 2 088.0 | 2 620.8 | 3 874.8 | 1 584.7 | 1 834.9 | 1 167.4 | 63 590.4 |

Cores | 20 | 20 | 20 | 20 | 40 | 40 | 128 | 64 |

For low-lying states, the agreement between variational VHCI and HI-RRBPM (G) is very good. At comparable computational cost, the accuracy of variational VASCI and VHCI is similar. The PT2 correction to VHCI produces a significant improvement of low-lying energies, which now match HI-RRBPM (G) to an accuracy of about 1 cm^{−1}.

We calculated the maximum absolute error and RMSE with respect to HI-RRBPM for the first 25 eigenstates. We matched the states to Refs. 44 and 38 according to their mode assignments. Variational VHCI achieves marginally closer agreement to the large HI-RRBPM(G) calculation than VASCI. Perturbation theory provides a noticeable improvement over the variational result. Curiously, VHCI+PT2 produces a more accurate result for the smaller variational space, indicating that the calculation is probably not converged with respect to *ε*_{2}. Indeed, ideally the PT2 correction should be calculated with *ε*_{2} ≪ *ε*_{1}, but even *ε*_{2} = 0.2 cm^{−1} produces an accurate result for low-lying states. As a convergence test, we also performed a large VHCI calculation that optimizes only the ground state (not shown) and obtained a variational ZPE that is at least 4 cm^{−1} lower than the HI-RRBPM (G) answer; therefore, some of our VHCI results that are lower than HI-RRBPM (G) may actually be more accurate, which would explain some of the discrepancies seen in the comparison.

For the high-lying states, which are much harder to converge as we discussed above, we only present results for those states that match the mode assignment from HI-RRBPM (G) and we do not attempt to calculate error statistics. For states with matching assignments, we see good agreement between variational VHCI and HI-RRBPM (G), especially at *ε*_{1} = 1 cm^{−1}. The PT2 correction is not as helpful for high-lying states as it is for low-lying ones because the variational space of the low-lying states is tightly converged and additional corrections are well-captured by perturbation theory. In some cases, the PT2 correction worsens the agreement with HI-RRBPM (G), although this may be a result of an incorrect mode assignment inside a dense spectrum of excited states. We do not attempt any extrapolation because *ε*_{2} = 0.2 cm^{−1} is not small enough to converge the PT2 correction; a stochastic PT2 implementation^{28} would be clearly beneficial. Overall, we are confident that the eigenvalue spectrum obtained from our VHCI calculations is accurate enough to be useful in real-world applications. The comparison of the overall CPU times, presented at the bottom of Table IV, underscores the competitiveness of VHCI as a means of accurately solving the vibrational Schrödinger equation on high-dimensional anharmonic potentials.

## IV. CONCLUSIONS

We have introduced vibrational heat-bath configuration interaction based on the original principles of HCI^{26} but with adaptations for vibrational Hamiltonians. VHCI+PT2 performed well on our four test molecules, achieving quantitative accuracy for the 12-dimensional systems acetonitrile and ethylene while outperforming VASCI and vDMRG. VHCI also performed well on larger systems, especially for low-lying states. High-lying states are a challenge due to their highly multi-configurational character and dense energetic spacing. Convergence of high-lying states could be improved by implementing a state-specific algorithm.^{33–35}

In future work, the implementation of semistochastic PT2^{27–29} will be critical for studying even larger systems. Extension of VHCI to alternative representations of the PES will be needed in order to eliminate error due to truncation and allow direct comparison to experiment; the efficiency of the algorithm in this context remains to be seen. Additionally, work is in progress to use VSCF single-mode basis functions within VHCI, which should greatly reduce the size of the configuration space needed to converge the variational stage. Furthermore, VHCI can be straightforwardly extended to efficiently calculate spectroscopic intensities based on the dipole moment surface.^{35,52,53} Consideration of spectroscopic activity can also be used to target eigenstates more efficiently, as recently implemented in A-VCI.^{54} Finally, it will be interesting to apply VHCI to more strongly anharmonic systems, such as molecular clusters^{55,56} or floppy molecules,^{6,57,58} where truncated expansions of the PES will not be sufficient.

## SUPPLEMENTARY MATERIAL

See the supplementary material for VHCI energies and assignments of all states calculated for the four molecules considered in this work.

## ACKNOWLEDGMENTS

We thank Dr. James Brown and Dr. Marc Odunlami for providing us with accurate expanded potential energy surfaces and reference energies for acetonitrile and ethylene oxide and Professor Sandeep Sharma for helpful discussions. J.H.F. was supported, in part, by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1644869. We acknowledge computing resources from Columbia University’s Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR), Contract No. C090171, both awarded on April 15, 2010. The Flatiron Institute is a division of the Simons Foundation.

## DATA AVAILABILITY

The software and data that support the findings of this study are openly available in the supplementary material and in our GitHub repository at https://github.com/berkelbach-group/VHCI with Zenodo digital object identifier http://doi.org/10.5281/zenodo.4116070.

## REFERENCES

_{2}CN

_{2}

_{2}O)

_{6}

_{7}O

_{3}

^{+}and H

_{9}O

_{4}

^{+}(eigen) and comparison with experiment

^{+}

_{3}up to 35000 cm

^{−1}