\def\header{\it Computational Procedure \hfil Page \folio \hfil \today}

\medskip
\centerline{\bf V. COMPUTATIONAL PROCEDURE}

\smallskip
\noindent {\bf A. Overview}

In this section we will outline our computational procedure
for reactive scattering.
The overall strategy is given in Fig. (7).
This schematic indicates the four steps used to
procure the desired scattering matrices:
   \item{1.}An optimized mesh (or element definition)
      is constructed in program MESHER at a multitude of fixed hyperradii.
      This program was designed to minimize
      computational time, to improve the accuracy of the surface functions,
      and to automate the creation of input data for the surface function
      program.
   \item{2.}The optimized mesh is employed by program SURFUNC
      (a finite element code) to
      generate surface functions,
      and eigenvalues of the surface function
      Hamiltonian Eq. (16). The numerical values of these surface functions
      at each mesh point is essential for all of the subsequent steps.
      This step {\it currrently} requires significant computer time
      {\it and} disk storage.
   \item{3.}The programs, MATELM and OVERLAP (which are independent of one
      another and can be executed simultaneously) use the surface functions to
      build coupling matrix integrals of Eqs. (21)-(23) and
      overlap integrals of Eq. (20) respectively.
   \item{4.}The APH3D integrator (scattering program) calculates
       matrices and transition probabilities.
       The transition matrices are safeguarded
       on disk for subsequent analysis (to predict differential cross sections,
       integrated cross sections, resonances, Wigner cusps, etc.).

\noindent
Steps 1 through 3 are
independent of the total scattering energy.  Consequently, once all of this
information has been accumulated the scattering can be done very efficiently
at many energies.  It should be {\it emphasized} that steps 1 through 3 are exclusively
necessary for the APH region.  Also, a considerable amount of step 4 is also
energy independent and the energy independent portion is inscribed on disk.
Each of these steps will now be discussed in more detail.

\Figsize 2.75 truein
\Fig 7. Schematic diagram of the four-step computational
procedure used to procure S-matrices. \par

\smallskip
\noindent {\bf B. Meshing}

MESHER establishes an optimized mesh for the surface function
program.  SURFUNC requires as input mesh made up of discrete points
(unique $\theta,\chi_i$ pairs)
called {\it nodes} and boundary conditions that are placed on the mesh.

For the reactions considered in this proposal the potential is a function of three of the
APH coordinates, ($\rho,\theta,\chi_i$). The surface functions will be calculated at
a series of fixed $\rho$ values on the ($\theta,\chi_i$) surface.  The nodal
points are specified by unique $\theta,\chi_i$ coordinate pairs. A set of
$3\times 3$ nodes form an {\it element}. Accordingly, we are using
quadrilateral elements with 9 nodes per element.
Hundreds of nonoverlapping elements of different sizes
make a {\it mesh}.

The shape of the mesh (size of each of the elements) is determined by two
conflciting characteristics of any finite element code. First, the more
nodal points included in the mesh, more accurate will be the eigenfunctions
and eigenvalues. Second, the more nodal points in the mesh, the more computer
time will be required to calculate the desired surface functions. The finite
element code uses a polynomial [$(\theta,\chi)$ space] in each element,
as a basis
for the surface functions.  In regions where the amplitudes of the surface
functions are large it is desirable to have as fine a mesh as one can afford.
However,  in
the regions where the surface functions are nearly zero one can use fairly
large elements or a coarse grid. Ideally, the mesh would be very dense in the
strongly attractive regions, very sparse in the very repulsive classically
forbidden regions, and have a smooth transition between the two regions. To obtain
this desired result we start with a very sparse mesh and then subdivide it
many times, concentrating each time on the deeper part of the potential.

The meshing process can be divided into five parts. The first part
initializes variables, then makes the crudest mesh. Next, several crude divisions
are made to give the mesh the character of a sparse mesh in the
repulsive region and a finer mesh in the attractive regions. Then several, more
selective divisions are made which provide as many nodes as possible in the
deepest part of the wells. Once the third step assigns to the mesh the
desired number of nodal
points, the mesh is determined. Then nodal points and
elements are numbered and put into a form that the surface function program
understands.  Finally, boundary conditions are imposed on the mesh to get
surface functions of the desired symmetry. An optimized mesh has about
3,000 to 5,000 nodal points instead of the $\approx$50,000 nodal points
of a uniform mesh. This procedure allows one to obtain highly
accurate surface functions with minimal cost. A typical mesh for
$\HHH$ is shown in Fig. (8).

\Figsize 2.25 truein
\Fig 8. The position of the nodal points for the
Porter-Karplus potential energy surface for $\HHH$.
The nodal points are located in regions where the potential is most
attractive and is fairly sparse in the classically forbidden regions.
See Fig. (5)\par

The Porter-Karplus potential energy surface
for $\HHH$ with the superimposed
mesh is displayed in Fig. (9a) corresponding to $\rho=3.6210$.
The optimized mesh is conviently restricted to the lower
regions of the potential. Currently, the elements are restriced to be
rectangular in [$(\theta,\chi_i)$ space] and we {\bf propose} to
augment our optimization procedure by allowing the elements to deform
slightly to a superior distribution
of points in the classically allowed and weakly classically forbidden regions.

\Figsize 3.5 truein
\Fig 9. a) The Porter-Karplus interaction potential
for the $\HHH$ system at a hyperradius of 3.6210$a_0$.
b) The surface function $\Phi^0_{10}(\theta,\chi_i;\rho)$ at $\rho=3.6210a_0$.
c) The surface function $\Phi^0_{20}(\theta,\chi_i;\rho)$ at $\rho=3.6210a_0$.
d) The surface function $\Phi^0_{30}(\theta,\chi_i;\rho)$ at $\rho=3.6210a_0$.\par

\smallskip
\noindent {\bf C. Surface Functions}

The surface function Hamiltonian in Eq. (16)
is a partial differential equation in the two variables $\theta$ and $\chi_i$.
The finite element method is used to solve this eigenvalue equation for
the eigenfunctions [surface functions
$\Phi^J_{t\Lambda}(\theta,\chi_i,\rho_h)$] and eigenenergies $\Escr^J_{t\Lambda} (\rho_h)$.

The finite element method is a powerful numerical technique for solving
partial differential equations. The fundamental idea of this method
element is that surface functions which are continuous can be
approximated by a set of piecewise-continuous polynomials defined over a
finite number of subdomains called elements. A
major advantage
of the finite element method is that the size of the elements can be
varied across the domain. Hence small elements can be used in physically
accessible regions of space and large elements in physically
inaccessible regions of space. Two types of elements are
used to model a two-dimensional domain, the triangle and
quadrilateral. Originally triangular shaped elements were used almost
exclusively.
It is now generally accepted that the quadrilateral elements are more
efficient. In Fig. (8) we show the domain subdivided into a finite set
of elements. In each element the surface
function is approximated by an interpolating polynomial of the form
$$\Phi^J_{t\Lambda}(\theta,\chi_i)
=a_1^k+a_2^k\theta+a_3^k\chi_i+a_4^k\theta^2 +a_5^k\theta\chi_i+a_6^k
\chi_i^2+a_7^k\theta^2\chi_i+a_8^k\theta\chi_i^2+a_9^k\theta^2\chi_i^2
\eqno(32)$$
This equation has 9 coefficients, so the element must have 9
nodes to uniquely determine the
polynomial coefficients. These nodes are at the corners of the
quadrilateral, the midpoint of each of the sides, and the center of the
quadrilateral.

The finite element method is based on the variational equation
$$\int_0^{\pi/2}d\theta \int_0^{2\pi} d\chi_i\delta\Phi^J_{t\Lambda}
\left[T_h+V+{15h^2\over 8\mu\rho^2}\right]\Phi_{t\Lambda}=0,\eqno(33)$$
where $\delta\Phi^J_{t\Lambda}$ is the variation of the surface function.
Using Green's first identity and the fact that the subdomains or elements
are pairwise disjoint, this equation can be written
$$\sum_{k=1}^N\int_{(\Delta\theta)_k} d\theta\int_{(\Delta\chi_i)_k} d\chi_i
\left\{\left[ 4\sin 2\theta{\partial\Phi^J_{t\Lambda} \over \partial
\theta}\right]^2
+\left[{2\cos\theta\over\sin\theta}{\partial\Phi^J_{t\Lambda} \over \partial \chi_i} \right]^2
+\left[V+{15\hbar^2\over 8\mu}{\partial\Phi^J_{t\Lambda} \over \partial
\theta} \right]^2\right\}=0\eqno(34)$$
where the integral is over the $k^{th}$ element. The integrals over the
domain of each element are evaluated {\it via} 3-point Gauss-Legendre
quadratures in the $\theta$ and $\chi_i$ integrations.

Since the surface functions are continuous, the nodal points for
adjoining elements are equal. To facilitate this procedure a
``connectivity matrix'' is used to
relate the local node numbers $j$, of the n$^{th}$ element to the
global node numbers $k$
$$M(n,j)=k.\eqno(35)$$
Note that node numbering is critical, since it
determines the bandwidth of the matrices which are to be diagonalized.
We have tried to be prudent about the node numbering but have not tried
to optimize this numbering. We {\bf propose} to optimize
this numbering procedure in the future.

 Eq. (34) can be written as a generalized eigenvalue equation
$$[H-\Escr U]\Phi=0\eqno(36)$$
where
$$H_{j^{\prime}j}=\sum_{i^{\prime}=1}^{N_p}\sum_{i=1}^{N_p}C^n_{i^\prime
j^\prime}\Hscr_{i_\prime i} C^n_{ij},\eqno(37)$$
$$U_{j^{\prime}j}=\sum_{i^{\prime}=1}^{N_p}\sum_{i=1}^{N_p}C^n_{i^\prime
j^\prime}\Uscr_{i_\prime i} C^n_{ij}\eqno(38)$$
and the connection coefficient is
$$C^n_{ij}=\delta_{j,M(n,i)}\eqno(39)$$
The $j^{th}$ element of the eigenvectors in (36) corresponds to the
value of the surface function at the $j^{th}$ node. The eigenvalues are the
eigenenergies of Eq. (16). We use the subspace iteration method to
solve for the eigenvalues and eigenvectors.

For a zero interaction potential the surface function Hamiltonian
can be solved analytically
for hyperspherical harmonics. For $J=0$, the hyperspherical harmonics
are
$$\Phi^J_{t\Lambda}=(1-\cos\theta)^{\vert m_t\vert/2}P_{n_t}^{(\vert m_t \vert,0)}(\cos\theta)
\cases \sin m_t\chi_i\cr\cos m_t\chi_i\endcases,\eqno(40)$$
where $P_{n_t}^{(\vert m_t \vert,0)}(\cos\theta)$ is a Jacobi Polynomial. The
associated eigenvalues are
$$\Escr^J_{t\Lambda}={15\over 4}+\lambda(\lambda+4),\eqno(41)$$
where $\lambda=4n_t+\vert m_t\vert$.
These analytic solutions have been compared
with the numerical solutions to check for accuracy and rate of convergence.
We found that the finite element method converges
rapidly to exact results as the number of elements for all of
the solutions is increased.

For collision energies below the three body breakup threshold the interaction
potential drastically restricts the region over which the surface functions
may have an appreciable amplitude.
We also calculated surface functions corresponding to the
Porter-Karplus potential for $\HHH$.
As pointed out by Kuppermann another advantage of hyperspherical coordinates
is that  mappings\ref{21-23} of the potential energy surface
are easily visualized.
The Porter-Karplus interaction potential energy surface
and first 3 surface functions are shown in Fig. (10) for $\rho_h=6.02a_0$.

\Figsize 3.5 truein
\Fig 10. Same as Fig. (9) except with $\rho$ fixed at
$\rho=7.68180a_0$. \par

The eigenenergies
as a function of $\rho$ are plotted in Fig. (11).
For large values of $\rho$, the eigenenergies are approaching the
rotational-vibrational eigenenergies
of the isolated diatom H$_2$ as one would expect.

\Figsize 3 truein
\Fig 11. Plot of the surface function eigenvalues
in eV as a function of hyperradius $\rho$ in $a_0$.\par

The shape of the surface functions drastically changes shape as a
function of $\rho$ as illustrated Fig. (12). This clearly indicates
that no one basis is optimum {\it i.e.} using only one basis would
require solving a multitude of coupled equations.

\Figsize 3.5 truein
\Fig 12. Surface function $\Phi^0_{10}$
for the $\DHH$ isotopic
substitution reaction using the Porter Karplus potential energy surface
at 4 different values of the hyperradius $\rho$.
a) $\rho=2.97540a_0$;
b) $\rho=4.42720a_0$;
c) $\rho=6.58880a_0$;
d) $\rho=9.80720a_0$. \par

The lowest-lying surface functions for $\DHH$ are shown in Figs. (13).

\Figsize 4 truein
\Fig 13. Lowest lying surface functions
for the $\DHH$ isotopic substitution reaction using the Porter-Karplus
potential energy surface at a hyperradius of $\rho=6.58880a_0$.
a) The surface function $\Phi^0_{10}(\theta,\chi_i;\rho)$.
b) The surface function $\Phi^0_{20}(\theta,\chi_i;\rho)$.
c) The surface function $\Phi^0_{30}(\theta,\chi_i;\rho)$.
d) The surface function $\Phi^0_{40}(\theta,\chi_i;\rho)$.
e) The surface function $\Phi^0_{50}(\theta,\chi_i;\rho)$.
f) The surface function $\Phi^0_{60}(\theta,\chi_i;\rho)$.\par


Calculation of the surface functions, $\Phi$, is currently the
expensive step of our procedure. Each surface function step requires
only 6 seconds of Cray-XMP time, but huge numbers of them are
needed.  We {\bf propose} to optimize and vectorize our surface function
code.  We predict that we can compose a finite element code that will
execute $\approx 10$ times as fast.
We also intend to experiment with alternative numerical methods in
a quest for a proficient procedure for calculating surface
functions.




\smallskip
\noindent {\bf D. Matrix Elements}

Once the surfaces functions have been determined at the nodal points the
value of each surface function can be evaluated at {\it any} point,
because the finite
element code uses low-order piecewise continuous polynomials as a
basis.
Now it is easy to construct coupling matrix integrals (matrix
elements).
For each element in the mesh a low-order Gauss-Legendre quadrature
is used to evaluate the contribution to the total integral. These
partial contributions are then summed to give the matrix coupling integral
between states $\ket{i}$ and $\ket{f}$,
$$\bra{i}\hat\Oscr\ket{f}=\sum_{k=1}^N\bra{i}\hat\Oscr\ket{f}_k$$
where the subscript $K$ indicates an integral over the disjoint
subdomains defined by the $k^{th}$ element. Since the coordinate
representations of the states $\ket{i}$ and $\ket{f}$ are simple
polynomials,  one can construct coupling matrix integrals in this
manner for any operator, including differential operators such as the
coriolis coupling operator given in Eqs. (13) and (14).

\smallskip
\noindent {\bf E. Overlap Integrals}

One might think that evaluating the overlap integrals is even easier
than calculating the coupling matrix integrals. However, the overlap
integrals are integrals between surface fuctions which have been
evaluated at different values of $\rho$, and since the mesh
changes as $\rho$ changes, the elements for the different surface
functions are not identical. Since we have simply subdivided our
elements, we have the mesh for both functions. To evaluate the overlap
integrals we first check to see if the elements are the same size.
If they are identical for both values of $\rho$, we simply evaluate
their partial contribution to the overlap integral. If the two
elements are not identical, we know thaw the smaller element was
produced by subdividing the larger element. Accordingly we use the
smallest elements in evaluating the contribution to the partial
overlap integral. \smallskip
\noindent {\bf F. Scattering}

The CC equations of reactive scattering can be written in matrix
notation as $$\left[I{d^2\over d\rho^2}+K^2(\rho)-U_{eff}(\rho)
\right]\psi=0,\eqno(24)$$
where $K^2$ is the kinetic energy matrix
$$\left[K^2\right]_{t\Lambda,t^{\prime}\Lambda^{\prime}}=
\delta_{tt^{\prime}}\delta_{\Lambda\Lambda^{\prime}}{2\mu\over \hbar^2}
\left[E-\Escr^J_{t\Lambda}(\rho)\right],\eqno(25)$$
and $U_{eff}(\rho)$ is the effective potential energy matrix
$$\eqalign{
\left[U_{eff}(\rho)\right]=
&\left[{\rho_j^2 \over \rho^2}\Escr_{t \Lambda}^J
(\rho_j)
\delta_{tt^{\prime}} + \langle \Phi_{t \Lambda}^J (\theta, \chi_i
\rho_j) \vert V(\rho, \theta, \chi_i)-{\rho_j^2 \over \rho^2}
V(\rho_j, \theta,\chi) \vert
\Phi_{t^{\prime} \Lambda}^J (\rho_j) \rangle \right]
\Psi_{t^{\prime}\Lambda}^{Jn}\cr
&+{\mu \over 2 \hbar^2}\sum_{t^{\prime}}
\biggl[\lambda_+(J, \Lambda) \lambda_+ (J, \Lambda+1) \langle
\Phi_{t \Lambda}^J \vert (A-B) \vert \Phi_{t^{\prime} \Lambda+2}^J \rangle
\psi_{t^{\prime} \Lambda+2}^{Jn}\cr
&+ \lambda_-(J, \Lambda) \lambda_-(J, \Lambda-1) \langle \Phi_{t \Lambda}^J
\vert (A-B) \vert \Phi_{t^{\prime} \Lambda-2}^J \rangle
\psi_{t^{\prime} \Lambda-2}^{Jn}\biggr].\cr}\eqno(26)$$
As $/rho$ becoms large the effective potential energy matrix
approaches zero and $K^2$ approaches a constant matrix. It should be
emphasized that we are using a sector-adiabatic basis, {\it i.e.},
different surface functions are employed from sector-to-sector, but
the basis is constant within each sector.

\Aside We must clarify some terminology
before proceeding futher. Regions are substantial areas of
configuration space. We have three distinct regions: APH, Delves and
Jacobi (See Fig. 1). Frame transformations are use to transform the CC
equations from one region to another. The APH and Delves regions are
subdivided in sectors. The basis for expanding the wavefunction
changes from sector-to-sector but is constant within each sector.
A single basis is used for the Jacobi region, consequently it has
only one sector. Each sector in the regions is subdivided into
intervals. At the midpoint of each interval the coupling matrix is
diagonalized. Finally, the intervals are subdivided to form steps.
This is schematically represented in Fig. (14).\par

\Figsize 2.75 truein
\Fig 14. Schematic diagram of the procedure used to
solve the CC equations [See Eq. (24)]. The log-derivative
region is subdivided into sectors (surface functions are calculated at
the midpoint of each sector) and then the sectors are subdivided into
steps. The equations are propagated from step to using Eqs. (29-31). The
VIVAS region is subdivided into sectors and then subdivided into
intervals and finally subdivided into steps. The uncoupled equations are
propagated from step to step using piecewise analytic functions. \par


The integrator section of this APH3D program is a modified version
(more than 11,000 lines of computer code have been added) of the
PJ-Hybrid
integrator.\ref{32} With the PJ-Hybrid integrator one can determine
the scattering S-matrix by using either the variable-interval
variable-step method\ref{32a} (VIVAS), or the logarithmic-derivative
method\ref{32d}, or, preferably, a combination of both\ref{32b,32c}.

The logarithmic derivative method is essentially a stabilized
wavefunction-following method and is most efficient in regions
where the de Broglie wavelength of the scattering particle is long or
where the potential is rapidly varying. In this method, the
second-order differential equations [Eq. (24)] are transformed to a
set of non-linear first-order differential equations (Ricatti
equations)
$$Y^{\prime}(\rho)+K^2(\rho)-U_{eff}(\rho)+Y^2(\rho)=0.\eqno(27)$$
The logarithmic derivative
$$Y(\rho)=\left[\psi^J_{t\Lambda}\right]
\left[{d\over d\rho}\psi^J_{t\Lambda}\right]^{-1}\eqno(28)$$
is propagated numerically in each sector. The recurrence relation
being $$[I+Y(\rho_n)]=[2I-U_n]-[I+hY(\rho_n)]^{-1}.\eqno(29)$$
Initially we have the condition
$$[I+Y(\rho_0)]=[I-U_0]+hY(\rho_0),\eqno(30)$$
where
$$U_n=\cases h^2[K^2-U_{eff}(\rho_n)]/3, & \text{if $n=0,N$}\cr
2h^2[K^2-U_{eff}(\rho_n)]/3, & \text{if $n=2,4,\dots N-2$} \cr
8I-48\left\{6+h^2[K^2-U_{eff}(\rho_n)]\right\}^{-1},& \text{if $n=1,3,
\dots N-1$}\endcases\eqno(31)$$
At the end of each sector the logarithmic derivative is transformed to
the new local sector basis using the overlap integrals Eq. (20).

The VIVAS integrator is a potential following method that is most
efficient in regions where the potential is slowly varying. This
method is closely related to the Gordon, Magnus or R-matrix methods.
In the VIVAS method,  sectors are subdivided into {\it intervals}. At
the midpoint of each interval, the potential interaction is
diagonalized; this step defines a fixed basis for propagating the
solutions across the interval. At the midpoint and only at that point
the CC equations are diagonal and hence completely uncoupled. Within
each interval first-order perturbative corrections to the zeroth-order
solutions are calculated. Since the diagonal elements usually change
faster across an intercal than the off-diagonal elements, the interval
is further subdivided into a series of {\it steps}. Within each step
analytic approximations to the solutions of the diagonal interaction
matrix are generated. At the end of each interval the solution is
stabilized by the R-matrix (inverse logarithmic derivative) evaluation
at the interval boundary. The R-matrix is transformed into the local
basis for the next interval, and propagation is continued to the
sector boundary [See Fig. (14)]. Accuracy is maintained by controlling
the interval and step size dynamically. The process of diagonalizing
the interaction potential  is independent of the total energy $E$, so
the eigenvalues and transformation matrices are saved for subsequent
calculations. At the sector boundaries the R-matrix is again
transformed into the local sector basis. Hence overlap integrals of
the surface functions for adjoining sectors must be evaluated.


We have developed a computational procedure that is extremely efficient
after accumulating extensive energy-independent information. Indeed,
we are able evaluate S-matrices in 15 to 60 seconds (Cray-XMP time)
per total scattering energy, even though we have made no made no
attempt to optimize or vectorize the computer code. We {\bf propose}
to streamline this scheme even further.

Results for $\HHH$ are shown in Figs. (15) through (17). These results
will not be discussed in detail in this proposal. However, from Fig. (15)
on sees that our results agree with other published results to within
mutual uncertainties. In Figs. (16) and (17) we present results at higher
energies than heretofore possible. These latter results show a second
resonace at roughly 1.4 eV.

\Figsize 3 truein
\Fig 15. Total reaction probabilities for specific
transitions $H+H_2(\nu_i,j_i)\longrightarrow H_2(\nu_f,j_f)$ for
$J=0$ and $(\nu_i,j_i)=(0,0)$ plotted vs. total energy in eV.  The curves
that are highest near the center of the figure are for $(\nu_f,j_f)=(0,1)$;
the middle curves are for $(0,2)$; and the lower curves are for $(0,0)$.
The present results for the three transitions are the open triangles, asterisks,
and open diamonds, respectively; the results of Kuppermann and Hipes are
the open hexagons, the open squares, and the plus signs; and the results of
Schatz and Kuppermann are the dark squares. \par

\Figsize 3 truein
\Fig 16. a) Total reaction probabilities summed over all
final states for $\HHH$.  The open circles are out of state
$(\nu_i,j_i)=(1,0)$; the open squares are out of state $(1,2)$; and
the open triangles are out of state $(1,4)$.
b) Total reaction probabilities summed over all
final states for $\HHH$.  The open circles are out of state
$(\nu_i,j_i)=(0,0)$; the open squares are out of state $(0,2)$;
the open triangles are out of state $(0,4)$; and the open diamonds
are out of state $(0,6)$. \par

\Figsize 3 truein
\Fig 17. a) Total reaction probabilities for specific
transitions $H+H_2(\nu_i,j_i)\longrightarrow H_2(\nu_f,j_f)$ for
$J=0$ and $(\nu_i,j_i)=(0,0)$ plotted vs. total energy in eV.
The open circles to final state
$(\nu_i,j_i)=(0,0)$; the open squares are to final state $(0,1)$;
the open triangles are to final state $(0,2)$.
b) Total reaction probabilities for specific
transitions $H+H_2(\nu_i,j_i)\longrightarrow H_2(\nu_f,j_f)$ for
$J=0$ and $(\nu_i,j_i)=(0,0)$ plotted vs. total energy in eV.
The open circles to final state
$(\nu_i,j_i)=(0,0)$; the open squares are to final state $(0,1)$;
the open triangles are to final state $(0,2)$; and the open diamonds
are to final state $(0,3)$. \par


\smallskip
\noindent {\bf G. In summary}

We have designed an efficient numerical method for studying reactive
scattering that can produce exact (and approximate) results
for systems and energies which have heretofore been impossible to
study. We conclude that our APH$_{3D}$ reactive scattering
calculations are accurate and viable and ascertain that this along
with other first-class research signals the dawning of a brighter day
in reactive scattering. However our approach has several numerical
advantages over other existing methods.
\item{1.} We have demonstrated that our mesher produces an optimum
nonuniform mesh that assists rapid and accurate determination of
surface functions. The uniform meshes currently being used by other
researchers would produce up to 50,000 nodal points at each $\rho$
value. We are able to produce more accurate surface
functions with less computer time than other existing approaches.
\item{2.} We are not restricted to symmetric mass combinations.
\item{3.} In our method, the most time-consuming numerical
calculations are independent of the total energy. Accordingly, we
are able map out cross sections over extensive energy ranges.
This is necessary to analyze resonances, unambiguously compare with
cross section experiments and produce accurate rate constants.


We have demonstrated that our method is viable and that we agree with
with other published results for $\HHH$ and $\DHH$\ref{7} however, it
is now necessary to critically test, optimize, vectorize and to find the
limitations of the current computer code.
We {\bf propose} to improve the computational proceedure by implementing the
following:
\item{1.} Allowing the elements to deform slightly to give a superior
distribution of points in the classically allowed and weakly forbidden
regions.
\item{2.} Rewrite or optimize the surface function computer code so that
it will execute more efficiently and take advantage of the vectorizability
of the Cray-XMP machines. We estimate that we can streamline this
code to execute 10 times
faster. We will also experiment with other numerical techniques such as
using distributed gaussians for solving the surface function Hamiltonian.
\item{3.} Include photodissociation processes.
\item{4.} Include collision induced dissociation processes.
\item{5.} Include effects of the conical intersections.
\item{6.} Include electronically nonadiabatic processes.
\par
We estimate that 200 Cray hours will be required for code development.

