Vibrational heat-bath configuration interaction

Expressed in terms of mass-weighted normal mode coordinates $Q_{i}$ with frequencies $\omega_{i}$, the nuclear Hamiltonian is given by

$$H=\frac{1}{2}\sum_{i=1}^{D}\left[-\frac{\partial^{2}}{\partial Q_{i}^{2}}+\omega_{i}^{2}Q_{i}^{2}\right]+V_{\mathrm{an}}(Q_{1},Q_{2},\ldots,Q_{D}),$$

(1)

where $D=3N-6$ is the number of normal mode degrees of freedom and $V_{\mathrm{an}}$ is the anharmonic part of the electronic ground state potential energy surface (PES). We assume that the anharmonic part of the PES has a normal mode expansion

$$\begin{split}V_{\mathrm{an}}(\bm{Q})&=\sum_{i\leq j\leq k}V_{ijk}Q_{i}Q_{j}Q_{k}+\sum_{i\leq j\leq k\leq l}V_{ijkl}Q_{i}Q_{j}Q_{k}Q_{l}+\ldots\\ &=\sum_{i\leq j\leq k}W_{ijk}\bar{Q}_{i}\bar{Q}_{j}\bar{Q}_{k}+\sum_{i\leq j\leq k\leq l}W_{ijkl}\bar{Q}_{i}\bar{Q}_{j}\bar{Q}_{k}\bar{Q}_{l}+\ldots\\ \end{split}$$

(2)

where the anharmonic force constants are partial derivatives of the PES, e.g. $V_{ijk}=(\partial^{3}V_{\mathrm{an}}/\partial Q_{i}\partial Q_{j}\partial Q_{k})_{\bm{Q}=0}$. In the second line of Eq. (2), we define $\bar{Q}_{i}=(a_{i}^{\dagger}+a_{i})$ and

$$W_{ijk\ldots}=\frac{V_{ijk\ldots}}{\sqrt{2^{p}\omega_{i}\omega_{j}\omega_{k}\ldots}},$$

(3)

where $p$ is the order of the anharmonicity. In practice the anharmonic PES expansion is truncated, commonly at fourth or sixth order. In this work, 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 $|\bm{n}\rangle=|n_{1},n_{2},\ldots,n_{D}\rangle$ formed from the harmonic oscillator orbitals $\phi_{n_{i}}(Q_{i})=\langle Q_{i}|n_{i}\rangle$, leading to the wavefunction

$$\ket{\Psi}=\sum_{{\bm{n}}\in{\mathcal{V}}}c_{\bm{n}}\ket{{\bm{n}}}.$$

(4)

In VCI, the variational space ${\mathcal{V}}$ is commonly truncated by 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 ${\mathcal{V}}$ by criteria other than 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 $\ket{{\bm{m}}}$ is added to ${\mathcal{V}}$ if $|H_{{\bm{m}}{\bm{n}}}c_{\bm{n}}|\geq\varepsilon_{1}$, where $\varepsilon_{1}$ is a user-defined convergence parameter. In electronic structure, the matrix elements $H_{{\bm{m}}{\bm{n}}}$ 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 because the modes can have occupation numbers greater than one. For example, consider an anharmonic PES with cubic and quartic terms. The Hamiltonian matrix elements between product states $\ket{{\bm{m}}}$ and $\ket{{\bm{n}}}$ that differ in their occupancy by one quantum in mode $i$, by two quanta in modes $i$ and $j$ (including $i=j$), and so on, are given by

	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i)}$	$\displaystyle=\sum_{j}W_{ijj}\langle{\bm{m}}|\bar{Q}_{i}\bar{Q}_{j}^{2}|{\bm{n}}\rangle$		(5a)
	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i,j)}$	$\displaystyle=\sum_{k}W_{ijkk}\langle{\bm{m}}|\bar{Q}_{i}\bar{Q}_{j}\bar{Q}_{k}^{2}|{\bm{n}}\rangle$		(5b)
	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i,j,k)}$	$\displaystyle=W_{ijk}\langle{\bm{m}}|\bar{Q}_{i}\bar{Q}_{j}\bar{Q}_{k}|{\bm{n}}\rangle$		(5c)
	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i,j,k,l)}$	$\displaystyle=W_{ijkl}\langle{\bm{m}}|\bar{Q}_{i}\bar{Q}_{j}\bar{Q}_{k}\bar{Q}_{l}|{\bm{n}}\rangle.$		(5d)

The factors $\langle{\bm{m}}|\bar{Q}_{i}\ldots|{\bm{n}}\rangle$ that are not zero are products of terms of the form $\sqrt{n_{i}}$ 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 $\ket{{\bm{m}}}$ and $\ket{{\bm{n}}}$ differ in the occupancy of three modes $(i,j,k)$, then

	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i+,j+,k+)}$	$\displaystyle=W_{ijk}\sqrt{m_{i}m_{j}m_{k}}$		(6a)
	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i+,j+,k-)}$	$\displaystyle=W_{ijk}\sqrt{m_{i}m_{j}n_{k}}$		(6b)
	$\displaystyle H_{{\bm{m}}{\bm{n}}}^{(i+,j-,k-)}$	$\displaystyle=W_{ijk}\sqrt{m_{i}n_{j}n_{k}}$		(6c)

and so on, where the notation indicates that the occupation of a given mode in $\ket{{\bm{m}}}$ is bigger or smaller than in $\ket{{\bm{n}}}$. Therefore, in addition to the identity of the modes with different occupations, the VCI matrix elements 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 $\langle{\bm{m}}|\bar{Q}_{i}\ldots|{\bm{n}}\rangle$ are products of terms of the form $\sqrt{n_{i}}$, 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}=\sum_{j}(2+\delta_{ij})W_{ijj}$ and $W_{ij}=\sum_{k}(2+\delta_{jk}+\delta_{ij}\delta_{jk})W_{ijkk}$, which approximately account for single and double mode excitations. Then, given an approximate wavefunction of the form (4) expanded over some space ${\mathcal{V}}$, we propose to add state $\ket{{\bm{m}}}$ to the variational space if $\ket{{\bm{m}}}$ and $\ket{{\bm{n}}}$ differ in their occupancy by one quantum in mode $i$, by two quanta in modes $i$ and $j$ (including $i=j$), and so on, and

$$\left|W_{ijk\ldots}c_{\bm{n}}\right|\geq\varepsilon_{1}.$$

(7)

The above criterion can be implemented very efficiently by pre-sorting the matrix elements $W_{ijk\ldots}$, 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:

1. 1.

   Sort the list of $W_{ijk\ldots}$, largest to smallest.

2. 2.

   Initialize the space ${\mathcal{V}}$ according to a small total number of quanta.

3. 3.

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

4. 4.

   Add states to ${\mathcal{V}}$ via the criterion (7), as follows. For each state $\ket{{\bm{n}}}$, traverse the sorted list of $W_{ijk\ldots}$.

   1. (a)

      If $|W_{i}c_{\bm{n}}|\geq{\varepsilon}_{1}$, add all states $\ket{{\bm{m}}}$ that are included in $Q_{i}\ket{{\bm{n}}}$.

   2. (b)

      If $|W_{ij}c_{\bm{n}}|\geq{\varepsilon}_{1}$, add all states $\ket{{\bm{m}}}$ that are included in $Q_{i}Q_{j}\ket{{\bm{n}}}$.

   3. (c)

      If $|W_{ijk\ldots}c_{\bm{n}}|\geq{\varepsilon}_{1}$, add all states $\ket{{\bm{m}}}$ that are included in $Q_{i}Q_{j}Q_{k}\ldots\ket{{\bm{n}}}$. At anharmonic order $p$, there are $O(2^{p})$ such states to add.

   4. (d)

      If $|W_{ijk\ldots}c_{\bm{n}}|<{\varepsilon}_{1}$, break and go on to the next $\ket{{\bm{n}}}$.

   Return to step 3 until the calculation is converged.

A number of possible convergence criteria can be devised; here, following the original HCI paper [24], 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_{{\bm{m}}{\bm{n}}}$ are zero due to the properties of harmonic oscillators; the states $\ket{{\bm{m}}}$ and $\ket{{\bm{n}}}$ 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 [1, 2, 3, 4], which is why we intentionally work in the noninteracting normal mode basis.

Crucially, the presorting of scaled and/or summed anharmonic force constants $(W_{i},W_{ij},W_{ijk\ldots})$ combined with the criterion (7) means that most states $\ket{{\bm{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. 25, 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_{\mathrm{s}}$ lowest in energy). Then we perform the addition step 4 for each of those states, combining all of their added basis states $\ket{{\bm{m}}}$ in the updated ${\mathcal{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 $N_{\mathrm{s}}$.

The above procedure can be applied to a PES expanded to arbitrary order in the normal mode coordinates. 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.4; although we could define new effective force constants for the screening procedure, e.g. $W_{ijk}^{\prime}=W_{ijk}+\sum_{k}W_{ijkll}$, instead we make the approximation $\sum_{l}W_{ijkll}\approx\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.

To the variational energy of each state $E_{\mathrm{var}}$ we add the second-order perturbation theory (PT2) correction

$$\Delta E_{2}\approx\sum_{\bm{m}}^{({\varepsilon}_{2})}\frac{\left(\sum_{\bm{n}}H_{{\bm{m}}{\bm{n}}}c_{\bm{n}}\right)^{2}}{E_{\mathrm{var}}-H_{{\bm{m}}{\bm{m}}}}.$$

(8)

In exact PT2, the summation over ${\bm{m}}$ includes all basis states that are connected to the variational space ${\mathcal{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 [24]. For VHCI, we again use criterion (7), with a cutoff $\varepsilon_{2}<\varepsilon_{1}$, to determine whether basis state $\ket{{\bm{m}}}$ should be included in the perturbative space. When $\varepsilon_{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 [26], in order to study the vibrational structure of even larger systems.

All simulations besides those on naphthalene were performed on a 4-core (8-thread) Intel Core i7-6700 3.4 GHz desktop 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 [38] with extensive modifications. Our VHCI code is available on GitHub [39]. The Hamiltonian matrix is stored in a sparse format using the Eigen linear algebra library [40]. The eigenvalues and eigenvectors of interest are calculated with the Lanczos algorithm as implemented in the SPECTRA linear algebra library [41]. We verified our code by comparing to the literature and to results obtained with the PyVCI software package [34].

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.

We first present results for acetonitrile, CH₃CN, a 6-atom, 12-dimensional system that has become one of the canonical benchmarks for new methods in vibrational structure theory [15, 42, 17, 43, 36]. We use the quartic PES described in Refs. 44, 15, 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 [14, 37]. All PESs used in this study are available as input files on our GitHub repository [39].

We calculated the energies of the first 70 states of acetonitrile, of which the first 20 are reported in Table 1. Data for all 70 states can be found in the supplementary material. We compare our results to those obtained using A-VCI [37], which we consider to be numerically exact (similar to those obtained with the nonproduct quadrature approach [14]); for reference, the A-VCI results are obtained with approximately 2.5 million 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 Ref. 37. 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 [36], which is a selected CI technique that is similar in spirit to VHCI.

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 ${\varepsilon}_{1}$ that controls the number of variational states $N_{\mathcal{V}}$. Using ${\varepsilon}_{1}=1$ cm^-1 results in a variational space with $N_{\mathcal{V}}=29\,900$ basis states, which enables comparison to the largest reported VASCI calculation, with $N_{\mathcal{V}}=30\,038$. 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 ($\varepsilon_{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 ${\varepsilon}_{1}=1$ cm^-1 calculation reported here took less than 3 minutes for the variational stage and less than 30 minutes for the full PT2 correction (on an 8-core desktop CPU).

In Table 1, we also present results of a larger calculation with ${\varepsilon}_{1}=0.1$ cm^-1, resulting in $N_{\mathcal{V}}=125~{}038$. 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 minutes on the same 8-core desktop CPU. This larger calculation still benefits from minor corrections with PT2, which takes about three hours 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.

Next, we study ethylene, C₂H₄, 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. 45, which was calculated entirely using CCSD(T)/cc-pVQZ, and we use the PyPES software package [46] 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 sixteen states are shown in Table 2 and those for the additional 84 states can be found in the supplementary material.

Through testing, we confirmed that our VHCI+PT2 results are converged with ${\varepsilon}_{1}=0.75$ cm^-1 for the quartic case (with full PT2) and ${\varepsilon}_{1}=0.5$ cm^-1 for the sixth-order case (with ${\varepsilon}_{2}=0.01$ cm^-1), resulting in $N_{\mathcal{V}}=$ 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 [45, 34, 47].

In Table 2, we also compare our results to those obtained with the vibrational density-matrix renormalization group (vDMRG) [17] using the same fourth- and sixth-order truncated PES. 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, $N_{\mathcal{V}}$. We see that VHCI (obtained with $N_{\mathrm{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 simply due to slight differences in the PES. As a check on the latter, we have also included in Table 2 the VCI excitation energies previously reported by the authors of the PyPES software package [34] 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 function of 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 2, we also show results obtained in Ref. 45 using a pruned VCI basis with up to 13 vibrational quanta for the untruncated PES. In its current form, VHCI requires the truncated normal mode expansion (2) and thus comparison to the results obtained for an untruncated PES is important for assessing the future potential of VHCI. 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-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.

Next, we study ethylene oxide, C₂H₄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.[48], 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. 16, 42, 43, 36, which is available along with our source code on our GitHub repository [39]. Like acetonitrile, ethylene oxide represents a well-studied system that is suitable for benchmarking new approaches to solving the vibrational structure problem. In Table 3, we present results for the ten lowest and ten highest states from a VHCI calculation targeting $N_{\mathrm{s}}=200$ states. We performed the calculations with two different variational cutoff values ${\varepsilon}_{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 ${\varepsilon}_{2}=0.01$ cm^-1.

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 ${\varepsilon}_{1}$, for a variety of values of the perturbative screening parameter ${\varepsilon}_{2}$. We see that the curves with ${\varepsilon}_{2}=1$ cm^-1, 0.1 cm^-1, and 0.01 cm^-1 agree very closely for all ${\varepsilon}_{1}$, with ${\varepsilon}_{2}=0.1$ cm^-1 and 0.01 cm^-1 being indistinguishable, confirming convergence with respect to ${\varepsilon}_{2}$. In Table 3 and Fig. 2, we compare our results to A-VCI calculations [16] that were obtained from an optimized variational space containing over 7 million 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 ${\varepsilon}_{1}$ becomes small. Addition of the PT2 correction leads to more rapid convergence with respect to ${\varepsilon}_{1}$. For example, for the ground state we see that results obtained with ${\varepsilon}_{1}=20$ cm^-1 and converged PT2 gives an answer that is closer to exact than that with ${\varepsilon}_{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 ${\varepsilon}_{1}\leq 1$ cm^-1. Table 3 shows the band assignment (obtained from the variational calculation) following the mode labeling convention of Ref. 16 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 ${\varepsilon}_{1}=1$ cm^-1 and the A-VCI reference; several of the other ten highest excitations are also omitted from the VASCI results [36]. Because the incorrect band assignments of high-lying states prevents 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-2 cm^-1 and maximum AE below 5 cm^-1; the addition of PT2 yields accuracy better than 1 cm^-1 for both ${\varepsilon}_{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 $N_{\mathcal{V}}=150\,000$, but the closely-spaced states 47 and 48 have the wrong energetic ordering for VASCI and for VHCI with ${\varepsilon}_{1}=2$ cm^-1. For ${\varepsilon}_{1}=1$ cm^-1, we obtain all 200 states with the correct assignments, enabling directe 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 calculation took less than 8 hours on an 8-core desktop CPU. As a rough comparison to other methods, we also present the timings for all calculations in Table 3.

Finally, following standard HCI practice in 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 $\Delta E_{2}$ seeking to extrapolate to $\Delta E_{2}=0$, following Ref. 25. 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 ${\varepsilon}_{1}$ changes is challenging. We performed extrapolation for the first 50 states using a linear fit to results obtained with ${\varepsilon}_{1}=1$ cm^-1 and 2 cm^-1, as shown graphically in 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 perfoming, although other schemes could be considered in the future. In Table 2, 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.

As a final test of VHCI, we consider naphthalene, C₁₀H₈, 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. [49], 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) [43]. 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. 36 with 1 million and 1.5 million basis states, obtained using 25 000 and 15 000 core states, respectively. Following Refs. 43, 36, we calculate the 128 lowest states of naphthalene, of which the ten lowest and ten highest are presented in Table 4; results for all states can be found in the supplementary material. We show results with VHCI variational cutoff values ${\varepsilon}_{1}=1.5$ cm^-1 and 1 cm^-1, producing approximately 1.3 million and 2.3 million basis states, respectively. Although a fully converged PT2 correction is intractable, we calculate an approximate PT2 correction with heat-bath screening; we used ${\varepsilon}_{2}=0.2$ cm^-1, producing a perturbative space containing approximately 15 million states.

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 Ref. 43 and Ref. 36 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 ${\varepsilon}_{2}$. Indeed, ideally the PT2 correction should be calculated with ${\varepsilon}_{2}\ll{\varepsilon}_{1}$, but even ${\varepsilon}_{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 ${\varepsilon}_{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 ${\varepsilon}_{2}=0.2$ cm^-1 is not small enough to converge the PT2 correction; a stochastic PT2 implementation [26] 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. Comparison of the overall CPU times, presented at the bottom of Table 4, underscores the competitiveness of VHCI as a means of accurately solving the vibrational Schrödinger equation on high-dimensional anharmonic potentials.