\documentclass[10pt]{article}

\usepackage{color}

\setlength{\oddsidemargin}{0.0in}
\setlength{\textwidth}{6.5in}
\setlength{\topmargin}{-0.5in}
\setlength{\textheight}{9in}

%===================================================
\usepackage{hhline}
\usepackage{dcolumn}
\newcolumntype{C}{>{$}c<{$}}
\newcolumntype{L}{>{$}l<{$}}
\newcolumntype{R}{>{$}r<{$}}

%================================================
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage[notref,notcite,draft]{showkeys}
\newcommand{\labelsec}{\label}
\newcommand{\labeleq}{\label}
\newcommand{\labelfig}{\label}
\newcommand{\labeltbl}{\label}
\newcommand{\Sec}[1] {Sec. {\ref{#1}} }
\newcommand{\Tab}[1] {Table {\ref{#1}} }
\newcommand{\Eq}[1] {({\ref{#1}})}
\newcommand{\Fig}[1] {Fig. {\ref{#1}} }
\newcommand{\APHphit}{\overline{\Phi}^p_{\tau\Lambda}(\gamma,\chi;\rho_\xi)}
\newcommand{\APHphi}{\Phi^p_{\tau\Lambda}(\gamma,\chi;\rho_\xi)}
\newcommand{\APHphicut}{\Psi^{p}_{\tau\Lambda}(\gamma,\chi;\rho_\xi)}
\newcommand{\I}[2]{
${\begin{pmatrix}
  #1 1 & 0 \\
  0 & #2 1
\end{pmatrix}}$}

\newcommand{\M}[4]{
${\begin{pmatrix}
  #1 \frac{1}{2} & #2 \frac{\sqrt{3}}{2} \\
  #3 \frac{\sqrt{3}}{2} & #4 \frac{1}{2}
\end{pmatrix}}$}

\newcommand{\Rchi}[1]{
$\frac{#1 \pi}{3}+\chi$ }

\newcommand{\Schi}[1]{
$\frac{#1 \pi}{3}-\chi$ }

%=================================================

\begin{document}
\title{Quantum reactive scattering in three dimensions using hyperspherical (APH) coordinates: Periodic Distributed Approximating Functionals (PDAF) method for surface functions}
\author{
Keming Zhang,
Gregory A. Parker\footnote{Supported under National Science Foundation Grant CHE-9710383 and PHY-0100704.},\\
       Department of Physics and Astronomy\\
       University of Oklahoma,
       Norman, OK 73019\\
Donald J. Kouri\footnote{Supported under National Science Foundation Grant CHE-0074311 and R. A. Welsch Foundation Grant E-0608.},\\
       Department of Chemistry and Department of Physics\\
       University of Houston,
       Houston, TX 77204-5641\\
David K. Hoffman\footnote{The Ames Laboratory is supported by the Department of Energy under Contract No 2-7405-ENG82},\\
       Department of Chemistry and Ames Laboratory,\\
       Iowa State University,
       Ames, IA 50011\\
and\\
Srinivasan S. Iyengar\\
       The Henry Eyring Center for Theoretical Chemistry and Department of Chemistry\\
       The University of Utah,
       Salt Lake City, UT 84112}

\date{\today}
\maketitle
%ABSTRACT

Periodic distributed approximating functionals (PDAFs) are proposed and used to obtain a coordinate representation for the \underline Adiabatically Adjusting \underline Principal Axis \underline Hyperspherical (APH) coordinates kinetic energy operator. The approach is tested and accurate results for adiabatic surface functions for the reaction $\text{F+H}_2\rightarrow \text{HF+H}$ are calculated and compared to those of some existing methods.

%INTRODUCTION
\section{Introduction}

In the hyperspherical coordinates formulation of triatomic reactive scattering, the total wave function is represented using the Wigner D-functions for the three Euler angles describing the spatial orientation of the three particle plane, and basis functions for the two internal hyperspherical angles. The dependence on the hyperradius, $\rho$ , is then determined by propagating the set of coupled channel(CC) differential equations from small $\rho$, where the solutions must be regular, to large $\rho$ where they can be projected onto the arrangement channels to determine the scattering matrix. In the methods using hyperspherical coordinates that treat all particles symmetrically\cite{Launay1991, Lepetit1991, PackParker1987, PackParker1989, Bacic1990, Kress1990, Launay1989, Linderberg1989}, one obtains ``surface functions,'' the basis functions of the two hyperangles, by solving a two dimensional (2D) Schr\"odinger equation. This equation, which is discussed in more detail later in this paper, depends parametrically on $\rho$ and must be solved at many values of $\rho$. In addition, a large number of surface functions must be obtained at each $\rho$, and hence efficient computational procedures are essential.

The first accurate 3D reactive scattering calculations employing hyperspherical coordinates used finite element methods (FEMs)\cite{PackParker1987, PackParker1989, Bacic1990, Kuppermann1989, Linderberg1989, Lagana1991, Parker1989} to solve for the surface functions. Although these FEMs give fairly accurate results, they are inefficient and not robust. Another commonly used method involves the use of the discrete variable representation (DVR)\cite{Bacic1990, Kress1990, Kress1989, Kress1990A, Kress1991}, which is most efficient at small $\rho$ where the surface functions are delocalized. At larger values of $\rho$, where the surface functions are highly localized, the DVR points still cover the whole space, making it much less efficient. In a few cases the DVR is even more expensive than the FEM because of the need for many grid points in a small, localized region. Other methods such as the finite basis representation (FBR)\cite{Launay1991, Lepetit1991, Launay1989} of Launay and LeDourneuf, and the method of Wolniewicz and Hinze[29] are also efficient only at large $\rho$. The analytic basis method (ABM) uses primitive basis functions centered in the arrangement channels, and provides a very compact representation and thus is quite effective at large $\rho$, but is inefficient and tends to diverge at small $\rho$.

In this paper, we present the Periodic Distributed Approximating Functional (PDAF) method, an approach that is efficient at both small and large values of $\rho$. The PDAF method is similar to the Symmetry-Adapted Hermite-Distributed Approximating Functional (SA-HDAF) approach of Iyengar, Parker, Kouri and Hoffman \cite{Iyengar1999}, but differs importantly by involving only real symmetric matrices. In addition, in the present approach only the surface functions are obtained, while in the SA-HDAF approach \cite{Iyengar1999} the full 3D wavefunction was obtained directly, by using an iterative procedure. The symmetry adaptation of the PDAF is carried  out here in a similar, but simpler, fashion as compared Ref. \cite{Iyengar1999}. Here we employ, the Distributed Approximating Functional (DAF) concept, but the PDAF differs from all previous DAF's. However, like other DAF's it is both accurate and efficient. A more detailed exposition of the connections between the PDAF and the Christoffel-Darbeaux formula for DAFs based on orthogonal polynomials \cite{vardaf} will follow \cite{PDAF-vardaf}. The sequential diagonalization-truncation technique\cite{WithinDVR} is employed to project the large-size Hamiltonian matrix into a smaller matrix using a projection matrix which is obtained by solving a one-dimensional eigensystem, thus significantly reducingthe memory requirements and the computation time.

This paper is organized as follows. In \Sec{sec:PDAF} we introduce the PDAFs and derive their formulae. In \Sec{sec:HReduction} the ro-vibrational triatomic Hamiltonian in the APH coordinates system is presented and the symmetrization and reductions of the Hamiltonian are illustrated. The PDAF approach is then tested in \Sec{sec:Test}. Surface functions for the $\text{FH}_2$ scattering system are computed, and the eigen-energies and the matrix elements are calculated and compared to those of existing methods (FEM, ABM and DVR) in \Sec{sec:CalcFH2}. In \Sec{sec:Conclusion} we present our conclusions.

%SECTION 2
\section{Periodic Distributed Approximating Function (PDAF)}\labelsec{sec:PDAF}
From the definition of the Dirac Delta function we know that
\begin{equation}
    f(x)=\int ^{\infty}_{-\infty} \delta(x-x')f(x')dx'
\end{equation}
for any continuous function $f(x)$. Consider periodic functions, $f_p(x)$, for convenience scaled to have period $2\pi$.
Expressing the integration range as an infinite sum of segments of length $2\pi$ we can write
\begin{equation}
    f_p(x)=\sum^\infty_{m=-\infty} \int ^{2\pi}_{0} \delta(x-x'-2m\pi)f_p(x'+2m\pi)dx'.
\end{equation}
Interchanging the integral and sum, and using the periodicity of $f_p(x)$, we have
\begin{equation}
    f_p(x)=\int ^{2\pi}_{0} \sum^\infty_{m=-\infty}  \delta(x-x'-2m\pi)f_p(x')dx'
\end{equation}
or
\begin{equation}\labeleq{eq:PDAFint0}
    f_p(x)=\int ^{2\pi}_{0} \delta_p(x-x')f_p(x')dx'
\end{equation}
where
\begin{equation}
    \delta_p(x-x')=\sum^\infty_{m=-\infty}  \delta(x-x'-2m\pi)
\end{equation}
is the periodic delta function. Since $\delta(x)$ is an even
function, it's easy to see that $\delta_p(x)$ is even and is also
periodic with a period of $2\pi$. The $k^{\text{th}}$ derivative of
$f_p(x)$ is then given by
\begin{equation}\labeleq{eq:PDAFint}
    f^{(k)}_p(x)=\int ^{2\pi}_{0} \delta_p^{(k)}(x-x')f_p(x')dx'.
\end{equation}

Expanding $f_p(x)$ and $\delta_p(x-x^\prime)$ in a Fourier series,
\begin{equation}
    \delta_p(x-x^\prime)=\frac{a_0}{2}+\sum^{\infty}_{n=1}a_n\cos{n(x-x^\prime)}
\end{equation}
\begin{equation}
    f_p(x)=\frac{b_0}{2}+\sum^{\infty}_{n=1} \left( b_n\cos{n x}+c_n\sin{nx} \right),
\end{equation}
and using the identities
\begin{eqnarray}
    \int^{2\pi}_{0}\cos{m(x-x')}\cos{nx'}dx'&=&\left\{
    \begin{array}{ll}
      \delta_{m n}\pi\cos{nx}&   n\neq 0\\
      \delta_{m n}2\pi        &   n=0
    \end{array}
     \right.\\
    \int^{2\pi}_{0}\cos{m(x-x')}\sin{nx'}dx'&=&\delta_{m n}\pi\sin{nx}
\end{eqnarray}
we obtain 
\begin{equation}
    a_n=\frac{1}{\pi},\;\;\;\;\;n=0,1,2,\ldots
\end{equation}
Thus, we obtain
\begin{equation}\labeleq{eq:PDAFexact}
    \delta_{p}(x-x^\prime)=\frac{1}{\pi}\left(  \frac{1}{2}+\sum_{n=1}^{\infty}{\cos{n(x-x^\prime)}} \right),
\end{equation}
for the periodic delta function.

We define the partial sum of the periodic delta function
to be the continuous Periodic Distributed Approximating Function (PDAF)
\begin{equation}\labeleq{eq:PDAF}
    \delta_{p,M}(x-x^\prime)\equiv\frac{1}{\pi}\left[  \frac{1}{2}+\sum_{n=1}^{M}{\cos{n(x-x^\prime)}} \right]
    =\frac{1}{2\pi}\left\lbrace\frac{\cos[(M-1)(x-x^\prime)]-\cos[M(x-x^\prime)]}{1-\cos(x-x^\prime)}\right\rbrace
\end{equation}
which is the basic result used in this work. The last expression is obtained by noting that $\delta_{p,M}(x-x^\prime)$ is the real part of a sum of exponentials which can be written as a geometric sum and thus done analytically. By definition, as $M$ increases the PDAF approaches the
periodic delta function, {\it i.e.}
\begin{equation}
\delta_{p}(x-x^\prime)=\lim_{M\rightarrow\infty}\delta_{p,M}(x-x^\prime).
\end{equation}

%\marginpar{\tiny Don Kouri may like to add a paragraph or two here, after Eq. (14).}

One may obtain the fully discretized PDAF by approximating the integral over $x'$ in \Eq{eq:PDAFint} using a trapezoidal
quadrature
\begin{equation}\label{eq:PDAFintd}
    f^{(k)}_{p,M}(x)=\sum^N_{j=1}\delta^{(k)}_{p,M}(x-x_j)f_p(x_j)\Delta x
\end{equation}
where $N$ is the number of grid points, $\Delta x =\frac{2\pi}{N}$ and $x_j=(j-\frac{1}{2})\Delta x$. 
This choice of quadrature points is particularly useful if we are solving differential equations with singular points at the two ends ($x=0$ and $x=2\pi$). 
As the grid points are fixed, $\delta_{p,M}^{(k)}(x)$ acts like the discretized $k^{th}$ derivative operator.
One can also discretize $x$ using the same grid points as used in the numerical quadrature.
If we treat $f_{p,M}(x)$ as a column vector $[f_{j}=f(x_j)]$, we have
\begin{equation}
    f_{i,M}^{(k)}=\sum^{N}_{i=1}D^{(k)}_{ij}f_{j}
\end{equation}
where
\begin{eqnarray} \label{eq:PDAF1}
D^{(k)}_{ij}=\delta_{p,M}^{(k)}(x_i-x_j)\Delta x=
\begin{cases}
(-1)^{\frac{k+1}{2}}\frac{\Delta x}{\pi}\sum^{M}_{n=1}n^k \sin{n[(x_i-x_j)]x}    &\text{$k$ odd }\\ \\
\frac{\Delta x}{2\pi}\delta_{0,k}+(-1)^{\frac{k  }{2}}\frac{1}{\pi}\sum^{M}_{n=1}n^k \cos{n[(x_i-x_j)]}    &\text{$k$ even}
\end{cases}
\end{eqnarray}
is the $k^{th}$ derivative operator in matrix form.
The differential operators, $D^{(k)}$ are periodic Toeplitz matrices which depend on only one parameter $M$.
The column vector $f_{i,M}^{(k)}$ is the PDAF approximate to the $k^{th}$ derivative of $f_p(x)$, {\it i.e.}
\begin{equation}
f_{i,M}^{(k)}\approx \frac{d^k}{dx^k}f_p(x)\Big\vert_{x=x_i}.
\end{equation}
As $M$ increases, the PDAF gives an increasingly accurate representation of the Dirac delta function. However, for small $M$ the
integrand of \Eq{eq:PDAFint} is smoother and the trapezoidal rule is more accurate. 
Hence, we seek an optimal compromise
value of $M$. Since,
the integrand involves a Fourier function, the theory of Gaussian quadrature
suggests a relation of $M=N/2$.
Through numerical experimentation (see \Sec{sec:Test}) we find that $M=N/2$,
gives the most accurate results, where $N$ is the number of grid points in the $2\pi$ range.
In this paper, we will use $D^{(0)}$, $D^{(1)}$ and $D^{(2)}$,
consequently we need $\delta_{p,M}^{(0)}$, $\delta_{p,M}^{(1)}$ and $\delta_{p,M}^{(2)}$.
They can be obtained by simply differentiating \Eq{eq:PDAF}. 

We see that $\delta^{(k)}_{p,M}(x)$ are even(odd) periodic functions when $k$ is even(odd). The functions,
$\delta^{(0)}_{p,M}(x)$, $\delta^{(1)}_{p,M}(x)$, and $\delta^{(2)}_{p,M}(x)$ which are the ones used in this paper,
are shown in \Fig{fig:PDAFs}. In both panels of \Fig{fig:PDAFs} $N=20$, in the left panel $M=5$ and in the right panel $M=N/2=10$. One can readily see the symmetry of each PDAF.
Comparing the left and right panels in \Fig{fig:PDAFs} we see that increasing $M$ makes the
PDAFs in the right panel sharply peaked and a better approximate representation of the periodic Dirac Delta function and its derivatives.

The equations obtained in this section for the PDAF, namely \Eq{eq:PDAF},
\Eq{eq:D} and \Eq{eq:PDAF1}, may be compared with those obtained using the
SA-HDAF \cite{Iyengar1999} for functions with periodic symmetry, namely Eq. (5)
and associated expressions in Ref. \cite{Iyengar1999}. It is clear that both
representations are Toeplitz-symmetric, however, the expression in
Ref. \cite{Iyengar1999} has the additional characteristic of being banded, which
is not the case in the present formulation. However, as we will see further in
this paper, this does not present any problems in the current implementation and
in fact helps in enhancing the accuracy of the procedure introduced here.

\begin{figure}
 \begin{center}
 \includegraphics[width=4in,height=3in]{graph/PDAF.eps}
 \end{center}
 \caption{$\delta^{(k)}_{p,M}(x)$ for different $M$, where the number of grid points ($N=20$) is used in both panels.  In the left panel $M=5$ and in the left panel and $M=N/2=10$.  Note $\delta_{p,M}^{(0)}$ and $\delta_{p,M}^{(2)}$ are symmetric about $x=0$, while $\delta_{p,M}^{(1)}$ is anti-symmetric. Increasing $M$ makes the PDAFs sharply peaked and a better representation of the respective delta function and its derivatives.}\labelfig{fig:PDAFs}
\end{figure}





\section{Symmetrization and Reductions of the Triatomic Hamiltonian in APH Coordinates}\labelsec{sec:HReduction}
\subsection{APH Hamiltonian}
The detailed reactive scattering theory formulated in adiabatically adjusting principle axes hyperspherical (APH) coordinates has been presented
previously\cite{PackParker1987}, and we repeat only the essentials here.
In this approach, one needs sector adiabatic basis functions $\Phi^p_{\tau\Lambda}$ of the APH hyperangles,
and in this work we choose surface functions, $\Phi^p_{\tau\Lambda}$ defined by the equation
\begin{equation}\labeleq{eq:Schrodinger}
    H\Phi^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)=\mathcal{E}^p_{\tau\Lambda}(\rho_\xi)\Phi^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)
\end{equation}
where $H$ is a portion of the full Hamiltonian omitting parts of the orbital angular momentum. Thus we take
\begin{equation}\labeleq{eq:H}
    H= -\frac{\hbar^2}{2\mu\rho_\xi^2}\left[\frac{4}{\sin{2\theta}} \frac{\partial}{\partial\theta}
         \sin{2\theta}\frac{\partial}{\partial\theta}+\frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\chi^2}\right]
     +\frac{15\hbar^2}{8\mu\rho_\xi^2}
       +C\hbar^2\Lambda^2
       +V(\rho_\xi,\theta,\chi)
\end{equation}
The first term in the Hamiltonian $H$ is the ``hyperspherical'' part of the kinetic energy operator, and
\begin{equation}
C=\frac{1}{\mu\rho_\xi^2(1-\sin\theta)}
\end{equation}
is part of the centrifugal potential. The potential energy $V$ used here is the complete potential energy surface (PES) of Brown \textit{et al}.\cite{FH2_T5A}, and the $\mathcal{E}^p_{\tau,\Lambda}(\rho_\xi)$ are the eigenenergies of $H$ at $\rho_\xi$.
The variable $\theta$ is the APH bending angle; its range is $0\leq\theta\leq\pi/2$, with $\pi/2$ describing linear configurations and 0 describing triangular symmetric top configurations. $\chi$ is the APH kinematic angle measured from the ``incident'' arrangement channel; it measures motion between arrangement channels, and its range is $-\pi\leq\chi\leq\pi$. The angles $\theta$ and $\chi$ cover the upper half of the surface of an internal coordinate sphere which we loosely call the ``hypersphere''. (More precisely, the surface of the hypersphere is the 5D space covered by $\theta$, $\chi$, and the three Euler angles which describe the orientation of the principal axes in space. In fact, the hypersphere may be decomposed into two commuting subgroups, $O(2)$ and $O(3)$, and this aspect has been exploited in computing surface functions with fixed total angular momentum. \cite{PackParker1987}) As one can see from \Eq{eq:Schrodinger}, the surface functions $\Phi^p_{\tau,\Lambda}$ and eigenenergies $\mathcal{E}^p_{\tau,\Lambda}$
depend parametrically on the hyperradius $\rho_\xi$. They are needed at a set of $\rho$ values
$\{\rho_\xi\}, \xi=1,2,...,n_\rho$, and are used as adiabatic basis functions for expanding the full
wavefunction in each sector where
${(\rho_\xi+\rho_{\xi-1})}/{2}\le\rho\le{(\rho_{\xi+1}+\rho_\xi)}/{2}$ for sector $\xi$.
As shown elsewhere\cite{PackParker1987} this adiabatic-by-sector (or sector-Adiabatic) expansion of the wavefunction gives rise to a set of
coupled second order differential equations. These Coupled-Channel (CC) equations must be numerically
integrated from a small value of the hyperradius $\rho$ where the full wavefunction is zero to a
large value of the hyperradius where asymptotic boundary conditions are applied.
The term, $15\hbar^2/8\mu\rho^2_\xi$, in \Eq{eq:H} results from removing the first derivative terms
from the coupled channel (CC) equations; it is a constant in this equation and can be folded
into the $\mathcal{E}^p_{\tau,\Lambda}$ if desired. 
In the numerical calculations we use atomic units throughout and therefore $\hbar=1$. 
The quantity $\mu$ is the three-body reduced mass of the system arising from using mass-scaled coordinates.
The three quantum numbers ($\Lambda, p, \tau$) labeling $\mathcal{E}^p_{\tau,\Lambda}$ and $\Phi^p_{\tau,\Lambda}$ are:
$\Lambda$ is the component of the total angular momentum along the APH body-frame (BF) $z$-axis
(the axis of least inertia of the three-body system),
$p$ is the parity quantum number with $p=0$ or $1$, the parity of
$\Phi^p_{\tau,\Lambda}$ under the parity transformation $\chi\rightarrow\chi\pm\pi$ is $(-1)^p$,
and $\tau=1,2,...,n_\Phi$, which indexes the solutions in order of increasing energy.\\

Equation \Eq{eq:H} differs from equation (164) of Reference\cite{PackParker1987} slightly
because it omits a rotational term of the form $\frac{1}{2}(A+B)\hbar^2[J(J+1)-\Lambda^2]$ as mentioned above.
As pointed out by Launay and LeDourneuf \cite{Launay1991},
this gives surface functions $\Phi^p_{\tau,\Lambda}$ which are independent of the total angular momentum $J$,
so  many fewer surface functions must be calculated.
The omitted term is easily included in the CC equations along with the remaining Coriolis and asymmetric top terms.
This surface function basis is expected to produce rapid convergence of the CC expansion to the
exact solution provided triangular symmetric top $(\theta=0)$ configurations are unimportant,
which is the case for many reactions.\\

The full wavefunction must be continuous and regular everywhere. This requires that
$\Phi^p_{\tau,\Lambda}$  must also be a continuous function of $\chi$ at $-\pi$ and $\pi$ and regular everywhere.
For systems with two or three identical atoms, the surface functions have other symmetries in addition to the parity, $p$, already defined and these symmetries will be exploited. The surface functions are real and normalized according to
\begin{equation}\label{eq:norm}
\int^{\pi}_{-\pi}d\chi\int_0^{\pi/2}  \Phi^{p'}_{\tau'\Lambda}(\theta,\chi;\rho_\xi)\Phi^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)\sin{2\theta} d\theta
=\delta_{\tau^\prime \tau}\delta_{p^\prime p}.
\end{equation}
\\

We now define
\begin{eqnarray}
    H_\theta&=&-\frac{4}{\sin{2\theta}} \frac{\partial}{\partial\theta} \sin 2\theta\frac{\partial}{\partial\theta}\labeleq{eq:Htheta} \\
    H_\chi&=&-\frac{\partial^2}{\partial\chi^2}  \\
    H_V&=&\frac{15\hbar^2}{8\mu\rho_\xi^2}+C\hbar^2\Lambda^2+V(\rho_\xi,\theta,\chi)
\end{eqnarray}
so that the surface function Hamiltonian, \Eq{eq:H}, is
\begin{equation}\labeleq{eq:Hsimp}
H=\frac{\hbar^2}{2\mu\rho_\xi^2}\left[H_\theta+\frac{1}{\sin^2\theta}H_\chi\right]+H_V.
\end{equation}
The discretized Hamiltonian is a matrix operator. 
In the rest of this paper, we use the term Hamiltonian to refer to this matrix.

\subsection{Symmetrization of APH Hamiltonian}\label{sub:symAPH}
The Hamiltonian $H$ in \Eq{eq:Hsimp} is real but it is not symmetric, because $H_\theta$ is not symmetric. 
Therefore, if we use this form, we have to solve a non-symmetric matrix eigensystem requiring a large amount of memory and CPU time. 
It is therefore critical in the computation of APH surface functions to symmetrize $H$. 
The matrix $H_\theta$ is the non-symmetric part of $H$, and it is not periodic. 
We first extend $H_\theta$ to a periodic form and symmetrize it. Then the result is extended to symmetrize $H$.\\

Although $H_\theta$ initially defined on $0\leq\theta\leq\frac{\pi}{2}$, 
we note that it is invariant under the transformation $\theta\rightarrow\pi-\theta$ and hence we can expand the domain
of $\psi(\theta)$ to the full real space by defining,
\begin{eqnarray}
  \psi(\theta)&=&\psi(\pi-\theta), \;\;\;\frac{\pi}{2}\leq\theta\leq\pi \label{eq:A1-2}\\
  \psi(\theta+k\pi)&=&\psi(\theta), \;\;\;0\leq\theta\leq\pi\;\mathrm{and}\; k=0,\pm1,\pm2,... \label{eq:A1-3}
\end{eqnarray}
Now $\psi(\theta)$ is a periodic function with a periodicity of $\pi$. Note that
\begin{eqnarray}
  H_\theta(\theta)&=&H(\pi-\theta)\\
  H_\theta(\theta+k\pi)&=&H(\theta), \;\;\; k=0,\pm1,\pm2,...
\end{eqnarray}
One can see that the extended $\psi(\theta)$ satisfies \Eq{eq:eigen-theta} for any $\theta$, i.e. the domain of $H_\theta$ is also extended to the full real space although it keeps its original form.\\

We comment that similar expressions for $H_\theta$ were obtained in Ref.
\cite{Iyengar1999} by using a symmetry adaptation procedure,
wherein $H_\theta$ was projected on the right side onto a function
belonging to the $A_1$ irreducible representation of the point
group $C_{2v}$, and on the left side onto a function belonging to
the $A_2$ irreducible representation of the point group $C_{2v}$.
The choice of the $A_1$ irreducible representation conforms with
the symmetry requirements presented in Eqs. \Eq{eq:A1-1},
\Eq{eq:A1-2} and \Eq{eq:A1-3} above. The $A_2$ irreducible
representation was chosen due to the fact that the first
derivative of an $A_1$ function is an $A_2$ function in $C_{2v}$.
The resultant representation for $H_\theta$  in Ref.
\cite{Iyengar1999} is, however, not symmetric.

Here, to obtain a symmetric representation for $H_\theta$, we transform it in two steps. First we introduce a continuous mapping function,
\begin{equation}
    \theta(\gamma)=\frac{\pi}{4}[1-\cos(\gamma+2k\pi)]-k\pi, \;\;\;k=0,\pm1,\pm2,...
\end{equation}
We note that the first derivative of this mapping is also continuous, but its second derivative is discontinuous at $\theta=k\pi/2$. This problem is handled by the mapping function itself: if
we use a uniform grid for $\gamma$, it is equivalent to using a non-uniform grid for $\theta$ which is dense near the singularities ($\theta=k\pi/2$) and sparse otherwise, as shown in \Fig{fig:theta-gamma}. Since more grid points are used in the $\theta=k\pi/2$ region, high accuracy can be obtained despite the discontinuity.\\

Substituting $\theta$ with $\gamma$ in \Eq{eq:Htheta}, we obtain,
\begin{equation}
    H_\theta(\gamma)=-\frac{64}{\pi^2} \frac{1}{\sin 2\theta(\gamma)\sin \gamma} \frac{\partial}{\partial \gamma}
    \frac{\sin 2\theta(\gamma)}{\sin \gamma}\frac{\partial}{\partial \gamma}
\end{equation}
$H_\theta(\gamma)$ is very similar to $H_\theta(\theta)$ except that it has a periodicity of $2\pi$.
\begin{figure}
 \begin{center}
 \includegraphics[width=4in,height=3in]{graph/theta-gamma.eps}
 \end{center}
 \caption{Grid mapping function $\theta(\gamma)$. The uniform grid in $\theta$ corresponds to the non-uniform grid for $\theta$ which is dense near the singularities( $\theta=k\pi/2$) and sparse otherwise.}\labelfig{fig:theta-gamma}
\end{figure}\\

In the second step, we define the transformation function
\begin{equation}
    T=\sqrt{\sin 2\theta \sin \gamma}
\end{equation}
We also observe that
\begin{eqnarray}
    T^\dag&=&T\\
    T^{-1}&=&\frac{1}{\sqrt {\sin 2\theta \sin \gamma} }
\end{eqnarray}
and it is important to note that the singularities in $T^{-1}$ are avoided by our choice of quadrature (see \Sec{sec:PDAF}). Applying $T^\dag$ to \Eq{eq:eigen-theta} and inserting the identity $T^{-1}T$ between $H_\theta$ and $\psi(\theta)$, we obtain
\begin{equation}\labeleq{eq:SymThetaEigenEquatioin}
    H_\gamma T \psi(\theta)=\lambda T \psi(\theta)
\end{equation}
where
\begin{eqnarray}
    H_\gamma&=&T^\dag H_\theta(\gamma) T^{-1} \nonumber\\
    &=&-\frac{64}{\pi^2}\frac{1}{\sqrt{\sin 2\theta(\gamma)\sin\gamma}}\frac{\partial}{\partial \gamma}\frac{\sin 2\theta(\gamma)}{\sin \gamma}\frac{\partial}{\partial \gamma}        \frac{1}{\sqrt{\sin 2\theta(\gamma)\sin \gamma}}
\end{eqnarray}
is a symmetric Hamiltonian.\\

The eigenfunction of $H_\gamma$ is
\begin{equation}\label{eq:psitilde}
  \tilde{\psi}(\gamma)=T(\gamma)\psi(\theta(\gamma))
\end{equation}
which is always zero at $\gamma=0, \pi,$ or $2\pi$. It is not hard to see that extending $\psi(\theta)$ from $0\leq\theta\leq\frac{\pi}{2}$ to $0\leq\theta\leq\pi$ introduces undesired degenerate eigenvalues. If we reduce matrix $H_\gamma$ from $0\leq\theta\leq2\pi$ to range $0\leq\theta\leq\pi$ using the symmetry of $\psi$, we can eliminate those undesired degenerate eigenvalues. We use a grid similar to the one used for PDAF;
\begin{equation}\label{}
  \gamma_j=(j-\frac{1}{2})\Delta\gamma, \;\;\;j=1, 2N_\gamma
\end{equation}
where $2N_\gamma$ is the number of gird points, and $\Delta\gamma=\frac{2\pi}{2N_\gamma}=\frac{\pi}{N_\gamma}$ is the distance between two consecutive grid points. For simplicity, we use a subscript $j$ to denote a term evaluated at $\gamma_j$, and we also define the following terms,
\begin{eqnarray}
    t_j&=&\frac{1}{\sqrt{\sin 2\theta(\gamma_j)\sin\gamma_j}}\\
    s_j&=&\frac{\sin 2\theta(\gamma_j)}{\sin \gamma_j}
\end{eqnarray}
Applying the derivative operator $D^{(1)}$ in \Eq{eq:D} (substituting $x$ with $\gamma$), the eigen-equation \Eq{eq:SymThetaEigenEquatioin} can be written as
\begin{eqnarray}
  \lambda\tilde{\psi}_i&=&\sum_{j=1}^{2N}\sum_{k=1}^{2N}-\frac{64}{\pi^2}t_i D^{(1)}_{ik}s_k D^{(1)}_{kj}t_j \tilde{\psi}_j \nonumber\\
    &=&(\sum_{j=1}^{N}+\sum_{j=N+1}^{2N})\sum_{k=1}^{2N}-\frac{64}{\pi^2}t_i D^{(1)}_{ik}s_k D^{(1)}_{kj}t_j \tilde{\psi}_j \nonumber\\
    &=&\sum_{j=1}^{N}\sum_{k=1}^{2N}-\frac{64}{\pi^2}t_i  D^{(1)}_{ik}s_k (D^{(1)}_{kj}+ D^{(1)}_{k\;2N-j+1})t_j \tilde{\psi}_j \nonumber\\
\label{eq:eigenHrgamma}&=&\sum_{j=1}^{N}{H^r_\gamma}_{ij} \tilde{\psi}_j
\end{eqnarray}
where
\begin{equation}\label{eq:Hrgamma}
    {H^r_\gamma}_{ij}=\sum_{k=1}^{2N}-\frac{64}{\pi^2}t_i  D^{(1)}_{ik}s_k (D^{(1)}_{kj}+ D^{(1)}_{k\;2N-j+1})t_j
\end{equation}
In the above derivation, $\tilde{\psi}_j=\tilde{\psi}_{2N-j+1}$, $t_j=t_{2N-j+1}$, $s_j=s_{2N-j+1}$ and the periodicity of $D^{(1)}$ is employed. One should note that in \Eq{eq:eigenHrgamma} $1\leq i\leq N_\gamma$ and $1\leq j\leq N_\gamma$, which means $0\leq\gamma\leq\pi$, so $H^r_\gamma$ is reduced to the original domain. $H^r_\gamma$ is the desired symmetric matrix form for $H_\theta$. It is, however, interesting to note that the
bracketted quantity, $(D^{(1)}_{kj}+ D^{(1)}_{k\;2N-j+1})$, in \Eq{eq:Hrgamma}
is not an anti-symmetric matrix, in contrast to Ref. \cite{Iyengar1999}. The full
$H^r_\gamma$ presented in \Eq{eq:Hrgamma} is symmetric by contstruction.

Similarly we can transform and discretize $H$ to obtain its symmetric form. Using \Eq{eq:Hsimp} the Schr\"odinger equation \Eq{eq:Schrodinger} reads
\begin{equation}\labeleq{}
\left\{\frac{\hbar^2}{2\mu\rho_\xi^2}\left[H_\theta+\frac{1}{\sin^2\theta}H_\chi\right]+H_V\right\}
\Phi^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)=\mathcal{E}^p_{\tau\Lambda}(\rho_\xi)\Phi^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)
\end{equation}

We multiply by $T^\dag$ from the left on both sides of the above equation, and insert the identity $T^{-1}T$, to obtain
\begin{equation}\labeleq{eq:sym-Schrodinger}
H^s\overline{\Phi}^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)=\mathcal{E}^p_{\tau\Lambda}(\rho_\xi)\overline{\Phi}^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)
\end{equation}
where
\begin{equation}\labeleq{eq:sym-H}
H^s=\frac{\hbar^2}{2\mu\rho_\xi^2}\left[H^r_\gamma+\frac{1}{\sin^2\theta}H_\chi\right]+H_V
\end{equation}
and
\begin{equation}
\overline{\Phi}^p_{\tau\Lambda}(\gamma,\chi;\rho_\xi)=T\Phi^p_{\tau\Lambda}(\theta,\chi;\rho_\xi)
\end{equation}
The independence of the variables $\gamma$(or $\theta$) and $\chi$ is used in the derivation. After discretization of the APH Hamiltonian, $H^s$ is symmetric, because $H^r_\gamma$ is symmetric, $H_V$ is diagonal, and $H_\chi$ is obviously symmetric.

\subsection{Reduction of the $H_\chi$ using point group symmetry}
If the Hamiltonian $H^s$ of a system commutes with a point group $G$, i.e. the potential $H_V$ is symmetric under the operations of $G$, we can reduce the size of the Hamiltonian matrix by transforming the Hamiltonian into the irreducible representations of $G$.

Let $G=\{R\}$, where $\{R\}$ is a set of symmetry operations such as rotations and reflections. The order of $G$ is $h$, i.e, $G$ contains $h$ symmetry operations, $R$. Suppose $H_\chi$ satisfies the following eigen-equation,
\begin{equation}\label{eq:eigen-Hchi}
  H_\chi \phi^{[z,\kappa]}(\chi)= \lambda \phi^{[z,\kappa]}(\chi)
\end{equation}
where $\phi^{[z,\kappa]}$ transform according the $\kappa$-th column of the $z$-th irreducible representation of $G$, i.e. $\phi^{[z,\kappa]}$ is the $\kappa$-th basis function for the $z$-th irreducible representation, and $\lambda$ is the corresponding eigenvalue. We will discretize \Eq{eq:eigen-Hchi} using a uniform
\begin{equation}\label{eq:chischeme}
  \chi_j=(j-1/2)\Delta\chi,\;\;\;j=1, hN_\chi
\end{equation}
where $hN_\chi$ is the number of grid points, $\Delta\chi=\frac{2\pi}{hN_\chi}$ is the spacing between consecutive grid points, and clearly $N_\chi$ is the number of grid points in interval $[0, \frac{2\pi}{h}]$. Then \Eq{eq:eigen-Hchi} can be written as,
\begin{equation}\label{eq:discHchi}
  \sum_{j=1}^{hN_\chi}{H_\chi}(\chi_i,\chi_j)\phi^{[z,\kappa]}(\chi_j)= \lambda \phi^{[z,\kappa]}({\chi_i}), \;\;\;\kappa=1,l_z
\end{equation}
where $l_z$ is the dimension of the $z$-th irreducible representation. The full range of $\chi$ can be generated by applying each symmetry operation $R$ on the range $[1, N\chi]$. We can then write \Eq{eq:discHchi} as
\begin{equation}\label{eq:RdiscHchi}
  \sum_{j=1}^{N_\chi}\sum_{R}{H_\chi}(\chi_i,R\chi_j)\phi^{[z,\kappa]}(R\chi_j)= \lambda \phi^{[z,\kappa]}({\chi_i}), \;\;\;\kappa=1,l_z
\end{equation}
From the well-known relations \cite{Weissbluth1978},
\begin{equation}\label{eq:PR-1}
  \phi^{[z,\kappa]}(R^{-1}\chi)=\sum_{\kappa'} \phi^{[z,\kappa']}(\chi)\Gamma^{[j]}_{\kappa',\kappa}(R)
\end{equation}
where $\Gamma^{[z]}_{\kappa',\kappa}(R)$ is the $(\kappa',\kappa)$-th element of the $z$-th irreducible representation matrix of $R$, we have,
\begin{eqnarray}\label{eq:PR}
  \phi^{[z,\kappa]}(R\chi)&=& \sum_{\kappa'}          \phi^{[z,\kappa']}(\chi)\Gamma^{[z]}_{\kappa',\kappa}(R^{-1})\nonumber\\
&=& \sum_{\kappa'}\phi^{[z,\kappa']}(\chi)\Gamma^{[z]\dag}_{\kappa',\kappa}(R)\nonumber\\
&=& \sum_{\kappa'}\phi^{[z,\kappa']}(\chi)\Gamma^{[z]*}_{\kappa,\kappa'}(R)
\end{eqnarray}
The following obvious properties of $\Gamma^{[z]}_{\kappa',\kappa}(R)$ are used in the above derivation,
\begin{equation}\label{eq:GammaProperty}
    \Gamma^{[z]}(R^{-1})={{\Gamma^{[z]}}^\dag(R)=\Gamma^{[z]}}^{-1}(R)
\end{equation}
Substituting \Eq{eq:PR} into \Eq{eq:RdiscHchi}, we obtain,
\begin{equation}\label{eq:red-Hchi}
  \sum_{j=1}^{N_\chi}\sum_{R}{H_\chi}(\chi_i,R\chi_j)  \sum_{\kappa'}\phi^{[z,\kappa']}(\chi_j)\Gamma^{[j]*}_{\kappa,\kappa'}(R)
    = \lambda \phi^{[z,\kappa]}({\chi_i}), \;\;\;\kappa=1,l_z
\end{equation}
The expression \Eq{eq:red-Hchi} is a set of coupled equations which can be written concisely in a matrix form. For triatomic reactions, the relevant irreducible representations are either one-dimensional or two-dimensional. For the one-dimensional case, $\kappa=\kappa'=1$, so we simply omit them. Hence we have,
\begin{equation}\label{eq:red-1D-Hchi}
  \sum_{j=1}^{N_\chi} H^r_\chi(\chi_i,\chi_j)\phi^{[z]}(\chi_j)
    = \lambda \phi^{[z]}({\chi_i})
\end{equation}
where
\begin{equation}\label{eq:redHchi-1D}
  H^r_\chi(\chi_i,\chi_j)=\sum_{R}{H_\chi}(\chi_i,R\chi_j)\Gamma^{[z]*}(R)
\end{equation}
is the reduced $H_\chi$ in a one-dimensional irreducible representation.\\

For the two-dimensional case, we have,
\begin{equation}\label{eq:red-2D-Hchi}
  \sum_{j=1}^{N_\chi} H^r_\chi(\chi_i,\chi_j)
   \begin{pmatrix}
    \phi^{[z,1](\chi_j)}\\
    \phi^{[z,2](\chi_j)}\
  \end{pmatrix}
    = \lambda \begin{pmatrix}
    \phi^{[z,1](\chi_i)}\\
    \phi^{[z,2](\chi_i)}\
  \end{pmatrix}
\end{equation}
where
\begin{equation}\label{eq:redHchi-2D}
  H^r_\chi(\chi_i,\chi_j)=\sum_{R}H_\chi(\chi_i,R\chi_j)\Gamma^{[z]*}(R)
\end{equation}
is the reduced $H_\chi$ in a two-dimensional irreducible representation. It has the same form as in the one-dimensional case, but it is a $2\times2$ matrix. $H_\chi^r$ can be proved to be symmetric if the associated irreducible representation matrices are real (see Appendix \ref{app:symproof}).  There are three types of triatomic interactions:

\begin{description}
  \item[ABC] All atoms are distinct.  The point group corresponding to this type is $C_2$
  \item[AAB] Two atoms are identical. The point group corresponding to this type is $C_{2v}$
  \item[AAA] All atoms are identical. The point group corresponding to this type is $C_{6v}$
\end{description}
For all the above groups in each type, the irreducible representation matrices are all real, so $H_\chi^r$ is symmetric. We provide the representation matrices in Appendix \ref{app:group} for convenience.  \\

One should note that if we confine $\chi$ to $[0,\frac{2\pi}{h}]$, $H_\chi^r$ is a $N_\chi\times N_\chi$ matrix in an one-dimensional irreducible representation and a $2N_\chi\times 2N_\chi$ matrix in a two-dimensional irreducible representation. When the eigensystem of \Eq{eq:red-1D-Hchi} or  \Eq{eq:red-2D-Hchi} is solved, eigenfunctions in the range $[0,\frac{2\pi}{h}]$ are obtained. One can employ \Eq{eq:PR-1} to compute the eigenfunctions on the full range of $\chi$.\\

One can see that using $H^r_\chi$ instead of $H_\chi$ doesn't affect the symmetrization of $H$ in \Sec{sub:symAPH}. Hence we shall use the following definition of $H^s$ in the rest of this paper
\begin{equation}\labeleq{eq:sym-H-group}
H^s=\frac{\hbar^2}{2\mu\rho_\xi^2}\left[H^r_\gamma+\frac{1}{\sin^2\theta}H^r_\chi\right]+H_V
\end{equation}

The derivation presented here may be contrasted with the symmetry-adaptation
procedure provided in Ref. \cite{Iyengar1999}. In Ref. \cite{Iyengar1999}
the symmetry-adaptation
was carried out by projection of the Hamiltonian (or the derivative operator)
onto two different sets of projectors, one on each side of the
Hamiltonian matrix. This is useful for cases where the Hamiltonian (or the
derivative operator) changes the symmetry of the function it acts on
(for example, as in the case of the $d/dx$ operator). In the present case,
the symmetry-adaptation introduced in this section is only used to adapt
the $\chi$-part of the Hamiltonian, which is totally symmetric and hence
does not change the symmetry of the functions it acts on (unlike the operator
corresponding to the $\theta$-part).

\subsection{Reduction of the Hamiltonian using projection}\label{sub:red-prj}

Suppose we discretize $H^s$ using $N_\chi$ grid points in $\chi$ coordinate and $N_\gamma$ in $\gamma$. The size of the Hamiltonian matrix ($N_\gamma N_\chi \times N_\gamma N_\chi)$ to be diagonalized is very large for systems of physical interest when highly accurate eigenvalues and eigenfunctions are needed.
Since, we wish to obtain only a few of the lowest eigenstates to high accuracy, we consider here a projection technique to reduce the size of the Hamiltonian matrix, which leads to reduction in the CPU time and memory requirements. First we find a nearly complete basis for the desired lowest eigenvectors, and then we project the Hamiltonian matrix to this basis: $H^{\text{cut}}=\widetilde{P}H^sP$, where $P$ is the projection matrix and $\widetilde P$ is the transpose of $P$. The matrix $H^{\text{cut}}$ is small. We solve the eigensystem of $H^{\text{cut}}$ and map it back to the original basis to obtain the approximate eigensystem of the $H$. Details are given below.\\

Noting that $H^r_\gamma$ depends only on $\gamma$ and $\frac{\hbar^2}{2\mu\rho_\xi^2}$ is a constant, and introducing the identity matrices $I_\gamma$ and $I_\chi$ in the $\gamma$ and $\chi$ spaces respectively, we can rewrite the Hamiltonian $H^s$ in \Eq{eq:sym-H-group} in the tensor product form,
\begin{eqnarray}
    H^s&=&\frac{\hbar^2}{2\mu\rho_\xi^2} H^r_\gamma\otimes I_\chi + (\frac{\hbar^2}{2\mu\rho_\xi^2} \frac{1}{\sin^2{\theta}}I_\gamma\otimes H^r_\chi + H_V I_\gamma\otimes I_\chi)
\end{eqnarray}
Now $H^s$ is divided into two terms
\begin{eqnarray}
H_1^s&=&\frac{\hbar^2}{2\mu\rho_\xi^2} H_\gamma\otimes I_\chi\\
H_2^s&=&\frac{\hbar^2}{2\mu\rho_\xi^2} \frac{1}{\sin^2{\theta}}I_\gamma\otimes H_\chi + H_V I_\gamma\otimes I_\chi
\end{eqnarray}
The first term, $H_1^s$, is independent of $\chi$, while the second term depends on both $\gamma$ and $\chi$. Since we are only interested in the lowest eigenvalues and eigenfunctions of $H^s$, and $H_1^s$ is independent of $\chi$, we can expect to get accurate results using only a subeigenspace of $H_2^s$ instead. We solve the eigensystem of $H_2^s$ for each fixed $\gamma_j$,
\begin{equation}
    H_2^s(\gamma_j)\alpha^{\gamma_j}_k=\lambda^{\gamma_j}_k\alpha^{\gamma_j}_k    \;\;\;\;j=1,N_\gamma
\end{equation}
where $N_\gamma$ is the number of the grid points used in $\gamma$ and $(\lambda^{\gamma_j}_k, \alpha^{\gamma_j}_k)$ is the $k$-th eigenpair for the given $\gamma_j$. Then we discard the eigenvectors with large eigenvalues. The rest of the eigenvectors consist of a nearly complete basis for the lowest eigenvectors of $H^s$. In our implementation, we sort all the eigenvectors obtained above in ascending order of eigenvalues, and choose only a number of the eigenvectors corresponding to the smallest eigenvalues. After normalization we construct the projection matrix,
\begin{equation}
P=
\begin{pmatrix}
    P_1  &   &   &  \\
    &   P_2  &   &  \\
    &   &   \ddots   & \\
    &   &   &   P_{N_\gamma}\\
\end{pmatrix}
\end{equation}
where
\begin{equation}
    P_i=(\alpha^{\gamma_i}_1, \alpha^{\gamma_i}_2,\ldots,\alpha^{\gamma_i}_{m_{\gamma_i}})
\end{equation}

Because $P$ is an orthonormal, nearly complete basis for the lowest eigenvectors, $\APHphit$,  we have
\begin{equation}
    P \widetilde{P}\APHphit\approx\APHphit
\end{equation}
Substituting $\APHphit$ in \Eq{eq:sym-Schrodinger} and multiplying by $\widetilde{P}$ from the left we get,
\begin{equation}\label{}
  H^{\text{cut}} \APHphicut
    \approx\mathcal{E}^p_{\tau\Lambda}(\rho_\xi)\APHphicut
\end{equation}
where
\begin{equation}
 H^{\text{cut}}=\widetilde{P} H P
\end{equation}
and
\begin{equation}
 \APHphicut=\widetilde{P}\APHphit
\end{equation}
We can see that $P$ is a $N_\gamma N_\chi\times N^{\text{cut}}_\chi$ matrix, where
\begin{equation}
    N^{\text{cut}}=\sum^{N_\gamma}_{i=1}m_{\gamma_i}
\end{equation}
so $H^{\text{cut}}$ is a $N^{\text{cut}}\times N^{\text{cut}}$ matrix. Its size is much smaller than the size of $H^s$ if $N^{\text{cut}}$ is much smaller than $N_\gamma N_\chi$, which is the case if we want only a few eigenvalues. $N^{\text{cut}}$ should be taken to be as small as possible, consistent with the desired level of convergence.\\

Noting that $H^s$ is symmetric, we have
\begin{equation}\label{}
    \widetilde{H}^{\text{cut}}=\widetilde{P}\widetilde{H}^sP=\widetilde{P}H^sP=H^{\text{cut}}
\end{equation}
Thus $H^{\text{cut}}$ is symmetric.




\section{Numerical Tests}\labelsec{sec:Test}
Numerical tests are carried out to study three aspects of the approach. First, we test the behavior with respect to parameter $M$ in the PDAF $\delta_{p,M}^{(k)}$. Next we examine the $H^r_\chi$ using different irreducible representations of various groups. Finally, the accuracy of $H^r_\gamma$ is tested.
\subsection{Test of PDAFs}
We tested PDAFs using different periodic functions, all of them showing similar results. Here we present only the test for the eigenvalues of $H_\chi$. Using the PDAF matrix defined in \Eq{eq:D}, we obtain the discrete presentation of $H_\chi$,
\begin{eqnarray}
    {H_\chi}_{ij}&=&-\sum_{k=1}^{N_\chi}D^{(2)}_{ij}\nonumber\\
                &=&-\sum_{k=1}^{N_\chi}{\delta^{(2)}_{p,M}}(\chi_i-\chi_j)\Delta\chi
\end{eqnarray}
where $N_\chi$ is the number of grid points used in $(0, 2\pi)$, $\Delta\chi=2\pi/N_\chi$ and $\chi_j=(j-\frac{1}{2})\Delta\chi$. We know that the correct eigenvalues of $H_\chi$ should be $0, 1, 1, 4, 4, \cdots, k^2, k^2, \cdots$, so it's easy to check the difference between the computed eigenvalues and the exact ones. We set the criteria of $10^{-6}$ as the maximum tolerable error and then count how many good eigenvalues can be obtained. A typical result is shown in \Fig{fig:TestM}, in which we set $N_\chi=32$ and various $M$ from 1 to 35. One can see that when $M=N_\chi/2$, almost all eigenvalues are good. This suggests the optimal choice of $M$.
\begin{figure}
 \begin{center}
 \includegraphics[width=4in,height=3in]{graph/TestM.eps}
 \end{center}
 \caption{The number of good eigenvalues reaches maximum when $M=N_\chi/2$ ($N_\chi=32$), where $M$ is the parameter in PDAF $\delta^{(k)}_{p,M}(x)$, and $N_\chi$ is the number of gird points in range $(0, 2\pi)$.}\label{fig:TestM}
\end{figure}

\subsection{Test of $H^r_\chi$}
Tests of $H^r_\chi$ were done for all the aforementioned irreducible representations of the relevant point groups. All of them give good results. Here we present only one typical result computed using the $C_{2v}$ group, since only $C_{2v}$ is involved in the $FH_2$ reaction calculation.\\

Similar to $H_\chi$, we can obtain the PDAF representation of $H^r_\chi$. Note that $C_{2v}$ has only one-dimensional irreducible representations, so we use \Eq{eq:redHchi-1D} for $H^r_\chi$. We obtain the discrete representation of $H^r_\chi$ as,
\begin{equation}\label{}
{H_\chi}_{ij}=-\sum_{k=1}^{N_\chi}\sum_{R}{\delta^{(2)}_{p,M}}(\chi_i-R\chi_j)\Gamma^{[z]}\Delta\chi
\end{equation}
Note that the grid scheme is defined in \Eq{eq:chischeme}. We set $N_\chi=8$, and the eigenvalues $\lambda_k$ for each of the irreducible representations are computed, and the $\sqrt{\lambda_k}$ are shown in table \ref{tbl:eigenHchiC2v} for easy comparison. One can see that if we combine the results from all the irreducible representations, we will obtain the (approximate) eigenvalues $0, 1, 1, 4, 4,\cdots,k^2,k^2,\cdots$, which are obtained from $H_\chi$. We also note that the degenerate eigenvalues in $H_\chi$ are no longer degenerate in $H^r_\chi$.

\begin{table}[htp!]
\begin{tabular}{l|cccc}
\hline
\hline
$k$ &        $A_1$       &       $B_1$      &        $A_2$     &        $B_2$    \\
\hline
1  &   0.00000052456066 & 2.00000000000009 & 0.99999999999989 & 1.00000000000007\\
2  &   2.00000000000007 & 3.99999999999997 & 3.00000000000004 & 3.00000000000002\\
3  &   3.99999999999999 & 5.99999999999999 & 5.00000000000000 & 5.00000000000000\\
4  &   5.99999999999999 & 7.99999999999999 & 7.00000000000000 & 7.00000000000000\\
5  &   7.99999999999999 & 10.0000000000000 & 9.00000000000001 & 9.00000000000001\\
6  &   10.0000000000000 & 12.0000000000000 & 11.0000000000000 & 11.0000000000000\\
7  &   12.0000000000000 & 14.0000000000000 & 13.0000000000000 & 13.0000000000000\\
8  &   14.0000000000000 & 22.6274169979695 & 15.0000000000000 & 15.0000000000000\\
\hline
\end{tabular}\\
%\parbox{73mm}{
\caption{$\sqrt{\lambda_k}$, square roots of eigenvalues of $H^r_\chi$ in irreducible representations of group $C_{2v}$\label{tbl:eigenHchiC2v}}
%}
\end{table}

\subsection{Test of $H^r_\gamma$}
$H^r_\gamma$ in \Eq{eq:Hrgamma} has the same eigenvalues as $H_\theta$ in \Eq{eq:Htheta}. If we rewrite $H_\theta$ as
\begin{equation}
    H_\theta=-\frac{8}{\sin{2\theta}} \frac{\partial}{\partial 2\theta} \sin 2\theta\frac{\partial}{\partial 2\theta}
\end{equation}
we can see that the eigen-equation of $H_\theta$, \Eq{eq:eigen-theta}, is a Legendre differential equation if we write the eigenvalue as
\begin{equation}
    \lambda_l=8l(l+1)\;\;l=0,1,...
\end{equation}
Thus $H^r_\gamma$ has eigenvalues of $8l(l+1)$. To evaluate the accuracy of the computed eigenvalue, we define $S_l(\lambda_l)$ as follows,
\begin{equation}\label{}
  S_l(\lambda_l)=
  \begin{cases}
    -\lg|\lambda_l| & \text{l=0}, \\
    -\lg\frac{\lambda_l-8l(l+1)}{8l(l+1)} & l>0
  \end{cases}
\end{equation}
where $\lambda_l$ is the $l$-th computed eigenvalue, and $S_l(\lambda_l)$ gives the approximate number of significant digits. \Fig{fig:AccurHgamma} shows the results when we set $N_\gamma=16,32,48 $ and $64$. One can see that for each $N_\gamma$ there are about $N_\gamma/3$ eigenvalues of high accuracy. This is accurate enough for the $FH_2$ calculation, because typically the order of the Hamiltonian matrix is about $1000$ but we require less than 300 eigenvalues.\\

It is worth mentioning that we also tested the non-symmetric form, $H_\theta$, using the PDAF presentation. Although we can get as highly accurate results as we get from $H_\chi$ or $H^r_\chi$, we still choose the symmetric form, because the non-symmetric matrix costs more time and memory to diagonalize. Moreover, when we add the potential, the accuracy of the non-symmetric form will decrease to the accuracy of the symmetric $H_\gamma^r$.
\begin{figure}
 \begin{center}
 \includegraphics[width=4in,height=3in]{graph/AccurHgamma.eps}
 \end{center}
 \caption{Accuracy (number of significant digits) in eigenvalues of $H^r_\gamma$. For all $N_\gamma$, the smallest $N_\gamma/3$ eigenvalues are of high accuracy.}\label{fig:AccurHgamma}
\end{figure}

\section{Calculations for the $FH_2$ System} \labelsec{sec:CalcFH2}
In this section we report the results of PDAF calculations of surface functions and the matrix elements, and compare them with the results of the DVR and ABM methods. The system chosen as a nontrivial example is the $\text{F+H}_2\rightarrow \text{HF+H}$ reaction; its treatment requires generation of a large basis of surface functions. The potential energy surface (PES) used is that of Brown \textit{et al}.\cite{FH2_T5A} commonly called the T5A surface, and we choose the zero of energy to be at the bottom of the asymptotic HF potential wells. This PES has been used in many calculations\cite{Launay1991,Zhang1991,Bacic1990, Kress1989, Yu1989,Manolopoulos1990,Neuhauser1990} on this reaction, and plots of the PES and surface functions showing their appearance in APH coordinates have also been published. Arrangement 1 or i(initial) is taken to be $\text{F+H}_2$  reactants.\\

The calculations are for $\Lambda=0$ and even parity ($p=0$) and include all functions connecting to the even $j$ rotational states of the $\text{F+H}_2$ arrangement. Because of the symmetry due to the identical H atoms, this only requires including in the PDAF calculations the $A_1$ irreducible representation, in which the surface functions are even under reflection about $\chi=0$.\\
\subsection{Grid Size and Mapping}
The grid sizes (the number of grid points),$N_\chi$, $N_\gamma$, and $N^{\text{cut}}$ are determined by the convergence test. We present only the results here. To attain five significant digits in the lowest 100 eigenvalues for any $\rho_\xi$,

\begin{eqnarray}
  \labeleq{eq:Ngamma}  N_\gamma&=&\lfloor 25.852+5.85119\rho_\xi-0.130102\rho_\xi^2-0.0042517*\rho_\xi^3 \rfloor\\
  \labeleq{eq:Ncut}    N^{\text{cut}}&=&\lfloor 1428.06+29.2007\rho_\xi-6.37755\rho_\xi^2-0.119048*\rho_\xi^3 \rfloor\\
  \labeleq{eq:Nchi}    N_\chi&=&\lfloor 4\pi\max(N_\gamma)\rfloor/h
\end{eqnarray}
where $h$ is the order of the associated symmetry group. To get $N_\chi$, we use the maximum of $N_\gamma$ because it makes $N_\chi$ identical for all $\rho_\xi$, thus making it easy to compute the overlap matrices. Also, using the maximum $N_\gamma$ doesn't significantly increase the computation time, because we reduce the matrix size according to $N^{\text{cut}}$, as discussed in \Sec{sub:red-prj}.\\

Although $N_\chi$ is chosen to be the same for all $\rho_\xi$, $N_\gamma$ is different. Hence, before we compute the overlap matrix, we have to map the wave functions to a uniform grid, to reduce the computation time for the overlap matrix.

Let $N^u_\gamma$ be the number of uniform grid points in $\gamma$; $N_\gamma$ is the original number of grid points. Applying \Eq{eq:PDAFintd} to the wave function
\begin{eqnarray}\labeleq{eq:gridmapping}
  \Phi^p_{\tau\Lambda}(\gamma_j,\chi;\rho_\xi)&=&\sum ^{2N_\gamma}_{k=1} \delta_{p,M}(\gamma_j-\gamma_k)\Phi^p_{\tau\Lambda}(\gamma_k,\chi;\rho_\xi)\Delta\gamma\nonumber\\
&=&\sum ^{N_\gamma}_{k=1}\left\{ \delta_{p,M}(\gamma_j-\gamma_k)\Phi^p_{\tau\Lambda}(\gamma_k,\chi;\rho_\xi)+
\delta_{p,M}(\gamma_j-2\pi+\gamma_k)\Phi^p_{\tau\Lambda}(2\pi-\gamma_k,\chi;\rho_\xi)\right\}\Delta\gamma \nonumber\\
&=&\sum ^{N_\gamma}_{k=1}\left[\delta_{p,M}(\gamma_j-\gamma_k)+\delta_{p,M}(\gamma_j-2\pi+\gamma_k)\right]\Phi^p_{\tau\Lambda}(\gamma_k,\chi;\rho_\xi)\Delta\gamma
\end{eqnarray}
where $\Delta\gamma=\pi/N_\gamma$, and $\gamma_j$,$\gamma_k$ represent the coordinates in the above two grid schemes respectively,
\begin{eqnarray}
  \gamma_j&=&(j-\frac{1}{2})\Delta\gamma^u\\
  \gamma_k&=&(k-\frac{1}{2})\Delta\gamma
\end{eqnarray}
where $\Delta\gamma^u=\pi/N^u_\gamma$\\
\subsection{Eigenvalues}
The atomic masses used are:
\begin{eqnarray}
\text{mass of F}&=&18.9984032\;\text{a.u.}\\
\text{mass of H}&=&1.00782503\;\text{a.u.}
\end{eqnarray}
The calculations were performed at five representative $\rho$ values ranging from the smallest to the largest values that were needed in our reactive scattering calculations\cite{Bacic1990,Kress1989}. The precise values of $\rho$ chosen have no particular significance, and some were chosen simply because the convergence of the DVR method had already been studied at these values.\\

The results are given in \Tab{tbl:E2.2} through \Tab{tbl:E9.0} for the five values of $\rho$ chosen. In each table, the energy eigenvalues of the highest ten important or open surface function states are given in eV. The omitted lower eigenvalues always agree to more significant figures than those shown. Also shown is $\overline{\mathcal{E}}(n)$, the average of the first $n$ eigenvalues. This gives a convenient measure of the overall agreement of the methods.
%Table
\begin{table}[htp!]
\begin{tabular}{ccccc}
\hline
\hline
$\tau$&PDAF&ABM&DVR&FEM\\
\hline
1&6.9176 &6.9177 &6.9177 &6.9177\\
2&7.0341 &7.0341 &7.0341 &7.0342\\
3&7.1299 &7.1298 &7.1296 &7.1299\\
4&7.3044 &7.3042 &7.3038 &7.3043\\
5&7.4934 &7.4932 &7.4927 &7.4933\\
6&7.7007 &7.7007 &7.7008 &7.7009\\
7&7.9227 &7.9219 &7.9202 &7.9222\\
8&8.0786 &8.0789 &8.0786 &8.0791\\
9&8.1548 &8.1524 &8.1475 &8.1529\\
10&8.2495 &8.2502 &8.2494 &8.2503\\
$\overline{\mathcal{E}}(10)$&7.5986 &7.5983 &7.5974 &7.5985\\
\hline
\end{tabular}\\
\parbox{73mm}{
\caption{PDAF, ABM, DVR, and FEM surface function energies $\mathcal{E}_\tau$ and average energies in eV at $\rho=2.2\;a_0$. \label{tbl:E2.2}}}
\end{table}

%Table
\begin{table}[htp!]
\begin{tabular}{ccccc}
\hline
\hline
$\tau$&PDAF&ABM&DVR&FEM\\
\hline
11&1.9845 &1.9845 &1.9845 &1.9849\\
12&1.9975 &1.9975 &1.9975 &1.9979\\
13&2.0514 &2.0515 &2.0515 &2.0519\\
14&2.1090 &2.1090 &2.1090 &2.1093\\
15&2.1275 &2.1275 &2.1275 &2.1281\\
16&2.1352 &2.1352 &2.1353 &2.1357\\
17&2.1993 &2.1994 &2.1994 &2.1998\\
18&2.2578 &2.2578 &2.2578 &2.2584\\
19&2.2730 &2.2730 &2.2730 &2.2737\\
20&2.3197 &2.3197 &2.3197 &2.3204\\
$\overline{\mathcal{E}}(20)$&1.8551 &1.8551 &1.8551 &1.8554\\
\hline
\end{tabular}\\
\parbox{73mm}{
\caption{PDAF, ABM, DVR, and FEM surface function energies $\mathcal{E}_\tau$ and average energies in eV at $\rho=3.0375828\;a_0$. \label{tbl:E3.0}}}
\end{table}

%Table
\begin{table}[htp!]
\begin{tabular}{ccccc}
\hline
\hline
$\tau$&PDAF&ABM&DVR&FEM\\
\hline
91&2.0699 &2.0701 &2.0701 &2.0720\\
92&2.0721 &2.0734 &2.0727 &2.0760\\
93&2.0864 &2.0896 &2.0869 &2.0916\\
94&2.0941 &2.0941 &2.0943 &2.0961\\
95&2.1087 &2.1107 &2.1093 &2.1132\\
96&2.1215 &2.1232 &2.1220 &2.1267\\
97&2.1487 &2.1509 &2.1495 &2.1549\\
98&2.1643 &2.1650 &2.1645 &2.1672\\
99&2.1793 &2.1806 &2.1799 &2.1833\\
100&2.1892 &2.1892 &2.1893 &2.1918\\
$\overline{\mathcal{E}}(100)$&1.3974 &1.3977 &1.3975 &1.3985\\
\hline
\end{tabular}\\
\parbox{73mm}{
\caption{PDAF, ABM, DVR, and FEM surface function energies $\mathcal{E}_\tau$ and average energies in eV at $\rho=4.9747966\;a_0$. \label{tbl:E4.9}}}
\end{table}

%Table
\begin{table}[htp!]
\begin{tabular}{ccccc}
\hline
\hline
$\tau$&PDAF&ABM&DVR&FEM\\
\hline
91&2.0559 &2.0564 &2.0569 &2.0679\\
92&2.0749 &2.0754 &2.0760 &2.0928\\
93&2.0980 &2.0981 &2.0981 &2.1023\\
94&2.0985 &2.0990 &2.0994 &2.1141\\
95&2.1098 &2.1098 &2.1099 &2.1148\\
96&2.1194 &2.1196 &2.1198 &2.1278\\
97&2.1268 &2.1273 &2.1278 &2.1443\\
98&2.1538 &2.1542 &2.1543 &2.1537\\
99&2.1551 &2.1551 &2.1551 &2.1608\\
100&2.1597 &2.1602 &2.1608 &2.1621\\
$\overline{\mathcal{E}}(100)$&1.3788 &1.3788 &1.3790 &1.3811\\
\hline
\end{tabular}\\
\parbox{73mm}{
\caption{PDAF, ABM, DVR, and FEM surface function energies $\mathcal{E}_\tau$ and average energies in eV at $\rho=7.2989993\;a_0$.\label{tbl:E7.2}}}
\end{table}

%Table
\begin{table}[htp!]
\begin{tabular}{ccccc}
\hline
\hline
$\tau$&PDAF&ABM&DVR&FEM\\
\hline
91&2.0747 &2.0747 &2.0748 &2.0805\\
92&2.0774 &2.0779 &2.0782 &2.0863\\
93&2.0957 &2.0964 &2.0974 &2.1051\\
94&2.1010 &2.1010 &2.1009 &2.1065\\
95&2.1153 &2.1153 &2.1154 &2.1233\\
96&2.1189 &2.1194 &2.1197 &2.1284\\
97&2.1250 &2.1252 &2.1253 &2.1315\\
98&2.1461 &2.1470 &2.1482 &2.1541\\
99&2.1542 &2.1545 &2.1546 &2.1562\\
100&2.1558 &2.1559 &2.1560 &2.1564\\
$\overline{\mathcal{E}}(100)$&1.3791 &1.3791 &1.3792 &1.3812\\
\hline
\end{tabular}\\
\parbox{73mm}{
\caption{PDAF, ABM, DVR, and FEM surface function energies $\mathcal{E}_\tau$ and average energies in eV at $\rho=9.0 a_0$. \label{tbl:E9.0}}}
\end{table}

\subsection{Matrix Elements}
The APH surface functions $\APHphi$ are ``sector adiabatic''; that is, they change from sector to sector, but are independent of $\rho$ within a sector. Thus, when the APH wave function is substituted into the Schr\"odinger equation, the resulting exact coupled channel or close coupling equations are of the form\cite{PackParker1987}
\begin{eqnarray}\label{eq:3D-eigen}
    \left[\frac{\partial^2}{\partial\rho^2}+\frac{2\mu E}{\hbar^2}\right]\psi^{Jpn}_{\tau\Lambda}  (\rho)=\frac{2\mu}{\hbar^2}\sum_{\tau'A'}{\langle \Phi^{Jp}_{\tau\Lambda}\hat{D}^{Jp}_{\Lambda M} |H_i|   \Phi^{Jp}_{\tau'\Lambda'} \hat{D}^{Jp}_{\Lambda' M}\rangle \psi^{Jpn}_{\tau'\Lambda'}{\rho} }
\end{eqnarray}
The matrix elements are obtained in Reference \cite{PackParker1987} as,
\begin{eqnarray}\labeleq{eq:MatElements}
            \langle \Phi^{Jp}_{\tau\Lambda}\hat{D}^{Jp}_{\Lambda M}|H_i|
            \Phi^{Jp}_{\tau'\Lambda'}\hat{D}^{Jp}_{\Lambda' M}\rangle  \nonumber
            &=&\frac{\rho_\xi^2}{\rho^2}\mathcal E_{\tau\Lambda}(\rho_\xi)\delta_{\tau\tau'}\delta_                 {\Lambda\Lambda'} \nonumber \\
            &&+\delta_{\Lambda\Lambda'}\langle  \Phi^{Jp}_{\tau\Lambda}\left| V(\rho,\theta,\chi)-            \frac{\rho_\xi^2}{\rho^2}V(\rho_\xi,\theta,\chi) \right|\Phi^{Jp}_{\tau'\Lambda}                    \rangle  \nonumber \\
            &&+ \langle  \Phi^{Jp}_{\tau\Lambda}\hat{D}^{Jp}_{\Lambda M} \left| \frac{A-B}{2}(J_x^2-J_y^2) + T_c \right| \Phi^{Jp}_{\tau'\Lambda'} \hat D^{Jp}_{\Lambda' M} \rangle
\end{eqnarray}
where $\rho_\xi$ denotes the $\xi^{th}$ hyper-radius sampled from interval $[\rho_{\min},\rho_{\max}]$. The $\rho_\xi$ are given by\\
\begin{equation}
    \rho_\xi=[\rho_{\min}+(\xi-1)\Delta\rho_1](1+\Delta\rho_2)^{\xi-1}
\end{equation}
This algorithm spaces the sector centers logarithmically. Given a sector with sector center $\rho_\xi$, we evaluate matrix elements at the three rho values given by
\begin{eqnarray}
    \rho_1&=&\frac{\rho_{\xi-1}+\rho_{\xi}}{2}\\
    \rho_2&=&\rho_\xi\\
    \rho_3&=&\frac{\rho_{\xi}+\rho_{\xi+1}}{2}
\end{eqnarray}
All the matrix elements in \Eq{eq:MatElements} are independent of $E$, so that they can be calculated once, stored, and used at many scattering energies.\\

The wave functions $\Phi^p_{\tau,\Lambda}(\gamma,\chi;\rho_\xi)$ should be normalized before computation of the matrix elements. The normalization factor $N$ can be calculated in the desired irreducible representation easily according to \Eq{eq:norm} as shown below,
\begin{eqnarray}\labeleq{eq:Norm}
1&=&N\int^{2\pi}_0d\chi\int^{\pi/2}_0d\theta\sin{2\theta}{\Phi^p_{\tau\Lambda}}^2(\theta,\chi;\rho_\xi)\nonumber\\
&=&N\int^{2\pi}_0d\chi\int^{\pi}_0d\gamma\frac{\pi}{4}\sin{\gamma}\sin{2\theta}{\Phi^p_{\tau\Lambda}}^2(\gamma,\chi;\rho_\xi)\nonumber\\
   &=&N\sum^{h N_\chi}_{i=1}\sum^{N_\gamma}_{j=1} \Delta\chi \Delta\gamma \frac{\pi}{4}\sin{\gamma_j}\sin{2\theta(\gamma_j)}{\Phi^p_{\tau\Lambda}}^2(\gamma_j,\chi_i;\rho_\xi)\nonumber\\
   &=&N\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1} \sum_{R} \Delta\chi \Delta\gamma \frac{\pi}{4}\sin{\gamma_j}\sin{2\theta(\gamma_j)}\Phi^p_{\tau\Lambda}(\gamma_j,R^{-1}\chi_i;\rho_\xi)\Phi^p_{\tau\Lambda}(\gamma_j,R^{-1}\chi_i;\rho_\xi)\nonumber\\
   &=&N\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1}\sum_{R}\Delta\chi \Delta\gamma \frac{\pi}{4}\sin{\gamma_j}\sin{2\theta(\gamma_j)}{\Gamma^{[z]}}^2(R){\Phi^p_{\tau\Lambda}}^2(\gamma_j,\chi_i;\rho_\xi)\nonumber\\
&=&h N\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1}\Delta\chi \Delta\gamma \frac{\pi}{4}\sin{\gamma_j}\sin{2\theta(\gamma_j)}{\Phi^p_{\tau\Lambda}}^2(\gamma_j,\chi_i;\rho_\xi)\nonumber
\end{eqnarray}
Thus
\begin{equation}
N=\left\{\frac{1}{4}h\pi\Delta\chi\Delta\gamma\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1} \sin{\gamma_j}\sin{2\theta(\gamma_j)}{\Phi^p_{\tau\Lambda}}^2(\gamma_j,\chi_i;\rho_\xi) \right\}^{-\frac{1}{2}}
\end{equation}
After $\APHphi$ is normalized, we compute the matrix elements. The first term on the right-hand side of \Eq{eq:MatElements} is just a local internal energy which together with the $E$-term on the left-hand side of \Eq{eq:3D-eigen} determines a local wave number.

The second term of \Eq{eq:3D-eigen} is often called a potential matrix element. It is small on the sector and can be evaluated with the same quadratures used in getting the surface functions. Similar to the evaluation of the normalization factor, we obtain for the potential matrix elements,
\begin{eqnarray}\labeleq{eq:Potmat}
\langle \Phi^{Jp}_{\tau\Lambda}(\rho_\xi)\mid V(\rho)-\frac{\rho_\xi^2}{\rho^2}V(\rho_\xi) \mid \Phi^{Jp}_{\tau'\Lambda} (\rho_\xi)\rangle \nonumber
&=& \frac{1}{4}h\pi\Delta\chi\Delta\gamma\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1} \sin{\gamma_j}\sin{2\theta(\gamma_j)}\times \nonumber\\
&&\;\;\Phi^p_{\tau\Lambda}(\gamma_j,\chi_i;\rho_\xi)\left[V(\rho)-\frac{\rho_\xi^2}{\rho^2}V(\rho_\xi)\right]
\Phi^p_{\tau'\Lambda}(\gamma_j,\chi_i;\rho_\xi)
\end{eqnarray}

The Coriolis term can be simplified as
\begin{eqnarray}
    \langle \Phi^{Jp}_{\tau\Lambda}\hat D^{Jp}_{\Lambda M}\vert T_c\vert \Phi^{Jp}_{\tau'\Lambda'}\hat D^{Jp}_{\Lambda'M}\rangle \nonumber
    &=&\frac{-\hbar^2}{2\mu\rho^2}\left< \Phi^{p}_{\tau\Lambda}\left|\frac{\cos\theta}{\sin^2\theta}\frac{\partial}{\partial\chi} \right|\Phi^p_{\tau'\Lambda'} \right>  \nonumber\\
&&\times[(1+\delta_{\Lambda0})(1+\delta_{\Lambda'0})]^{-1/2}\nonumber \\
&&\times[\lambda_+(J,\Lambda)\delta_{\Lambda',\Lambda+1}-\lambda_-(J,\Lambda)\delta_{\Lambda',\Lambda-1}\nonumber \\
&&+\lambda_-(J,\Lambda)(-1)^{J+\Lambda+p}\delta_{\Lambda,1-\Lambda}]\label{eq:coriolis}
\end{eqnarray}
where
\begin{equation}
    \lambda_\pm(J,\Lambda)=[(J\pm\Lambda+1)(J\mp\Lambda)]^{\frac{1}{2}}
\end{equation}
and
\begin{eqnarray}
   \left< \Phi^{p}_{\tau\Lambda}\left|\frac{\cos\theta}{\sin^2\theta}\frac{\partial}{\partial\chi} \right|\Phi^p_{\tau'\Lambda'} \right>  \nonumber
&=& \frac{1}{2}h\pi\Delta\chi^2\Delta\gamma\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1}\sum^{N_\chi}_{k=1}
 \sin{\gamma_j}\frac{\cos^2\theta(\gamma_j)}{\sin\theta(\gamma_j)}\times \nonumber\\
&&\;\;\Phi^p_{\tau\Lambda}(\gamma_j,\chi_i;\rho_\xi)
    \delta^{(1)}_{p,M}(\chi_i-\chi_k)\Phi^p_{\tau'\Lambda'}(\gamma_j,\chi_k;\rho_\xi)
\end{eqnarray}

It should be noted that the last term in the bracket in \Eq{eq:coriolis} is only nonzero when $\Lambda=0$ or 1 and also that because the $\rho$ dependence of the operator factors out, the matrix elements over the $\APHphi$ only need to be evaluated once on each sector. 
They are readily evaluated using the PDAF $\APHphi$ and their quadrature points since the PDAF code generates the derivatives of $\APHphi$ directly.\\

The asymmetric top terms of \Eq{eq:MatElements} can be explicitly written as
\begin{eqnarray}
 \langle \tau\Lambda\vert \frac{1}{2}(A-B)(J_x^2-J_y^2)\vert \tau'\Lambda'\rangle \nonumber
 &=&\frac{1}{4}\hbar^2\langle \Phi^p_{\tau\Lambda}\vert A-B\vert \Phi^p_{\tau'\Lambda'}\rangle [(1+\delta_{\Lambda0})(1+\delta{\Lambda'0})]^{-1/2}\nonumber \\
    &&\times[\lambda_+(J,\Lambda)\lambda_+(J,\Lambda+1)\delta_{\Lambda',\Lambda+2}\nonumber \\
    &&+\lambda_-(J,\Lambda)\lambda_-(J,\Lambda-1)\delta_{\Lambda',\Lambda-2}\nonumber \\
    &&+(-1)^{J+\Lambda+p}\lambda_-(J,\Lambda)\lambda_-(J,\Lambda-1)\delta_{\Lambda',2-\Lambda}]\label{eq:asymtop}
\end{eqnarray}
where
\begin{eqnarray}
    A=\frac{1}{\mu\rho_\xi^2(1+\sin^2\theta)} \\
    B=\frac{1}{2\mu\rho_\xi^2\sin^2\theta}
\end{eqnarray}

\begin{eqnarray}
\langle \Phi^p_{\tau\Lambda}\vert A-B\vert \Phi^p_{\tau'\Lambda'}\rangle\nonumber
&&= \frac{1}{4}h\pi\Delta\chi\Delta\gamma\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1} \sin{\gamma_j}\sin{2\theta(\gamma_j)}\times \nonumber\\
&&\;\;\Phi^p_{\tau\Lambda}(\gamma_j,\chi_i;\rho_\xi)(A-B)
\Phi^p_{\tau'\Lambda'}(\gamma_j,\chi_i;\rho_\xi)
\end{eqnarray}

The third term in the bracket in \Eq{eq:asymtop} is always zero if $|\Lambda-\Lambda'|>2$, and the $\rho$ dependence of $A-B$ again factors out, giving the same simplifications and allowing the evaluation of the integrals by the same methods as for the Coriolis terms.

At the boundaries between sectors, the $R$ matrix is transformed by an orthogonal transformation, which requires calculation of the overlap matrix elements. The formula for the overlap matrix elements is given by ,
\begin{eqnarray}\labeleq{eq:Overlap}
&&\langle \Phi^{Jp}_{\tau\Lambda}(\rho_\xi)\mid \Phi^{Jp}_{\tau'\Lambda} (\rho_{\xi'})\rangle \nonumber\\
&&= \frac{1}{4}h\pi\Delta\chi\Delta\gamma^u\sum^{N_\chi}_{i=1}\sum^{N_\gamma}_{j=1} \frac{\pi}{4}\sin{\gamma_j}\sin{2\theta(\gamma_j)}\times \nonumber\\
&&\;\;\Phi^{Jp}_{\tau\Lambda}(\gamma_j,\chi_i;\rho_\xi)\Phi^{Jp}_{\tau'\Lambda}(\gamma_j,\chi_i\rho_{\xi'})
\end{eqnarray}

We compare the computational efficiency of the PDAF, DVR, ABM, and FEM methods on the same computer (PIII 866MHz), by computing eigenvalues, potential matrix elements, and overlap matrices at 100 $\rho_\xi$'s. $\rho_\xi$ starts from $2.0\;a_0$ and ends at $9.0\;a_0$ and is evenly spaced in between. The computation times (CPU time) are shown in \Fig{fig:Time-Rho}. 
The FEM method takes much more time than the others, thus we do not include it in \Fig{fig:Time-Rho}. We see that the PDAF method is much faster than the DVR over most of the range of $\rho_\xi$. If $\rho_\xi$ is very small (less than $2.2\;a_0$), the DVR is the most efficient method. Although the ABM works slightly faster than the PDAF method, it turns out that the ABM diverges when $\rho_\xi$ is small.

\begin{figure}
 \begin{center}
  \includegraphics[width=4in,height=3in]{graph/Time-Rho.eps}
 \end{center}
 \caption{CPU Time of PDAF, DVR and ABM}\labelfig{fig:Time-Rho}
\end{figure}

We compare only the potential matrix elements and the overlap matrix elements here, computed using the ABM, DVR and PDAF. The comparison of the potential matrix elements is shown in \Fig{fig:SPotmat} and one can see that all the three methods agree with each other very well when $\rho>3.3 a_0$. However, when $\rho<3.3$, the ABM result diverges from the PDAF and DVR significantly.
The comparison of the overlap matrix elements is shown in \Fig{fig:SOverlap}. 
One can see clearly from the figure that the PDAF gives results very close to the DVR at small $\rho$ ($rho<3.4 a_0$), and it also agrees with the ABM very well at large $\rho$ ($rho>4.4 a_0$). 
The ABM gives very different results from the PDAF and DVR at small $\rho$, and DVR gives very different results from PDAF and ABM at large $\rho$. 
We know that the ABM is very accurate at large $\rho$ and the DVR at small $\rho$, and thus the comparisons show that the PDAF is accurate both for small and large $\rho$.
\begin{figure}
 \begin{center}
 \includegraphics[width=4in,height=3in]{graph/SPotmat.eps}
 \end{center}
 \caption{Comparison of the Sum of potential matrix elements using PDAF, ABM and DVR}\labelfig{fig:SPotmat}
\end{figure}

\begin{figure}
 \begin{center}
 \includegraphics[width=4in,height=3in]{graph/SOverlap.eps}
 \end{center}
 \caption{Comparison of the Sum of the overlap matrix elements using PDAF, ABM and DVR.}\labelfig{fig:SOverlap}
\end{figure}

\section{Conclusion}\labelsec{sec:Conclusion}
In this paper we have presented a periodic distributed approximating functional (PDAF) method for calculating the surface function basis needed in hyperspherical formulations of reactive scattering theory. PDAF functions are introduced and shown to be capable of providing an accurate, efficient representation of the derivative operators.

Test calculations on the $\text{F}+\text{H}_2$ system with the T5A PES, comparing the PDAF, ABM, DVR and FEM methods showed that the FEM is always the least efficient of the four, and the ABM is the most efficient method for large $\rho$ but is not accurate at small $\rho$. On the other hand, the DVR is the most efficient method for small $\rho$ but is not accurate for large $\rho$. 
The PDAF is efficient at both large $\rho$ like the ABM and small $\rho$ like the DVR and it is accurate for all $\rho$. As a result, we find that PDAF is the best method for surface function calculations in hyperspherical reactive scattering calculations.


\section{Acknowledgments}
We would like to thank the National Science Foundation NSF CHE-9710383
and NSF PHY-0100704 for support. One of us (GAP) would also like to
thank the members of Group T-12 at the Los Alamos National
Laboratory (LANL) for their helpfulness and
hospitality extended during summer visits.

\clearpage

\appendix

\section{$H^r_\chi$ is symmetric}\label{app:symproof}
\begin{equation}
  H^r_\chi(\chi_i,\chi_j)=\sum_{R}{H_\chi}(\chi_i,R\chi_j)\Gamma^{[\tau]*}(R)
\end{equation}
To prove that $H^r_\chi$ is symmetric, it is sufficient to show
\begin{equation} \label{eq:sym}
{H^r_\chi}_{\kappa,\kappa'}(\chi_i,\chi_j)={H^r_\chi}_{\kappa',\kappa}(\chi_j,\chi_i)
\end{equation}
where
\begin{equation} \label{eq:Hele}
{H^r_\chi}_{\kappa,\kappa'}(\chi_i,\chi_j)=\sum_{R}{H_\chi}(\chi_i,R\chi_j)\Gamma_{\kappa,\kappa'}^{[z]*}(R)
\end{equation}

We know that $H_\chi$ is symmetric, so
\begin{eqnarray}\label{eq:Hrr}
{H^r_\chi}_{\kappa',\kappa}(\chi_j,\chi_i)
&=&\sum_{R}{H_\chi}(\chi_j,R\chi_i)\Gamma_{\kappa',\kappa}^{[z]*}(R)\nonumber\\
&=&\sum_{R}{H_\chi}(R\chi_i,\chi_j)\Gamma_{\kappa',\kappa}^{[z]*}(R)\nonumber\\
&=&\sum_{R}{H_\chi}(R^{-1}\chi_i,\chi_j)\Gamma_{\kappa',\kappa}^{[z]*}(R^{-1})\nonumber\\
&=&\sum_{R}{H_\chi}(R^{-1}\chi_i,\chi_j)\Gamma_{\kappa',\kappa}^{[z]*{-1}}(R)\nonumber\\
&=&\sum_{R}{H_\chi}(R^{-1}\chi_i,\chi_j)\Gamma_{\kappa',\kappa}^{[z]*\dag}(R)\nonumber\\
&=&\sum_{R}{H_\chi}(R^{-1}\chi_i,\chi_j)\Gamma_{\kappa,\kappa'}^{[z]}(R)
\end{eqnarray}
Because $H_\chi(\chi_i,\chi_j)$ ($DAF^{(2)}(\chi_i-\chi_j)$) depends on $|\chi_i-\chi_j|$, and $R$ is a length-preserving operation, we have
\begin{eqnarray}\label{eq:Hchir}
H_\chi(R^{-1}\chi_j,\chi_i)
&=&H_\chi(R^{-1}\chi_j,\chi_i)\nonumber\\
&=&H_\chi(R R^{-1}\chi_j,R\chi_i)\nonumber\\
&=&H_\chi(\chi_j,R\chi_i)
\end{eqnarray}
Combining (\ref{eq:Hrr}) and (\ref{eq:Hchir}) gives
\begin{equation} {H^r_\chi}_{\kappa,\kappa'}(\chi_i,\chi_j)=\sum_{R}{H_\chi}(\chi_i,R\chi_j)\Gamma_{\kappa,\kappa'}^{[z]}(R)
\end{equation}
Therefore, if $\Gamma^{[z]}$ is real, (\ref{eq:sym}) is true, and $H_\chi^r$ is symmetric.


\section{Irreducible representations of group $C_2$, $C_{2v}$ and $C_{6v}$}\label{app:group}
The irreducible representation matrices for some point groups frequently used in APH surface function 
computations are given in Tables \ref{tbl:C2}, \ref{tbl:C2v} and \ref{tbl:C6v}. 
The first column in each table gives the names of the irreducible representations. 
The second row gives the transformation when a symmetry operation $R$ acts on $\chi$. The parity $p$ is also given in each table.

%C2
\begin{table}[htp!]
\begin{tabular}{l|cc|r}
\hline
\hline
$C_2$&$E$&$C_2$&$p$\\
$R\chi$&$\chi$&$\pi+\chi$\\
\hline
$A$&1&1&0\\
$B$&1&-1&1\\
\hline
\end{tabular}\\
\parbox{73mm}{\caption{Irreducible representations of group $C_2$\labeltbl{tbl:C2}}}
\end{table}

%C2v
\begin{table}[htp!]
\begin{tabular}{l|cccc|r}
\hline
\hline
$C_{2v}$&$E$&$C_2$&$\sigma_v$&$\sigma_v'$&$p$\\
$R\chi$&$\chi$&$\pi+\chi$&$2\pi-\chi$&$\pi-\chi$&\\
\hline
$A_1$&1&1&1&1&0\\
$A_2$&1&1&-1&\-1&0\\
$B_1$&1&-1&1&-1&1\\
$B_2$&1&-1&-1&1&1\\
\hline
\end{tabular}\\
\parbox{73mm}{\caption{Irreducible representations of group $C_{2v}$\labeltbl{tbl:C2v}}}
\end{table}

%C6v
\begin{table}[htp!]
\begin{tabular}{l|cccccc|r}
\hline
\hline
$C_{6v}$&$E$&$C_6$&$C_6^5$&$C_3$&$C_3^2$&$C_2$&$p$\\
$R\chi$&$\chi$&\Rchi{}&\Rchi{5}&\Rchi{2}&\Rchi{4}&$\pi+\chi$&\\
\hline
$A_1$&1&  1&  1&  1&  1&  1&0\\
$A_2$&1&  1&  1&  1&  1&  1&0\\
$B_1$&1& -1& -1&  1&  1& -1&1\\
$B_2$&1& -1& -1&  1&  1& -1&1\\
$E_1$&\I{}{}&\M{}{-}{}{}&\M{}{}{-}{}&\M{-}{-}{}{-}&\M{-}{}{-}{-}& \I{-}{-}&1\\
$E_2$&\I{}{}&\M{-}{-}{}{-}&\M{-}{}{-}{-}&\M{-}{}{-}{-}&\M{-}{-}{}{-}& \I{}{}&0\\
\hline
$R$& $\sigma_v$ & $\sigma_v'$ & $\sigma_v''$ & $\sigma_d$ & $\sigma_d'$ & $\sigma_d''$ & $p$\\
$R\chi$ &$2\pi-\chi$&\Schi{2}&\Schi{4}&\Schi{}&$\pi-\chi$&\Schi{5}&\\
\hline
$A_1$& 1&  1&  1&  1&  1&  1&0\\
$A_2$&-1& -1& -1& -1& -1& -1&0\\
$B_1$& 1&  1&  1& -1& -1& -1&1\\
$B_2$&-1& -1& -1&  1&  1&  1&1\\
$E_1$&\I{}{-}&\M{-}{}{}{}&\M{-}{-}{-}{}&\M{}{}{}{-}&\I{-}{}& \M{}{-}{-}{-}&1\\
$E_2$&\I{}{-}&\M{-}{-}{-}{}&\M{-}{}{}{}&\M{-}{}{}{}&\I{}{-}& \M{-}{-}{-}{}&0\\
\hline
\end{tabular}\\
\parbox{73mm}{\caption{Irreducible representations of group $C_{6v}$ \labeltbl{tbl:C6v}}}
\end{table}



\begin{thebibliography}{References}
\bibitem{Launay1991} J. M. Launay and M. LeDourneuf, Chem. Phys. Lett. \textbf{169}, 473 (1990); J. M. Launay, Theor. Chim. Acta \textbf{79}, 183 (1991).

\bibitem{Lepetit1991}B. Lepetit and J. M. Launay, J. Chem. Phys. \textbf{95}, 5159 (1991).
\bibitem{PackParker1987} R. T Pack, G. A. Parker, J. Chem. Phys. \textbf{87}, 3888 (1987).
\bibitem{PackParker1989} R. T Pack, G. A. Parker, J. Chem. Phys. \textbf{90}, 3551 (1989).
\bibitem{Bacic1990}Z. Ba\u{c}i\'{c}, J. D. Kress, G. A. Parker, and R. T Pack, J. Chem. Phys. \textbf{92}, 2344 (1990).

\bibitem{Kress1990}J. D. Kress, Z. Ba\u{c}i\'{c}, G. A. Parker, and R. T Pack,
J. Chem. Phys. \textbf{94}, 8055 (1990).

\bibitem{Launay1989} J. M. Launay and M. LeDourneuf, Chem. Phys. Lett. \textbf{163}, 178 (1989).

\bibitem{Linderberg1989} J. Linderberg, S. B. Padkjaer, Y. Ohrn, and B. Vessal, J. Chem. Phys. \textbf{90}, 6254 (1989).

\bibitem{Kuppermann1989}A. Kuppermann and P. G. Hipes, J. Chem. Phys. \textbf{84}, 5962 (1986); S. A. Cucaro, P. G. Hipes, and A. Kuppermann, Chem. Phys. Lett. \textbf{154}, 155(1989), \textbf{157}, 440 (1989).

\bibitem{Lagana1991}(a) A. Lagan\`a, R. T Pack, and G. A. Parker, Faraday Discuss. Chem. Soc. \textbf{84}, 409 (1987); (b) \textbf{91}, 386 (1991).

\bibitem{Parker1989} G. A. Parker, R. T Pack, A. Lagan\`a, B. J. Archer, J. D. Kress,a and Z. Ba\u{c}i\'{c}, in textit{Super Computer Algorithms for Reactivity, Dynamics and Kinetics of Small Molecules}, edited by A. Lagan\`a (Kluwer, Holland, 1989), p. 105.

\bibitem{Iyengar1999} S. S. Iyengar, G. A. Parker, D. J. Kouri, D. K. Hoffman, J. Chem. Phys. \textbf{110}, 10283 (1999).

\bibitem{vardaf} D.K.~Hoffman, T.L.~Marchioro II, M.~Arnold, Y.~Huang, W.~Zhu
and D.J.~Kouri, {\it J. Math. Chem.} {\bf 20} (1996) 117-140.

\bibitem{PDAF-vardaf} D.J.~Kouri, D.K.~Hoffman, K. Zhang, G. A. Parker, and
S. S. Iyengar.

\bibitem{WithinDVR}See Ref. \cite{Bacic1990} and therein.

\bibitem{Kress1989}J. D. Kress, Z. Ba\u{c}i\'{c}, G. A. Parker, and R. T Pack, Chem. Phys. Lett. \textbf{157}, 484 (1989).

\bibitem{Kress1990A}J. D. Kress, R. T Pack, G. A. Parker,
Chem. Phys. Lett. \textbf{170}, 306 (1990).

\bibitem{Kress1991}J. D. Kress, Chem. Phys. Lett. \textbf{179}, 510 (1991).

\bibitem{Wolniewicz1989}L. Wolniewicz, J. Chem. Phys. \textbf{90}, 371 (1989);

\bibitem{FH2_T5A} F. B. Brown, R. Steckler, D. W. Schwenke, D. G. Truhlar, and B. C. Garrett, J. Chem. Phys. \textbf{82}, 188 (1985); R. Steckler, D. G. Truhlar, and B. C. Garrett, \textit{ibid}. \textbf{82}, 5499 (1985).

\bibitem{Zhang1991}J. Z. H. Zhang, Chem. Phys. Lett. \textbf{181}, 63 (1991).
\bibitem{Yu1989}C. H. Yu, D. J. Kouri, M. Zhao, D. G. Truhlar, and D. W, Schwenke, Chem. Phys. Lett. \textbf{157}, 491 (1989); Int. J. Quantum Chem. Symp. \textbf{23}, 45 (1989).

\bibitem{Manolopoulos1990} D. E. Manolopoulos, M. D'Mello, and R. E. Wyatt, J. Chem. Phys. \textbf{93}, 403 (1990).
\bibitem{Neuhauser1990} D. Neuhauser, M. Baer, R. S. Judson, and D. J. Kouri, Chem. Phys. Lett. \textbf{169}, 372 (1990); D. Neuhauser, R. S. Judson, R. L. Jaffe, M. Baer, and D. J. Kouri, \textit{ibid}. \textbf{176}, 546 (1991).
\bibitem{Weissbluth1978} M. Weissbluth, Atoms and Molecules, Academic Press, INC. (London) Ltd, 1978.
\end{thebibliography}


% Set the ending of a LaTeX document
\end{document}

