Chinese Journal of Chemical Physics  2017, Vol. 30 Issue (6): 761-770

#### The article information

Xue-ming Li, Zhi-gang Sun

An Exact Propagator for Solving the Triatomic Reactive Schrödinger Equation

Chinese Journal of Chemical Physics, 2017, 30(6): 761-770

http://dx.doi.org/10.1063/1674-0068/30/cjcp1711220

### Article history

Accepted on: December 25, 2017
An Exact Propagator for Solving the Triatomic Reactive Schrödinger Equation
Xue-ming Lia,b, Zhi-gang Suna
Dated: Received on November 15, 2017; Accepted on December 25, 2017
a. State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical Computational Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian 116023, China;
b. University of Chinese Academy of Sciences, Beijing 100049, China
*Author to whom correspondence should be addressed. Zhi-gang Sun, E-mail:zsun@dicp.ac.cn
Abstract: The exact short time propagator, in a form similar to the Crank-Nicholson method but in the spirit of spectrally transformed Hamiltonian, was proposed to solve the triatomic reactive time-dependent schrödinger equation. This new propagator is exact and unconditionally convergent for calculating reactive scattering processes with large time step sizes. In order to improve the computational efficiency, the spectral difference method was applied. This resulted the Hamiltonian with elements confined in a narrow diagonal band. In contrast to our previous theoretical work, the discrete variable representation was applied and resulted in full Hamiltonian matrix. As examples, the collision energy-dependent probability of the triatomic H+H2 and O+O2 reaction are calculated. The numerical results demonstrate that this new propagator is numerically accurate and capable of propagating the wave packet with large time steps. However, the efficiency and accuracy of this new propagator strongly depend on the mathematical method for solving the involved linear equations and the choice of preconditioner.
Key words: Time-dependent wavepacket method    Spectral difference    Spectrally transformed Hamiltonian    Exact short time propagator    Reactive scattering
Ⅰ. INTRODUCTION

The application of time-dependent quantum wave packet (TDWP) method for studying reactive scattering processes becomes a more and more popular computational tool [1-5]. Usually the reactive scattering processes are simulated using the TDWP method following a similar procedure: starting with a well defined initial wave packet, then choose a suitable grid representation and an efficient time propagator to evolve the wavepackt. Finally, we are able to extract the scattering attributes from the evolving wave packet. Thus, we roughly have two aspects to improve the numerical efficiency of the TDWP method. One is to improve the numerical efficiency of grid points, and the other is to improve the numerical efficiency of the time propagator.

In the past decades, there have been much effort on developing an efficient grid representation method. Since the generalized discrete variable representation (DVR) [6] and fast Fourier transformation (FFT) [7] were introduced into the chemical dynamics field in the 1980s, many kinds of DVR methods based upon classical polynomials have been proposed for further improving the numerical efficiency [8-13]. With these advances, now accurate rovibrational eigenvalues of 12D molecular potential energy surface without any approximation can be found [14], and molecular reactive scattering dynamics can be simulated including nine degrees of freedom [4]. However, DVR methods essentially are spectral and global methods which are utilized with classical polynomials, thus they usually produce a full matrix for a single degree of freedom.

Different from the DVR method with full Hamiltonian matrix, the Hamiltonian matrix constructed by the local methods, such as the finite element and finite difference method [15], can generate banded matrices. At the same time, the potential energy retains the diagonal form. This can save much computational effort. However, low order finite element and finite difference method converge slowly with increasing grid points, unlike the exponential convergence of the DVR method and spectral method [16, 17]. The spectral elements [18-20] and the spectral difference [21-27] are put forward to improve the numerical efficiency. The spectral difference methods can be either realised by using the distributed approximating functional (DAF) approach [24-26] or through the classical higher order finite difference formula. The spectral difference method accomplished by the sum acceleration method based upon Sinc DVR method [21, 22, 28] essentially is the same as the classical higher order finite difference method. The spectral difference method improves the convergence of the grid points substantially, and eventually leads to spectrally exponentially convergence in some cases. At the same time, the Hamiltonian matrix keeps banded matrix with increasing orders, except that with increasing width of the band. As we will see in the following, the efficiency of the grid points is nearly the same as that of the grid points in DVR when the spectral difference method with order high enough is applied. Thus the computational effort for evaluating the action of the Hamiltonian operator on the wave function reduces much.

To evolve the wave function, the exponential evolution operator (e$^{-i\hat{H}t/\hbar}$ where $\hat{H}$ is the Hamiltonian) must be approximated. In the past decades, there is much interest in investigating how to efficiently propagate the wave function but with high accuracy. The time propagation can be realized by expanding the exponential operator in terms of orthogonal polynomials such as Chebyshev polynomials [29] or Faber polynomials [30, 31], and so forth [28, 32-34]. Although the adaptation of orthogonal polynomials leads to a high accuracy, the required iteration numbers have a direct ratio with the complication of the system. An alternative and popular propagator is based on the split operator (SO) [35-38] algorithm. One of the interesting properties of the SO is that it conserves the norm of the wave packet, even when an unphysically large time step is used. As a result, the propagation is exceedingly stable. Currently, most reactive scattering reactions such as H+HBr [3], N+H$_2$ [39], O+O$_2$ [40] are calculated by the SO propagator.

Besides the above mentioned methods, there are some other popular propagators. In 1998, Gray and Balint-Kurti first introduced the Chebyshev real wave packet (CRWP) operator to deal with the reactive scattering calculation [41]. The interesting thing of the CRWP method is that only real arithmetic is carried out, which not only improves the numerical efficiency but also reduces memory requirements [42]. The intriguing aspect of the CRWP approach is that it uses the concept of the transformed Hamiltonian $H$, in which the Hamiltonian $H$ is modified by $f(H)$, where the choice for $f$ is a $arccos$ maping [43]. In this way, the Chebyshev polynomials form a complete set and the propagation in the Chebyshev order ($k$) domain can be easily and accurately realized by the three-term Chebyshev recursion relation. Since its introduction, it has been used to study many reactive scattering problems, including Br+H$_2$ [44], O+O$_2$ [5], F+HCl [45] etc. The short-time iterative Lanczos method also has popular application where the spectral range of the Hamiltonian is large [46, 47]. Three-dimension reactive scattering processes, such as H+O$_2$ [48], F+H$_2$ [49] and O+HCl$_2$ [50] reaction, have been studied by the short-time iterative Lanczos method. Sun et al. [51] have demonstrated that the second-order SO method is more efficient when the wavepacket is propagated on the flat potential energy surface (PES), while for the PES with a deep potential, the CRWP method can perform more efficient for the long-lived resonance states. However, for these reactions with deep wells to obtain fully converged scattering information, the CRWP method also requires a large number of iterations. Therefore, it is highly desirable to explore other propagation schemes in order to achieve a better efficiency on studying such complex formation reactions.

Based on the concept of spectrally transformed Hamiltonian, Chen-Guo introduced a new propagator in a form of the Crank-Nicholson method [52] for solving the vibrational state [53, 54]. The numerical examples revealed that this new short time propagator has machine round off accuracy, but with arbitrarily larger time step. In 2011, Sun and Yang moved a step forward by applying the spectrally transformed propagator in the form of the Crank-Nicolson method (named as SCN propagator hereafter) for solving the one-dimension quantum scattering process passing over a barrier [55]. The result revealed again that the SCN propagator is exact with machine accuracy and the most impressive thing is that it can propagate the wave packet with large time-step $\Delta{t}$=100 a.u. without losing accuracy, in contrast with the second-order SO propagator with time-step $\Delta{t}$=10 a.u. and the traditional CN propagator with time-step small as $\Delta{t}$=2 a.u. to obtain the results with error of about 1%. We note that due to the unconditional stability and unitary, especially the good ability for treating the presence of Coulomb singularity, the Crank-Nicholson method has been widely used for solving the electronic TDSE involving the Coulomb potential [55-57]. However, because of the lower efficiency of the lower order finite difference method and its large dispersive error, the CN propagator has not arisen much attention in the reactive scattering field.

In this work, the SCN propagator was applied to real triatomic reactive processes. In order to improve the numerical efficiency, the spectral difference method was adopted. At the same time, we applied a suitable preconditioner to accelerate the convergence for solving the linear equation arising for propagating the wave function with the SCN propagator. The reactive scattering reactions, H+H$_2$ and O+O$_2$, were taken as the numerical example for illustrating the properties of the SCN propagator.

Ⅱ. THEORY A. Hamiltonian and basis for triatomic reactive scattering

For the triatomic reactive scattering in the reactant coordinates in the body-fixed (BF) frame, the Hamiltonian for a given total angular momentum $J$ can be written as

 $\begin{eqnarray} \hat{H} &=& -\frac{\hbar^2}{2\mu_{R_\alpha}}\frac{\partial^2}{\partial R _\alpha^2}-\frac{\hbar^2}{2\mu_{r_\alpha}}\frac{\partial^2}{\partial r_\alpha^2}+ \frac{(\hat{J}-\hat{j})^2}{2\mu_{r_\alpha}R_\alpha^2}+\nonumber\\ &&{}\frac{\hat{j}^2}{2\mu_{r_\alpha}r_\alpha^2}+ V \label{har} \end{eqnarray}$ (1)

where $\hat{J}$ is the total angular momentum operator, $\hat{j}$ is the rotational angular momentum operator of BC, $\mu_{R_\alpha}$ is the reduced mass between the center of mass of A and BC, and $\mu_{r_\alpha}$ is the reduced mass of BC. The Hamiltonian and wave packet are discretized in the same mixed discretized representation. For the two radial degrees of freedom ($R$ and $r$ direction), the spectral difference method and the Sinc DVR method were applied for evaluating the kinetic energy, and the spherical harmonic basis $y_{jK}(\theta_{\alpha})$ is used for representing the angular degree of freedom, where DVR was adopted for switching between the DVR and finite basis representation for evaluating the action of potential and rotational energy on the wave function, respectively.

The time-dependent wave packet in the BF frame can be expressed as [58]

 $\begin{eqnarray} \Psi^{JM\epsilon}(\vec{R}_\alpha, \vec{r}_\alpha) = \sum\limits_{K_\alpha}\bar{D}^{J\epsilon\ast}_{MK_\alpha}(\vec{\Omega}_\alpha)\psi_\alpha(R_\alpha, r_\alpha, \theta^{K_\alpha}_\alpha;K_\alpha)\nonumber\\ \end{eqnarray}$ (2)

where $\bar{D}^{J\epsilon\ast}_{MK_\alpha}(\vec{\Omega}_\alpha)$ is the parity-adapted normalized rotation matrix, depending only on the Euler angles $\vec{\Omega}_\alpha$,

 $\bar{D}^{J\epsilon\ast}_{MK_\alpha}(\vec{\Omega}_\alpha) = \sqrt{\frac{2J+1}{8\pi(1+\delta_{K_\alpha0)}}}\cdot\nonumber\\ \;\;\;\;\;\;\left[D^{J\ast}_{MK_\alpha}(\vec{\Omega}_\alpha)+ \epsilon(-1)^{J+K_\alpha}D^{J\ast}_{M-K_\alpha}(\vec{\Omega}_\alpha)\right]\quad\quad$ (3)

where $\epsilon$=$(-1)^{j+l}$ is the parity of the system, with $l$ is the orbital angular momentum quantum number, being $K_\alpha$ is the projection of the total angular momentum $J$ on the BF $z$ axis and ${D}^{J\epsilon}_{MK_\alpha}(\vec{\Omega}_\alpha)$ is the Wigner rotation matrix. The wave function $\psi_\alpha(R_\alpha$, $r_\alpha$, $\theta^{K_\alpha}_\alpha$; $K_\alpha$) only depends on the internal coordinates ($R_\alpha$, $r_\alpha$, $\theta^{K_\alpha}_\alpha$) and $K_\alpha$, and can be expanded as

 $\begin{eqnarray} \psi_\alpha(R_\alpha, r_\alpha, \theta^{K_\alpha}_\alpha;K_\alpha) &=& \sum\limits_{n, m, j}\Big\{F^{K_\alpha}_{nmj} u_n(R_\alpha)\phi_m(r_\alpha)\cdot\nonumber\\&&{}y_{jK_\alpha}(\theta_\alpha)\Big\} \end{eqnarray}$ (4)

where $n$ and $m$ are the translational basis labels, $u_n(R_\alpha)$ and $\phi_m(r_\alpha)$ are the translational basis functions for $R$ and $r$ respectively, and $y_{jK_\alpha}$=$\sqrt{\displaystyle\frac{(2j+1)}{4\pi}}d_{K_{\alpha 0}}^j$ are spherical harmonics. $d_{K_{\alpha 0}}^j$ is a reduced Wigner rotational matrix [59] with $K_v$=0.

B. Construction of the initial wave packet

In a wave packet calculation, before its propagation, we first need to construct an initial wave packet. For a triatomic reaction scattering, the initial wave packet in space-fixed (SF) representation ($v_0$, $j_0$, $l_0$) can be constructed simply as

 $\begin{eqnarray} \psi_{\alpha v_{0} j_{0}l_{0}}^{JM\epsilon}(t=0)=G(R_\alpha)\phi_{v_{0}j_{0}} (r_{\alpha})|JMj_{0}l_0\epsilon\rangle \label{init1} \end{eqnarray}$ (5)

where $|JMj_{0}l_0\epsilon\rangle$ is total angular momentum eigenfunction in the SF representation with parity of system $\epsilon$=$(-1)^{j_0+l_0}$, $\phi_{v_{0}j_{0}(r_{\alpha})}$ is the rovibrational eigenfunction for molecule BC, and $G(R_\alpha)$ is Gaussian wave packet

 $\begin{eqnarray} G(R_\alpha)=\left(\frac{2}{\pi {\sigma}^2}\right)^{1/4}\exp\left[-\frac{(R_\alpha-R_{\alpha}^c)}{\sigma}^2 -ik_c R_\alpha\right] \end{eqnarray}$ (6)

In order to propagate the wave packet in the BF representation, one should transform $|JMj_{0}l_0\epsilon\rangle$ in Eq.(5) to its BF representation counterpart as [60]

 $\begin{eqnarray} |JMj_{0}l_0\epsilon\rangle &=&\sum\limits_{K_{\alpha}\ge0}C_{l_{0}K_{\alpha}}^ {Jj_{0}\epsilon}|JMj_{0}K_{\alpha}\rangle\nonumber\\ &=&\sum\limits_{K_{\alpha}\ge0}C_{l_{0}K_{\alpha}}^ {Jj_{0}\epsilon}D^{J\epsilon\ast}_{MK_\alpha}(\Omega_\alpha)y_{j_{0}K_{\alpha}} \end{eqnarray}$ (7)

where $C_{l_{0}K}^{Jj_{0}\epsilon}$ is the parity-adapted orthogonal transform matrix between the SF and BF representations [61-64] and given by

 $\begin{eqnarray} C_{lK}^{Jj\epsilon}=\sqrt{\frac{2l+1}{2J+1}}\sqrt{(2-\delta_{K, 0})} \langle jKl_0|JK\rangle \end{eqnarray}$ (8)

where $\langle\cdots\rangle$ is the Clebsch-Gordan coefficient.

C. The spectral difference method for radial degrees of freedom

The spectral difference method was realized by the weighted cardinal basis functions using the DAF method [27]:

 $\begin{eqnarray} \chi_{n}(x, \Delta x, \sigma) = C_{n}(x, \Delta x)\omega_{n}(x, \sigma) \end{eqnarray}$ (9)

with $\omega_{n}(x, \sigma)$ defined as a Gaussian function form,

 $\begin{eqnarray} \omega_{n}(x, \sigma) = \rm{e}^{-(x-x_{n})^2/2\sigma^2} \end{eqnarray}$ (10)

and $C_{n}(x, \Delta x)$ is a sinc function [65, 66],

 $\begin{eqnarray} C_{n}(x, \Delta x)={\rm sinc}[(x-x_{n})/\Delta x] \end{eqnarray}$ (11)

By decaying to zero, the weight Gaussian function $\omega_n$($x$, $\sigma$) can make derivatives of $\chi_n$($x$, $\Delta x$, $\sigma$) at $x_m$ have no coupling with the grid point at $x_n$ at far distance. The parameter $\sigma$ controls the rate of decay. Wei et al. [26] have shown the $\sigma$ should be chosen such that points lying just beyond the bandwidth have a negligible contribution to the derivative. So by adjusting $\sigma$ appropriately, we can just ignore the terms of the derivative outside the bandwidth. As the terms in the differential equation without derivatives can translate into diagonal matrices, we can get the sparser Hamiltonian matrix by the spectral difference method than by the DVR method. The Hamiltonian matrix can be obtained as follows:

 $\begin{eqnarray} T_{m, n} &=& \Delta''_{|m-n|}\omega_{n}(x, \sigma)+2\Delta'_{|m-n|}\omega'_{n}(x, \Delta x)+\nonumber\\ &&{}\delta_{m, n}\omega''_{n}(x, \sigma) \end{eqnarray}$ (12)

where $\omega'_{n}(x, \Delta x)$ and $\omega''_{n}(x, \Delta x)$ is the first and second derivatives of the Gaussian weighted function respectively. And $\Delta_{|m-n|}'$ and $\Delta_{|m-n|}''$ are the first and second derivatives of sinc function given by

 $\Delta _{|m - n|}' = \left\{ {\begin{array}{*{20}{c}} {0,\;\;\;\;\;\;\;\;\;\;\;\;\; |m - n| = 0}\\ {\displaystyle\frac{{{{( - 1)}^{|m - n| + 1}}}}{{\Delta x|m - n|}}, |m - n| \ne 0} \end{array}} \right.$ (13)
 $\Delta _{|m - n|}'' = \left\{ {\begin{array}{*{20}{c}} {\displaystyle\frac{{ - {\pi ^2}}}{{3\Delta {x^2}}}, }\;\;\;\;\;\;\;\;\;\;{|m - n| = 0}\\ {\displaystyle\frac{{2{{( - 1)}^{|m - n| + 1}}}}{{\Delta {x^2}|m - n{|^2}}}, }{|m - n| \ne 0} \end{array}} \right.$ (14)

It should be mentioned that the accuracy of the matrix depends on the width of the parameter $\sigma$ in the Gaussian function when the width is not large enough. This is similar to the application of classical high order finite-difference method. When limited grid points are applied, higher order method would lead to results of higher accuracy before obtaining converged numerical results. The effects of different widths of $\sigma$ on the numerical accuracy would be discussed in Section III.

D. Wave packet propagation using the SCN approach

The evolution of the wave packet under the time-independent Hamiltonian can be written as

 $\begin{eqnarray} |\Psi(r, \Delta{t})\rangle={\rm exp}(-iH\Delta{t})|\Psi(r, 0)\rangle \end{eqnarray}$ (15)

Usually, the propagator ${\rm exp}(-iH\Delta{t})$ in a short-time form has low order accuracy. For example, in the traditional CN scheme, an approximate expression of the propagator is written as the following formula:

 $\begin{eqnarray} \exp(-iH\Delta t)\approx\frac{1-iH\Delta{t}/2}{1+iH\Delta{t}/2}+O(\Delta{t}^{3}) \label{cn} \end{eqnarray}$ (16)

This is unitary and of second order accuracy, so relative small time-step $\Delta{t}$ has to be adopted because of its large dispersive error.

Given a spectrally transformed Hamiltonian $H'$ which has the relationship with the original Hamiltonian $H$ as

 $\begin{eqnarray} H'=2{\rm arctan}(\Delta{t}H/2)/\Delta{t} \end{eqnarray}$ (17)

The spectrally transformed Hamiltonian $H'$ has the same eigenfunctions $|\varphi_{n}(r)\rangle$ with the original Hamiltonian $H$ but with different eigenvalues. However, these eigenvalues $E'_{n}$ have the one-to-one relationship with the original eigenvalues $E_{n}$ with an analytic form:

 $\begin{eqnarray} &&\begin{array}{l} H\left| {{\varphi _n}(r)} \right\rangle = {E_n}\left| {{\varphi _n}(r)} \right\rangle \\ H'\left| {{\varphi _n}'(r)} \right\rangle = E{'_n}\left| {\varphi'_n(r)} \right\rangle \end{array} \label{EE} \end{eqnarray}$ (18)
 $\begin{array}{l} \left| {{\varphi _n}^\prime (r)} \right\rangle = \left| {{\varphi _n}(r)} \right\rangle {\rm{ }}\\ E' = 2\arctan (\Delta tE/2)/\Delta t \end{array}$ (19)

With the spectrally transformed Hamiltonian, we obtain the new time-dependent Schrödinger equation

 $\begin{eqnarray} i\hbar\frac{\partial}{\partial{t}}|\Psi(r, t)\rangle=H'|\Psi(r, t)\rangle \end{eqnarray}$ (20)

For this new TDSE with above spectrally transformed Hamiltonian, we obtain a new propagator based on the Hamiltonian $H'$ with the time-step size $\Delta{t}$:

 $\begin{eqnarray} \exp(-iH'\Delta{t})\equiv\frac{1-i\tan H'\Delta{t}/2}{1+i\tan H'\Delta{t}/2}\equiv\frac{1-iH\Delta{t}/2}{1+iH\Delta{t}/2} \end{eqnarray}$ (21)

The last term above is exactly the same as the right side of Eq.(16). Now this new propagator is an absolutely exact propagator for the TDSE with the spectrally transformed Hamiltonian $H'$. The most wonderful characteristic of the proposed SCN propagator is that it is without any approximation and is unconditionally stable and accurate with huge time-step to propagate a proper wave packet.

E. Calculation of the reactive scattering probability by the SCN propagator

For a reactive scattering process, usually the reaction probability can be obtained from the TDSE by the flux formalism as [67]:

 $\begin{eqnarray} P(E)&=&|S_{fi}(E)|^2\nonumber\\ &=&\frac{2\pi}{m}{\rm Im}\left[A_{fi}(R_{\beta}, E)^{*} \frac{\rm{d}}{\rm{d}R_{\beta}}A_{fi}(R_{\beta}, E)\right]\bigg|_{R_{\beta}=R_L}\quad\quad \end{eqnarray}$ (22)

where time-independent energy-resolved $A_{fi}(R_{\beta}, E)$ are calculated by the Fourier transform of the wavefunction with respect to time $t$.

 $\begin{eqnarray} A_{fi}(R_{\beta}, E)=\frac{1}{2\pi a_{\alpha i}(E)} \int_{-\infty}^{\infty}{\rm{d}}t {\rm{e}}^{iEt} A_{fi}(R_{\beta}, t) \end{eqnarray}$ (23)

The meaning of the relevant symbol one may refer Ref.[67]. While using the SCN propagator, the reaction probability can be calculated as

 $\begin{eqnarray} P(E')&=&|S_{fi}'(E')|^2\nonumber\\ &=&\frac{2\pi}{m}f(\theta) {\rm Im}[A_{fi}(R_{\beta}, E')^{*}]\cdot\nonumber\\ &&{}\hspace{-0.12cm} \frac{\rm{d}}{\rm{d}R_{\beta}}A_{fi}(R_{\beta}, E')|_{R_{\beta}=R_L} \end{eqnarray}$ (24)

where the time-independent energy-resolved $A_{fi}$($R_{\beta}$, $E'$) are calculated by the Fourier transform with respect to time $t$ again.

 $\begin{eqnarray} A_{fi}(R_{\beta}, E')=\frac{1}{2\pi a_{\alpha i}(E')} \int_{-\infty}^{\infty}{\rm{d}}t {\rm{e}}^{iE't} A_{fi}(R_{\beta}, t) \end{eqnarray}$ (25)

According to the relation between $E'$ and $E$ given by Eq.(19), we can get the factor $f(\theta)$

 $\begin{eqnarray} f(\theta)=\left(\frac{{\rm{d}}E'}{{\rm{d}}E}\right)^2=\left(\frac{1}{1+\Delta{t}E/2}\right)^2 \end{eqnarray}$ (26)

The derivation of $f(\theta)$ is similar to that in the CRWP approach. In the CRWP method, because of the introducing of the spectrally transformed Hamiltonian [41]:

 $\begin{eqnarray} H'=2\arccos(\Delta{t}H/2)/\Delta{t} \end{eqnarray}$ (27)

brings the factor

 $\begin{eqnarray} f(\theta)=\left(\frac{{\rm{d}}E'}{{\rm{d}}E}\right)^2=\left(\frac{1}{1-\sqrt{\Delta{t}E/2}}\right)^2 \end{eqnarray}$ (28)

of the reaction probability respect to the original reaction probabilities.

We estimate the error from the total reaction probabilities and the error is defined as

 $\begin{eqnarray} \sigma=\frac{1}{M}\sum\limits_{k=1}^{M}|P(E_k)-P^0(E_k)|/P^0(E_k) \label{error} \end{eqnarray}$ (29)

where $M$ is the number of the collision energies $E_k$ and $P^0(E_k)$ is the converged results.

F. Preconditioner for solving the linear equation

Using the SCN propagator, the most important issue is that how to efficiently solve the equation

 $\begin{eqnarray} \Psi(r, t+\Delta{t})=\frac{1-iH\Delta{t}/2}{1+iH\Delta{t}/2} \Psi(r, t) \label{qiuni} \end{eqnarray}$ (30)

For avoiding evaluation of inverse of the Hamiltonian matrix, Eq.(30) can be written as

 $\begin{eqnarray} (1+iH\Delta{t}/2)\Psi(t+\Delta{t})=(1-iH\Delta{t}/2)\Psi(t) \end{eqnarray}$ (31)

which now is a linear equation of type $\hat A\chi$=$\varphi$. According to Eq.(1), the operator $\hat{A}$ has the form:

 $\begin{eqnarray} \hat A=1+i(\hat T_R+\hat T_r+\hat T_j+V)\Delta{t}/2 \end{eqnarray}$ (32)

where $V$ is the PES of the triatomic system, $\hat T_j$=[($\hat{J}$$-$$\hat{j})^2$)/2$\mu_{R_\alpha}R_\alpha^2]$+$(\hat{j}^2/2\mu_{r_\alpha}r_\alpha^2)$ is the angular momentum operator. $\hat T_R$=$-(1/2\mu_{R_\alpha})(\partial^2/\partial{R_\alpha^2}$) and $\hat T_r$=$-(1/2\mu_{r_\alpha})(\partial^2/\partial{r_\alpha^2})$ are the radial kinetic operators of the triatomic system. There are many methods for solving the linear equation such as the BCG [68], CGS [69] and QMR [70] method. In this present work, TFQMR [71] method is used to solve the linear equation.

The iterative solvers for a linear equation can be accelerated by applying a suitable preconditioner, which means that the solution of the equation $\hat A\chi(x, t)$=$\varphi(x, t)$ can be efficiently obtained by solving the equation $\hat A_0^{-1}\hat A\chi(x, t)$=$\hat A_0^{-1}\varphi(x)$. The key properties of the preconditioner $\hat A_0^{-1}$ is that it should be a good approximation of $\hat A$, that is $\hat A_0^{-1}\hat A$$\approx$$I$ ($I$ means unit matrix). The preconditioner $\hat A_0^{-1}$ we used is the kinetic energy in the basis representation of the operators and potential energy as a constant, similar to that proposed by Peskin et al. [72]. The preconditioner can be implemented with two steps, the first is the kinetic preconditioner. In this step, the $\hat A_0$ has the form:

 $\begin{eqnarray} \hat{A_0}&=&1+i[\hat T_R+\hat T_r+(\hat{J}-\hat{j})^2/2\mu_{R_\alpha}{R_{0}^2}+\nonumber\\ &&{}\hat{j}^2/2\mu_{r_\alpha}{r_{0}^2}+V_0]\Delta{t}/2 \end{eqnarray}$ (33)

where $r_0$, $V_0$, $R_0$ are constant. So the equation becomes:

 $\hat{B}\chi(x, t)=\hat{A}_0^{-1}\varphi(x, t)$ (34)
 $\hat{B}=\hat{A}_0^{-1}\hat{A}=\hat{I}+\hat{A}_0^{-1}(\hat{A}-\hat{A_0})$ (35)

The action of the operator $\hat{B}$ on a function $\chi(x, t)$ can be realized first to multiply the local operator $(\hat{A}-\hat{A_0})$ to the function $\chi(x, t)$ on the grid representation. Then by transforming the function $\chi(x, t)$ into the momentum representation $\chi(p_x, p_t)$, we can evaluate the action of $\hat{A}_0^{-1}$ on the wave function by a simple scalar multiplication. Finally, by transforming back the $\chi$($p_x$, $p_t$) into the original space representation $\chi$($x$, $t$), we can continue the propagation for the next time step. Since both the local and the differential operators are diagonal in its respective representation, the major computational task becomes the representation switch of the wave function, which is the same as the nonpreconditioned $\hat{A}$ used. Therefore, in our calculation, the introduction of the preconditioner involves little numerical effort.

The second step of the proconditioner is the diagonal part of the Hamiltonian matrix in space representation. After the first step is done, the matrix representation of $\hat{B}$ is diagonally dominant. The diagonal part of $\hat B$ in the grid representation $\chi(x, t)$ is:

 $\begin{eqnarray} {\rm diag}(\hat B)=I+{\rm diag}(\hat A_0^{-1})(\hat A-\hat A_0) \end{eqnarray}$ (36)

where ${\rm diag}(\hat A_0^{-1})$ in the grid representation is a constant given by the average value of Eq.(33). In this step, the inverse of the diagonal part of $\hat{B}$, ${\rm diag}(\hat{B})^{-1}$, multiply on both sides of the equation, and the equation eventually becomes

 $\begin{eqnarray} {\rm diag}(\hat{B})^{-1}\hat{B}\chi(x, t)={\rm diag}(\hat{B})^{-1}\hat{A}_0^{-1}\varphi(x, t) \end{eqnarray}$ (37)

where the left side of the equation is

 $\begin{eqnarray} {\rm diag}(\hat{B})^{-1}\hat{B}&=&[\hat{I}+{\rm diag}(\hat{A_0^{-1}})(\hat{A}-\hat{A_0)}]\cdot\nonumber\\ &&{}[\hat{I}+\hat{A_0^{-1}}(\hat{A}-\hat{A_0)})] \end{eqnarray}$ (38)

The advantage of diagonal-preconditioner is that it only requires an extra vector multiplication and is thus of negligible extra computational task [72]. However, it accelerates the convergence speed for solving the linear equation.

G. Absorbing potential

Finally, we use a damping function in the following form to avoid the wave function reaching the grid boundary along $r$,

 $\begin{eqnarray} D(r) = \left\{ \begin{array}{ll} \exp\left[-\Delta_{t}C_{a}\left(\displaystyle\frac{r-r_a}{r_b-r_a}\right)\right], r_a\leq r \leq r_b \\ \exp\left[-\Delta_{t}C_{b}\left(\displaystyle\frac{r-r_b}{r_{\rm end}-r_b}\right)\right] +\exp[-\Delta_{t}C_{a}], \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;r_b\leq r \leq r_{\rm end} \end{array} \right. \end{eqnarray}$ (39)

The first part $r_a$$\leq$$r$$\leq$$r_b$ is designed for the product with low kinetic energies but the second part $r_b$$\leq$$r $$\leq$$r_{\rm end}$ is used for the product with relative higher kinetic energies.

Ⅲ. RESULTS AND DISCUSSION A. H+H$_2$ (${v_0=0, j_0=0}$) reaction

The H+H$_2$ reactive scattering has been extensively studied in the past decades [73, 74]. And it is also a benchmark reaction for the development of new theoretical methods. With advance in experimental and theoretical method, an excellent agreement between theory and experiment has been accomplished for this reaction on the BKMP2 PES [75]. In this study, the BKMP2 PES is employed in the following calculation for numerical illustrations.

The total reaction probabilities will be calculated by three different methods. One applies the SCN with the sinc DVR method for radial degrees of freedom (SCN method), one applies the SCN with spectral difference for radial degrees of freedom (SCNDAF method) and one applies the second order SO with the Sinc DVR method for radial degrees of freedom (SO method).

The comparison of the convergence behaviours of the total reaction probabilities for the H+H$_2$ reaction with total angular momenta $J$=0 with three methods, as a function of time step, is plotted in FIG. 1. For SCNDAF, we choose $\sigma$=0.4 a.u., which can give converged results. The error is calculated by Eq.(29). The number $M$ of the measured collision energies $E_k$ is 61 and the measured energy $E_k$ evenly distributes in the range of [0.3, 0.9] eV. From this figure, we can see that even with large time step $\Delta{t}$=90 a.u., the SCN and SCNDAF propagators can give a converged result with error about 1%. However, the SO method only with time step small as $\Delta{t}$=10 a.u. can give results with such accuracy.

 FIG. 1 The convergence behaviors of the total reaction probabilities as a function of time step of the SCN, the SCNDAF, and the SO method for the H+H$_2$ reactive scattering with $J$=0.

FIG. 2 shows the error of the reaction probabilities as a function of the Gaussian parameter $\sigma$ in the spectral difference with different threshold for abandoning the matrix elements far away from the diagonal. The time step $\Delta{t}$=90 a.u. was applied. The parameter $\sigma$ determines the width of the matrix band thus is vital to the numerical efficiency. The result indicates that the estimated optimal $\sigma$=0.4 a.u. and the threshold $10^{-5}$ is an optimal choice for this reaction. For $\sigma$=0.4, in this reaction with 70 grid points for every degree of freedom, the band width is about 9 elements when the elements of values less than 1.0$\times$10$^{-5}$ are eliminated. Thus much of the computational effort was saved, as compared to the DVR method.

 FIG. 2 The error of the total reaction probabilities as a function of the Gaussian parameter $\sigma$ in the spectral difference for different cutoff threshold $\varepsilon$ with time step $\Delta{t}$=90 a.u. for the H+H$_2$ reactive scattering with $J$=0.

The number of iterations required for this reaction as a function of the time step is plotted in FIG. 3 using the SCN propagator. The red line represents the iterating numbers without the preconditioner. The green line represents the iterating numbers with the preconditioner using the Sinc DVR for the two radial degrees of freedom. The blue line represents the iterating numbers with preconditioner using the spectral difference method. FIG. 3 demonstrates the importance of the preconditioner to the efficiency for solving Eq.(30). FIG. 4 presents total reaction probabilities as a function of collision energy for the H+H$_2$ reaction with total angular momenta $J$=0 calculated by the SCNDAF propagator with large time step $\Delta{t}$=90 a.u., compared with the results obtaining using the SO method with time step as $\Delta{t}$=2.0 a.u. We can see the results calculated by the SCNDAF propagator agree well with the SO method even when the SCNDAF propagator with so large time steps.

 FIG. 3 The number of iterations required for solving the linear equation as a function of the time step for the H+H$_2$ reactive scattering with $J$=0.
 FIG. 4 Total reaction probabilities calculated by the SCNDAF method with large time step $\Delta{t}$=90.0 a.u., comparing with the results calculated by the second-order SO propagator of $\Delta{t}$=2.0 a.u. for the H+H$_2$ reactive scattering with $J$=0.

Therefore, for the H+H$_2$ reaction, the SCNDAF propagator with large time steps is capable of giving converged results. The relative efficiency of the SCNDAF and second-order SO propagator can be estimated from the effective time-step. Since with larger time-step, each time-step requires more iterations for the SCNDAF. An optimal choice would be $\Delta{t}$=90 a.u. From FIG. 3 every time step requires the number of iterations is about 7.8, thus the effective time step is $90/7.8$$\approx11.5 a.u. While for the second-order SO method, the time step is about 20 a.u. for obtaining converged results [51]. We may conclude that the efficiency of this new SCNDAF propagator is less effective than the SO method. On the other side, the SCNDAF propagator can explore the sparsity of the band structure of the matrix, unlike the SO method. On balance, the efficiency of the SCNDAF propagator is comparable with the second order SO method. Anyway, the SCNDAF propagator is less effective than the CRWP method for this simple reaction [51]. B. O+O_2 (v_0=0, j_0=0) exchange reaction The O+O_2 exchange reaction is a prerequisite for the formation of ozone in Earth's atmosphere [76]. There is much interest in the O_3 complex-forming reaction [77]. The state-to-state differential and integral cross section for the O+O_2 isotope exchange reaction at collision energies relevant to atmospheric conditions by quantum wave packet scattering calculations have been reported [78]. The comparison of the convergence behaviors of the total reaction probabilities as a function of collision energy for the O+O_2 reaction with total angular momenta J=0 for three methods is plotted in FIG. 5. The error is calculated by Eq.(29). The number M of the collision energies E_k is 700 which evenly distributes in the range of [0.02, 0.16] eV. From the figure, we can see that even with large time step \Delta{t}=100 a.u., the SCN propagator can give converged results with error less than 1%. In contrast, the SO propagator only with time step small as \Delta{t}=3 a.u. can give results with such accuracy.  FIG. 5 The convergence behaviors of the total reaction probabilities as a function of time step of the SCN, the SCNDAF, and the SO method for the O+O_2 reactive scattering with J=0. The left panel of FIG. 6 shows the error of the reaction probabilities as a function of the Gaussian parameter \sigma for the spectral difference with different cutoff thresholds for R degree of freedom with time step \Delta{t}=100 a.u. The plots suggest that \sigma=0.3 and the cutoff threshold as 10^{-6} are optimal choices. Total 255 grid points are used for this reaction for R degree of freedom. With the matrix elements of values less than 1.0\times10^{-6} removed, the band width of the kinetic matrix of R degree of the freedom is about 19 elements. Thus the matrix is much sparser than the original one and save much computational effort for evaluating the action of Hamiltonian matrix on the wave function.  FIG. 6 Left: The error of the reaction probabilities as a function of the Gaussian parameter \sigma in the spectral difference for different cutoff threshold \varepsilon for R degree of freedom with time step \Delta{t}=100 a.u. for the O+O_2 reactive scattering with total angular momentum J=0. Right: The error of the reaction probabilities as a function of the Gaussian parameter \sigma in the spectral difference for different cutoff threshold \varepsilon for r degree of freedom with time step \Delta{t}=100 a.u. for the O+O_2 reactive scattering with total angular momentum J=0. The right panel of FIG. 6 shows the error of the reaction probabilities as a function of the Gaussian parameter \sigma for different cutoff threshold for r degree of freedom with time step \Delta{t}=100 a.u. The plots indicate that \sigma=0.15 and the cutoff threshold 10^{-6} are optimal choices. Total 89 grid points of this reaction for r degree of freedom are used. With cutoff threshold as 1.0\times10^{-6}, the kinetic matrix for r degree of freedom has band width only about 12 elements. Thus the matrix is much sparser than the original one and saves much computational effort for evaluating the action of Hamiltonian matrix on the wave function. Since there are deep wells for the O+O_2 reaction, which require denser grid points than for the simple H+H_2 reaction. Therefore the spectral difference is of particular advantage for such kind of reaction, and saves much computational effort as comparing with the DVR method. The number of iterations required for solving the linear equation of this reaction as a function of the time step are plotted in FIG. 7. The red line represents the iteration numbers as a function of time step without preconditioner using the DVR method. The green line represents the iteration numbers as a function of time step with the preconditioner using the DVR method. The blue line represents the iteration numbers as a function of time-step with preconditioner using the spectral difference method. FIG. 7 demonstrates the accelerated convergence for adoption of preconditioner for solving Eq.(30). However, the acceleration is not so obvious as that for the H+H_2 reaction.  FIG. 7 The number of iterations required for solving the linear equation for the O+O_2 reaction as a function of the time step. FIG. 8 presents total reaction probabilities for the O+O_2 reaction with total angular momenta J=0 calculated by the SCNDAF propagator using the spectral difference with large time step \Delta{t} as 100 a.u., comparing with that using the SO method with \Delta{t}=3.0 a.u. We can see the results agree well with each other even SCNDAF propagator with large time step.  FIG. 8 Total reaction probabilities for the O+O_2 reaction with total angular momenta J=0 calculated by the SCNDAF method with large time step \Delta{t}=100 a.u., comparing with the second-order SO method with time step \Delta{t}=3.0 a.u. The reaction O+O_2 has long-lived resonance states in a deep potential well, which requires long time propagation. We may expect that the efficiency of this SCNDAF propagator is better than the second-order SO method. Such as with \Delta{t}=100 a.u. as from FIG. 7, the effective time step with the precontidioner is 100/14$$\approx$7.15 a.u. While for the second-order SO method, the time step is 3.0 a.u. for producing the results with error about 1%. Anyway, as comparing the effective time-step 7.0$\times$2 a.u. for the CRWP method, the efficiency of the SCNDAF method is not appealing.

Finally, we should note that the above effective time-step is estimated based upon the number of iterations for solving the linear equation. In practical calculations, the inverse of the total Hamiltonian matrix for a multi-dimensional reactive process usually is impossible to prepare in advance for the propagation, unlike that in the 1D case, one has to evaluate the inverse of the Hamiltonian matrix for each time step, as did in the current work. Each iteration may require more than one action of the Hamiltonian on the wave function, which depends on the mathematical method for solving the linear equation. Such as for the TFQMR [71] method applied in the present work, each iteration requires four times the action of the Hamiltonian on the wave function. As a result, presently the SCN method is not so appealing as the popular SO and CRWP propagators, unless highly efficient preconditioner can be found.

On one hand, the SCN has no restriction on the spectral range of the Hamiltonian, unlike the CRWP method, which can similarly exploit the band structure of the Hamiltonian matrix. On the other hand, although the split operator method has no restriction on the spectral range of the Hamiltonian matrix, it cannot exploit the band-structure of the matrix produced by the DAF method. Thus, in a calculation with many degrees of freedom or involving potential with singularity, such the Coulomb potential, which usually have large/huge spectral range, the SCN method may have particular advantages, as comparing with the Chebyshev and split operator method. This deserve exploration in the future.

Ⅳ. DISCUSS AND CONCLUSION

In this work, the SCN method, which has the same form as the classic Crank-Nicholson method, but in a spirit of spectrally transformed Hamiltonian, was applied for studying triatomic reactive scattering. This propagator is of short time type, but is exact to machine round off accuracy.

In order to improve the numerical efficiency of the SCN propagator, the spectral difference method was applied to construct of the Hamiltonian matrix of the kinetic operator. This resulted in the matrix with band structure thus saving the computational effort for evaluating the action of the Hamiltonian matrix on the wave function. At the same time, optimal preconditioner was applied for solving the linear equation involving in the SCN propagator to accelerate the numerical convergence.

Numerical illustration with the triatomic H+H$_2$ and O+O$_2$ reaction demonstrate that the SCN propagator can use huge time step to evolve the scattering wave function, without losing any accuracy. The preconditioner accelerate the convergence significantly, and the spectral difference method saves much numerical effort as comparing with the DVR method. Anyway, the current implementation of the SCN propagator is not competing with the current popular split operator and the CRWP propagator, which have been demonstrated being powerful for simulating quantum reactive scattering processes.

Ⅴ. ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China (No.21222308, No.21103187, and No.21133006), the Chinese Academy of Sciences, and the Key Research Program of the Chinese Academy of Sciences.

Reference
 [1] D. H. Zhang, and J. Z. H. Zhang, J. Chem. Phys. 99 , 5615 (1993). DOI:10.1063/1.465954 [2] D. H. Zhang, and J. Z. H. Zhang, J. Chem. Phys. 100 , 2697 (1994). DOI:10.1063/1.466464 [3] B. Fu, and D. H. Zhang, J. Phys. Chem. A 111 , 9516 (2007). DOI:10.1021/jp073811z [4] M. Yang, J. Chem. Phys. 129 , 064315 (2008). DOI:10.1063/1.2967854 [5] Z. Sun, S. Y. Lee, H. Guo, and D. H. Zhang, J. Chem. Phys. 130 , 174102 (2009). DOI:10.1063/1.3126363 [6] J. C. Light, I. P. Hamilton, and J. V. Lill, J. Chem. Phys. 82 , 1400 (1985). DOI:10.1063/1.448462 [7] R. Kosloff, J. Phys. Chem. 92 , 2087 (1988). DOI:10.1021/j100319a003 [8] J. C. Light, and T. Carrington, Adv. Chem. Phys. 114 , 263 (2007). [9] Z. Bacic, and J. C. Light, Annu. Rev. Phys. Chem 40 , 469 (1989). DOI:10.1146/annurev.pc.40.100189.002345 [10] B. I. Schneider, Phys. Rev. A 55 , 3417 (1997). DOI:10.1103/PhysRevA.55.3417 [11] S. Bauch, D. Hochstuhl, K. Balzer, and M. Bonitz, J. Phys.:Conf. Ser. 220 , 012020 (2010). DOI:10.1088/1742-6596/220/1/012020 [12] V. Szalay, J. Chem. Phys. 99 , 1978 (1993). DOI:10.1063/1.465258 [13] K. M. Dunseath, J. M. Launay, M. Terao-Dunseath, and L. Mouret, J. Phys. B:At. Mol. Opt. Phys. 35 , 3539 (2002). DOI:10.1088/0953-4075/35/16/313 [14] X. G. Wang, T. C. Carrington Jr., J. Chem. Phys. 129 , 234102 (2008). [15] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Lecture Notes in Engineering, Vol. 49. New York: Springer (1989). [16] C. Canuto, M. Y. Hussaini, and A. Quarteroni, Spectral Methods in Fluid Dynamics. New York: Springer (1988). [17] D. Gottlieb, and S. A. Orszag, Numerical Analysis of Spectral Methods:Theory and Applications, Philadelphia. PA: SIAM (1997). [18] P. F. Batcho, Phys. Rev. A 57 , 4246 (1998). DOI:10.1103/PhysRevA.57.4246 [19] A. T. Patera, J. Comput. Phys 54 , 468 (1984). DOI:10.1016/0021-9991(84)90128-1 [20] Y. Maday, and A. T. Patera, State of the Art Surveys in Computational Mechanics, A. Noor Ed.. New York: American Society for Mechanical Engineers (1987). [21] J. P. Boyd, Appl. Numer. Math. 7 , 287 (1991). DOI:10.1016/0168-9274(91)90065-8 [22] J. P. Boyd, J. Comp. Phys. 103 , 243 (1992). DOI:10.1016/0021-9991(92)90399-J [23] J. P. Boyd, Comput. Methods Appl. Mech. Eng. 116 , 1 (1994). DOI:10.1016/S0045-7825(94)80003-0 [24] G. W. Wei, D. S. Zhang, D. J. Kouri, and D. K. Hoffmann, Phys. Rev. Lett. 79 , 775 (1997). DOI:10.1103/PhysRevLett.79.775 [25] D. K. Hoffman, G. W. Wei, D. S. Zhang, and D. J. Kouri, Phys. Rev. E 57 , 6152 (1998). [26] G. W. Wei, S. C. Althorpe, D. S. Zhang, D. J. Kouri, and D. K. Hoffman, Phys. Rev. A 57 , 3309 (1998). DOI:10.1103/PhysRevA.57.3309 [27] D. A. Mazziotti, Chem. Phys. Lett. 299 , 473 (1999). DOI:10.1016/S0009-2614(98)01324-4 [28] A. Vijay, and H. Metiu, J. Chem. Phys. 116 , 60 (2002). DOI:10.1063/1.1425824 [29] H. Tal-Ezer, and R. Kosloff, J. Chem. Phys. 81 , 3967 (1984). DOI:10.1063/1.448136 [30] Y. Huang, D. J. Kouri, and D. K. Hoffman, J. Chem. Phys. 101 , 10493 (1994). DOI:10.1063/1.468481 [31] D. Q. Xie, R. Q. Chen, and H. Guo, J. Chem. Phys. 112 , 5263 (2000). DOI:10.1063/1.481096 [32] A. Vijay, R. E. Wyatt, and G. D. Billing, J. Chem. Phys. 111 , 10794 (1999). DOI:10.1063/1.480483 [33] X. G. Hu, Phys. Rev. E 59 , 2471 (1999). DOI:10.1103/PhysRevE.59.2471 [34] R. Kosloff, and Annu, Rev. Phys. Chem. 45 , 145 (1994). DOI:10.1146/annurev.pc.45.100194.001045 [35] M. D. Feit, J. A. Fleck Jr, and A. Steiger, J. Comput. Phys. 47 , 412 (1982). DOI:10.1016/0021-9991(82)90091-2 [36] J. A. Fleck Jr, J. R. Morris, and M. D. Feit, Ann. Phys. 10 , 129 (1976). [37] Z. G. Sun, D. H. Zhang, and W. T. Yang, Phys. Chem. Chem. Phys. 14 , 1817 (2012). [38] W. T. Li, D. H. Zhang, and Z. G. Sun, J. Phys. Chem. A 14 , 9801 (2014). [39] T. S. Chu, K. L. Han, and A. J. C. Varanda, J. Phys. Chem. A 110 , 1666 (2006). DOI:10.1021/jp054572n [40] K. L. Yeh, D. Xie, and D. H. Zhang, J. Phys. Chem. A 107 , 7215 (2003). DOI:10.1021/jp034471u [41] S. K. Gray, G. G. Balint-Kurti, J. Chem. Phys. 108 , 950 (1998). [42] G. J. Kroes, and D. Neuhauser, J. Chem. Phys. 105 , 8690 (1996). DOI:10.1063/1.472650 [43] S. K. Gray, G. G. Balint-Kurti, J. Chem. Phys. 108 , 3 (1998). [44] W. L. Quan, P. Y. Tang, and B. Y. Tang, Int. J. Quantum. Chem. 107 , 657 (2007). DOI:10.1002/(ISSN)1097-461X [45] P. Defazio, and C. Petrongolo, J. Phys. Chem. A 113 , 4208 (2009). DOI:10.1021/jp8106414 [46] C. L. Leforestier, R. H. Bisseling, and R. Kosloff, J. Comput. Phys. 94 , 59 (1991). DOI:10.1016/0021-9991(91)90137-A [47] T. J. Park, and J. C. Light, J. Chem. Phys. 85 , 5870 (1986). DOI:10.1063/1.451548 [48] H. Zhang, and S. C. Smith, J. Chem. Phys. 117 , 5714 (2002). [49] N. Markovic, and G. D. Billing, Mol. Phys. 98 , 1771 (2000). [50] B. Poirier, J. Chem. Phys. 108 , 5216 (1998). DOI:10.1063/1.475958 [51] Z. Sun, S. Y. L. X. Xu, and D. H. Zhang, J. Chem. Phys. 108 , 950 (2009). [52] J. Crank, and P. Nicolson, Proc. Cambridge Philos. Soc 43 , 50 (1947). DOI:10.1017/S0305004100023197 [53] R. Q. Chen, and H. Guo, Chem. Phys. Lett. 252 , 201 (1996). DOI:10.1016/0009-2614(96)00147-9 [54] R. Chen, and H. Guo, Comput. Phys. Commun. 119 , 19 (1999). DOI:10.1016/S0010-4655(98)00179-9 [55] Z. Sun, and W. Yang, J. Chem. Phys. 134 , 041101 (2011). DOI:10.1063/1.3549570 [56] T. Cvitas, and S. C. Althorpe, Comput. Phys. Commun. 177 , 357 (2007). DOI:10.1016/j.cpc.2007.05.002 [57] G. Sun, and C. W. Trueman, Electron. Lett. 39 , 595 (2003). DOI:10.1049/el:20030416 [58] B. Fu, and D. H. Zhang, J. Phys. Chem. A 111 , 9516 (2007). DOI:10.1021/jp073811z [59] D. M. Brink, and G. R. Satchler, Angular Momentum(Clarendon, 2nd Edn. England: Oxford (1968). [60] J. Z. H. Zhang, J. Chem. Phys. 94 , 6047 (1991). DOI:10.1063/1.460442 [61] S. Y. Lin, and H. Guo, Phys. Rev. A 74 , 022703 (2006). DOI:10.1103/PhysRevA.74.022703 [62] G. C. Schatz, and A. Kuppermann, J. Chem. Phys. 65 , 4642 (1976). DOI:10.1063/1.432918 [63] R. N. Zare, Angular Momentum. New York: Wiley (1988). [64] J. Z. H. Zhang, and W. H. Miller, J. Chem. Phys. 91 , 1528 (1989). DOI:10.1063/1.457650 [65] J. Lund, and K. L. Bowers, Sinc Methods for Quadrature and Differential Equations. Philadelphia, PA: SIAM (1992). [66] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions. New York: Springer (1993). [67] J. Dai, and J. Z. H. Zhang, J. Chem. Phys. 100 , 6898 (1996). DOI:10.1021/jp9536662 [68] R. Fletcher, Lecture Notes in Mathematics, 506 , 73 (1976). DOI:10.1007/BFb0080109 [69] P. Sonneveld, SIAM J. Sci. Comput 10 , 36 (1989). DOI:10.1137/0910004 [70] R. W. Freund, and N. Nachtigal, Numer. Math 60 , 315 (1991). DOI:10.1007/BF01385726 [71] R. W. Freund, SIAM J. Sci. Comput 14 , 470 (1993). DOI:10.1137/0914029 [72] U. Peskin, W. H. Miller, and A. Edlund, J. Chem. Phys. 103 , 10030 (1995). DOI:10.1063/1.469906 [73] W. H. Miller, Annu. Rev. Phys. Chem. 41 , 245 (1990). DOI:10.1146/annurev.pc.41.100190.001333 [74] F. J. Aoiz, L. Banares, and V. J. Herrero, Int. Rev. Phys. Chem 24 , 119 (2005). DOI:10.1080/01442350500195659 [75] A. I. Boothroyd, W. J. Keogh, P. G. Martin, and M. R. Peterson, J. Chem. Phys 104 , 7139 (1996). DOI:10.1063/1.471430 [76] R. P. Wayne, Chemistry in Atomospheres. Oxford: Oxford University Press (1991). [77] A. L. V. V. Wyngarden, K. A. Mar, K. A. Boering, J. J. Lin, Y. T. Lee, S. Y. Lin, H. Guo, and G. Lendvay, J. Am. Chem. Soc. 129 , 2866 (2007). DOI:10.1021/ja0668163 [78] Z. Sun, L. Liu, S. Y. Lin, Reinhard, H. Guo, and D. H. Zhang, PNAS 107 , 555 (2010). DOI:10.1073/pnas.0911356107

a. 中国科学院大连化学物理研究所, 分子反应动力学国家重点实验室和理论与计算化学中心, 大连 116023;
b. 中国科学院大学, 北京 100049