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.

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 CIPSI37 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.

Expressed in terms of mass-weighted normal mode coordinates Qi with frequencies ωi, the nuclear Hamiltonian is given by

H=12i=1D2Qi2+ωi2Qi2+Van(Q1,Q2,,QD),
(1)

where D = 3N − 6 is the number of normal mode degrees of freedom and Van 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⟩ = |n1, n2, …, nD⟩ formed from the harmonic oscillator orbitals ϕni(Qi)=Qi|ni, leading to the wavefunction

|Ψ=nVcn|n.
(2)

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 |Hmncn| ≥ ε1, where ε1 is a user-defined convergence parameter. In electroinic structure theory, the matrix elements Hmn 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 Hmn matrix elements are expensive. In order to devise an analogous VHCI screening criterion, we consider an approximate Hamiltonian H̃ obtained from a truncated normal mode expansion of the anharmonic part of the PES,

Van(Q)=ijkVijkQiQjQk+ijklVijklQiQjQkQl+=ijkWijkQ¯iQ¯jQ¯k+ijklWijklQ¯iQ¯jQ¯kQ¯l+,
(3)

where the anharmonic force constants are partial derivatives of the PES, e.g., Vijk=(3Van/QiQjQk)Q=0. In the second line of Eq. (3), we define Q¯i=(ai+ai) and

Wijk=Vijk2pωiωjωk,
(4)

where p is the order of the anharmonicity. We emphasize that while we limit ourselves to using H̃ 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 |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

H̃mn(i)=jWijjm|Q¯iQ¯j2|n,
(5a)
H̃mn(i,j)=kWijkkm|Q¯iQ¯jQ¯k2|n,
(5b)
H̃mn(i,j,k)=Wijkm|Q¯iQ¯jQ¯k|n,
(5c)
H̃mn(i,j,k,l)=Wijklm|Q¯iQ¯jQ¯kQ¯l|n.
(5d)

The factors m|Q¯i|n 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 |n⟩ differ in the occupancy of three modes (i, j, k), then

H̃mn(i+,j+,k+)=Wijkmimjmk,
(6a)
H̃mn(i+,j+,k)=Wijkmimjnk,
(6b)
H̃mn(i+,j,k)=Wijkminjnk,
(6c)

and so on, where the notation indicates that the occupation of a given mode in |m⟩ is bigger or smaller than in |n⟩. Therefore, in addition to the 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 m|Q¯i|n 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 Wi = ∑j(2 + δij)Wijj and Wij = ∑k(2 + δjk + δijδjk)Wijkk, 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 |m⟩ and |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

Wijkcnε1.
(7)

The above criterion can be implemented very efficiently by pre-sorting the matrix elements Wijk…, of which there are at most a polynomial number, e.g., O(D4) for a quartic PES.

The ground-state VHCI algorithm is performed as follows:

  1. Sort the list of Wijk, largest to smallest.

  2. Initialize the space V according to a small total number of quanta.

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

  4. Add states to V via criterion (7) as follows. For each state |n⟩, traverse the sorted list of Wijk.

    • If |Wicn| ≥ ε1, add all states |m⟩ that are included in Qi|n⟩.

    • If |Wijcn| ≥ ε1, add all states |m⟩ that are included in QiQj|n⟩.

    • If |Wijk, …, cn| ≥ ε1, add all states |m⟩ that are included in QiQjQk, …, |n⟩. At anharmonic order p, there are O(2p) such states to add.

    • If |Wijk, …, cn| < ε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̃mn are zero due to the properties of harmonic oscillators; the states |m⟩ and |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–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 (Wi, Wij, Wijk, …, ) combined with criterion (7) means that most states |mare 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 Ns lowest in energy). Then, we perform the addition step 4 for each of those states, combining all of their added basis states |m⟩ 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 Ns.

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 = Wijk + ∑kWijkll, instead we make the approximation ∑lWijkll ≈ maxlWijkll. 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 Evar, we add the second-order perturbation theory (PT2) correction

ΔE2m(ε2)nHmncn2EvarHmm.
(8)

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 |m⟩ should be included in the perturbative space. When ε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.

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 package39 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̃ 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.

We first present results for acetonitrile, CH3CN, 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 approach16); for reference, the A-VCI results are obtained with ∼2.5 × 106 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.

TABLE I.

The 20 lowest-energy states from a calculation of the first 70 eigenstates of acetonitrile, with excitation energies given relative to the ZPE. Mode assignments are given based on the character of the basis state with the largest absolute CI coefficient using mode-numbering from Refs. 16 and 38. We also include the total vibrational angular momentum superscript and symmetry label for each mode. VHCI results are shown for two different values of ε1, which determines the number of variational basis states, NV, and are shown with and without the PT2 correction to the energies (“Var” and “Full PT2,” respectively). We compare to exact reference values from A-VCI18 and to VASCI with full PT2 correction and NV comparable to our ε1 = 1.0 cm−1 case.38 The maximum absolute error (Max. AE) and root mean squared error (RMSE) relative to A-VCI were calculated across all 70 states. All energies are in cm−1.

Methodε1 = 1.0ε1 = 0.1VASCI38 A-VCI18 
NV29 900125 03830 0382 488 511
Assign.Sym.VarFull PT2VarFull PT2Full PT2Var
ZPE  9 837.43 9 837.41 9 837.41 9 837.41 9 837.41 9 837.41 
ω11±1 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ω11±2 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ω110 A1 723.90 723.84 723.83 723.83 723.89 723.83 
ω4±1 A1 900.70 900.66 900.66 900.66 900.68 900.66 
ω9±1 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ω11±3 A1 1 086.65 1 086.56 1 086.56 1 086.55 1 086.64 1 086.55 
3ω11±3 A2 1 086.66 1 086.56 1 086.56 1 086.55 1 086.64 1 086.55 
3ω11±1 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 
ω4+ω11±1 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 A1 1 389.10 1 388.99 1 388.98 1 388.97 1 388.99 1 388.97 
ω9±1+ω11±1 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 
ω9±1+ω111 A2 1 395.08 1 394.93 1 394.91 1 394.90 1 395.00 1 394.90 
ω9±1+ω111 A1 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.1VASCI38 A-VCI18 
NV29 900125 03830 0382 488 511
Assign.Sym.VarFull PT2VarFull PT2Full PT2Var
ZPE  9 837.43 9 837.41 9 837.41 9 837.41 9 837.41 9 837.41 
ω11±1 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ω11±2 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ω110 A1 723.90 723.84 723.83 723.83 723.89 723.83 
ω4±1 A1 900.70 900.66 900.66 900.66 900.68 900.66 
ω9±1 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ω11±3 A1 1 086.65 1 086.56 1 086.56 1 086.55 1 086.64 1 086.55 
3ω11±3 A2 1 086.66 1 086.56 1 086.56 1 086.55 1 086.64 1 086.55 
3ω11±1 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 
ω4+ω11±1 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 A1 1 389.10 1 388.99 1 388.98 1 388.97 1 388.99 1 388.97 
ω9±1+ω11±1 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 
ω9±1+ω111 A2 1 395.08 1 394.93 1 394.91 1 394.90 1 395.00 1 394.90 
ω9±1+ω111 A1 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=29900 basis states, which enables comparison to the largest reported VASCI calculation, with NV=30038. 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=125038. 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.

Next, we study ethylene, C2H4, 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 package48 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.

TABLE II.

The 16 lowest-energy states from a calculation of the first 100 eigenstates of ethylene, with excitation energies given relative to the ZPE. Calculations were performed using fourth- and sixth-order truncations of the PES normal mode expansion. Both VHCI+PT2 calculations were converged with respect to variational and perturbative basis size. We compare to vDMRG19 for both truncations and to VCI36 with up to eight quanta per product state for the sixth-order truncation. We also compare to a calculation from Ref. 47 that used an untruncated version of the PES solved using VCI with a pruned basis set containing up to 13 quanta per product state. Assignments are given using the mode-numbering convention of Ref. 47. All energies are in cm−1.

PESFourth-orderSixth-orderUntruncated
MethodVHCI+Full PT2VHCI+PT2 (ε2 = 0.01)VCI(PyVCI)48 Pruned VCI47 
Assign.ε1 = 0.75vDMRG19 ε1 = 0.5vDMRG19 Ntot = 8Ntot = 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 
PESFourth-orderSixth-orderUntruncated
MethodVHCI+Full PT2VHCI+PT2 (ε2 = 0.01)VCI(PyVCI)48 Pruned VCI47 
Assign.ε1 = 0.75vDMRG19 ε1 = 0.5vDMRG19 Ntot = 8Ntot = 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 Ns = 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 package36 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.

FIG. 1.

The variational ground-state energy of ethylene as a function of the number of basis states NV. We compare the ground state energies of conventional VCI (red solid circles) and VHCI (blue solid squares) using a fourth-order (top) and sixth-order (bottom) truncation of the potential. VHCI calculations were performed by optimizing for the ground state (Ns = 1), and the VCI space was truncated by limiting the total number of quanta per product state. The black dashed lines represent the converged VHCI ground-state energy, which we consider to be exact. The orange dashed lines are the vDMRG values from Ref. 19. We also present the sparsity of the Hamiltonian matrix (right axis) as a function of NV for conventional VCI (red open circles) and VHCI (blue open squares).

FIG. 1.

The variational ground-state energy of ethylene as a function of the number of basis states NV. We compare the ground state energies of conventional VCI (red solid circles) and VHCI (blue solid squares) using a fourth-order (top) and sixth-order (bottom) truncation of the potential. VHCI calculations were performed by optimizing for the ground state (Ns = 1), and the VCI space was truncated by limiting the total number of quanta per product state. The black dashed lines represent the converged VHCI ground-state energy, which we consider to be exact. The orange dashed lines are the vDMRG values from Ref. 19. We also present the sparsity of the Hamiltonian matrix (right axis) as a function of NV for conventional VCI (red open circles) and VHCI (blue open squares).

Close modal

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.

Next, we study ethylene oxide, C2H4O, 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 Ns = 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.

TABLE III.

Ten lowest- and ten highest-energy states from the first 200 eigenstates of ethylene oxide, with excitation energies given relative to the ZPE. VHCI results are shown for two different values of ε1 with no PT2 correction (Var) and with PT2 at ε2 = 0.01 cm−1. We also show the first ten extrapolated energies obtained by a linear fit of Etot vs ΔE2, as shown in Fig. 2. We compare all results to numerically exact reference values from A-VCI18 and to VASCI with full PT2 correction and NV roughly comparable to our ε1 = 2 cm−1 case.38 We report the maximum absolute error (Max. AE) and root mean squared error (RMSE) of the first 50 states for all methods. We show the same error metrics for all 200 states for VHCI with ε1 = 1 cm−1, which is the only case for which all mode assignments match the exact result. Assignments use the mode-numbering convention of Ref. 18. All energies are in cm−1.

Methodε1 = 2.0ε1 = 1.0Linear ΔE2VASCI38 A-VCI18 
NV132 163259 070Extrap.150 0007 118 214
Assign.Varε2 = 0.01Varε2 = 0.01ε2 = 0.01Full PT2Var
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 … 24 
Methodε1 = 2.0ε1 = 1.0Linear ΔE2VASCI38 A-VCI18 
NV132 163259 070Extrap.150 0007 118 214
Assign.Varε2 = 0.01Varε2 = 0.01ε2 = 0.01Full PT2Var
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 … 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 calculations18 that were obtained from an optimized variational space containing over 7 × 106 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=150000, 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.

FIG. 2.

VHCI energy of the ground state (top row), the 50th state (middle row), and the 200th state (bottom row) of ethylene oxide. The left-hand column shows the energy as a function of the variational heat-bath cutoff parameter ε1 and includes the variational energy Evar (blue squares) and the total energy Etot = Evar + ΔE2 for various values of the perturbative heat-bath screening parameter ε2. The black dashed lines are the converged A-VCI values for each state, which can be considered to be exact. In the right-hand column, we present the total energy Etot (red circles) with converged PT2 (ε2 = 0.01 cm−1) as a function of the perturbative correction ΔE2. Note that we use a smaller y-range because the energy varies less with respect to ε1 when the PT2 stage is converged. We also show the linear fits (blue lines) of the most accurate two points (ε1 = 1 cm−1 and 2 cm−1) for the ground state and state 50, which were used to find the extrapolated energies presented in Table III.

FIG. 2.

VHCI energy of the ground state (top row), the 50th state (middle row), and the 200th state (bottom row) of ethylene oxide. The left-hand column shows the energy as a function of the variational heat-bath cutoff parameter ε1 and includes the variational energy Evar (blue squares) and the total energy Etot = Evar + ΔE2 for various values of the perturbative heat-bath screening parameter ε2. The black dashed lines are the converged A-VCI values for each state, which can be considered to be exact. In the right-hand column, we present the total energy Etot (red circles) with converged PT2 (ε2 = 0.01 cm−1) as a function of the perturbative correction ΔE2. Note that we use a smaller y-range because the energy varies less with respect to ε1 when the PT2 stage is converged. We also show the linear fits (blue lines) of the most accurate two points (ε1 = 1 cm−1 and 2 cm−1) for the ground state and state 50, which were used to find the extrapolated energies presented in Table III.

Close modal

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 ΔE2 seeking to extrapolate to ΔE2 = 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.

As a final test of VHCI, we consider naphthalene, C10H8, 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 × 106 and 1.5 × 106 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 × 106 and 2.3 × 106 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 × 106 states.

TABLE IV.

Ten lowest- and ten highest-energy states from a calculation of the first 128 eigenstates of naphthalene, with excitation energies given relative to the ZPE. VHCI results are shown for two different values of ε1 with no PT2 correction (Var) and with PT2 at ε2 = 0.2 cm−1. We compare to reference values from the smallest and largest HI-RRBPM calculations in Ref. 44 (“A” and “G,” respectively) as well as to VASCI38 with no PT2 correction and NV similar to our ε1 = 1.5 cm−1 case. Assignments are given using the mode-numbering convention of Ref. 51. All energies are in cm−1.

Methodε1 = 1.5ε1 = 1.0VASCI38 HI-RRBPM44 
NV1 322 3342 270 6721 000 0001 500 000AG
Assign.Varε2 = 0.2Varε2 = 0.2VarVarVarVar
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.0VASCI38 HI-RRBPM44 
NV1 322 3342 270 6721 000 0001 500 000AG
Assign.Varε2 = 0.2Varε2 = 0.2VarVarVarVar
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 implementation28 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.

We have introduced vibrational heat-bath configuration interaction based on the original principles of HCI26 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 PT227–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 clusters55,56 or floppy molecules,6,57,58 where truncated expansions of the PES will not be sufficient.

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

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.

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.

1.
R. B.
Gerber
,
B.
Brauer
,
S. K.
Gregurick
, and
G. M.
Chaban
, “
Calculation of anharmonic vibrational spectroscopy of small biological molecules
,”
PhysChemComm
5
,
142
(
2002
).
2.
G. M.
Chaban
,
J. O.
Jung
, and
R. B.
Gerber
, “
Ab initio calculation of anharmonic vibrational states of polyatomic systems: Electronic structure combined with vibrational self-consistent field
,”
J. Chem. Phys.
111
,
1823
1829
(
1999
).
3.
S.
Carter
,
J. M.
Bowman
, and
N. C.
Handy
, “
Extensions and tests of ‘multimode’: A code to obtain accurate vibration/rotation energies of many-mode molecules
,”
Theor. Chem. Acc.
100
,
191
198
(
1998
).
4.
S.
Carter
,
S. J.
Culik
, and
J. M.
Bowman
, “
Vibrational self-consistent field method for many-mode systems: A new approach and application to the vibrations of CO adsorbed on Cu(100)
,”
J. Chem. Phys.
107
,
10458
10469
(
1997
).
5.
J. M.
Bowman
,
S.
Carter
, and
X.
Huang
, “
Multimode: A code to calculate rovibrational energies of polyatomic molecules
,”
Int. Rev. Phys. Chem.
22
,
533
549
(
2003
).
6.
T. K.
Roy
and
R. B.
Gerber
, “
Vibrational self-consistent field calculations for spectroscopy of biological molecules: New algorithmic developments and applications
,”
Phys. Chem. Chem. Phys.
15
,
9468
(
2013
).
7.
V.
Barone
, “
Anharmonic vibrational properties by a fully automated second-order perturbative approach
,”
J. Chem. Phys.
122
,
014108
(
2005
).
8.
O.
Christiansen
, “
Møller–Plesset perturbation theory for vibrational wave functions
,”
J. Chem. Phys.
119
,
5773
5781
(
2003
).
9.
E. L.
Sibert
, “
Theoretical studies of vibrationally excited polyatomic molecules using canonical Van Vleck perturbation theory
,”
J. Chem. Phys.
88
,
4378
4390
(
1987
).
10.
S.
Banik
,
S.
Pal
, and
M. D.
Prasad
, “
Calculation of vibrational energy of molecule using coupled cluster linear response theory in bosonic representation: Convergence studies
,”
J. Chem. Phys.
129
,
134111
(
2008
).
11.
P.
Seidler
and
O.
Christiansen
, “
Automatic derivation and evaluation of vibrational coupled cluster theory equations
,”
J. Chem. Phys.
131
,
234109
(
2009
).
12.
O.
Christiansen
, “
Vibrational coupled cluster theory
,”
J. Chem. Phys.
120
,
2149
2159
(
2004
).
13.
J. M.
Bowman
,
K.
Christoffel
, and
F.
Tobin
, “
Application of SCF-SI theory to vibrational motion in polyatomic molecules
,”
J. Phys. Chem.
83
,
905
920
(
1979
).
14.
T. C.
Thompson
and
D. G.
Truhlar
, “
SCF CI calculations for vibrational eigenvalues and wavefunctions of systems exhibiting fermi resonance
,”
Chem. Phys. Lett.
75
,
87
90
(
1980
).
15.
K. M.
Christoffel
and
J. M.
Bowman
, “
Investigations of self-consistent field, SCF CI and virtual state configuration interaction vibrational energies for a model three-mode system
,”
Chem. Phys. Lett.
85
,
220
224
(
1982
).
16.
G.
Avila
and
T.
Carrington
, “
Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D
,”
J. Chem. Phys.
134
,
054126
(
2011
).
17.
A.
Leclerc
and
T.
Carrington
, “
Calculating vibrational spectra with sum of product basis functions without storing full-dimensional vectors or matrices
,”
J. Chem. Phys.
140
,
174111
(
2014
).
18.
M.
Odunlami
,
V.
Le Bris
,
D.
Bégué
,
I.
Baraille
, and
O.
Coulaud
, “
A-VCI: A flexible method to efficiently compute vibrational spectra
,”
J. Chem. Phys.
146
,
214108
(
2017
).
19.
A.
Baiardi
,
C. J.
Stein
,
V.
Barone
, and
M.
Reiher
, “
Vibrational density matrix renormalization group
,”
J. Chem. Theory Comput.
13
,
3764
3777
(
2017
).
20.
B.
Huron
,
J. P.
Malrieu
, and
P.
Rancurel
, “
Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions
,”
J. Chem. Phys.
58
,
5745
5759
(
1973
).
21.
R. J.
Buenker
and
S. D.
Peyerimhoff
, “
Individualized configuration selection in CI calculations with subsequent energy extrapolation
,”
Theor. Chim. Acta
35
,
33
58
(
1974
).
22.
R. J.
Harrison
, “
Approximating full configuration interaction with selected configuration interaction and perturbation theory
,”
J. Chem. Phys.
94
,
5021
5031
(
1991
).
23.
J. B.
Schriber
and
F. A.
Evangelista
, “
Communication: An adaptive configuration interaction approach for strongly correlated electrons with tunable accuracy
,”
J. Chem. Phys.
144
,
161106
(
2016
).
24.
N. M.
Tubman
,
J.
Lee
,
T. Y.
Takeshita
,
M.
Head-Gordon
, and
K. B.
Whaley
, “
A deterministic alternative to the full configuration interaction quantum Monte Carlo method
,”
J. Chem. Phys.
145
,
044112
(
2016
).
25.
N. M.
Tubman
,
C. D.
Freeman
,
D. S.
Levine
,
D.
Hait
,
M.
Head-Gordon
, and
K. B.
Whaley
, “
Modern approaches to exact diagonalization and selected configuration interaction with the adaptive sampling CI method
,”
J. Chem. Theory Comput.
16
,
2139
2159
(
2020
).
26.
A. A.
Holmes
,
N. M.
Tubman
, and
C. J.
Umrigar
, “
Heat-bath configuration interaction: An efficient selected configuration interaction algorithm inspired by heat-bath sampling
,”
J. Chem. Theory Comput.
12
,
3674
3680
(
2016
).
27.
A. A.
Holmes
,
C. J.
Umrigar
, and
S.
Sharma
, “
Excited states using semistochastic heat-bath configuration interaction
,”
J. Chem. Phys.
147
,
164111
(
2017
).
28.
S.
Sharma
,
A. A.
Holmes
,
G.
Jeanmairet
,
A.
Alavi
, and
C. J.
Umrigar
, “
Semistochastic heat-bath configuration interaction method: Selected configuration interaction with semistochastic perturbation theory
,”
J. Chem. Theory Comput.
13
,
1595
1604
(
2017
).
29.
J.
Li
,
M.
Otten
,
A. A.
Holmes
,
S.
Sharma
, and
C. J.
Umrigar
, “
Fast semistochastic heat-bath configuration interaction
,”
J. Chem. Phys.
149
,
214110
(
2018
).
30.
K. T.
Williams
,
Y.
Yao
,
J.
Li
,
L.
Chen
,
H.
Shi
,
M.
Motta
,
C.
Niu
,
U.
Ray
,
S.
Guo
,
R. J.
Anderson
,
J.
Li
,
L. N.
Tran
,
C. N.
Yeh
,
B.
Mussard
,
S.
Sharma
,
F.
Bruneval
,
M.
Van Schilfgaarde
,
G. H.
Booth
,
G. K. L.
Chan
,
S.
Zhang
,
E.
Gull
,
D.
Zgid
,
A.
Millis
,
C. J.
Umrigar
, and
L. K.
Wagner
, “
Direct comparison of many-body methods for realistic electronic Hamiltonians
,”
Phys. Rev. X
10
,
011041
(
2020
).
31.
Y.
Yao
,
E.
Giner
,
J.
Li
,
J.
Toulouse
, and
C. J.
Umrigar
, “
Almost exact energies for the Gaussian-2 set with the semistochastic heat-bath configuration interaction method
,”
J. Chem. Phys.
153
,
124117
(
2020
).
32.
J. J.
Eriksen
,
T. A.
Anderson
,
J. E.
Deustua
,
K.
Ghanem
,
D.
Hait
,
M. R.
Hoffmann
,
S.
Lee
,
D. S.
Levine
,
I.
Magoulas
,
J.
Shen
,
N. M.
Tubman
,
K. B.
Whaley
,
E.
Xu
,
Y.
Yao
,
N.
Zhang
,
A.
Alavi
,
G. K.-L.
Chan
,
M.
Head-Gordon
,
W.
Liu
,
P.
Piecuch
,
S.
Sharma
,
S. L.
Ten-no
,
C. J.
Umrigar
, and
J.
Gauss
, “
The ground state electronic energy of benzene
,”
J. Phys. Chem. Lett.
11
,
8922
8929
(
2020
).
33.
G.
Rauhut
, “
Configuration selection as a route towards efficient vibrational configuration interaction calculations
,”
J. Chem. Phys.
127
,
184109
(
2007
).
34.
M.
Neff
and
G.
Rauhut
, “
Toward large scale vibrational configuration interaction calculations
,”
J. Chem. Phys.
131
,
124129
(
2009
).
35.
P.
Carbonnière
,
A.
Dargelos
, and
C.
Pouchan
, “
The VCI-P code: An iterative variation-perturbation scheme for efficient computations of anharmonic vibrational levels and IR intensities of polyatomic molecules
,”
Theor. Chem. Acc.
125
,
543
554
(
2010
).
36.
M.
Sibaev
and
D. L.
Crittenden
, “
Balancing accuracy and efficiency in selecting vibrational configuration interaction basis states using vibrational perturbation theory
,”
J. Chem. Phys.
145
,
064106
(
2016
).
37.
Y.
Scribano
and
D. M.
Benoit
, “
Iterative active-space selection for vibrational configuration interaction calculations using a reduced-coupling VSCF basis
,”
Chem. Phys. Lett.
458
,
384
387
(
2008
).
38.
E.
Lesko
,
M.
Ardiansyah
, and
K. R.
Brorsen
, “
Vibrational adaptive sampling configuration interaction
,”
J. Chem. Phys.
151
,
164103
(
2019
).
39.
E. G.
Kratz
, LOVCI: Ladder operator vibrational configuration iteraction, https://github.com/kratman/VibCI,
2016
.
40.
J. H.
Fetherolf
, VHCI: Vibrational heat-bath configuration interaction v0.1, https://github.com/berkelbach-group/VHCI,
2020
.
41.
G.
Guennebaud
,
B.
Jacob
 et al, Eigen v3, http://eigen.tuxfamily.org,
2010
.
42.
Y.
Qiu
, SPECTRA: Sparse eigenvalue computation toolkit as a redesigned ARPACK, https://spectralib.org/,
2015
.
43.
J.
Brown
and
T.
Carrington
, “
Using an expanding nondirect product harmonic basis with an iterative eigensolver to compute vibrational energy levels with as many as seven atoms
,”
J. Chem. Phys.
145
,
144104
(
2016
).
44.
P. S.
Thomas
,
T.
Carrington
,
J.
Agarwal
, and
H. F.
Schaefer
, “
Using an iterative eigensolver and intertwined rank reduction to compute vibrational spectra of molecules with more than a dozen atoms: Uracil and naphthalene
,”
J. Chem. Phys.
149
,
064108
(
2018
).
45.
D.
Begue
,
P.
Carbonniere
, and
C.
Pouchan
, “
Calculations of vibrational energy levels by using a hybrid ab initio and DFT quartic force field: Application to acetonitrile
,”
J. Phys. Chem. A
109
,
4611
4616
(
2005
).
46.
R.
Garnier
,
M.
Odunlami
,
V.
Le Bris
,
D.
Bégué
,
I.
Baraille
, and
O.
Coulaud
, “
Adaptive vibrational configuration interaction (A-VCI): A posteriori error estimation to efficiently compute anharmonic IR spectra
,”
J. Chem. Phys.
144
,
204123
(
2016
).
47.
T.
Delahaye
,
A.
Nikitin
,
M.
Rey
,
P. G.
Szalay
, and
V. G.
Tyuterev
, “
A new accurate ground-state potential energy surface of ethylene and predictions for rotational and vibrational energy levels
,”
J. Chem. Phys.
141
,
104301
(
2014
).
48.
M.
Sibaev
and
D. L.
Crittenden
, “
The PyPES library of high quality semi-global potential energy surfaces
,”
J. Comput. Chem.
36
,
2200
2207
(
2015
).
49.
M.
Sibaev
and
D. L.
Crittenden
, “
PyVCI: A flexible open-source code for calculating accurate molecular infrared spectra
,”
Comput. Phys. Commun.
203
,
290
297
(
2016
).
50.
D.
Bégué
,
N.
Gohaud
,
C.
Pouchan
,
P.
Cassam-Chenaï
, and
J.
Liévin
, “
A comparison of two methods for selecting vibrational configuration interaction spaces on a heptatomic system: Ethylene oxide
,”
J. Chem. Phys.
127
,
164115
(
2007
).
51.
E.
Cané
,
A.
Miani
, and
A.
Trombetti
, “
Anharmonic force fields of naphthalene-h8 and naphthalene-d8
,”
J. Phys. Chem. A
111
,
8218
8222
(
2007
).
52.
I.
Baraille
,
C.
Larrieu
,
A.
Dargelos
, and
M.
Chaillet
, “
Calculation of non-fundamental IR frequencies and intensities at the anharmonic level. I. The overtone, combination and difference bands of diazomethane, H2CN2
,”
Chem. Phys.
273
,
91
101
(
2001
).
53.
R.
Burcl
,
S.
Carter
, and
N. C.
Handy
, “
Infrared intensities from the MULTIMODE code
,”
Chem. Phys. Lett.
380
,
237
244
(
2003
).
54.
V.
Le Bris
,
M.
Odunlami
,
D.
Bégué
,
I.
Baraille
, and
O.
Coulaud
, “
Using computed infrared intensities for the reduction of vibrational configuration interaction bases
,”
Phys. Chem. Chem. Phys.
22
,
7021
7030
(
2020
).
55.
K.
Kim
,
K. D.
Jordan
, and
T. S.
Zwier
, “
Low-energy structures and vibrational frequencies of the water hexamer: Comparison with benzene-(H2O)6
,”
J. Am. Chem. Soc.
116
,
11568
11569
(
1994
).
56.
Q.
Yu
and
J. M.
Bowman
, “
Classical, thermostated ring polymer, and quantum VSCF/VCI calculations of IR spectra of H7O3+ and H9O4+ (eigen) and comparison with experiment
,”
J. Phys. Chem. A
123
,
1399
1409
(
2019
).
57.
Z.
Bačić
and
J. C.
Light
, “
Theoretical methods for rovibrational states of floppy molecules
,”
Annu. Rev. Phys. Chem.
40
,
469
498
(
1989
).
58.
M. J.
Bramley
,
J. W.
Tromp
,
T.
Carrington
, and
G. C.
Corey
, “
Efficient calculation of highly excited vibrational energy levels of floppy molecules: The band origins of H+3 up to 35000 cm−1
,”
J. Chem. Phys.
100
,
6175
6194
(
1994
).

Supplementary Material