Hybrid variation-perturbation method for calculating rovibrational energy levels of polyatomic molecules
A procedure for calculation of rotation-vibration states of medium
sized molecules is presented. It combines the advantages of
variational calculations and perturbation theory. The vibrational
problem is solved by diagonalizing a Hamiltonian matrix, which is
partitioned into two sub-blocks. The first, smaller sub-block
includes matrix elements with the largest contribution to the energy
levels targeted in the calculations. The second, larger sub-block
comprises those basis states which have little effect on these
energy levels. Numerical perturbation theory, implemented as a
Jacobi rotation, is used to compute the contributions from the
matrix elements of the second sub-block. Only the first sub-block
needs to be stored in memory and diagonalized. Calculations of the
vibrational-rotational energy levels also employ a partitioning of
the Hamiltonian matrix into sub-blocks, each of which corresponds
either to a single vibrational state or a set of resonating
vibrational states, with all associated rotational levels.
Physically, this partitioning is efficient when the Coriolis
coupling between different vibrational states is small. Numerical
perturbation theory is used to include the cross-contributions from
different vibrational states. Separate individual sub-blocks are
then diagonalized, replacing the diagonalization of a large
Hamiltonian matrix with a number of small matrix diagonalizations.
Numerical examples show that the proposed hybrid
variational-perturbation method greatly speeds up the variational
procedure without significant loss of precision for both
vibrational-rotational energy levels and transition intensities. The
hybrid scheme can be used for accurate nuclear motion calculations on
molecules with up to 15 atoms on currently available
olecular rotation; vibration; variational; perturbation theory; infra red spectra; nuclear motion
Acknowledgment This work was supported by the ERC under Advanced Investigator Project 267219.
10.1080/0026897YYxxxxxxxx \jvol00 \jnum00 \jyear2014
Approaches used to compute vibration-rotation energy levels and wave functions include methods based on perturbation theory, effective Hamiltonians and on the use of the variational principle. Each of these methods has its advantages and disadvantages when used for larger molecules such as those with more than ten atoms.
Historically, perturbation theory was the first method employed to treat the many-body anharmonic problem [1, 2, 3, 4, 5, 6, 7, 8, 9], see also the review by Klein . Second-order perturbation theory is the most common analytic treatment for estimating the vibrational energy levels and including contributions from terms in the Hamiltonian beyond the harmonic approximation (see, for example, Refs. [4, 11, 12, 13]) as well as infrared intensities . However, there are drawbacks inherent in this method that prevent extensive use of it in practice. In particular the results obtained depend on the specific form of the Hamiltonian and the distribution of degenerate oscillators. Polyatomic molecules often show quasi-degeneracies between vibrational energy levels and use of second-order perturbation theory can result in significant errors. Furthermore, in practice anharmonic effects in polyatomic molecules are sufficiently large that they are often not converged with a second-order treatment. As a result perturbation theory calculations usually overestimate the anharmonic corrections, sometimes by a factor of two; see the final columns of Tables 3 and 4 given below.
Variational methods for the calculation of anharmonic energy levels were developed independently by a number of authors [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. Early implementations in general computer codes focused on the use of basis functions for triatomic molecules [34, 35, 36, 37]. More recent developments have involved the increasing use of the discrete variable representation [38, 39, 40] and the extension of the work to polyatomic molecules [41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. The main advantage of this method is that it allows the almost exact calculations of the vibrational-rotational energy levels and wave functions for a given anharmonic potential . However, variational methods work much better for few-atom systems since the size of the basis set grows factorially with the number of vibrational degrees of freedom in the molecule; hence the size of the Hamiltonian matrices, which must be computed and diagonalized, also grows very rapidly. For this reason, variational calculations for polyatomic molecules are difficult even on modern computers and therefore are not routinely used for molecules with more than six atoms, when only small basis sets can be afforded. The computational errors associated with incomplete basis sets are normally significantly larger than the errors from perturbation approaches; this point will be discussed further in the following section.
The treatment of ro-vibration states adds extra complexity to the calculation. In this case a rotational basis set, usually described with analytical rigid-rotor functions, must be introduced. For high rotational quantum number, , this leads to large Hamiltonian matrices even for small polyatomic molecules. However a two-step variational procedure can be used to mitigate the effects of this problem  meaning that it has long been possible to compute rotational excitation up to dissociation for small molecules [60, 61].
There are a number of modifications of the variational method which facilitates (ro-)vibrational calculations on larger molecules. For example, the use of vibrational self-consistent field theory by Gerber, Ratner and others [62, 63, 64]. Bowman, Carter and Handy [65, 66] used vibrational configuration interaction (VCI) to reduce the size of the diagonalization for large molecules. If the molecule has separable degrees of freedom, a significant reduction in the time taken to diagonalize the Hamiltonian matrices can be obtained. Scribano and Benoit [67, 68] used a hybrid approach where a modification of the VCI method take into account the interactions of individual configurations using perturbation theory. Similarly the general code MULTIMODE  allows such treatments of larger systems  with approximations involving the degree of coupling between vibrational modes.
Here we propose a new hybrid variational-perturbation theory method based on a physically reasonable division of the large, full, variational Hamiltonian matrix into weakly interacting sub-blocks. Second-order perturbation theory is used to include cross-interaction effects between these sub-blocks, which are then diagonalized separately. This allows one to replace a single diagonalization of a large Hamiltonian matrix by a series of diagonalizations of much smaller sub-matrices. We show that our hybrid method can greatly accelerate the variational procedure by eliminating diagonalization of large matrices without significant loss of precision in the computed vibrational-rotational energy levels and the intensity of transitions between them. This hybrid scheme is able to perform calculations for large molecules containing up to 15 atoms on currently available computers. Finally we note that a different version of the hybrid variational-perturbation approach has recently been proposed by Fabri et al  and implemented in the general Eckart-Watson Hamiltonian code DEWE . Their hybrid scheme uses a traditional single-state ro-vibrational perturbation method based on the variationally computed vibrational () eigenfunctions as a zero-order solution, where the vibrational problem is solved variationally and the ro-vibrational energies are derived using the second order perturbation theory expressions. This is different from our two-step hybrid scheme, where the perturbation method is an integral part of calculations. We believe that this is the key to extending variational calculations to larger molecules.
The method we propose is not dependent on the precise form of the Hamiltonian. Instead the requirement of the Hamiltonian is that it should result in a diagonally dominant matrix representation. For semi-rigid molecules this criterion is naturally satisfied the Eckart-Watson Hamiltonian  as well as number of related Hamiltonians [74, 75, 76]. However the method is unlikely to perform so well for Hamiltonians whose emphasis is not on making producing diagonally dominant matrices such those expressed in polyspherical coordinates . In this work, all calculations are performed using ANGMOL , a variational program for calculating ro-vibrational spectra of general polyatomic molecules. ANGMOL uses a Hamiltonian expressed in curvilinear internal coordinates and an Eckart embedding which is discussed briefly in the Appendix and extensively elsewhere .
2 A variational method for vibrational energy levels
To start we consider the nature of the error arising from use of the variational method for calculating vibrational () energy levels of polyatomic molecules. In the variational Rayleigh-Ritz procedure, the wave function, , is generally approximated by a finite expansion
in terms of some appropriate basis functions . Mathematically, this procedure is equivalent to finding the eigenvalues and eigenvectors of the Hamiltonian matrix whose elements are given by:
Eigenvalues of this matrix give values of the vibrational energy levels, , and the eigenvectors, , give the wave function, , when combined with the basis functions. For ease of use and better convergence of the basis functions, the are generally chosen to form a complete orthonormal set. Exact values of the energy levels and eigenfunctions are only guaranteed with an infinite number of basis functions.
The accuracy of the calculated energies depend on the number of basis functions used. As the basis set is extended, the calculated energy levels tend (mostly smoothly and monotonically) to their theoretical limit. The variational limit for a group of (lowest) energy levels is reached when the number of basis functions is sufficient for the chosen energy levels to differ from the exact value by less than some prescribed error. The physical meaning of the variational limit is that basis functions representing highly excited vibrational levels make only a small contribution to the eigenfunctions of the low-lying levels.
The number of basis functions needed to reach the variational limit depends on the number of energy levels to be calculated. In general, the greater the number of energy levels calculated, the larger the number of basis functions required. However, a good choice of basis functions also speeds the convergence: the closer the basis functions are to desired eigenfunctions, the fewer of them are required to achieve the variational limit. Mathematically, this means that the off-diagonal matrix elements () become small, the eigenvalues become close to the corresponding diagonal element (), and the corresponding eigenvectors tend towards being unit vectors.
Details of our method are give in the Appendix. We start by defining the vibrational basis as a product of one-mode functions
where is a either a Morse or a harmonic oscillator function, is the excitation number of the oscillator and is the number of vibrational degrees of freedom. Morse oscillator functions are commonly used for the X–H stretches and other stretching motion whose vibrations are strongly anharmonic, such as doubly-bonded CO. Harmonic oscillator functions are used for skeletal bonds in molecules with small anharmonicity, for example the CC bonds in hydrocarbons, and for all deformation (angular) oscillations. This form of the basis set helps to reduce the off-diagonal elements in the Hamiltonian matrix and thus has been found to give rapid convergence to the variational limit.
In variational procedures it is often considered desirable to use a complete, orthonormal basis. The bound states of the Morse oscillator form a variational basis of finite length basis and so are strictly speaking not longer complete. It is possible to develop a complete set based on the Morse oscillators ; however, our practical computation shows that rapid convergence can be achieved with a finite number of Morse oscillators allowing the variational limit to be reached.
Our implementation relies on the particular structure of the Hamiltonian matrix ordered by increasing polyad (total vibrational excitation) number,
where is some integer weighting which is often roughly proportional to the inverse of the frequency . For simplicity in this work we use for all . This gives the size of the basis set, in terms of the maximum polyad number, ,
which increases factorially with the number of vibrational degrees of freedom, .
With this definition of the basis set, the Hamiltonian matrix has a semi-banded structure, as shown in Fig. 1, and can be factorized into two regions: the region around the diagonal (Region A), containing large off-diagonal elements, and the outer region (Region B) with ‘small’ or zero off-diagonal elements. The width of the band, as discussed in the Appendix, depends on the form of the Hamiltonian and obviously on the precise definition of ‘small’. For example, with a potential function represented as a fourth-order polynomial in terms of the coordinate displacements and harmonic oscillator basis functions, the band only includes matrix elements with excitations whose basis functions are separated from the diagonal by . Use of the Morse oscillator basis functions gives a similar picture. Clearly distinguishable in Fig. 1 are eight off-diagonal bands corresponding to elements connected by , , , and . Thus for a fourth-order, anharmonic potential function, all matrix elements separated by lie in Region B which becomes increasingly large as either or increases. This rule has similarities with the well-known Slater’s rules for configuration interaction matrix elements in electronic structure calculations.
As an example we use the program ANGMOL to compute the vibrational band origins and corresponding band intensities for HO and HNO. In these calculations we employ semi-empirical potential energy functions represented as fourth-order polynomials which were obtained by fitting the corresponding potential parameters to experimentally derived energies [78, 79]. These empirical potential energy functions reproduce the experimental vibrational term values of the fundamental and first overtone states of HNO and HO with the root-means-squares (rms) errors less then 0.1 and 0.4 cm, respectively. Initial values of the potential parameters were obtained using ab initio calculations at the CCSD(T)/aug-cc-pVQZ level of theory. The corresponding dipole moments were represented in the form of a second-order polynomial fitted to the CCSD(T)/aug-cc-pVQZ ab initio values. Full details will be given elsewhere .
All calculations were performed using curvilinear vibrational coordinates, see Appendix. To represent the stretching coordinates we use Morse coordinates , where in the bond length displacements of the bond and is the standard Morse parameter. Valence angles were represented by , where and are the inter-bond angle and its equilibrium value, respectively. Basis functions were constructed as a direct product of the Morse oscillator functions for the stretches and harmonic oscillator functions for the angles. The harmonic part of the Hamiltonian is constructed to be diagonal for the bending part.
This form of variational basis functions has the advantage that the vibrational Hamiltonian matrix has relatively small off-diagonal elements, i.e. it is close to a diagonal form. This is due to (i) the property of the Morse oscillators which give a good description of the stretching modes and (ii) small bending anharmonic terms. Besides, in this basis the Hamiltonian matrix elements all have simple analytical forms and thus can be efficiently evaluated.
Tables 1 and 2 present calculated values of the wavenumbers and intensities of the vibrational transitions for HO and HNO as a function of the number of basis functions used in the variational calculations. In the following, the parameter () will be used to reference the polyad number for energy values targeted in the variational calculations and () to reference the corresponding number of vibrational states. The results for the water molecule suggest that to compute all vibrational levels associated with the polyad number to better than 5 cm, the basis set must include functions with at least up to . For example, to obtain all energies up to the second overtones () with this accuracy, must be at least 7. This is in accord with our assumption that for the potential energy function given as a fourth-order polynomial, the large matrix elements which belong in Region 1 are all associated with functions given by .
Calculating all vibrational term values to an accuracy better than 0.3 cm requires basis functions with . This means that the Hamiltonian matrix must include all the basis functions for which the difference in is up to 8.
For water, using the lowest possible basis, i.e. , gives errors of 28 cm for the fundamental band, 84 cm for the first overtone and 154 cm for , see Table 1. For HNO the situation even worse, with , the band centres of the fundamentals are not reproduced to within 500 cm.
It is important to note that the number of the target energy levels defined by grows rapidly with the number of degrees of freedom . For a triatomic molecule like water ( 3), to reach an accuracy better than 0.5 cm for the fundamental bands only (, ) requires 220 basis functions. For a pentatomic molecule, such as HNO ( 9), this requires 48 620 basis functions and for an 8-atom molecule such as CH ( = 18), is 4 686 825.
The vibrational energies of polyatomic molecules are often characterized by (accidental) resonances, such as Fermi resonances, resulting from strong mixing between vibrationally excited states. One consequence of this is that for HNO is increased to . Comparison of columns 2 of Tables 2 and Tables 4 shows that even when for 48 620 basis functions corresponding to , the convergence of the fundamental energy levels, i.e. , is only good to 4 cm.
For nitric acid the resonances can also give rise to large errors in band intensities for vibrational states involved in these resonances, especially when small basis sets are used. However, the sum of the bands intensities for all resonant states is less sensitive to the size of the basis set. We note, however, that for water even individual vibrational band intensities show relatively small dependence on the basis set size, see Table 1.
Tables 1 and 2 also show that the stretching energy levels converge faster, requiring less basis functions, than the bending modes. This is because the Morse oscillators are more compact than the harmonic oscillator basis functions and also because the stretching quanta of hydrogen bonds is normally about two times larger than that of the bending modes, i.e, it takes fewer stretching quanta to reach the same energy. However, the harmonic oscillator functions used for the bending modes also represent a reasonable choice giving small coupling matrix elements, which satisfy the requirements of perturbation theory.
3 A hybrid method for calculating vibrational energy levels
Variational calculations on many-atom systems rapidly become impractical. To address this problem we implement a mixed variational-perturbation theory approach. The idea is based on the observation that the calculated energies and wave functions of the lower-lying vibrational levels depend differently on different parts of the Hamiltonian matrix. As described above, the general Hamiltonian matrix can be factorised into two regions as illustrated in Fig. 2. The first region contains the zeroth-order contribution to the calculated vibrational energies as well as all other states strongly coupled to them. As noted above, these strongly-coupled levels are those for which the Hamiltonian matrix contains large off-diagonal elements with the states of interest. For a potential function represented as a fourth-order polynomial, the first block should contain couplings to the targeted states with all contributions with .
The second region of the matrix is generally much larger. This region includes all the remaining basis states that only have a small effect on the calculated energy levels in question. These basis states are ones for which the polyad numbers differ by more than 4 from those of the levels of interest, , because the corresponding off-diagonal elements are small.
The essence of our method is that instead of diagonalizing the whole Hamiltonian matrix, we diagonalize the smaller block 1 (see Fig. 2) with the matrix elements corrected by the second order perturbative contributions from Region 2 as described in detail below.
Assuming the size of the overall Hamiltonian matrix is defined by in connection with Eq. (5) and the size of block 1, , is given by
where is the number of the largest polyad included in block 1. To reach an accuracy of 1 cm or better, the following two conditions have to be met: and .
In the first step of our approach all the matrix elements from blocks 1 are computed along with the diagonal matrix elements of block 3. The second step involves computing the off-diagonal elements of block 2 and accounting for their effect on the matrix elements of the block 1 using one Jacobi rotation [81, 82]. Considering a contribution from the block 2 off-diagonal element , which couples the diagonal elements in block 1 and in block 3, see Fig. 2, these two diagonal elements are perturbatively adjusted using the Jacobi formula
Here is an initial (unperturbed) diagonal element and is an adjusted one. This formula automatically allows for cases where the unperturbed vibrational states are (quasi-)degenerate, i.e. . In the absence of the degeneracy, , we have , which is equivalent to the energy correction given by second-order perturbation theory (see e.g. ). Similarly, if there is degeneracy or quasi-degeneracy , we have and for degeneracy , we have . That is is always satisfied.
Strictly speaking, an off-diagonal element in block 2 should be included as perturbative contribution both to the relevant diagonal and off-diagonal elements from block 1. However, the changes to the off-diagonal are of the second-order leading to small contributions and are thus omitted here without significant loss of accuracy, as will be shown below.
In the final step of the algorithm, block 1 of the corrected Hamiltonian matrix is diagonalized. The resulting eigenfunctions are represented as expansions in terms of the basis functions associated with block 1 only, while the eigenvalues include contributions from all three regions of the matrix.
As in the examples above, we use HO and HNO to illustrate the method (see also Ref. , where VPT2 was used to study the IR spectroscopic properties of HNO). Tables 3 and 4 show similar results obtained using our hybrid method, where the vibrational term values calculated with the size of block 1, , increasing. We note that when , our hybrid method automatically becomes a full variational calculation.
The hybrid method provides accurate results for the vibrational energy levels even when the dimension of the final diagonalized block is much smaller than the dimension of the full Hamiltonian matrix . For example, to calculate the vibrational energy levels of water with an average accuracy of 1 cm, it is sufficient to explicitly include the basis functions with . That is, block 1 includes all the larger off-diagonal elements (at least ) and block 2 contains the smaller off-diagonal elements contributing to the target vibrational energies (at least ).
This reduction is especially valuable for larger molecules. For example, the rule for the eight-atomic molecule CH ( = 18) means that to calculate the fundamental band centres to within 1 cm, the dimension of block 1 is only ( 5), while the full size of the Hamiltonian matrix is ( 9). For HNO we have and ( 9).
However in case of a strong Fermi resonance between vibrational states the accuracy condition must be increased at least by 1 unit. For example, in order to obtain all vibrational term values HNO defined by the polyad number to within 4.0 cm, the basis set should include at least contributions, see Table 4.
Even with very small basis sets the hybrid method gives reasonable results. For example, with the minimum possible size defined by the error in the band centres of , , and of HO is 2 cm, 6 cm and 15 cm, respectively, which can be compared to the error of 28 cm, 84 cm and 154 cm, respectively, obtained using the purely variational procedure with the same basis set size. This and other examples collected in Tables 3 and 4 show that our hybrid method inherits the advantage of perturbation theory which performs reasonably well with small basis sets.
In the limit of , our variation-perturbation method turns into the pure perturbation calculation with the perturbation treated numerically rather than analytically. For comparison Tables 3 and 4 also give results obtained from perturbation theory using a single Jacobi rotation for each off-diagonal element.
Tables 3 and 4 also illustrate the performance of the hybrid approach for the vibrational band intensities of these two molecules. These results suggest that for small basis sets the hybrid intensities are slightly better than intensities obtained using a direct variational treatment of the same size. This is due to the reasonable quality achieved, even with small basis sets, of the corresponding energy levels and eigenfunctions, see above.
It should be noted that the hybrid method requires that the off-diagonal elements of the block 2 are small compared to the elements of block 1. This requirement can be satisfied if basis functions are close to the true eigenfunctions. Our tests show that the Morse oscillators used for the stretching coordinates and harmonic oscillators used for the other degrees of freedom satisfy this criterion.
Finally we note that it is not necessary to compute off-diagonal elements in block 3. This means that the total number of the off-diagonal elements calculated is , which is much less than the overall size of the Hamiltonian matrix .
4 Hybrid method for calculating ro-vibration energy levels
Inclusion of rotational motion into the problem adds an extra level of complexity, especially when highly excited rotational states are required. The dimension of the ro-vibrational problem is normally increased by the factor of , where is the total angular momentum quantum number. In the variational approach this leads to matrices which can be several order of magnitude larger than the corresponding pure vibrational () Hamiltonian matrices.
In order to make the vibrational part of the basis set more compact, we construct the ro-vibrational basis set as a direct product of the rigid-rotor Winston’s and a set of selected vibrational eigenfunctions of the problem obtained as described above. This approach is sometimes referenced in the literature as the contraction . Let us use as the number of vibrational eigenstates for this purpose. can be significantly reduced compared to the total number of vibrational states, , from the previous step as defined by the size of block 1. This is because not all of the vibrational states are equally important for the target vibrational bands . Indeed, because our expansion of the kinetic energy operator is truncated at the second order (see Appendix), many matrix elements with are either exactly zero (for the Harmonic oscillators) or very small.
As an example, consider the situation where we would like to compute a spectrum involving all ro-vibrational energies of HNO up to cm. This range corresponds to = 4 or = 715. From our experience of the HNO molecule, the corresponding ro-vibrational basis set has to include excitations at least up to cm to reach the convergence of 1 cm, which means including about = 8000 vibrational levels. Then the corresponding rotational contribution may include up to . This situation is common in the ExoMol project , which computes spectra of hot molecules. This leads to a set of ro-vibrational Hamiltonian matrices with the dimensions ranging from 8000 to = 1 608 000.
In the following we show how to exploit the advantages of the special structure of the ro-vibration Hamiltonian matrix in order to simplify this problem. For large polyatomic molecules, the ro-vibrational Hamiltonian matrix has a pronounced quasi-block character built around vibrational states. This is due to the fact, see Appendix, that the off-diagonal elements between two different vibrational states and always contribute substantially less to the calculated energy levels than the off-diagonal elements within a vibrational state, , where is the quantum number giving the projection of on the body-fixed z-axis () and runs over the vibrational basis set (). Physically, this corresponds to small Coriolis interactions between different vibrational states, which is achieved by the use of the Eckart embedding .
Figure 3 illustrates the structure of the ro-vibration matrix for the first four vibrational states of molecules HNO (). We can see that most of the off-diagonal elements are zero when reflecting the quadratic differential form of the kinetic energy operator. It can also be seen that the non-zero off-diagonal elements gather near the main diagonal and the matrix has an almost tridiagonal form. The near-diagonal elements are mostly associated with changes in the effective geometry of the rotating molecule through interactions with the corresponding vibrational state . The off-diagonal elements are associated with Coriolis coupling between vibrational states and .
Importantly, the contribution of the off-diagonal elements corresponding to the coupling between different vibrational states is much smaller than that of the off-diagonal elements belonging to the same vibrational state, because of the large energy separation between vibrational levels. This is illustrated in Fig. 3, where the actual contributions from the off-diagonal matrix elements to the final ro-vibrational energy levels for molecules HNO are shown. Typically contribution from off-diagonal elements corresponding to different vibrational states are very small, less than 0.001 cm.
This feature of the vibration-rotation Hamiltonian matrix is common for large molecules and underpins the ro-vibrational version of our hybrid approach, which is schematically represented in Fig. 4. Again we use the second-order perturbation theory (as defined by the Jacobi rotation) to transform the matrix to a block diagonal form built from rotational sub-matrices corresponding to different vibrational states . Thus the dimension of each degenerate rotational sub-block is only and we only consider sub-matrices that correspond to vibrational states. In the case of strong resonances between different vibrational states, these states can be combined into one, enlarged sub-matrix and treated together, i.e. for an -fold symmetry vibrational degeneracy, the dimension of the sub-block is . For details see the Appendix.
As above, we employ a single Jacobi rotation which we apply to the (complex) ro-vibrational Hermitian matrix. Our calculations show that, in contrast to the pure vibrational perturbation method, the best agreement with the variational solution is achieved when both the diagonal and off-diagonal elements are updated as given by
for the diagonal elements
and for the off-diagonal elements, where
Here and are the initial (unperturbed) matrix and perturbed matrix, respectively; runs from 1 to , and .
Equation (11) is a symmetrised version of the standard formula for the single Jacobi rotation with respect to the indices and . It should be noted that in the limit , Eq. (11) is identical with the corresponding expression for second-order perturbation theory:
The resulting block-diagonal form is then diagonalized for each sub-matrix separately, where indicates either a single vibrational state or a set of strongly interacting vibrational states. Thus our algorithm replaces the diagonalization of a huge  ro-vibrational matrix with a number of diagonalizations of much smaller [ or ] matrices, where is the number of resonance sub-states. It should be noted that in our algorithm, the ro-vibration wave functions, in contrast to those from an exact diagonalization, do not contain any contribution from other vibrational states. This leads to some errors in the subsequent calculation of the transition dipole moments and intensities of the ro-vibrational transitions . However, the relative error is rather small for the stronger and most important transitions because of the relatively small perturbation contribution from the off-diagonal () matrix elements . For other states, where intensity stealing is important , the resonant coupled vibrational states procedure can be used . Finally, we note that our procedure is both very quick and easily parallelized.
5 Transition intensities
Although the calculation of the ro-vibration energy levels and wave functions may be long, it is only an intermediate stage. The final step is the calculation of transition intensities which involves computing matrix elements of the dipole moment functions between the computed ro-vibration wave functions. For hot molecules, where a huge number of transition dipoles must be evaluated, this step becomes very time consuming and even prohibitively expensive .
Thus, even for HNO, the problem of calculating the vibration-rotation spectrum in the range of 0–4500 cm using a standard variational method without approximations is extremely difficult using existing computers. For larger molecules such calculations are currently almost impossible.
The separable form of our hybrid rotation-vibration wave functions, which are obtained in separate diagonalizations for each vibrational state, allows a speed-up, by several orders of magnitude, of the subsequent evaluation of transition intensities. This is because the wave functions are much more compact and represented by a number of independent parts. As we will show below, this approximation in case of HNO and HO, leads to relatively small errors in the overall shape and magnitude of the ro-vibrational spectra computed at room or higher temperatures.
Tables 5 and 6 show the pure rotational energy levels of HO and HNO for several values obtained using three different methods to find eigensolutions of the corresponding Hamiltonian matrices: full diagionalisation (F), the hybrid ro-vibrational approach described above (H) and a rigid rotor calculation represented by separate diagonalisations of the rotational sub-blocks with the off-diagonal couplings with other vibrational states neglected (R).
The results for water, Table 5, are of interest because the molecule exhibits relatively large off-diagonal elements and thus large centrifugal effects.
The difference illustrates the magnitude of the centrifugal and Coriolis contributions, while shows the quality of our hybrid method for different . For example, for the maximum value of is only 0.03 cm and is within 0.005 cm and for , is within 0.6 cm and is within 0.05 cm. At higher these residuals grow significantly, up to 161 and 44 cm, respectively, for .
For a given , both and increase with increasing rotational energy. For example, for for the lowest term value of 741.25 cm (), cm and cm, while for the highest term value 1805.18 cm () we obtain cm and cm.
Thus, we conclude that for a small molecule with large centrifugal effects, like water, the hybrid method does not provide accurate ro-vibrational energies. Indeed, it is well known that the contraction performs poorly for water  because of issues with linear geometry. However, for such systems purely variational calculations do not present a computational problem . For larger molecules it is generally not necessary to deal with issues associated with quasi-linearity.
For larger molecules with smaller rotational constants, centrifugal distortion effects are expected to decrease significantly. Therefore we also expect the hybrid method to perform better for such systems. Tables 6 illustrates this effect for HNO. For example, for for the maximal residuals we obtain cm and cm, which also correspond to the highest term value 1556.46 cm for the manifold (). For the rotational levels in the range , the average error of the hybrid method is within 0.02 cm, i.e. much better than that for water for the same energy range. The range is selected as it represents the states which are important for the room temperature ro-vibration spectrum of HNO. It should be noted that for large molecules with small rotational constants, like HNO, individual absorption lines are often not resolved because of the high density of the lines even at room temperature. For example, the spectrum of HNO contains lines in the region 0 – 4000 cm with intensities above 10 cmmolecules at 600 K. This suggests that calculations with the accuracy within 0.02 cm should be sufficient for most purposes.
As an example, Fig. 5 shows the absorption cross-sections calculated in the region of the and bands of HNO for 298 K using the full variational (‘Full’) calculation as well as the difference with the hybrid calculation. Here a Voigt line profile was used with parameters cm (a half width at half maximum (HWHM) of 0.153 cm), chosen to match spectra from the PNNL database . Further details of these spectra will be given elsewhere . For the central part of the band (), where the rotational lines are strong enough to be resolved, the two methods (full and hybrid) give almost identical results. At the edges of the bands () slight differences in the frequencies and intensities of the individual ro-vibration absorption lines appear. However the rotational structure is poorly resolved due to superposition of a large number of lines. Besides the intensity of the band diminishes rapidly with due to the decreasing population of the highly excited rotational levels, therefore the absolute difference is also negligible.
In this paper two hybrid variation-perturbation methods for computing vibration and ro-vibrational energies for large molecules are proposed. These methods yield significant speed-ups for the computation of ro-vibrational spectra, with only a minor loss of accuracy, at least for semi-rigid polyatomic molecules. We plan to use method to compute extensive line lists and cross sections for large molecules, including nitric acid , as part of the ExoMol project .
As an illustration of the efficiency of the hybrid method, calculations of the vibrational and ro-vibrational energies of a tri- and a penta-atomic molecule are presented. We can use these to project to even larger systems. Assuming that it is currently feasible to diagonalize a matrix of about 1 000 000 by 1 000 000 and using Eq. (5), we can expect to achieve an accuracy of 1 cm for all ro-vibrational term values for polyads for a given potential energy surface for a molecule containing up to 15 atoms using our hybrid procedure. For molecules containing up to 25 atoms, the expected accuracy is about 5 cm.
Our hybrid methods have several advantages:
They are relatively easy to implement as an extension to a pure variational computer program, which readily provides all required components. For example, all the calculations presented here were performed using ANGMOL , a variational program for ro-vibrational spectra of general polyatomic molecules. It required changes to less than 1% of the entire program to implement the hybrid methods.
They allow efficient parallelization and vectorization of the code.
They are very fast compared to full diagonalizaion. Not only is the diagonalization time greatly reduced, but also only a small proportion of the Hamiltonian matrices needs to be computed. Furthermore, the coupling elements can be evaluated on the fly and do not need to be stored in the memory.
Our method adopts several approximations arising from second-order perturbation theory and the rigid rotor model to full diagonalization approaching the variational limit. For the vibration-only problem, this approximation is easily controlled by altering the size of the block 1.
Our method has scope for further improvement. For example:
For the purely vibrational Hamiltonian matrix, one can change not only the diagonal but also the off-diagonal elements in block 1, as we do for the ro-vibrational Hamiltonian matrix. This will increase the accuracy at the expense of increased computer time.
One could use not only one but two or more consecutive Jacobi rotations. Again, this will increase the accuracy at the expense of increased computer time.
Simple perturbation theory corrections generally overestimate the effect of the perturbation. One could reduce this effect by multiplying the computed perturbations by some empirical coefficient . This will increase the accuracy without increasing the computer time. The problem with this approach is that the magnitude of depends on the number of degrees-of-freedom involved, the structure and size of the Hamiltonian matrix and value of perturbation adjustment.
Finally, we note that in our present implementation no advantage is taken of the pronounced polyad structure shown by many molecules. This could be achieved by a suitable choice of the coefficients in Eq. (4).
This work was supported by the ERC under Advanced Investigator Project 267219.
-  J.H. Van Vleck, Phys. Rev. 33, 467 (1929).
-  O.M. Jordahl, Phys. Rev. 45, 87 (1934).
-  W.H. Shaffer, H.H. Nielsen and L.H. Thomas, Phys. Rev. 56, 895 (1939).
-  H.H. Nielsen, Rev. Mod. Phys. 23, 90 (1951).
-  C. Bloch, Nucl. Phys. 6, 329 (1958).
-  C.T. G. Amat, H. H. Nielsen, Rotation-vibration of polyatomic molecules (Dekker, New York, 1971).
-  Jørgensen, F and T. Pedersen, Mol. Phys. 27, 33 (1974).
-  S. Califano, Vibrational States (John Wilay and Sons, New York, London, 1976).
-  V.G. Tyuterev and V.I. Perevalov, Chem. Phys. Lett. 74, 494 (1980).
-  D.J. Klein, J. Chem. Phys. 61, 786 (1974).
-  I.M. Mills, J. Mol. Spectrosc. 5, 334 (1961).
-  I.M. Mills, J. Mol. Spectrosc. 17, 164 (1965).
-  V. Barone, J. Chem. Phys. 122, 014108 (2005).
-  V. Barone, J. Bloino, C.A. Guido and F. Lipparini, Chem. Phys. Lett. 496, 157 (2010).
-  L.A. Gribov, Opt. Spectrosc 31, 456 (1971).
-  I. Suzuky, Bull. Chem Soc. Japan 44, 3277 (1971).
-  M.G. Bucknell, N.C. Handy and S.F. Boys, Mol. Phys. 28, 759 (1974).
-  M.G. Bucknell and N.C. Handy, Mol. Phys. 28, 777 (1974).
-  R.J. Whitehead and N.C. Handy, J. Mol. Spectrosc. 55, 356 (1975).
-  G.D. Carney and C.W. Kern, Intern. J. Quantum Chem. S9, 317 (1975).
-  R.J. Whitehead and N.C. Handy, J. Mol. Spectrosc. 59, 459 (1976).
-  G.D. Carney and R.N. Porter, J. Chem. Phys. 65, 3547 (1976).
-  G.D. Carney, S.R. Langhoff and L.A. Curtiss, J. Chem. Phys. 66, 3724 (1977).
-  G.D. Carney and R.N. Porter, Chem. Phys. Lett. 50, 327 (1977).
-  G.D. Carney, S. Giorgianni and K.N. Rao, J. Mol. Spectrosc. 80, 158 (1980).
-  G.D. Carney and R.N. Porter, Phys. Rev. Lett. 45, 537 (1980).
-  S. Carter and N.C. Handy, Mol. Phys. 47, 1445 (1982).
-  J. Tennyson and B.T. Sutcliffe, J. Chem. Phys. 77, 4061 (1982).
-  G. Brocks and J. Tennyson, J. Mol. Spectrosc. 99, 263 (1983).
-  S. Carter, N. Handy and B.T. Sutcliffe, Mol. Phys. 49, 745 (1983).
-  J. Tennyson and B.T. Sutcliffe, J. Chem. Phys. 79, 43 (1983).
-  D. Cropek and G.D. Carney, J. Chem. Phys. 80, 4280 (1984).
-  J. Tennyson and B.T. Sutcliffe, Mol. Phys. 51, 887 (1984).
-  J. Tennyson, Comput. Phys. Commun. 29, 307 (1983).
-  J. Tennyson, Comput. Phys. Commun. 42, 257 (1986).
-  J. Tennyson and S. Miller, Comput. Phys. Commun. 55, 149 (1989).
-  S. Carter, P. Rosmus, N.C. Handy, S. Miller, J. Tennyson and B.T. Sutcliffe, Comput. Phys. Commun. 55, 71 (1989).
-  J.C. Light, I.P. Hamilton and J.V. Lill, J. Chem. Phys. 82, 1400 (1985).
-  Z. Bačić and J.C. Light, Annu. Rev. Phys. Chem. 40, 469 (1989).
-  J.C. Light and T. Carrington Jr, Adv. Phys. Chem. 114, 263 (2000).
-  I.N. Kozin, M.M. Law, J. Tennyson and J.M. Hutson, Comput. Phys. Commun. 163, 117 (2004).
-  S.N. Yurchenko, W. Thiel and P. Jensen, J. Mol. Spectrosc. 245, 126 (2007).
-  P. Botschwina, Chem. Phys. 40, 33 (1979).
-  P. Botschwina, H. Haertner and W. Sawodny, Chem. Phys. Lett. 74, 156 (1980).
-  P. Botschwina, Chem. Phys. 68, 41 (1982).
-  P. Botschwina, Chem. Phys. 81, 73 (1983).
-  P. Botschwina and P. Sebald, J. Mol. Spectrosc. 100, 1 (1983).
-  Y.I. Ponomarev and A.I. Pavlyuchko, Optika I Spektroskopiya 44, 1262 (1978).
-  G.V. Yukhnevich, E.G. Kokhanova, A.I. Pavlyuchko and V.V. Volkov, J. Molec. Struct. (THEOCHEM) 23, 1 (1985).
-  A.I. Pavlyuchko and L.A. Gribov, Optika I Spektroskopiya 58, 1247 (1985).
-  A.I. Pavlyuchko, Optika I Spektroskopiya 67, 286 (1989).
-  A.I. Pavlyuchko, J. Struct. Chem. 36, 204 (1995).
-  A.I. Pavlyuchko, J. Struct. Chem. 36, 210 (1995).
-  A.I. Pavlyuchko, J. Struct. Chem. 38, 959 (1997).
-  A.A. Vigasin, F. Huisken, A.I. Pavlyuchko, L. Ramonat and E.G. Tarakanova, J. Mol. Spectrosc. 209, 81 (2001).
-  A.A. Vigasin, A.I. Pavlyuchko, Y. Jin and S. Ikawa, J. Mol. Struct. 742, 173 (2005).
-  A.I. Pavlyuchko and L.A. Gribov, J. Struct. Chem. 51, 444 (2010).
-  M. Pavanello, L. Adamowicz, A. Alijah, N.F. Zobov, I.I. Mizus, O.L. Polyansky, J. Tennyson, T. Szidarovszky, A.G. Császár, M. Berg, A. Petrignani and A. Wolf, Phys. Rev. Lett. 108, 023002 (2012).
-  J. Tennyson and B.T. Sutcliffe, Mol. Phys. 58, 1067 (1986).
-  S. Miller and J. Tennyson, Chem. Phys. Lett. 145, 117 (1988).
-  H.Y. Mussa and J. Tennyson, J. Chem. Phys. 109, 10885 (1998).
-  R.B. Gerber and M.A. Ratner, Chem. Phys. Lett. 68, 195 (1979).
-  J.M. Bowman, Acc. Chem. Res. 19, 202 (1986).
-  R.B. Gerber and M.A. Ratner, Adv. Chem. Phys. 70, 97 (1988).
-  S. Carter, J.M. Bowman and N.C. Handy, Theor. Chem. Acc. 100, 191 (1998).
-  N.C. Handy and S. Carter, Mol. Phys. 102, 2201 (2004).
-  Y. Scribano and D.M. Benoit, Chem. Phys. Lett. 458, 384 (2008).
-  Y. Scribano, D.M. Lauvergnat and D.M. Benoit, J. Chem. Phys. 133, 094103 (2010).
-  J.M. Bowman, S. Carter and X.C. Huang, Intern. J. Quantum Chem. 22, 533 (2003).
-  M.L. Senent, P. Palmieri, S. Carter and N.C. Handy, Chem. Phys. Lett. 354, 1 (2002).
-  C. Fabri, T. Furtenbacher and A.G. Császár, Mol. Phys. 112, 2462 (2014).
-  E. Mátyus, G. Czakó, B.T. Sutcliffe and A.G. Császár, J. Chem. Phys. 127, 084102 (2007).
-  J.K.G. Watson, Mol. Phys. 15, 479 (1968).
-  L.A. Gribov and A.I. Pavlyuchko, Variational Methods for Solving Anharmonic Problems in the Theory of Vibrational Spectra of Molecules (Nauka, Moscow, 1998; (in Russian)), (in Russian).
-  G.O. Sørensen, in Topics in Current Chemistry, edited by M. J. S. Dewar et al, Vol. 82 (Springer-Verlag, Dordrecht, 1979), p. 99.
-  S.N. Yurchenko, W. Thiel and P. Jensen, J. Mol. Spectrosc. 245, 126 (2007).
-  F. Gatti and C. Iung, Phys. Rep. 484, 1 (2009).
-  J. Tennyson, P.F. Bernath, L.R. Brown, A. Campargue, M.R. Carleer, A.G. Császár, L. Daumont, R.R. Gamache, J.T. Hodges, O.V. Naumenko, O.L. Polyansky, L.S. Rothmam, A.C. Vandaele, N.F. Zobov, A.R. Al Derzi, C. Fábri, A.Z. Fazliev, T. Furtenbacher, I.E. Gordon, L. Lodi and I.I. Mizus, J. Quant. Spectrosc. Radiat. Transf. 117, 29 (2013).
-  L.S. Rothman, I.E. Gordon, Y. Babikov, A. Barbe, D.C. Benner, P.F. Bernath, M. Birk, L. Bizzocchi, V. Boudon, L.R. Brown, A. Campargue, K. Chance, E.A. Cohen, L.H. Coudert, V.M. Devi, B.J. Drouin, A. Fayt, J.M. Flaud, R.R. Gamache, J.J. Harrison, J.M. Hartmann, C. Hill, J.T. Hodges, D. Jacquemart, A. Jolly, J. Lamouroux, R.J. Le Roy, G. Li, D.A. Long, O.M. Lyulin, C.J. Mackie, S.T. Massie, S. Mikhailenko, H.S.P. Müller, O.V. Naumenko, A.V. Nikitin, J. Orphal, V. Perevalov, A. Perrin, E.R. Polovtseva, C. Richard, M.A.H. Smith, E. Starikova, K. Sung, S. Tashkun, J. Tennyson, G.C. Toon, V.G. Tyuterev and G. Wagner, J. Quant. Spectrosc. Radiat. Transf. 130, 4 (2013).
-  A.I. Pavlyuchko, S.N. Yurchenko and J. Tennyson, J. Chem. Phys. (2014), ( submitted).
-  C.G.J. Jacobi, Crelle’s J. 30, 51 (1846), (in German).
-  J.H. Wilkinson, The Algebraic Eigenvalue Problem (Oxford University Press, Oxford, UK, 1965).
-  K.J. Feierabend, D.K. Havey, M.E. Varner, J.F. Stanton and V. Vaida, J. Chem. Phys. 124, 124323 (2006).
-  S.N. Yurchenko, R.J. Barber, A. Yachmenev, W. Thiel, P. Jensen and J. Tennyson, J. Phys. Chem. A 113, 11845 (2009).
-  J. Tennyson and S.N. Yurchenko, Mon. Not. R. Astron. Soc. 425, 21 (2012).
-  C. Eckart, Phys. Rev. 47, 552 (1935).
-  H.H. Nielsen, Intern. J. Quantum Chem. 1, 217 (1967).
-  S. Miller, J. Tennyson and B.T. Sutcliffe, J. Mol. Spectrosc. 141, 104 (1990).
-  S.N. Yurchenko and J. Tennyson, Mon. Not. R. Astron. Soc. 440, 1649 (2014).
-  J.C. Light, R.M. Whitnell, T.J. Pack and S.E. Choi, in Supercomputer Algorithms for Reactivity, Dynamics and Kinetics of Small Molecules, edited by A. Laganà, NATO ASI series C, Vol. 277 (Kluwer, Dordrecht, 1989), p. 187.
-  R.J. Barber, J. Tennyson, G.J. Harris and R.N. Tolchenov, Mon. Not. R. Astron. Soc. 368, 1087 (2006).
-  S.W. Sharpe, T.J. Johnson, R.L. Sams, P.M. Chu, G.C. Rhoderick and P.A. Johnson, Appl. Spectrosc. 58, 1452 (2004).
-  A.I. Pavlyuchko, A. Yachmenev, J. Tennyson and S.N. Yurchenko 2014, (work in progress).
|0 1 0||1595.0||1595.0||1595.0||1595.0||1595.0||1595.1||1595.2||1595.3||1598.4||1601.6||1626.1||1623.4|
|0 2 0||3152.0||3152.0||3152.1||3152.2||3152.5||3153.3||3153.5||3163.0||3173.4||3211.5||3236.0|
|1 0 0||3656.8||3656.8||3656.8||3656.8||3656.8||3656.9||3656.9||3657.1||3657.9||3659.3||3682.0||3672.4|
|0 0 1||3755.5||3755.5||3755.5||3755.5||3755.5||3755.5||3755.5||3755.6||3756.0||3756.4||3768.3||3767.4|
|0 3 0||4666.8||4667.1||4667.6||4668.9||4671.2||4671.7||4694.7||4716.8||4766.5||4820.7|
|1 1 0||5234.6||5234.6||5234.6||5234.7||5234.8||5235.2||5235.7||5239.7||5247.4||5297.2||5306.8|
|0 1 1||5331.2||5331.2||5331.2||5331.2||5331.3||5331.4||5331.6||5333.8||5337.2||5365.8||5383.7|
|1 2 0||6773.7||6773.9||6774.2||6774.8||6776.2||6777.6||6789.2||6808.7||6883.0||6914.0|
|0 2 1||6870.6||6870.6||6870.7||6870.9||6871.6||6872.0||6878.8||6887.2||6935.1||6962.3|
|2 0 0||7202.0||7202.0||7202.0||7202.0||7202.1||7202.3||7202.6||7203.8||7207.9||7241.0||7244.8|
|1 0 1||7250.1||7250.1||7250.1||7250.1||7250.2||7250.3||7250.5||7251.2||7253.8||7279.2||7290.7|
|0 0 2||7444.7||7444.7||7444.7||7444.7||7444.7||7444.8||7444.9||7445.4||7447.0||7462.7||7488.5|
|3 0 0||10599.3||10599.4||10599.4||10599.5||10599.8||10600.4||10601.9||10611.9||10649.9||10661.9|
|0 0 3||11032.7||11032.7||11032.8||11032.8||11032.9||11033.1||11033.5||11036.1||11054.7||11106.8|
|Transition||=20 , =1771 , /||PT|
|0 1 0||1595.0||1595.0||1595.0||1595.0||1594.9||1595.0||1595.0||1594.8||1594.8||1594.6||1594.7||1593.8||1593.7|
|0 2 0||3152.0||3152.0||3152.0||3151.9||3151.9||3152.1||3151.2||3152.1||3152.8||3147.8||3146.5||3146.6|
|1 0 0||3656.8||3656.8||3656.8||3656.8||3656.8||3656.8||3656.8||3656.8||3656.8||3657.0||3664.5||3660.2||3660.1|
|0 0 1||3755.5||3755.5||3755.5||3755.5||3755.5||3755.5||3755.4||3755.4||3755.1||3754.5||3751.1||3755.3||3756.6|
|0 3 0||4666.8||4666.6||4666.1||4666.5||4667.1||4663.6||4668.1||4672.2||4654.0||4651.7||4656.0|
|1 1 0||5234.6||5234.6||5234.6||5234.6||5234.6||5234.7||5234.6||5234.6||5235.9||5250.4||5248.5||5242.7|
|0 1 1||5331.2||5331.2||5331.2||5331.1||5331.1||5331.1||5331.0||5329.7||5327.8||5320.6||5325.6||5331.9|
|1 2 0||6773.7||6773.6||6773.6||6773.7||6774.1||6773.3||6774.1||6778.9||6794.0||6802.8||6786.5|
|0 2 1||6870.6||6870.6||6870.4||6870.4||6870.4||6869.8||6866.8||6863.8||6850.7||6850.6||6866.3|
|2 0 0||7202.0||7202.0||7202.0||7202.0||7202.0||7202.0||7202.0||7201.8||7202.5||7214.9||7206.2||7202.5|
|1 0 1||7250.1||7250.1||7250.1||7250.0||7250.0||7250.0||7250.0||7249.4||7249.0||7253.2||7252.4||7255.5|
|0 0 2||7444.7||7444.7||7444.7||7444.7||7444.6||7444.6||7444.6||7444.0||7443.0||7438.5||7449.4||7452.2|
|3 0 0||10599.3||10599.3||10599.2||10599.2||10599.2||10599.1||10598.4||10599.8||10613.4||10603.1||10604.9|
|0 0 3||11032.7||11032.7||11032.7||11032.7||11032.6||11032.5||11031.5||11029.8||11022.3||11046.0||11041.2|
|Transition||=14 , =817190 , /||PT|