Selected configuration interaction plus perturbation theory approaches have long been used to solve both the electronic and vibrational Schrödinger equations. In the last few years, many new selection algorithms have been developed for these approaches and applied to solve the electronic Schrödinger equation, but these algorithms have seen little to no use for solving the vibrational Schrödinger equation. Herein, we adapt one of the recently developed approaches, the adaptive sampling configuration interaction (ASCI) method, to calculate the vibrational excitations of molecules. The vibrational ASCI method has accuracy comparable to other high-accuracy approaches for solving the vibrational Schrödinger equation while requiring only modest computer resources. We study two different approaches for calculating excited states with vibrational ASCI and benchmark the method on acetonitrile and ethylene oxide. Finally, we demonstrate the applicability of the vibrational ASCI method to large systems by calculating the 128 lowest energy vibrational states of naphthalene, which has 48 vibrational degrees of freedom.

## I. INTRODUCTION

Computational approaches for calculating vibrational properties play an important role in the interpretation of experimental vibrational spectroscopy.^{1,2} Because of this importance, much effort has been devoted to developing methods for calculating vibrational properties that are highly accurate and can be applied to large chemical systems. These computational approaches normally solve the vibrational Schrödinger equation assuming the Born-Oppenheimer approximation, where the vibrating atoms of the molecule move on a single electronic potential energy surface. In practice, many of the approaches for solving the vibrational Schrödinger equation are similar to the approaches used for solving the electronic Schrödinger equation, such as vibrational self-consistent field,^{3–5} vibrational configuration interaction (VCI),^{3,6–9} vibrational perturbation theory (VPT),^{10–14} or vibrational coupled cluster theory.^{15–17}

When solving the vibrational Schrödinger equation using configuration interaction (CI)-based approaches, the infinite-dimensional Hilbert space is commonly discretized using a finite-basis representation where a vibrational configuration is written as a direct product of normal modes,^{18–21} but other finite-basis representations^{7,22–24} or discrete variable representations^{25,26} can be used. When using a VCI-based approach, the infinite-dimensional Hilbert space is also typically truncated at some maximum normal-mode excitation level, resulting in a finite-dimensional Hilbert space. Full VCI gives the exact answer in this configurational space, but similar to a full CI approach for solving the electronic Schrödinger equation, the number of configurations scales exponentially with respect to system size allowing it to be used for only the smallest systems and/or configurational spaces. In place of a full VCI expansion, a truncated VCI expansion can be used, but the Born-Oppenheimer potential energy surface is almost always anharmonic, which results in strong coupling among the normal modes, the need to include many configurations to describe each vibrational state accurately, and necessitates a truncated VCI expansion with large maximum excitation.^{27}

An alternative to pure VPT or VCI approaches is to combine the two in a unified approach of which there exist two main categories. In the first,^{12,28–34} VCI submatrices are constructed and diagonalized for sets of resonance states and then a perturbation correction is applied to the resulting vibrational states. As this relies on the *a priori* determination of resonance states, it can be difficult to correctly construct the submatrices. Alternatively, a VCI calculation can be performed in a configurational space including all states and then a perturbative correction applied. This requires a careful determination of the configurational space, as all strong couplings among states must be included, while still restricting the space to be small enough to allow the Hamiltonian matrix to be constructed and diagonalized.^{7,27,35–38} This determination is commonly implemented using a second-order VPT (VPT2)-based screening procedure, and methods of this type can be grouped under the name vibrational selected CI + perturbation theory (VSCI+PT). Even with such selection procedures, the size of the configurational space can grow large, resulting in computational and memory bottlenecks, particularly during the perturbative correction. In response to this issue, state-specific VSCI+PT approaches^{7,27,38} have been developed that calculate a single vibrational state of interest, but these methods can be computationally expensive for calculating an entire vibrational spectrum as a separate calculation must be performed for each vibrational state.

Other variational approaches for solving the vibrational Schrödinger equation exist that do not rely on a perturbative-based screening procedure that is seen in many VSCI+PT methods. These approaches are typically grouped into two classes: pruning approaches^{20,39–49} and contraction approaches.^{50–54} In pruning approaches, of which the VSCI+PT methods are a subset, the number of configurations in the configurational space is restricted using a pruning criterion. In contraction approaches, sub-Hamiltonians of the full-VCI Hamiltonian are constructed and diagonalized. Of the additional variational approaches, we highlight the adaptive VCI (A-VCI) method,^{55} which discretizes the Hamiltonian using nested bases, and the reduced-rank block power method (RRBPM)^{56–59} that factorizes the VCI matrix in canonical polyadic (CP) format^{60,61} as these two methods will later be used for benchmarking.

Of course, CI-based techniques have a long history of being used to solve the electronic Schrödinger equation with the same computational-expense related problems as VCI-based techniques for solving the vibrational Schrödinger equation.^{62} However, some of the most exciting developments in electronic structure theory in the past two decades have been new methods for finding an approximate full-CI solution such as the density matrix renormalization group (DMRG)^{63–66} or full CI quantum Monte Carlo (FCIQMC) methods.^{67,68} Using these methods, full-CI or near full-CI calculations on systems with active spaces of unprecedented size have been performed. Similar to the RRBPM family of methods for solving the vibrational Schrödinger equation, the DMRG method decomposes the full CI matrix using a tensor factorization. The DMRG method uses matrix product states (MPSs)^{66,69–71} rather than CP tensors. The MPS factorization is better at including strong correlation among degrees of freedom, and it has been hypothesized^{72} that the difficulty of including strong correlation using CP tensors might be why a CP tensor decomposition has not been used to solve the electronic Schrödinger equation. FCIQMC is a stochastic technique that samples the full CI space using walkers that propagate through the space of determinants.

There is also a long history of using selected CI plus perturbation theory methods (SCI+PT)^{73–75} to solve the electronic Schrödinger equation beginning with the CI using a perturbation selection done iteratively (CIPSI) method in 1973. Many of these SCI+PT methods designed for solving the electronic Schrödinger equation provide the initial inspiration for the VSCI+PT methods referenced previously. The last few years have seen a large renewal of interest in SCI+PT methods for solving the electronic Schrödinger equation^{76–87} and the development of powerful new selection schemes and ranking algorithms for finding the most important determinants to include in the configurational space as well as algorithmic improvements to speed up the perturbative energy correction. These new SCI+PT methods can perform high-accuracy approximate complete active space CI (CASCI) calculations with active spaces as large as 34 electrons in 152 orbitals on the Si_{2}H_{6} molecule^{78} with the adaptive sampling CI (ASCI) method and 28 electrons in 76 spatial orbitals with a relativistic Hamiltonian for the chromium dimer with the stochastic heat-bath CI method.^{83}

The adaption and application of the DMRG, FCIQMC, and recent SCI+PT methods to the solution of the vibrational Schrödinger equation has been slow. Although the DMRG was first used in quantum chemistry in the late 1990s,^{64,88–92} an implementation of the DMRG for calculating vibrational spectra (vDMRG) was not developed until 2017.^{93} FCIQMC^{67} was introduced in 2009, and to the best of our knowledge, no vibrational FCIQMC method has been derived. The original vDMRG study^{93} found excited states by performing a ground-state vDMRG calculation in a configurational space orthogonal to all lower energy states, which requires finding the *n*th excited state before the *n*+1th excited state can be found, which slows down the finding of higher-energy fundamentals of large molecules as the vibrational spectrum can be dense. More recently,^{72} vDMRG has employed shift-and-invert methods^{94,95} and the folded-spectrum operator^{96} to target excited states in specific energy regions at an increased computational cost relative to calculating the same number of excited states in the lowest energy region of the vibrational spectrum.

If a configurational space is well balanced between configurations that contribute to the ground and excited states, an SCI+PT method offers a straightforward way of finding excited states compared to DMRG- or FCIQMC-based approaches as diagonalization of the CI matrix gives the energy of all states. This conceptual simplicity likely explains in part the initial development and wide use of SCI+PT for solving the vibrational Schrödinger equation. While SCI+PT methods have been introduced in the last decade for solving the vibrational Schrödinger equation,^{35,36,38} these methods are inspired by the previous generation of SCI+PT-based approaches for solving the electronic Schrödinger equation. To the best of the authors’ knowledge, no VSCI+PT approach has made use of the powerful new selection algorithms that have been developed in the last few years.

In this study, we focus on adapting one of these recently introduced SCI+PT methods, the adaptive sampling configuration interaction (ASCI) method,^{77,78} to calculate vibrational excitations. The vibrational ASCI (VASCI) method is conceptually simple, easy to implement, and even with the rudimentary implementation presented here, is able to calculate vibrational excitations of medium and large molecules at a level of accuracy similar to that of the current best variational methods for solving the vibrational Schrödinger equation. The VASCI method only requires modest computer resources, which when combined with its high accuracy, demonstrates the power of these recently introduced SCI+PT methods for calculating vibrational excitations.

## II. METHODS

### A. Adaptive sampling configuration interaction

We briefly outline the main steps of the ASCI method.^{77} Although the ASCI method has seen recent algorithmic improvements centered on the use of sorting algorithms,^{78} the goal of the present study is to demonstrate proof of principle for the applicability of the ASCI method to the solution of the vibrational Schrödinger equation. Because of this, we will ignore many of the recent improvements to the ASCI method in our VASCI implementation.

In the ASCI method, the wave function is expanded as a linear combination of determinants, $Di$, as

The ASCI algorithm has two distinct stages: a variational stage and a perturbative stage. In the variational stage, the best set of determinants and their corresponding coefficients of the summation in Eq. (1) are found. During the variational stage, two determinant subspaces are utilized: a core space and a target space. The number of determinants in these subspaces is specified by the parameters *cdets* and *tdets*, respectively. The core space and initial wave function are initialized in a determinant space of size less than or equal to *cdets*. From this core space, all connected determinants, i.e., all determinants not in the core space with the nonzero Hamiltonian matrix element with a determinant in the core space, are found and then using the set of wave function coefficients {*C*_{i}} from the initial guess and the Hamiltonian matrix elements, an estimate of the wave function coefficient, *A*_{i}, for all core and connected determinants is found using an equation similar to the first-order perturbation estimate of the CI coefficient in Epstein-Nesbet perturbation theory,^{75}

From the set of coefficients {*A*_{i}}, the *tdets* determinants corresponding to the coefficients of largest absolute value are selected as the target space. The Hamiltonian is then constructed and diagonalized in the target space basis. The eigenpair with the lowest eigenvalue is defined as the ground state energy and ground state wave function. From this ground state wave function, the *cdets* determinants with coefficients of the largest absolute value from the set of {*C*_{i}} in Eq. (1) are chosen as the new core space. Using this core space, the ASCI procedure repeats starting with an estimation of the first-order perturbation estimate of the CI coefficients for all core and connected determinants using Eq. (2). The entire variational stage of the ASCI procedure is repeated until the change in the ground state energy between iterations has fallen below some predefined threshold, which was set to 10^{−5} hartree in the original study.

After the variational stage, a perturbative stage begins where the second-order Epstein-Nesbet perturbative correction to the energy^{75} is calculated as

where the sum runs over all determinants that are not in the final target space but are connected to the wave function, *ψ*_{elec}. The original implementation^{77} of the perturbative stage was computationally demanding and required a large amount of memory due to the need to store a partial sum for each determinant that is connected to a determinant in the final target space, which limited its applicability for large systems. Recent algorithmic improvements to the perturbative energy correction, which will not be implemented in this study, have largely eliminated the memory bottleneck, while also increasing parallel efficiency.^{78} Alternative approaches for evaluating Eq. (3) include stochastic and semistochastic methods,^{86} which both eliminate the memory bottleneck and greatly reduce the central processing unit (CPU) requirements compared to the original ASCI approach. We note that there is some debate in the literature about the need for stochastic approaches, as it has been claimed^{97} that an efficient deterministic algorithm can be faster than the current state of the art stochastic approaches by multiple orders of magnitude. At present, the issue of which approach is best remains unsettled.

The ASCI algorithm can be considered to be a refinement of the CIPSI algorithm.^{75} The key improvement of the ASCI algorithm is the use of a core space with *cdets* determinants, which both keeps the number of connected determinants that must be evaluated using Eq. (2) small and speeds up the calculation of Eq. (2). This is essential for ensuring a fast method as the number of connected determinants rapidly increases with system size. While the ranking algorithm in Eq. (2) has been shown to efficiently select the most important determinants to include in the target space, other recently developed selected CI methods have ranking algorithms based on an integral-driven search metric, whereby a determinant $Di$ is added to the target space if $HijCj$ is greater than a cutoff value, where *j* corresponds to a determinant $Dj$ in the core space. Such an approach is used in the heat-bath CI family of methods.^{76,83} The evaluation of the ranking criteria in integral-driven search methods is faster than the evaluation of Eq. (2) in the ASCI method. It has been argued that the increase in computational expense due to the use of Eq. (2) is offset by a better choice of determinants for the target space.^{97} This leads to a more compact representation of the wave function, thereby reducing the computational expense of the perturbative energy correction as a smaller number of determinants in the target space can be used to obtain a given level of accuracy. Much like whether a deterministic vs stochastic approach for computation of the perturbative energy correction is best, this issue remains unsettled.

### B. Vibrational adaptive sampling configuration interaction

The adaptation of the ASCI method to find the ground state of the vibrational Schrödinger equation requires minimal changes to the ASCI algorithm. Instead of the electronic Hamiltonian, we use an approximate Watson Hamiltonian. For this study, where we will use quartic force fields and neglect Coriolis effects, the Hamiltonian is written as

where *Q*_{i} is the *i*th normal mode coordinate and *F*’s are the force constants. All summations in Eq. (4) run over all normal modes, the number of which is given by *N*_{vib}. Additionally, rather than a basis of Slater determinants, in vibrational ASCI (VASCI), a basis of Hartree products of vibrational normal modes is used with the wave function written as a linear combination of these basis states,

where **n** is a vector giving the vibrational excitation level of each of the normal modes of the Hartree product, $\Phi nQ1,\u2026,QNvib$. $\Phi nQ1,\u2026,QNvib$ is defined as

where $\varphi ni$ is a harmonic oscillator basis function for mode *i* with excitation level *n*_{i} specified from vector **n**.

Using these definitions in VASCI, the wave function coefficient estimate corresponding to Eq. (2) in ASCI is

while the perturbative energy correction corresponding to Eq. (3) in ASCI is

where the sum in Eq. (8) is over all Hartree products not in the target space, but with the nonzero matrix element with a Hartree product in the target space.

The generation of all Hartree products connected to the current core space and the evaluation of the wave function coefficient estimates in Eq. (7) are more computationally demanding in VASCI than the generation of connected determinants and the evaluation of the wave function coefficient estimate using Eq. (2) in ASCI as the vibrational Hamiltonian couples states that differ by up to four quanta rather than by up to two spin orbitals as for the electronic Hamiltonian. Therefore, when generating the connected Hartree products for the wave function coefficient estimate, any Hartree products that differ by up to four quanta with a Hartree product in the current core space must be included in the connected Hartree product space. This can be computationally expensive. Consequently, we will explore truncating the generation of the connected Hartree product space such that only Hartree products that differ by up to two or three quanta are included. As shown in Sec. IV, this has a negligible effect on the accuracy of the VASCI method. Similarly, the generation of the space of connected Hartree products for the calculation of the perturbative energy correction in Eq. (8) should include all Hartree products that differ by four quanta or less with a Hartree product in the final target space. We will also explore reducing the size of the connected Hartree product space in Eq. (8) by truncating the generation of the connected Hartree products to all Hartree products that differ by two or three quanta from a Hartree product in the final target space. All other aspects of a VASCI calculation of the vibrational ground state wave function and energy proceed similarly to that of the original ASCI method. As the vibrational wave function is a linear combination of Hartree products rather than Slater determinants, we rename the variables *cde*ts and *tdets* in the ASCI method to *cprods* and *tprods* in the VASCI method for consistency.

While the vibrational ground state wave function provides important information about a system such as the zero point vibrational energy (ZPVE), the goal of a VCI calculation is typically to calculate the energies of excited states relative to the ground state. Diagonalization of a Hamiltonian, as in the ASCI method, gives the full energy eigenvalue spectrum and the original ASCI study used this spectrum to find the energies of excited electronic states.^{77} Unfortunately, as the selection algorithm of the original ASCI method focuses on finding the best set of determinants for the ground state, the calculated energies of the excited states were poor. Therefore, the adaption of the ASCI method to calculate the energies of vibrational excited states requires a modification of the original ASCI method. We note that another recently developed SCI+PT method, stochastic heat-bath CI (SHCI), has had more success in calculating electronic excited states by using a selection algorithm that includes determinants that are important for both ground and excited states in the configurational space.^{82,84}

We propose two modifications of the original ASCI algorithm to calculate excited states using the VASCI method. The first proposed method, called type I, is similar to the approach adopted for calculating excited states in the SHCI method.^{84} It begins with a specification of the number of states that will be obtained, *N*_{states}. The number of excited states obtained will be one less than *N*_{states}. The variational stage then proceeds. An initial guess for all *N*_{states} states is obtained from a truncated VCI calculation. All Hartree products in the truncated VCI calculation define the initial core space. From the coefficients of each of the vibrational states, a maximum coefficient vector, **C**_{max}, is constructed as

where *C*_{j,i} is the *i*th expansion coefficient in Eq. (5) of the *j*th vibrational state. The elements of **C**_{max} are then used to generate the wave function coefficient estimate for all core and connected Hartree products using Eq. (7). The Hartree products corresponding to the *tprods* estimated wave function coefficient estimates with the largest absolute value define the target space. The Hamiltonian is then constructed and diagonalized in the target space basis. The *N*_{states} eigenpairs with the lowest eigenvalue are taken to be the *N*_{states} vibrational wave functions and vibrational energies. From the *N*_{states} sets of wave function coefficients, a new maximum coefficient vector is found using Eq. (9). The Hartree products corresponding to the *cprods* elements of the maximum coefficient vector with the largest absolute value are taken to be the new core space and the variational stage repeats by estimating the wave function coefficients using Eq. (7) for all core and connected Hartree products. The variational stage of type I VASCI is repeated until the change in energy of every state between iterations is less than a user-defined cutoff value. Once the variational stage has converged, the perturbative stage occurs identically to a ground-state calculation, but now the perturbation energy of *N*_{states} states must found using Eq. (8). The final vibrational excitation energies are then calculated as the difference between the perturbation-corrected ground and excited state energies.

The second proposed modification of the ASCI algorithm to calculate excited states, called type II, begins similarly to type I VASCI with a specification of *N*_{states} and the generation of an initial guess for the wave function coefficients of the *N*_{states} states using truncated VCI. As in type I VASCI, the number of excited states obtained is one less than *N*_{states}. For type II VASCI, we define **C**_{k} as the vector of wave function coefficients from Eq. (5) of state *k.* For each state *k*, a set of wave function coefficient estimates, {*A*_{i}}_{k}, are obtained using Eq. (7) for all core and connected Hartree products. From each set of wave function coefficient estimates, the Hartree products corresponding to the *tprods* elements of $Aik$ with the largest absolute value are combined to form a total target space.

For a VASCI type II calculations with *N*_{states} states and a value of *tprods*, the maximum number of Hartree products in the total target space is *N*_{states} ^{*} *tprods*, but as some Hartree products will be included from multiple vibrational states, in practice, the total number of Hartree products in the total target space is less. The Hamiltonian is constructed and diagonalized in the total target space basis. The *N*_{states} eigenpairs with smallest eigenvalues are taken to be the *N*_{states} vibrational wave functions and vibrational energies. For each vibrational state, the set of Hartree products corresponding to the *cprod* wave function coefficients with the largest absolute value in Eq. (5) is found and these sets are combined to form a total core space. Similar to the construction of the total target space, the maximum size of the total core space is *N*_{states} ^{*} *cprods*, but due to Hartree products being included in the core space from more than one state, the size of the total core space is normally less. With the determination of the total core space, the variational stage of type II VASCI repeats starting with the determination of a set of wave function coefficient estimates of the core and connected Hartree products for each state using Eq. (7). The variational stage of the type II VASCI procedure ends when the change in energy of every state between iterations is less than a user-defined cutoff value. With the end of the variational stage, the perturbative stage begins and is performed identically to type I VASCI to obtain the final vibrational excitation energies.

The type II approach is similar to performing a VASCI calculation on each of the *N*_{states} vibrational states, but with the total core and total target spaces being the union of the core and target spaces of each of the individual states, respectively. Therefore, each state is approximately treated equally as the number of Hartree products included in the total core and total target spaces are the same for each vibrational state. The overlaps of the core and target spaces of the states change each iteration leading to the size of the total core and total target spaces changing each iteration. In the type II approach, it is difficult to estimate *a priori* what the maximum size of the total core and total target spaces will be in order to determine memory requirements and estimated run times for the calculations.

In type I VASCI, the number of direct products in the core and target spaces remains constant from iteration to iteration, which gives more control over memory requirements and making it easier to estimate the run time. For each of the *N*_{states} states, there is no minimum number of Hartree products estimated to be important for describing that state included in the target space. Although this could be seen as a drawback of type I VASCI, in practice, we hypothesize it to be beneficial as states that are highly localized will have fewer Hartree products devoted to their description in the target space, while delocalized vibrational states will have more, which could lead to a more balanced overall description for a given value of *tprods*.

The VASCI method shares some similarities to other methods for calculating the vibrational excitation energies of molecules.^{35,36,38} These similarities are to be expected, as the ASCI method is a refinement of the CIPSI method that provides the foundation for many VSCI+PT methods. The key improvement of the VASCI method is the restriction of the space for which the sum in Eq. (7) is calculated to Hartree products connected to one of the *cprods* Hartree products in the core space. The calculation of Eq. (7) is computationally expensive, and by restricting the number of Hartree products that need to be checked, VASCI calculations can be run with a large number of Hartree products in the target space.

## III. IMPLEMENTATION AND METHODOLOGY

All VASCI calculations were performed with an in-house Python code. For certain portions of performance critical code such as the calculation of the wave function coefficient estimates in Eq. (7) and the construction of the Hamiltonian matrix elements, Cython was used to increase computational efficiency. Parts of the Cython code used OpenMP to enable parallel computation. We note that this Python code is not optimized and large speed ups through further optimization of the code are likely to be possible. As will be seen, even with this nonoptimized code, the VASCI method is able to efficiently calculate vibrational excitations relative to other computational techniques. For all VASCI calculations, the Hamiltonian matrix is stored in sparse matrix format and diagonalized using the sparse diagonalization routines in SciPy. The SciPy sparse diagonalization routines are a wrapper to ARPACK,^{98} which uses the implicitly restarted Lanczos method^{99} to calculate the *N*_{states} lowest eigenpairs.

For VCI calculations that use a normal mode direct product basis, there are an infinite number of Hartree product basis functions due to the normal mode excitation level having no upper bound. Because of this, full VCI and many VSCI+PT calculations are performed in a restricted space where the total number of vibrational excitation quanta, $ntot=\u2211iNvibni$, where *n*_{i} is the excitation level of mode *i*, is less than or equal to an user-defined value. This is referred to as a polyad truncation in the literature.^{45,100} We note that other truncation scheme are possible.^{45,101} For all of the VASCI calculations, we will employ the restriction *n*_{tot} ≤ 10. In practice, increasing the restriction of *n*_{tot} ≤ 10 to a larger value presents no difficulty.

For all VASCI calculations, the variational stage is iterated until the change of energy relative to the previous iteration has fallen below 0.5 cm^{−1} for all states. The energy of the *i*th excited vibrational state, *ν*_{i}, is expressed relative to the ZPVE or ground state, *E*_{0}, as

where *E*_{i} is the absolute energy of the *i*th excited vibrational state. State assignment for each vibrational state was determined by finding the Hartree products that correspond to the coefficient with the largest absolute value in Eq. (5).

All VASCI calculations were performed using Intel Xeon Gold 5138 2.00 GHz CPUs.

## IV. COMPUTATIONAL RESULTS AND DISCUSSION

### A. Acetonitrile

The first test system is the acetonitrile molecule, which has 12 normal modes and is a standard benchmark for newly introduced methods for calculating vibrational excitations.^{24,59,93} To enable a comparison of the VASCI method to previous results, we use the quartic force field of Pouchan and co-workers^{102} for which force field constants are computed at the CCSD(T)/cc-pVTZ level for the harmonic constants and B3LYP/cc-pVTZ level for the cubic and quartic constants. Similar to previous studies, we calculate the 70 vibrational states with lowest energy for all VASCI calculations.

The first test of the VASCI method is on the relative accuracies of the type I and type II VASCI methods. As the maximum size of the total core and target spaces are unknown for type II VASCI beforehand, three type II VASCI calculations were performed first using values of *cprods* and *tprod*s for each of the individual states of 25 and 800, 50 and 2000, and 50 and 2500, respectively. These three type II VASCI calculations had a maximum size of the total core and total target spaces of 800 and 13 408, 926 and 24 723, and 1441 and 30 038, respectively. These values were then used as *cprods* and *tprods* for three type I VASCI calculations. For all of these calculations, an initial singles and doubles VCI (VCISD) guess was used and the connected Hartree product space was restricted to connected Hartree products singly or doubly excited from a Hartree product in the core or target space and not in the core or target space, respectively. A subset of these results from these calculations is found in Table I. Due to the small differences between the type I and type II VASCI calculations, we present all energies to 0.01 cm^{−1}. For the rest of the study, we will only present results to 0.1 cm^{−1}. The full results are found in the supplementary material.

cprods
. | 800 . | 926 . | 1441 . | . | |||
---|---|---|---|---|---|---|---|

tprods
. | 13 408 . | 24 723 . | 30 038 . | . | |||

Type . | I . | II . | I . | II . | I . | II . | Smolyak^{a}
. |

ZPVE | 9837.44 | 9837.42 | 9837.43 | 9837.42 | 9837.42 | 9837.41 | 9837.41 |

ν11 | 361.02 | 361.02 | 361.01 | 361.01 | 361.01 | 361.01 | 360.99 |

361.02 | 361.02 | 361.01 | 361.01 | 361.01 | 361.01 | 360.99 | |

2ν11 | 723.24 | 723.24 | 723.23 | 723.23 | 723.22 | 723.22 | 723.18 |

723.24 | 723.25 | 723.23 | 723.24 | 723.23 | 723.23 | 723.18 | |

2ν11 | 723.89 | 723.91 | 723.88 | 723.90 | 723.89 | 723.89 | 723.83 |

ν4 | 900.69 | 900.70 | 900.68 | 900.69 | 900.68 | 900.68 | 900.66 |

ν9 | 1034.17 | 1034.17 | 1034.14 | 1034.14 | 1034.14 | 1034.14 | 1034.13 |

1034.17 | 1034.17 | 1034.14 | 1034.14 | 1034.14 | 1034.14 | 1034.13 | |

3ν11 | 1086.68 | 1086.68 | 1086.67 | 1086.66 | 1086.66 | 1086.64 | 1086.55 |

3ν11 | 1086.68 | 1086.68 | 1086.67 | 1086.66 | 1086.66 | 1086.64 | 1086.55 |

3ν11 | 1087.92 | 1087.92 | 1087.91 | 1087.91 | 1087.89 | 1087.88 | 1087.78 |

1087.92 | 1087.92 | 1087.91 | 1087.90 | 1087.90 | 1087.88 | 1087.78 | |

ν4 + ν11 | 1259.88 | 1259.90 | 1259.87 | 1259.88 | 1259.86 | 1259.86 | 1259.88 |

1259.91 | 1259.91 | 1259.89 | 1259.89 | 1259.87 | 1259.86 | 1259.88 | |

ν3 | 1389.02 | 1389.03 | 1388.99 | 1389.00 | 1388.99 | 1388.99 | 1388.97 |

ν9 + ν11 | 1394.80 | 1394.83 | 1394.77 | 1394.79 | 1394.77 | 1394.76 | 1394.69 |

1394.83 | 1394.81 | 1394.78 | 1394.78 | 1394.79 | 1394.79 | 1394.69 | |

ν9 + ν11 | 1395.03 | 1395.04 | 1394.99 | 1394.99 | 1395.00 | 1395.00 | 1394.91 |

ν9 + ν11 | 1397.81 | 1397.81 | 1397.78 | 1397.79 | 1397.78 | 1397.77 | 1397.69 |

MAE | 0.55 | 0.93 | 0.53 | 0.54 | 0.49 | 0.50 | |

RMSE | 0.26 | 0.29 | 0.23 | 0.24 | 0.20 | 0.20 |

cprods
. | 800 . | 926 . | 1441 . | . | |||
---|---|---|---|---|---|---|---|

tprods
. | 13 408 . | 24 723 . | 30 038 . | . | |||

Type . | I . | II . | I . | II . | I . | II . | Smolyak^{a}
. |

ZPVE | 9837.44 | 9837.42 | 9837.43 | 9837.42 | 9837.42 | 9837.41 | 9837.41 |

ν11 | 361.02 | 361.02 | 361.01 | 361.01 | 361.01 | 361.01 | 360.99 |

361.02 | 361.02 | 361.01 | 361.01 | 361.01 | 361.01 | 360.99 | |

2ν11 | 723.24 | 723.24 | 723.23 | 723.23 | 723.22 | 723.22 | 723.18 |

723.24 | 723.25 | 723.23 | 723.24 | 723.23 | 723.23 | 723.18 | |

2ν11 | 723.89 | 723.91 | 723.88 | 723.90 | 723.89 | 723.89 | 723.83 |

ν4 | 900.69 | 900.70 | 900.68 | 900.69 | 900.68 | 900.68 | 900.66 |

ν9 | 1034.17 | 1034.17 | 1034.14 | 1034.14 | 1034.14 | 1034.14 | 1034.13 |

1034.17 | 1034.17 | 1034.14 | 1034.14 | 1034.14 | 1034.14 | 1034.13 | |

3ν11 | 1086.68 | 1086.68 | 1086.67 | 1086.66 | 1086.66 | 1086.64 | 1086.55 |

3ν11 | 1086.68 | 1086.68 | 1086.67 | 1086.66 | 1086.66 | 1086.64 | 1086.55 |

3ν11 | 1087.92 | 1087.92 | 1087.91 | 1087.91 | 1087.89 | 1087.88 | 1087.78 |

1087.92 | 1087.92 | 1087.91 | 1087.90 | 1087.90 | 1087.88 | 1087.78 | |

ν4 + ν11 | 1259.88 | 1259.90 | 1259.87 | 1259.88 | 1259.86 | 1259.86 | 1259.88 |

1259.91 | 1259.91 | 1259.89 | 1259.89 | 1259.87 | 1259.86 | 1259.88 | |

ν3 | 1389.02 | 1389.03 | 1388.99 | 1389.00 | 1388.99 | 1388.99 | 1388.97 |

ν9 + ν11 | 1394.80 | 1394.83 | 1394.77 | 1394.79 | 1394.77 | 1394.76 | 1394.69 |

1394.83 | 1394.81 | 1394.78 | 1394.78 | 1394.79 | 1394.79 | 1394.69 | |

ν9 + ν11 | 1395.03 | 1395.04 | 1394.99 | 1394.99 | 1395.00 | 1395.00 | 1394.91 |

ν9 + ν11 | 1397.81 | 1397.81 | 1397.78 | 1397.79 | 1397.78 | 1397.77 | 1397.69 |

MAE | 0.55 | 0.93 | 0.53 | 0.54 | 0.49 | 0.50 | |

RMSE | 0.26 | 0.29 | 0.23 | 0.24 | 0.20 | 0.20 |

^{a}

Reference 24.

Table I shows that both type I and type II VASCI accurately calculate the vibrational spectra of the lowest 70 states of acetonitrile with the most accurate calculations having a maximum absolute error (MAE) of 0.49 cm^{−1} and a root-mean-squared error (RMSE) of 0.20 cm^{−1} for type I VASCI, and a MAE of 0.50 cm^{−1} and a RMSE of 0.20 cm^{−1} for type II VASCI. For a given number of *cprods* and *tprods*, the type I and type II VASCI methods perform similarly. Therefore, because of the additional control over memory requirements and run time that a type I VASCI calculation possesses, we will use type I VASCI for the rest of the study. We will refer to type I VASCI as simply VASCI for the rest of the study when there is no ambiguity of context.

Additional benchmarking calculations of the VASCI method were performed where different maximum excitation levels of the truncated VCI expansion were used for the starting guess. Benchmarking calculations were also run where the connected Hartree product space for which the wave function coefficient estimates in Eq. (7) were generated or the perturbative energy sum was performed over in Eq. (8) included Hartree products that were not only doubly but also triply or quadruply excited from a Hartree product in the core or final target space. From these calculations, it was found that the starting guess and the maximum allowed excitation level of the connected Hartree products had very small effects on the accuracy of the VASCI calculations on acetonitrile. Tables of results from these calculations are found in the supplementary material. Therefore, for the rest of the study, we will start with the smallest truncated VCI expansion that still provides enough states that the initial diagonalization can obtain *N*_{states} states and when generating the space of connected Hartree products, restrict it to those that are singly or doubly excited from a Hartree product in the core or target space and not in the core or target space, respectively.

We also examined the improvement in the calculated vibrational excitation energies due to the perturbative energy correction in Eq. (8). A subset of these results is presented in Table II. The full set of results is found in the supplementary material. The perturbative energy correction is seen to have a large effect on the accuracy of the VASCI method with the MAE decreasing by at least a factor of three and the RMSE decreasing by at least a factor of five for all combinations of *cprods* and *tprods*. As expected, the perturbative energy correction is larger for systems with smaller values of *cprods* and *tprods*, as the variational energy for these smaller values is less accurate. For the largest values of *cprods* and *tprods* (840 and 30 000, respectively), the variational calculation gives good agreement with the reference values with a MAE of 1.8 cm^{−1} and a RMSE of 1.0 cm^{−1}, but as the perturbative correction always increases the accuracy in the results presented here, a VASCI calculation should make use of it if computationally feasible.

cprods
. | 210 . | 420 . | 420 . | 840 . | 840 . | . | |||||
---|---|---|---|---|---|---|---|---|---|---|---|

tprods
. | 10 000 . | 10 000 . | 20 000 . | 20 000 . | 30 000 . | . | |||||

. | Var . | Pert . | Var . | Pert . | Var . | Pert . | Var . | Pert . | Var . | Pert . | Smolyak^{a}
. |

ZPVE | 9844.6 | 9837.9 | 9837.8 | 9837.5 | 9837.7 | 9837.5 | 9837.7 | 9837.4 | 9837.7 | 9837.4 | 9837.4 |

ν11 | 361.1 | 361.0 | 361.1 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 |

362.1 | 361.1 | 361.1 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | |

2ν11 | 723.7 | 723.3 | 723.4 | 723.2 | 723.4 | 723.2 | 723.3 | 723.2 | 723.3 | 723.2 | 723.2 |

724.8 | 723.4 | 728.7 | 723.5 | 728.7 | 723.5 | 723.3 | 723.2 | 723.3 | 723.2 | 723.2 | |

2ν11 | 725.4 | 724.0 | 729.4 | 724.2 | 729.3 | 724.2 | 724.0 | 723.9 | 724.0 | 723.9 | 723.8 |

ν4 | 901.9 | 900.8 | 901.2 | 900.7 | 901.1 | 900.7 | 900.7 | 900.7 | 900.7 | 900.7 | 900.7 |

ν9 | 1036.2 | 1034.4 | 1034.5 | 1034.2 | 1034.5 | 1034.2 | 1034.3 | 1034.1 | 1034.3 | 1034.1 | 1034.1 |

1036.2 | 1034.4 | 1034.5 | 1034.2 | 1034.5 | 1034.2 | 1034.3 | 1034.1 | 1034.3 | 1034.1 | 1034.1 | |

ν11 | 1088.3 | 1086.8 | 1092.0 | 1086.9 | 1092.0 | 1086.9 | 1087.1 | 1086.7 | 1087.1 | 1086.7 | 1086.6 |

3ν11 | 1088.3 | 1086.8 | 1092.0 | 1086.9 | 1092.0 | 1086.9 | 1087.1 | 1086.7 | 1087.1 | 1086.7 | 1086.6 |

3ν11 | 1089.6 | 1088.0 | 1093.3 | 1088.2 | 1093.3 | 1088.2 | 1088.4 | 1087.9 | 1088.3 | 1087.9 | 1087.8 |

1089.6 | 1088.0 | 1093.3 | 1088.2 | 1093.3 | 1088.2 | 1088.4 | 1087.9 | 1088.3 | 1087.9 | 1087.8 | |

ν4 + ν11 | 1261.3 | 1260.0 | 1260.5 | 1259.9 | 1260.5 | 1259.9 | 1260.0 | 1259.9 | 1259.9 | 1259.9 | 1259.9 |

1262.3 | 1260.2 | 1260.7 | 1260.0 | 1260.7 | 1259.9 | 1260.3 | 1259.9 | 1260.2 | 1259.9 | 1259.9 | |

ν3 | 1391.9 | 1389.3 | 1395.1 | 1389.4 | 1395.0 | 1389.3 | 1389.3 | 1389.0 | 1389.3 | 1389.0 | 1389.0 |

ν9 + ν11 | 1398.0 | 1395.2 | 1401.2 | 1395.3 | 1401.1 | 1395.1 | 1395.2 | 1394.8 | 1395.1 | 1394.8 | 1394.7 |

1398.1 | 1395.1 | 1401.2 | 1395.3 | 1401.1 | 1395.1 | 1395.2 | 1394.8 | 1395.1 | 1394.8 | 1394.7 | |

ν9 + ν11 | 1398.3 | 1395.3 | 1401.4 | 1395.5 | 1401.3 | 1395.4 | 1395.4 | 1395.0 | 1395.4 | 1395.0 | 1394.9 |

ν9 + ν11 | 1400.8 | 1398.1 | 1404.0 | 1398.2 | 1403.9 | 1398.1 | 1398.2 | 1397.8 | 1398.1 | 1397.8 | 1397.7 |

MAE | 5.3 | 1.0 | 9.9 | 1.5 | 10.4 | 1.7 | 1.9 | 0.5 | 1.8 | 0.5 | |

RMSE | 3.6 | 0.6 | 6.4 | 0.6 | 6.6 | 0.6 | 1.1 | 0.2 | 1.0 | 0.2 |

cprods
. | 210 . | 420 . | 420 . | 840 . | 840 . | . | |||||
---|---|---|---|---|---|---|---|---|---|---|---|

tprods
. | 10 000 . | 10 000 . | 20 000 . | 20 000 . | 30 000 . | . | |||||

. | Var . | Pert . | Var . | Pert . | Var . | Pert . | Var . | Pert . | Var . | Pert . | Smolyak^{a}
. |

ZPVE | 9844.6 | 9837.9 | 9837.8 | 9837.5 | 9837.7 | 9837.5 | 9837.7 | 9837.4 | 9837.7 | 9837.4 | 9837.4 |

ν11 | 361.1 | 361.0 | 361.1 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 |

362.1 | 361.1 | 361.1 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | 361.0 | |

2ν11 | 723.7 | 723.3 | 723.4 | 723.2 | 723.4 | 723.2 | 723.3 | 723.2 | 723.3 | 723.2 | 723.2 |

724.8 | 723.4 | 728.7 | 723.5 | 728.7 | 723.5 | 723.3 | 723.2 | 723.3 | 723.2 | 723.2 | |

2ν11 | 725.4 | 724.0 | 729.4 | 724.2 | 729.3 | 724.2 | 724.0 | 723.9 | 724.0 | 723.9 | 723.8 |

ν4 | 901.9 | 900.8 | 901.2 | 900.7 | 901.1 | 900.7 | 900.7 | 900.7 | 900.7 | 900.7 | 900.7 |

ν9 | 1036.2 | 1034.4 | 1034.5 | 1034.2 | 1034.5 | 1034.2 | 1034.3 | 1034.1 | 1034.3 | 1034.1 | 1034.1 |

1036.2 | 1034.4 | 1034.5 | 1034.2 | 1034.5 | 1034.2 | 1034.3 | 1034.1 | 1034.3 | 1034.1 | 1034.1 | |

ν11 | 1088.3 | 1086.8 | 1092.0 | 1086.9 | 1092.0 | 1086.9 | 1087.1 | 1086.7 | 1087.1 | 1086.7 | 1086.6 |

3ν11 | 1088.3 | 1086.8 | 1092.0 | 1086.9 | 1092.0 | 1086.9 | 1087.1 | 1086.7 | 1087.1 | 1086.7 | 1086.6 |

3ν11 | 1089.6 | 1088.0 | 1093.3 | 1088.2 | 1093.3 | 1088.2 | 1088.4 | 1087.9 | 1088.3 | 1087.9 | 1087.8 |

1089.6 | 1088.0 | 1093.3 | 1088.2 | 1093.3 | 1088.2 | 1088.4 | 1087.9 | 1088.3 | 1087.9 | 1087.8 | |

ν4 + ν11 | 1261.3 | 1260.0 | 1260.5 | 1259.9 | 1260.5 | 1259.9 | 1260.0 | 1259.9 | 1259.9 | 1259.9 | 1259.9 |

1262.3 | 1260.2 | 1260.7 | 1260.0 | 1260.7 | 1259.9 | 1260.3 | 1259.9 | 1260.2 | 1259.9 | 1259.9 | |

ν3 | 1391.9 | 1389.3 | 1395.1 | 1389.4 | 1395.0 | 1389.3 | 1389.3 | 1389.0 | 1389.3 | 1389.0 | 1389.0 |

ν9 + ν11 | 1398.0 | 1395.2 | 1401.2 | 1395.3 | 1401.1 | 1395.1 | 1395.2 | 1394.8 | 1395.1 | 1394.8 | 1394.7 |

1398.1 | 1395.1 | 1401.2 | 1395.3 | 1401.1 | 1395.1 | 1395.2 | 1394.8 | 1395.1 | 1394.8 | 1394.7 | |

ν9 + ν11 | 1398.3 | 1395.3 | 1401.4 | 1395.5 | 1401.3 | 1395.4 | 1395.4 | 1395.0 | 1395.4 | 1395.0 | 1394.9 |

ν9 + ν11 | 1400.8 | 1398.1 | 1404.0 | 1398.2 | 1403.9 | 1398.1 | 1398.2 | 1397.8 | 1398.1 | 1397.8 | 1397.7 |

MAE | 5.3 | 1.0 | 9.9 | 1.5 | 10.4 | 1.7 | 1.9 | 0.5 | 1.8 | 0.5 | |

RMSE | 3.6 | 0.6 | 6.4 | 0.6 | 6.6 | 0.6 | 1.1 | 0.2 | 1.0 | 0.2 |

^{a}

Reference 24.

### B. Ethylene oxide

The second test system is ethylene oxide, which has 15 vibrational modes and has been previously well studied using other techniques for calculating vibrational excitations.^{7,38,55,57} For these VASCI calculations, the force field of Pouchan and co-workers^{103} was used with the missing force field values from the original study included, as outlined in Ref. 101. The force field has harmonic constants computed at the CCSD(T)/cc-pVTZ level and cubic and quartic anharmonic constants calculated at the B3LYP/6-31+G(d,p) level. To allow a comparison with the previous studies, the 200 vibrational states with lowest energy were calculated. We performed three VASCI calculations with different values of *cprods* and *tprods* and compared to the hierarchical intertwined (HI)-RRBPM^{57} and A-VCI^{55} methods. A subset of these results is presented in Table III. The full set of VASCI results is found in the supplementary material.

cprods
. | 5000 . | 5000 . | 10 000 . | . | . | . | |||
---|---|---|---|---|---|---|---|---|---|

tprods
. | 50 000 . | 100 000 . | 150 000 . | HI-RRBPM^{a}
. | . | ||||

. | Var . | Pert . | Var . | Pert . | Var . | Pert . | A . | D . | A-VCI^{b}
. |

ZPVE | 12 462.3 | 12 461.7 | 12 461.8 | 12 461.6 | 12 461.7 | 12 461.6 | 12 461.9 | 12 461.7 | 12 461.4 |

ν15 | 794.1 | 793.1 | 793.6 | 792.9 | 793.2 | 792.8 | 793.3 | 793.1 | 792.6 |

ν12 | 823.1 | 822.3 | 822.8 | 822.2 | 822.5 | 822.2 | 822.2 | 822.0 | 821.9 |

ν5 | 879.3 | 878.6 | 879.0 | 878.5 | 878.7 | 878.4 | 878.5 | 878.3 | 878.3 |

ν8 | 1018.9 | 1017.7 | 1018.4 | 1017.5 | 1017.9 | 1017.4 | 1017.8 | 1017.5 | 1017.1 |

ν4 | 1122.4 | 1121.1 | 1121.8 | 1120.9 | 1121.3 | 1120.8 | 1122.3 | 1121.9 | 1121.2 |

ν11 | 1125.3 | 1124.0 | 1124.8 | 1123.8 | 1124.2 | 1123.8 | 1124.6 | 1124.4 | 1123.6 |

ν14 | 1147.6 | 1146.2 | 1147.0 | 1146.0 | 1146.4 | 1146.0 | 1146.4 | 1146.2 | 1145.7 |

ν7 | 1149.9 | 1148.6 | 1149.3 | 1148.4 | 1148.8 | 1148.3 | 1148.6 | 1148.4 | 1148.0 |

ν3 | 1272.7 | 1271.6 | 1272.3 | 1271.4 | 1271.8 | 1271.3 | 1271.3 | 1270.9 | 1270.8 |

ν10 | 1469.6 | 1467.9 | 1468.9 | 1467.6 | 1468.1 | 1467.5 | 1468.1 | 1467.7 | 1467.3 |

ν5 + ν7 + ν14 | 3174.3 | 3166.8 | 3172.7 | 3165.8 | 3169.5 | 3165.2 | 3199.0 | 3169.4 | 3163.5 |

ν5 + 2ν7 | 3176.2 | 3169.6 | 3175.1 | 3168.8 | 3172.4 | 3168.3 | 3200.5 | 3171.5 | 3166.6 |

ν3 + ν11 + ν15 | 3182.6 | 3173.1 | 3179.4 | 3171.4 | 3176.9 | 3170.9 | 3202.7 | 3175.5 | 3168.0 |

ν3 + ν4 + ν15 | 3185.9 | 3176.2 | 3182.8 | 3174.5 | 3179.9 | 3173.9 | 3204.7 | 3179.5 | 3171.5 |

ν7 + 2ν8 | 3188.7 | 3180.1 | 3186.5 | 3178.7 | 3183.4 | 3178.0 | 3211.7 | 3182.1 | 3175.9 |

ν2 + ν5 + ν12 | 3192.6 | 3184.1 | 3189.8 | 3182.8 | 3186.8 | 3182.2 | 3213.1 | 3189.1 | 3180.2 |

ν12 + 3ν15 | 3207.9 | 3193.8 | 3203.5 | 3191.8 | 3197.6 | 3190.7 | N/A | N/A | 3184.9 |

ν3 + ν7 + ν15 | 3204.8 | 3195.0 | 3201.5 | 3193.5 | 3198.3 | 3192.7 | 3216.7 | 3199.0 | 3189.7 |

ν3 + ν11 + ν12 | 3211.1 | 3202.8 | 3208.7 | 3201.6 | 3205.8 | 3200.9 | 3218.9 | 3205.7 | 3199.2 |

ν3 + ν14 + ν15 | 3216.3 | 3208.5 | 3214.2 | 3207.3 | 3211.6 | 3206.7 | 3226.4 | 3209.2 | 3204.4 |

Core hours | 10.4 | 30.4 | 67.1 | 384.0 | 6681.6 | 1756.4 | |||

Cores | 1 | 1 | 2 | 32 | 32 | 24 |

cprods
. | 5000 . | 5000 . | 10 000 . | . | . | . | |||
---|---|---|---|---|---|---|---|---|---|

tprods
. | 50 000 . | 100 000 . | 150 000 . | HI-RRBPM^{a}
. | . | ||||

. | Var . | Pert . | Var . | Pert . | Var . | Pert . | A . | D . | A-VCI^{b}
. |

ZPVE | 12 462.3 | 12 461.7 | 12 461.8 | 12 461.6 | 12 461.7 | 12 461.6 | 12 461.9 | 12 461.7 | 12 461.4 |

ν15 | 794.1 | 793.1 | 793.6 | 792.9 | 793.2 | 792.8 | 793.3 | 793.1 | 792.6 |

ν12 | 823.1 | 822.3 | 822.8 | 822.2 | 822.5 | 822.2 | 822.2 | 822.0 | 821.9 |

ν5 | 879.3 | 878.6 | 879.0 | 878.5 | 878.7 | 878.4 | 878.5 | 878.3 | 878.3 |

ν8 | 1018.9 | 1017.7 | 1018.4 | 1017.5 | 1017.9 | 1017.4 | 1017.8 | 1017.5 | 1017.1 |

ν4 | 1122.4 | 1121.1 | 1121.8 | 1120.9 | 1121.3 | 1120.8 | 1122.3 | 1121.9 | 1121.2 |

ν11 | 1125.3 | 1124.0 | 1124.8 | 1123.8 | 1124.2 | 1123.8 | 1124.6 | 1124.4 | 1123.6 |

ν14 | 1147.6 | 1146.2 | 1147.0 | 1146.0 | 1146.4 | 1146.0 | 1146.4 | 1146.2 | 1145.7 |

ν7 | 1149.9 | 1148.6 | 1149.3 | 1148.4 | 1148.8 | 1148.3 | 1148.6 | 1148.4 | 1148.0 |

ν3 | 1272.7 | 1271.6 | 1272.3 | 1271.4 | 1271.8 | 1271.3 | 1271.3 | 1270.9 | 1270.8 |

ν10 | 1469.6 | 1467.9 | 1468.9 | 1467.6 | 1468.1 | 1467.5 | 1468.1 | 1467.7 | 1467.3 |

ν5 + ν7 + ν14 | 3174.3 | 3166.8 | 3172.7 | 3165.8 | 3169.5 | 3165.2 | 3199.0 | 3169.4 | 3163.5 |

ν5 + 2ν7 | 3176.2 | 3169.6 | 3175.1 | 3168.8 | 3172.4 | 3168.3 | 3200.5 | 3171.5 | 3166.6 |

ν3 + ν11 + ν15 | 3182.6 | 3173.1 | 3179.4 | 3171.4 | 3176.9 | 3170.9 | 3202.7 | 3175.5 | 3168.0 |

ν3 + ν4 + ν15 | 3185.9 | 3176.2 | 3182.8 | 3174.5 | 3179.9 | 3173.9 | 3204.7 | 3179.5 | 3171.5 |

ν7 + 2ν8 | 3188.7 | 3180.1 | 3186.5 | 3178.7 | 3183.4 | 3178.0 | 3211.7 | 3182.1 | 3175.9 |

ν2 + ν5 + ν12 | 3192.6 | 3184.1 | 3189.8 | 3182.8 | 3186.8 | 3182.2 | 3213.1 | 3189.1 | 3180.2 |

ν12 + 3ν15 | 3207.9 | 3193.8 | 3203.5 | 3191.8 | 3197.6 | 3190.7 | N/A | N/A | 3184.9 |

ν3 + ν7 + ν15 | 3204.8 | 3195.0 | 3201.5 | 3193.5 | 3198.3 | 3192.7 | 3216.7 | 3199.0 | 3189.7 |

ν3 + ν11 + ν12 | 3211.1 | 3202.8 | 3208.7 | 3201.6 | 3205.8 | 3200.9 | 3218.9 | 3205.7 | 3199.2 |

ν3 + ν14 + ν15 | 3216.3 | 3208.5 | 3214.2 | 3207.3 | 3211.6 | 3206.7 | 3226.4 | 3209.2 | 3204.4 |

Core hours | 10.4 | 30.4 | 67.1 | 384.0 | 6681.6 | 1756.4 | |||

Cores | 1 | 1 | 2 | 32 | 32 | 24 |

For the ten vibrational excitations of lowest energy, the VASCI, HI-RRBPM, and A-VCI methods all perform similarly. For all states, the energy of the A-VCI method is the lowest of the three methods. As the A-VCI method is variational, it will be considered the reference value. For the calculation with the smallest number of *cprods* and *tprods*, VASCI performs slightly worse than the most accurate HI-RRBPM calculation, labeled “D” in Table III. The VASCI calculation with the largest number of *cprods* and *tprods* performs slightly better than the most accurate HI-RRBPM calculation, while the remaining VASCI calculation does better for certain states and worse for other relative to the HI-RRBPM method.

For the ten highest energy vibrational excitations found by both the VASCI and A-VCI methods, the VASCI, HI-RRBPM, and A-VCI methods show more discrepancy, which is consistent with the generally observed trend in VCI related methods that states of high energy are more difficult to converge than states of low energy. Once again, for all states, the energy of the A-VCI method is the lowest of the three methods and will be considered the reference value for the other calculations. For these ten states, all VASCI calculations give more accurate results than the most accurate HI-RRBPM calculation. For the VASCI calculation with largest values of *cprods* and *tprods*, the variational vibrational excitation energies with no perturbative energy correction are of similar accuracy to that of the HI-RRBPM method. Even with the superior accuracy of the VASCI method relative to the HI-RRBPM method for states of higher energy, the VASCI method does perform worse for these states than those with lower energy with errors typically in the range of 2–3 cm^{−1} for the most accurate VASCI calculation.

The good performance of the VASCI method is more significant when the computational timings of the methods are compared. The most accurate and computationally expensive VASCI calculation took 67.1 core hours on two Intel Xeon Gold 5138 2.0 GHz CPUs, which is about a factor of 100 less than the most accurate and computationally expensive HI-RRBPM calculation, which took 6681.6 core hours on 32 AMD Opteron 6386 SE 2.8 GHz CPUs. The most computationally expensive VASCI calculation is over a factor of 25 less computationally expensive than the A-VCI method, which took 1756.4 core hours on 24 Haswell Intel Xeon E5-2680 2.8 GHz CPUs. Given the different computer architectures employed for these calculations, any comparison of timing information should be used cautiously, but even with that caution, the computational efficiency of the VASCI method appears to be high. When combined with the similar accuracy of the VASCI method relative to the HI-RRBPM and A-VCI, these results demonstrate the potential utility of the VASCI method for calculating vibrational excitations.

### C. Naphthalene

With the demonstration that the VASCI method can accurately calculate the vibrational excitations of molecules with up to 15 normal modes, we choose naphthalene as the final test system, which, along with other polycyclic aromatic hydrocarbons, has attracted significant recent interest as the infrared bands from these molecules are responsible for some of the strongest emission spectra in the observed interstellar medium.^{104} Naphthalene has 48 vibrational modes and represents a significantly tougher challenge for calculating accurate vibrational excitation energies than the two previous systems. Additionally, compared to the two previous systems, the vibrational excitations of naphthalene have been less well studied with the HI-RRBPM method only recently being applied to the system.^{56}

Following the previous HI-RRBPM study as a guide, we use the quartic force field of Trombetti and co-workers,^{105} for which all force constants were calculated at the B97-1/TZ2P level of theory. For all the VASCI calculations, the 128 vibrational states with lowest energy were found. At present, our VASCI implementation is unable to perform the perturbative energy correction due to our limited current computer resources, so we present only variational energies. As mentioned previously, our current VASCI implementation is nonoptimized and there exist recent improvements to the ASCI method or other related methods that could be leveraged to allow the perturbative energy correction to be performed on naphthalene. This will be the subject of future research. Even without the perturbative energy correction, the accuracy of the VASCI method is competitive with the HI-RRBPM technique for the naphthalene system.

A subset of the 128 calculated vibrational excitation energies is presented in Table IV. The full set is found in the supplementary material. For the ten excited states with smallest energy, the HI-RRBPM technique performs better than the VASCI method. For these modes, the accuracy of the VASCI method typically lies in between that of the HI-RRBPM calculations labeled “A” and “B,” but in some cases even the “A” results, which are the least accurate of the HI-RRBPM results, are more accurate than the best VASCI results. For these modes, the use of the VASCI method does not offer an increase in accuracy for a given number of core hours as the computational expense of the VASCI calculations also lies between that of the “A” and “B” HI-RRBPM calculations.

cprods
. | 15 000 . | 20 000 . | 25 000 . | HI-RRBPM^{a}
. | |||
---|---|---|---|---|---|---|---|

tprods
. | 1 500 000 . | 1 200 000 . | 1 000 000 . | A . | B . | D . | G . |

ZPVE | 31 773.6 | 31 774.3 | 31 774.4 | 31 782.2 | 31 768.7 | 31 767.1 | 31 766.0 |

ν48 | 166.3 | 166.3 | 166.4 | 165.8 | 165.8 | 165.1 | 164.6 |

ν13 | 179.9 | 179.7 | 179.9 | 184.9 | 179.2 | 178.6 | 178.2 |

2ν48 | 336.1 | 336.3 | 336.4 | 338.2 | 332.1 | 330.4 | 329.4 |

ν13 + ν48 | 349.4 | 349.4 | 349.5 | 365.7 | 345.2 | 343.1 | 342.0 |

ν24 | 363.8 | 361.6 | 361.4 | 372.9 | 357.9 | 355.2 | 354.4 |

2ν13 | 366.6 | 363.4 | 363.4 | 397.3 | 361.1 | 358.9 | 357.7 |

ν16 | 396.3 | 394.2 | 394.1 | 405.4 | 388.9 | 388.0 | 387.7 |

ν28 | 475.3 | 471.5 | 470.9 | 468.2 | 464.9 | 464.0 | 463.5 |

ν47 | 483.1 | 481.7 | 479.7 | 477.1 | 473.8 | 472.9 | 472.4 |

3ν48 | 501.8 | 502.7 | 502.5 | 506.6 | 500.3 | 497.1 | 495.5 |

ν13 + ν23 | 1 007.1 | 996.1 | 991.6 | 1 091.2 | 1010.7 | 990.3 | 975.0 |

ν28 + ν44 | 1 008.1 | 999.9 | 991.3 | 1 094.6 | 1 012.7 | 993.0 | 978.3 |

ν9 + ν28 | 1 010.7 | 1 004.1 | 997.7 | 1 098.7 | 1 017.7 | 996.9 | 981.3 |

ν24 + ν36 | 1 015.4 | 1 010.3 | 1 006.6 | 1 099.9 | 1 019.6 | 999.3 | 982.9 |

ν12 + ν24 | 1 015.7 | 1 007.3 | 1 002.2 | 1 100.6 | 1 024.0 | 1 005.5 | 985.1 |

v44 + ν47 | 1 013.3 | 1 006.8 | 1 001.8 | 1 106.1 | 1 027.3 | 1 006.3 | 987.9 |

ν12 + ν13 + ν48 | 996.9 | 990.5 | 989.3 | 1 114.3 | 1 028.0 | 1 007.5 | 988.1 |

2ν13 + ν36 | 1 008.4 | 1 002.7 | 1 000.9 | 1 116.8 | 1 030.1 | 1 010.4 | 989.8 |

ν24 + ν28 + ν48 | 1 020.6 | 1 010.0 | 1 004.3 | 1 120.5 | 1 033.0 | 1 010.8 | 992.3 |

ν12 + 2ν13 | 1 007.1 | 1 001.5 | 999.7 | 1 128.6 | 1 036.6 | 1 017.6 | 998.7 |

Core hours | 1 834.9 | 1 866.9 | 1 584.7 | 1 167.4 | 6 451.2 | 49 971.2 | 63 590.4 |

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

cprods
. | 15 000 . | 20 000 . | 25 000 . | HI-RRBPM^{a}
. | |||
---|---|---|---|---|---|---|---|

tprods
. | 1 500 000 . | 1 200 000 . | 1 000 000 . | A . | B . | D . | G . |

ZPVE | 31 773.6 | 31 774.3 | 31 774.4 | 31 782.2 | 31 768.7 | 31 767.1 | 31 766.0 |

ν48 | 166.3 | 166.3 | 166.4 | 165.8 | 165.8 | 165.1 | 164.6 |

ν13 | 179.9 | 179.7 | 179.9 | 184.9 | 179.2 | 178.6 | 178.2 |

2ν48 | 336.1 | 336.3 | 336.4 | 338.2 | 332.1 | 330.4 | 329.4 |

ν13 + ν48 | 349.4 | 349.4 | 349.5 | 365.7 | 345.2 | 343.1 | 342.0 |

ν24 | 363.8 | 361.6 | 361.4 | 372.9 | 357.9 | 355.2 | 354.4 |

2ν13 | 366.6 | 363.4 | 363.4 | 397.3 | 361.1 | 358.9 | 357.7 |

ν16 | 396.3 | 394.2 | 394.1 | 405.4 | 388.9 | 388.0 | 387.7 |

ν28 | 475.3 | 471.5 | 470.9 | 468.2 | 464.9 | 464.0 | 463.5 |

ν47 | 483.1 | 481.7 | 479.7 | 477.1 | 473.8 | 472.9 | 472.4 |

3ν48 | 501.8 | 502.7 | 502.5 | 506.6 | 500.3 | 497.1 | 495.5 |

ν13 + ν23 | 1 007.1 | 996.1 | 991.6 | 1 091.2 | 1010.7 | 990.3 | 975.0 |

ν28 + ν44 | 1 008.1 | 999.9 | 991.3 | 1 094.6 | 1 012.7 | 993.0 | 978.3 |

ν9 + ν28 | 1 010.7 | 1 004.1 | 997.7 | 1 098.7 | 1 017.7 | 996.9 | 981.3 |

ν24 + ν36 | 1 015.4 | 1 010.3 | 1 006.6 | 1 099.9 | 1 019.6 | 999.3 | 982.9 |

ν12 + ν24 | 1 015.7 | 1 007.3 | 1 002.2 | 1 100.6 | 1 024.0 | 1 005.5 | 985.1 |

v44 + ν47 | 1 013.3 | 1 006.8 | 1 001.8 | 1 106.1 | 1 027.3 | 1 006.3 | 987.9 |

ν12 + ν13 + ν48 | 996.9 | 990.5 | 989.3 | 1 114.3 | 1 028.0 | 1 007.5 | 988.1 |

2ν13 + ν36 | 1 008.4 | 1 002.7 | 1 000.9 | 1 116.8 | 1 030.1 | 1 010.4 | 989.8 |

ν24 + ν28 + ν48 | 1 020.6 | 1 010.0 | 1 004.3 | 1 120.5 | 1 033.0 | 1 010.8 | 992.3 |

ν12 + 2ν13 | 1 007.1 | 1 001.5 | 999.7 | 1 128.6 | 1 036.6 | 1 017.6 | 998.7 |

Core hours | 1 834.9 | 1 866.9 | 1 584.7 | 1 167.4 | 6 451.2 | 49 971.2 | 63 590.4 |

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

^{a}

Reference 56.

For the ten excited states of the calculated 128 with highest energy, the performance of the VASCI method relative to the HI-RRBPM technique is better than for the lowest energy excited states. The performance of the most accurate VASCI calculation is most similar to that of the HI-RRBPM calculation labeled “D,” but we note that for two of the ten excitation energies shown here, the VASCI calculated frequency is over 15 cm^{−1} lower than “D” calculation. For both of these states, the difference between the VASCI calculation and the most accurate HI-RRBPM calculation, labeled “G,” is less than 1 cm^{−1}. The efficiency of the VASCI method for calculating the states of higher energy is also demonstrated as the most accurate VASCI calculation is a factor of 30 and 40 less computationally expensive than the “D” and “G” HI-RRBPM calculations, respectively. The HI-RRBPM calculations were performed on 2.5 GHz Intel E7-8867 (v3) CPUs. From these results, the VASCI method appears to be a promising technique for calculating the vibrational excitation energies of large systems, especially to obtain the vibrational excitation energy of states in higher energy regions of the vibrational spectrum.

## V. CONCLUSIONS

The VASCI method has been derived and implemented by adapting the ASCI method to calculate vibrational excitation energies. For the three test systems presented here, the VASCI method offers a good balance of accuracy and computational expense compared to the other state of the art methods for calculating vibrational spectra such as the HI-RRBPM and A-VCI methods. This is especially evident for calculations on larger systems such as the naphthalene molecule. ASCI is just one of many recently introduced SCI+PT methods for calculating electronic energies, and it remains to be seen if other methods such as SHCI^{76,84,86} or adaptive configuration interaction^{79,87} may be better suited than ASCI for solving the vibrational Schrödinger equation.

At present, the VASCI code is implemented using Python and Cython. As mentioned previously, the code is not well optimized and significant speed ups are likely to be achieved either by adapting recent algorithmic improvements for calculating the electronic energy with SCI+PT methods or by porting the code to C or C++. As an example, the generation of Hartree products connected to the core or target space is done using a recursive function, which is known in Python to be slow, especially relative to a pure C or C++ approach. These improvements, especially that of a faster perturbative correction, should allow the full VASCI method to be applied to larger systems than in this study.

At the same time, for many of the electronic SCI+PT algorithmic improvements, it is not immediately obvious how to adapt them as they rely on using a bit string to describe the occupancy of the spin orbitals as the occupancy of each spin orbital is zero or one. Each vibrational degree of freedom can have multiple quanta in it leading to a nonbinary occupancy, which makes using a bit string to describe the occupancy of a vibrational product of normal modes in the same manner as the electronic SCI+PT algorithms difficult. These issues are left to be the subject of future studies. Even without any improvements, the present implementation of the VASCI method is a promising technique for calculating the vibrational excitations of medium-to-large sized molecules.

## SUPPLEMENTARY MATERIAL

See the supplementary material for all VASCI vibrational excitation energies for the acetonitrile, ethylene oxide, and naphthalene systems.

## ACKNOWLEDGMENTS

K.R.B. thanks the University of Missouri for startup funding. The computations for this work were performed on the high performance computing infrastructure provided by Research Computing Support Services and in part by the National Science Foundation under Grant No. CNS-1429294 at the University of Missouri, Columbia, MO.