The article information
 Xueming Li, Zhigang Sun
 李学明, 孙志刚
 An Exact Propagator for Solving the Triatomic Reactive Schrödinger Equation^{†}
 一种精确求解三原子反应散射的传播子
 Chinese Journal of Chemical Physics, 2017, 30(6): 761770
 化学物理学报, 2017, 30(6): 761770
 http://dx.doi.org/10.1063/16740068/30/cjcp1711220

Article history
 Received on: November 15, 2017
 Accepted on: December 25, 2017
b. University of Chinese Academy of Sciences, Beijing 100049, China
The application of timedependent quantum wave packet (TDWP) method for studying reactive scattering processes becomes a more and more popular computational tool [15]. 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 [813]. 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 [1820] and the spectral difference [2127] are put forward to improve the numerical efficiency. The spectral difference methods can be either realised by using the distributed approximating functional (DAF) approach [2426] 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
Besides the above mentioned methods, there are some other popular propagators. In 1998, Gray and BalintKurti 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
Based on the concept of spectrally transformed Hamiltonian, ChenGuo introduced a new propagator in a form of the CrankNicholson 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 CrankNicolson method (named as SCN propagator hereafter) for solving the onedimension 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 timestep
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
For the triatomic reactive scattering in the reactant coordinates in the bodyfixed (BF) frame, the Hamiltonian for a given total angular momentum
$ \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
The timedependent 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) = \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}_{MK_\alpha}(\vec{\Omega}_\alpha)\right]\quad\quad $  (3) 
where
$ \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
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 spacefixed (SF) representation (
$ \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
$ \begin{eqnarray} G(R_\alpha)=\left(\frac{2}{\pi {\sigma}^2}\right)^{1/4}\exp\left[\frac{(R_\alphaR_{\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
$ \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
$ \begin{eqnarray} C_{lK}^{Jj\epsilon}=\sqrt{\frac{2l+1}{2J+1}}\sqrt{(2\delta_{K, 0})} \langle jKl_0JK\rangle \end{eqnarray} $  (8) 
where
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
$ \begin{eqnarray} \omega_{n}(x, \sigma) = \rm{e}^{(xx_{n})^2/2\sigma^2} \end{eqnarray} $  (10) 
and
$ \begin{eqnarray} C_{n}(x, \Delta x)={\rm sinc}[(xx_{n})/\Delta x] \end{eqnarray} $  (11) 
By decaying to zero, the weight Gaussian function
$ \begin{eqnarray} T_{m, n} &=& \Delta''_{mn}\omega_{n}(x, \sigma)+2\Delta'_{mn}\omega'_{n}(x, \Delta x)+\nonumber\\ &&{}\delta_{m, n}\omega''_{n}(x, \sigma) \end{eqnarray} $  (12) 
where
$ \Delta _{m  n}' = \left\{ {\begin{array}{*{20}{c}} {0,\;\;\;\;\;\;\;\;\;\;\;\;\; m  n = 0}\\ {\displaystyle\frac{{{{(  1)}^{m  n + 1}}}}{{\Delta xm  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
The evolution of the wave packet under the timeindependent 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
$ \begin{eqnarray} \exp(iH\Delta t)\approx\frac{1iH\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 timestep
Given a spectrally transformed Hamiltonian
$ \begin{eqnarray} H'=2{\rm arctan}(\Delta{t}H/2)/\Delta{t} \end{eqnarray} $  (17) 
The spectrally transformed Hamiltonian
$ \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 timedependent 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
$ \begin{eqnarray} \exp(iH'\Delta{t})\equiv\frac{1i\tan H'\Delta{t}/2}{1+i\tan H'\Delta{t}/2}\equiv\frac{1iH\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
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 timeindependent energyresolved
$ \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 timeindependent energyresolved
$ \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
$ \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
$ \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
Using the SCN propagator, the most important issue is that how to efficiently solve the equation
$ \begin{eqnarray} \Psi(r, t+\Delta{t})=\frac{1iH\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})=(1iH\Delta{t}/2)\Psi(t) \end{eqnarray} $  (31) 
which now is a linear equation of type
$ \begin{eqnarray} \hat A=1+i(\hat T_R+\hat T_r+\hat T_j+V)\Delta{t}/2 \end{eqnarray} $  (32) 
where
The iterative solvers for a linear equation can be accelerated by applying a suitable preconditioner, which means that the solution of the equation
$ \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
$ \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
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
$ \begin{eqnarray} {\rm diag}(\hat B)=I+{\rm diag}(\hat A_0^{1})(\hat A\hat A_0) \end{eqnarray} $  (36) 
where
$ \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 diagonalpreconditioner 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 potentialFinally, we use a damping function in the following form to avoid the wave function reaching the grid boundary along
$ \begin{eqnarray} D(r) = \left\{ \begin{array}{ll} \exp\left[\Delta_{t}C_{a}\left(\displaystyle\frac{rr_a}{r_br_a}\right)\right], r_a\leq r \leq r_b \\ \exp\left[\Delta_{t}C_{b}\left(\displaystyle\frac{rr_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
The H+H
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
FIG. 2 shows the error of the reaction probabilities as a function of the Gaussian parameter
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
Therefore, for the H+H
The O+O
The comparison of the convergence behaviors of the total reaction probabilities as a function of collision energy for the O+O
The left panel of FIG. 6 shows the error of the reaction probabilities as a function of the Gaussian parameter
The right panel of FIG. 6 shows the error of the reaction probabilities as a function of the Gaussian parameter
Since there are deep wells for the O+O
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 timestep 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
FIG. 8 presents total reaction probabilities for the O+O
The reaction O+O
Finally, we should note that the above effective timestep 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 multidimensional 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 bandstructure 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 CONCLUSIONIn this work, the SCN method, which has the same form as the classic CrankNicholson 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
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.
[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/17426596/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. TeraoDunseath, and L. Mouret, J. Phys. B:At. Mol. Opt. Phys. 35 , 3539 (2002). DOI:10.1088/09534075/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/00219991(84)901281 
[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/01689274(91)900658 
[22]  J. P. Boyd, J. Comp. Phys. 103 , 243 (1992). DOI:10.1016/00219991(92)90399J 
[23]  J. P. Boyd, Comput. Methods Appl. Mech. Eng. 116 , 1 (1994). DOI:10.1016/S00457825(94)800030 
[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/S00092614(98)013244 
[28]  A. Vijay, and H. Metiu, J. Chem. Phys. 116 , 60 (2002). DOI:10.1063/1.1425824 
[29]  H. TalEzer, 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/00219991(82)900912 
[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. BalintKurti, 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. BalintKurti, 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)1097461X 
[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/00219991(91)90137A 
[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/00092614(96)001479 
[54]  R. Chen, and H. Guo, Comput. Phys. Commun. 119 , 19 (1999). DOI:10.1016/S00104655(98)001799 
[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 
b. 中国科学院大学, 北京 100049