The triatomic Hamiltonian in Jacobi coordinates for the H+ O system
in the BF frame is given by:
Using the basis set expansion, given in Eq. (3), we derived the equations of motion; these are given in detail in Ref th:anth1998. The equations of motion are tridiagonal in the projection quantum number : all terms are diagonal in except for the the Coriolis terms, which couple to .
In our previous calculations on this system,[1,2]
we used a symplectic integrator to propagate a complex
wave packet according to the time dependent time-dependent Schrödinger
equation. The large calculations, however, are computationally more
demanding, not only because of the greater number of states, but also
because of the larger spread in the energies involved. Thus, the time-step
required by the symplectic method is reduced and consequently the number
of time-steps needed to converge the calculations is increased. In order
to reduce our overall computation time, we employ the real wave packet
propagation scheme of Gray and Balint-Kurti. In this
method, the Hamiltonian in Eq. (4) is represented by a finite,
real matrix H and the wave packet is represented by a real finite vector . Similar methods have been developed by
Kouri and coworkers,[62,63],
Mandelshtam and Taylor,[64,65],
Kroes and Neuhauser,
and Chen and Guo.[67,68]
The method is based upon the following ideas:
[61,69] (i) reaction probabilities may
be inferred from knowledge of the real part of the wave packet alone,
(ii) reaction probabilities, pertinent to a Hamiltonian operator,
may be inferred from a wave packet which has evolved under a modified
time-dependent Schrödinger equation in which the Hamiltonian operator has
been replaced by a function of itself, , (iii) it is possible
to choose such that the act of propagating forward one time
step is accomplished by a simple Chebychev iteration involving a single
evaluation of the Hamiltonian on a real vector:
The real wave packet method results in a far more efficient calculation than we could have achieved using the symplectic integrator. For all , the total propagation time is 30,000 a.u. ( 0.72 ps); for the real wave packet calculations this is an ``effective time''. The symplectic method required a reduction in the size of the time-step as increased from 10 a.u. for to 4 a.u. for . Thus, the calculation required 18,000 complex evaluations of (where is the vector representation of ) whereas 45,000 such evaluations were required for (i.e. 36,000 and 90,000 real matrix-vector multiplications). From trial calculations, we estimate that the symplectic approach would have required more than evaluations for . In sharp contrast, the real wave packet approach required only 37,500 real matrix-vector multiplications for and for . Thus, we achieve approximately an order-of-magnitude improvement over the method used in our previous work. Furthermore, we were able to use the same effective time step  for all except for , where we had to reduce it somewhat.
To reduce the amount of memory that the calculations require, we transform from the associated Legendre basis or finite basis representation (FBR) to a grid or discrete variable representation (DVR) in to compute the action of the potential energy operator, . In the DVR representation, the potential energy matrix is diagonal, whereas the FBR representation requires a reduced potential matrix with elements given by . Just as the symmetry of the O molecule allows us to use only odd states in the Legendre expansion, it also allows us to use only half of the Gauss-Legendre grid points in the DVR representation.[72,73] Using this FBR-DVR approach we reduce memory usage in our calculations by 100-150 MBytes.
Because has to be transformed from FBR to DVR and back in each iteration, the evaluation of can be less efficient in CPU time than in the straight FBR method. Initially, we implemented the DVR-FBR transformation as a series of matrix-vector multiplications, where the stride in the BLAS[74,75,76] routine is the total number of DVR points divided by the number of angular points (Note that this stride will be less than , because of the point-selection mechanism we are using). However, if is the transformation matrix, where rotational basis functions label rows and angular points columns, we can rewrite the transformation as , where is simply rewritten as a matrix where angular points label columns and points and points label rows. In the resulting matrix the rows are labeled by -points and points and the columns by rotational basis functions. Of course in a FORTRAN program this matrix is easily used as a vector in the calculation with the rotational energy operator, which gives us and . The reverse transformation to the DVR is then obtained through , which gives us and subsequently . This turns out to make the current evaluation of at least as efficient as the one used in our earlier papers. Note that, because we now use a grid-based method for instead of a basis, we can, in principle, apply a similar algorithm as used for and to remove grid points that do not play a role in the calculation. This would allow us to create a grid tailor-made for a specific reactive system. This algorithm was not used in the calculations presented in this paper, however.
We use a recently developed flux-based approach for extracting the energy-dependent reaction probabilities, , from the real wave packets, where is the total energy in the system.