Iterative quantum computations of HO2 bound states and resonances for J = 4 and 5

Bound and resonance states of HO2 have been calculated by both the complex Lanczos homogeneous filter diagonalisation (LHFD) method1,2 and the real Chebyshev filter diagonalisation method3,4 for non-zero total angular momentum J = 4 and 5.

For bound states, the agreement between the two methods is quite satisfactory; for resonances while the energies are in good agreement, the widths are only in general agreement.

The relative performances of the two iterative FD methods have also been discussed in terms of efficiency as well as convergence behaviour for such a computationally challenging problem.

A helicity quantum number Ω assignment (within the helicity conserving approximation) is performed and the results indicate that Coriolis coupling becomes more important as J increases and the helicity conserving approximation is not a good one for the HO2 resonance states.


Exact calculations for total angular momentum J > 0 are essential in fully understanding quantum reaction dynamics.

However, these J > 0 calculations are still very challenging even for triatomic reactions, especially when dealing with complex-forming systems.

The major reason for this situation is the so called ‘angular momentum catastrophe’:5 many J > 0 calculations have to be performed, and the size of the Hamiltonian matrix increases linearly with J.

For these non-zero J calculations, it is apparently impractical to employ conventional direct diagonalisation methods due to the requirement of a significant computer core memory.

On the other hand, iterative methods are well suited to solve this large-scale eigenvalue problem.

Among different iterative methods, two of the most powerful methods which have emerged are the real Chebyshev filter diagonalisation method (RCFD)3,4,6 and complex Lanczos filter diagonalisation method.1,2,7–9

The basic idea of filter diagonalisation10 is to act with operator functions f(EiĤ) onto a generic initial wavepacket to generate filtered states at a sequence of chosen filter energies Ei.

Eigenstates in the energy window can then be obtained by diagonalising a small-sized Hamiltonian matrix in a subspace spanned by the filtered states.

However, a significant difficulty presents itself in that a large storage requirement is needed for storing the filtered states in the primary representation in earlier FD versions.

Later, Wall and Neuhauser proposed an auto-correlation function based filter diagonalisation scheme,11 in which eigen-values in the energy windows can be obtained from a single sequence of auto-correlation functions without the need to explicitly construct the filtered states.

This essentially eliminates the memory bottleneck.

Mandelshtam and Taylor3 combined this scheme with their modified real Chebyshev propagation approach, and implemented a box-like low storage filter diagonalisation formalism.

The most important advantage associated with this approach is that one can employ a real algorithm with a single, extended Chebyshev vector recursion.

This leads to substantially greater efficiency in comparison with step-by-step propagation of a complex wavepacket.

Very recently, Li and Guo proposed a very efficient scheme of doubling Chebyshev correlation functions for calculating narrow resonances,12 with numerical validation of their proposal for two molecular systems.

Meanwhile, Neumaier and Mandelshtam derived a pseudo-time Schrödinger equation and provided a rigorous proof that an exact doubling formula exists for damped Chebyshev propagation.13

Another attractive approach is to use iterative methods based primarily upon the Lanczos algorithm of tri-diagonalisation.14,15

Lanczos-based FD methods have been introduced by Yu and Smith7–9 and by Chen and Guo.16

These methods are especially useful when one only considers a small section of the entire spectrum as in most applications.

In addition, they have the desirable property of avoiding most of the ghost eigenvalues that will appear in a normal Lanczos diagonalisation.

A key feature of Smith et al.’s Lanczos representation FD methods8,17 is that the entire calculation is carried out within the Lanczos subspace.

This has the major advantage that one does not incur the memory cost of having to explicitly store the filtered states in the primary representation.

Consequently a single Lanczos subspace serves to analyse the entire spectrum, since the tridiagonal representation of the Hamiltonian makes it a straightforward exercise to generate a new set of filtered states in the next energy window by solving QMR or MINRES equations implicitly within the Lanczos subspace.

Very recently, we have developed a simpler and more efficient Lanczos homogeneous filter diagonalization (LHFD) algorithm based on a very simple homogeneous filtering recursion within the Lanczos representation.1

Using the LHFD method, J = 1–3 bound and resonance calculations on HO2 have been performed in a previous paper.18

In that paper we compared our calculations with Wu and Hayes’s results19 for low lying bound states for J = 1–3, but for resonances no other calculations have been available except our own Lanczos results.

For J = 4 and 5, to our knowledge, no calculations have been published for either resonance or bound states, so we will use the two most powerful FD methods presently available (namely, the Chebyshev and Lanczos approaches) to perform the J = 4 and 5 calculations.

The computational tasks are demanding for a system with a deep well such as this, in particular for the purpose of resolving individual resonances.

It will also be of some interests to compare the relative efficiency of Lanczos-based and Chebyshev-based approaches to the computation of molecular bound states and resonances for more challenging non-zero J calculations.

Some such comparative studies have appeared for bound state problems,20–24 and for resonance calculations.25,26

In the above comparative studies,23,24,26 it was found that in RCFD the convergence behaviour can be accurately predicted by the spectral range and the physical density of states.

For resonance calculations, the Chebyshev algorithm required roughly the same number of iterations as the Lanczos method if doubling schemes12,26 for computing the autocorrelation functions are employed, but due to its use of real algebra the Chebyshev FD methods may gain some advantages.

In the present work we present a comparison of Lanczos-based and Chebyshev-based FD methods for the much more demanding problem of the resonance calculations of the HO2 radical for J = 4 and 5.

Our objective is to examine how the relative efficiency of the Lanczos and Chebyshev iterative approaches scales with the order of difficulty of the calculation.

It is now well known1,9 that convergence of the very narrow resonances requires ca. 40 000 complex-symmetric Lanczos iterations for J = 0 HO2 case, and it will be interesting to see the convergence behaviour for higher J values of HO2 for the two iterative methods.

Due to the computational constraints, approximate quantum methods such as adiabatic rotation (AR),27J-shifting28 and helicity conserving (HC)29 approximations are commonly used for non zero J calculations.

The key issue in these approximations is whether a reasonably good quantum number Ω associated with the projection of total angular momentum on a body fixed axis exists.

If the substates Ω of the wave function for J > 0 are heavily coupled, the Coriolis coupling between the states cannot be ignored and any attempts to assign the helicity quantum number Ω will fail.

Instead of explicitly computation of bound and resonance states within these approximations, we will examine this issue by a helicity quantum number Ω assignment for both bound and resonance states.

The rest of this article is organised as follows.

In section 2 we briefly describe the theoretical methods needed to characterize both bound and resonance states for non-zero total angular momentum.

In section 3 we shall present the results of J = 4 and 5 bound and resonance calculations and give discussions concerning the relative performances of the two iterative FD methods and the helicity Ω assignment.

Section 4 concludes.


In general, we treat the three internal Jacobi coordinates (R,r,γ) in discrete variable representation (DVR), while the three Eulerian angles (θ,ϕ,ψ) are described in a basis set.30–32

This procedure is very efficient because the potential part of the Hamiltonian matrix is diagonal, which can reduce the memory requirement substantially.

Using symmetry-adapted symmetric top eigenfunctions to expand the total wave function, one can get the Hamiltonian matrix explicitly represented in DVR (for more details the reader is referred to ref. 18):withIn eqn. (1), we have used Ω-dependent DVR for γ coordinate, which is obtained by either diagonalizing the coordinate, operator (x = cos γ) matrix or by a Gauss–Jacobi quadrature scheme.

We have compared the two DVR schemes, and the DVR points as well as the transformation matrix T from the two methods are nearly the same.

We calculate at the outset and then store the neighbouring Coriolis coupling matrices [see the last two terms in eqn. (1) for the details of the matrix elements] for each Ω component.

The memory requirement for the coupling matrices is not large, and whenever they are needed in the iterations, we use them to perform the Hamiltonian matrix vector multiplications directly within the DVR.

For R and r coordinates, we have used potential optimised DVR.33

In the two iterative methods, the basic propagation is quite similar in principle, which is a three-term recursion.

In Lanczos iteration, we use the basic Lanczos algorithm for complex-symmetric matrices34 to project the non-Hermitian absorbing potential augmented Hamiltonian into a Krylov subspace.

We need to store the two complex vectors in Lanczos iterations for later FD analysis to extract physical information such as energies and widths.

Similarly, in Chebyshev propagation, a real modified Chebyshev polynomial propagation proposed by Mandelshtam and Taylor35–37 has been employed.

Again, we will store the correlation functions in the Chebyshev propagation for later FD analysis.

After setting up tridiagonal Lanczos subspace, we perform LHFD inside the subspace to extract the bound and resonance information.

The key issue in LHFD is to solve the homogeneous linear system(EjTM)|ϕ(Ej)〉 = 0to generate filtered state ϕ(Ej) at different filter energies Ej.

Here a backward three-term substitution recursion is employed.18

Then we construct the overlap matrix with elements Sjj = (ϕ(Ej)|ϕ(Ej)) and subspace Hamiltonian matrix with elements Wjj = (ϕ(Ej)|TM|ϕ(Ej)).

The complex energies, {ε} can be obtained through solving the small-size generalized complex-symmetric eigenvalue problem WB = SBε.

The above procedures can be repeated for any chosen energy windows.

Several FD versions based on the Chebyshev recursion exist, such as Wall and Neuhauser’s version,11 Mandelshtam and Taylor’s version3 and Chen and Guo’s version.6

We note that Guo et al.’s version is essentially identical to Mandelshtam and Taylor’s version using a box filter to compute resonances.

In this paper, we will employ Mandelshtam’s real Chebyshev FD version.

For more details about RCFD, the reader is referred to refs. 3,4.


The triatomic HO2 Hamiltonian matrix was set up in terms of reactant Jacobi coordinates, and the HO2 DMBE IV PES38 was employed as we did for previous J = 0–3 bound and resonance calculations.1,2,18,25

For the two radial coordinates, a potential-optimized discrete variable representation33 (PODVR) was utilized to reduce the size of the Hamiltonian matrix and the details can be found in ref. 18 For the γ variable, Ω-dependent symmetry-adapted DVR functions were employed to take account of the odd O−O exchange parity.

The resulting direct product basis set was further contracted by discarding those points whose potential energies were higher than the cut-off energy Vcutoff = 2.0 eV, resulting in the final basis size of approximately 110700(J + 1) for even spectroscopic symmetry and approximately 110700J for odd spectroscopic symmetry.

The details of the absorbing or damping potential is similar to those used in .ref. 18

No stabilisation procedure39 has been attempted due to the huge computational resources required, and this will have some effects on resonance widths (see below).

The bound state energies as well as the resonance energies and widths for two chosen energy windows for J = 4 and 5 have been computed.

The first energy window is for the lowest bound state energies from −0.08 eV to 0.92 eV.

Here the zero energy point is referred to as the ground state energy of HO2 for J = 0, which is −2.015861 eV relative to the H + O2 dissociation limit.

This energy window is relatively easy to calculate and a Lanczos subspace size M = 5000–10 000 is enough to converge all the energies in this window.

In RCFD, Chebyshev iteration number needed is roughly n = 20 000–30 000.

In Tables 1–4 we have listed the selected lowest bound state energies for each symmetry of J = 4 and 5 from both LHFD and RCFD methods for comparison.

Inspection of the energies shows that the agreement between them is quite satisfactory and 4 digits of relative accuracy have been achieved for most of the energies.

The second energy window we have chosen is close to and above dissociation threshold, namely, the highest lying bound state energies and lowest resonance energies and widths from 2.10 eV to 2.18 eV.

Since these are the first calculated results, the convergence has been carefully tested and in Fig. 1 we have plotted the convergence behaviour for one resonance at E = 2.1227 eV for J = 5 (even symmetry) case from the two iterative FD methods.

From Fig. 1(a) one can see that M > 1 60 000 Lanczos iterations can well converge the resonances and in all our calculations we have used a larger M = 2 20 000 Lanczos subspace size for this energy range.

Similarly, Fig. 1(b) indicates that M > 5 00 000 Chebyshev iterations are needed to converge the resonance and in all our calculations we have used a larger M = 6 00 000 Chebyshev iteration number for the second energy window.

These results show that the convergence is very slow for resonances and almost monotonic.

The difficulty of convergence for HO2 system is related to both the changes of spectral range and the physical densities of states as J increases.

By comparison with previous J = 0–3 results, we also found that the iteration number needed to converge the resonances increases as J increases.

For example, in RCFD method, 1 00 000 iterations are enough to converge most of the resonances for J = 0 case; similarly, in LHFD method, 40 000 iterations are enough to converge most of the resonances for J = 0 case.

It is understandable because there are more energy levels to converge for high J values and we believe that this behaviour is general for higher J values.

In Fig. 2 we have plotted the resonance widths versus energies for J = 4–5 for both even and odd symmetries from the two methods.

Inspection of the resonances in these plots, one can see that while for resonance energies the agreement is good, for resonance widths the differences are relatively large.

In principle, the complex resonance energies should be the same for the two methods if a perfect absorbing potential and a perfect damping operator are respectively used.

However, it is very challenging to obtain the perfect absorbing or damping parameters using stabilisation method39 or optimizing method26 for the HO2 system, especially for J > 0 calculations and we do not attempt this in the present work.

For LHFD, we did only 6 test calculations by varying absorbing strength parameter V0 in J = 5 even symmetry case.

The results indicate that the best V0 should fall in between 0.02 eV and 0.06 eV and we take V0 = 0.02 eV in this paper.

For RCFD, the damping potential parameters are based upon previous J = 0 calculations.

Since approximate absorbing boundary conditions have been imposed here via two quite different methods it will be apparent that the resonance widths reported in this paper for either method are simply estimations of the true resonance widths.

Having said this, observation of the plots from Fig. 2(a) to Fig. 2(d) indicates that the agreement for most resonances between the two methods is not too bad, albeit not perfect.

For the HO2 system, other investigators have also pointed out40,41 that the resonances widths in particular are sensitive to many details of the calculations (e.g., basis set size, grid size, cut off energy), which adds to the complexity of accurately calculating the resonance widths.

In principle, a stabilisation procedure39 or optimizing approach26 can be employed to determine more accurate resonance widths in the future.

Here it is interesting to note that although the resonance widths do not agree very well, the general variation trends with energy are still similar for the two methods.

In Tables 1–4 an Ω assignment has been given for the low lying bound states, supposing that the helicity conserving or adiabatic rotation approximation holds for most of them (because there exist near degeneracy for the same Ω components from both symmetries, it is possible to assign them by comparing the calculated energies from even and odd symmetries).

The purpose of Ω assignment is to investigate the importance of the Coriolis coupling for this system.

If this assignment is successful, then helicity conserving calculations or even much simpler adiabatic rotation approximations should be accurate, which will save quite a lot of computational time.

We note that two bound states at about 0.318641 eV from the J = 5 (even symmetry) case in Table 3 are too degenerate to be distinguished and we have put them together.

In J = 5 (odd symmetry) case, the two bound states are still distinguishable.

For the high-lying bound states as well as for the resonances, we have failed to assign them unambiguously.

For example, we have analysed the high lying bound state energies near the dissociation threshold from J = 4–5 calculations for both even and odd spectroscopic symmetries, respectively (the results for the high lying bound state energies are not shown here, and they can be acquired from us upon request).

While only several of them can be assigned tentatively, most of them cannot be assigned with confidence.

The indication is that the mixing of different Ω components is so strong that Ω is no longer a good quantum number.

Of course, the difficulties in assignment also arise from the fact that the spacings between these high lying bound states and resonance states are becoming smaller and smaller.

For overlapping resonances, the assignment is further complicated by the mixing of different resonance states, i.e., at one resonance energy, other neighbouring resonances might interfere with this ‘main’ resonance.

For this system, it seems that HC calculations or adiabatic rotation approximation can give reasonably accurate results only for low bound state energies.

This observation is in consistent with the previously reported J > 0 total reaction probability calculations for this system, which show that for HO2 the Coriolis coupling is important and cannot be ignored.42

Interestingly, this situation is in contrast to the H2O and HOCl system, for which HC or AR is a good approximation.22,43,44

Finally, some more comparisons will be made for the two seemingly similar iterative methods.

The Chebyshev method is time dependent in essence, and the time can be thought of as the iteration number n (in Chebyshev propagation the system evolves under the function of the Hamiltonian).

The way to add the damping function in the iteration is similar to the implementation of damping at the edges of a TD wave packet.

In contrast, the Lanczos method is time-independent in essence, and Lanczos recursion is to transform the Hamiltonian from primary representation into the tridiagonal Lanczos representation.

The complex absorbing potential is included in the Hamiltonian, which is quite different from the implementation of damping in the Chebyshev iteration.

In both iterative methods, there is similar convergence behaviour, i.e., the eigenvalues are converged progressively from bound states to resonances, but Lanczos method converges eigenvalues in the low/upper part more quickly than Chebyshev method.

This can be seen from the convergence behaviour in the energy window 1.

For high lying bound states and resonances such as those in energy window 2, the two methods are becoming comparable.

We have not implemented the doubling schemes12,26 for the autocorrelation functions in this paper, however, one notes that if this is done then the Chebyshev FD method will require only slightly more iterations than Lanczos method.

The basic matrix-vector multiplication procedures of the two methods are very similar.

In our Lanczos code, the matrix-vector multiplications are realised through multiplication of a real matrix onto a complex vector with a slight additional overhead for the diagonal complex absorbing potential.

In the Chebyshev algorithm both the Hamiltonian matrix and the vector are real.

Of course, the operation of real-complex multiplication such as storing and retrieving and writing will take more time than that of real-real multiplication.

Thus considering its real algebra and doubling scheme, the Chebyshev FD method may be marginally more efficient than the Lanczos method in cpu time.

This is still similar to previous investigations for J = 0 calculations,25,45 even though the order of difficulty of the problem is increased markedly in this application.

Of course, the real Chebyshev method will be more efficient in terms of memory, due to real vectors being employed in the recursion.

Here we note that we have used Lanczos method to accurately estimate Hmin and Hmax used in our Chebyshev calculations, i.e., the spectral range and the physical density of states are the same for the two methods.


In this paper two iterative filter diagonalisation techniques have been applied to HO2J > 0 calculations and converged bound state energies as well as resonance energies and widths have been obtained for J = 4–5.

Regarding the computer memory requirement and the efficiency, the two iterative methods are roughly comparable to one another for resonance calculations, with RCFD having a marginal advantage particularly in so far as memory is concerned.

This behaviour will not change even though the order of difficulty of the problem is increased significantly from J = 0 to higher J values.

For bound state and resonance energies, agreement between the two methods is quite satisfactory.

For resonance widths, while the general trend with energy is similar for the two methods, the details of the fine structures are different.

We associate this difference with the different ways in which approximate absorbing boundary conditions are imposed by the two methods.

While in principal using a stabilisation procedure or finding a perfect absorbing potential for all resonances should resolve these differences, such approaches still represent a major computational challenge for long-lived complex-forming systems, especially for the non-zero total angular momentum values explored in this paper.

Currently, we are developing alternative methods to compute resonances without employing complex absorbing potentials or damping functions, based upon Lanczos arbitrary boundary inhomogeneity method46–48.