next up previous
Next: Parallel Method Up: Theory and Computational details Previous: Coordinates and Basis set

Hamiltonian, Propagation, and Analysis

The triatomic Hamiltonian in Jacobi coordinates for the H+ O$_2$ system in the BF frame is given by:

$\displaystyle \hat{H} = -\frac{\hbar^2}{2\mu_R}\frac{\partial^2}{\partial R^2} ...
...{J}-\hat{j})^2}{2\mu_R R^2} + \frac{\hat{j}^2}{2\mu_r r^2} + V
(R,r,\vartheta),$     (4)

where $\mu_R$ and $\mu_r$ are the reduced mass of the H + O$_2$ collision complex and the reduced mass of O$_2$, respectively. $V(R,r,\vartheta)$ is the intermolecular potential. In our case, we use the DMBE IV surface.[48]

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 $\Omega$: all terms are diagonal in $\Omega$ except for the the Coriolis terms, which couple $\Omega$ to $\Omega\pm1$.

In our previous calculations on this system,[1,2] we used a symplectic integrator[60] to propagate a complex wave packet according to the time dependent time-dependent Schrödinger equation. The large $J$ 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.[61] In this method, the Hamiltonian in Eq. (4) is represented by a finite, real matrix H and the wave packet $\Psi$ is represented by a real finite vector ${\bf q}$. Similar methods have been developed by Kouri and coworkers,[62,63], Mandelshtam and Taylor,[64,65], Kroes and Neuhauser,[66] 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, $f(\hat{H})$, (iii) it is possible to choose $f(\hat{H})$ 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:

{\bf q_{k+1} } = {\bf A} \cdot (- {\bf A q_{k-1}} + 2{\bf H_s q_k}),
\end{displaymath} (5)

where ${\bf H_s}= a_s{\bf H} + b_s$. $a_s$ and $b_s$ are chosen in such a way that all eigenvalues of $\mathbf H_s$ lie between $-1$ and 1. The matrix $\mathbf A$ ensures that the wave packet is absorbed at the edge of the grid to avoid any unphysical reflections. $\bf A$ is defined as the matrix representation of the operators $A^R_\lambda$ and $A^r_\nu$ from paper I.

The real wave packet method results in a far more efficient calculation than we could have achieved using the symplectic integrator. For all $J$, the total propagation time is 30,000 a.u. ($\approx$ 0.72 ps); for the real wave packet calculations this is an ``effective time''.[70] The symplectic method required a reduction in the size of the time-step as $J$ increased from 10 a.u. for $J=0$ to 4 a.u. for $J = 10$. Thus, the $J=0$ calculation required 18,000 complex evaluations of ${\bf H\Psi}$ (where $\bf\Psi$ is the vector representation of $\Psi$) whereas 45,000 such evaluations were required for $J = 10$ (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 $180,000$ ${\bf H\Psi}$ evaluations for $J=35$. In sharp contrast, the real wave packet approach required only 37,500 real matrix-vector multiplications $\bf Hq$ for $J = 15-25$ and $42,500$ for $J=35$. 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 [61] for all $J$ except for $J=35$, 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 $\vartheta$ to compute the action of the potential energy operator, ${\bf V}$.[71] In the DVR representation, the potential energy matrix ${\bf V}$ is diagonal, whereas the FBR representation requires a reduced potential matrix ${\bf V^\Omega_{j,j'}}$ with elements given by $<\Theta^\Omega_j\vert
V\vert \Theta^\Omega_{j'}>$. Just as the symmetry of the O$_2$ molecule allows us to use only odd $j$ 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 $\bf q$ has to be transformed from FBR to DVR and back in each iteration, the evaluation of $\bf Hq$ can be less efficient in CPU time than in the straight FBR method. Initially, we implemented the DVR-FBR transformation as a series of $n_R\times n_r$ 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 $n_R\times n_r$, because of the point-selection mechanism we are using). However, if ${\bf U}$ is the transformation matrix, where rotational basis functions label rows and angular points columns, we can rewrite the transformation as ${\bf QU}^T$, where $\bf Q$ is simply $\bf q$ rewritten as a matrix where angular points label columns and $R$ points and $r$ points label rows. In the resulting matrix $\bf X$ the rows are labeled by $R$-points and $r$ points and the columns by rotational basis functions. Of course in a FORTRAN program this matrix $\bf X$ is easily used as a vector $\bf x$ in the calculation with the rotational energy operator, which gives us $\bf x'$ and $\bf X'$. The reverse transformation to the DVR is then obtained through $\bf QX'$, which gives us $\bf C'$ and subsequently $\bf c'$. This turns out to make the current evaluation of $\bf Hq$ at least as efficient as the one used in our earlier papers. Note that, because we now use a grid-based method for $\vartheta$ instead of a basis, we can, in principle, apply a similar algorithm as used for $R$ and $r$ 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.[77] 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, $P_{\nu_i,j_i,\Omega_i}^{J,p}
(E)$, from the real wave packets[69], where $E$ is the total energy in the system.

next up previous
Next: Parallel Method Up: Theory and Computational details Previous: Coordinates and Basis set
Anthony J. H. M. Meijer 2000-10-05