The article information
 DequanYu, HeHe, GunnarNyman, ZhigangSun
 Hermiticity of Hamiltonian Matrix using the Fourier Basis Sets in BondBondAngle and Radau Coordinates
 Chinese Journal of Chemical Physics, 2016, 29(1): 112122
 http://dx.doi.org/10.1063/16740068/29/cjcp1507141

Article history
 Received on July 4, 2015
 Accepted on July 13, 2015
b. School of Physics and Electronic Technology, Liaoning Normal University, Dalian 116029, China;
c. Department of Chemistry, Physical Chemistry, Göteborg University, SE412 96 Göteborg, Sweden;
d. Center for Advanced Chemical Physics and 2011 Frontier Center for Quantum Science and Technology, University of Science and Technology of China, Hefei 230026, China
Quantum dynamics calculations with molecular Hamiltonians have witnessed enormous progress. In the past years,a number of numerical techniques have been proposed to reduce the integration time for constructing the Hamiltonian matrix for a chosen basis set and to improve the accuracy of the numerical results ^{[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]}. This makes it possible to accurately simulate complex molecular dynamics and to calculate rovibrational spectra ^{[13, 14, 15, 16, 17, 18, 19]}. However,exact quantum calculations beyond triatomic and tetraatomic molecules are still a formidable task. Improving the computational methods will be of continued interest in the theoretical field.
Quantum dynamics calculations solving the timedependent molecular Schrödinger equation have proven to be invaluable for understanding dynamical process like photodissociation and femtosecond realtime experiments ^{[20, 21, 22]}. In an accurate triatomic molecular calculation,the Hamiltonian used is often described in Jacobi ^{[23]},Radau ^{[24, 25]},hyperspherical ^{[26, 27]},or bondbond angle coordinates ^{[25, 28]}. Fourier series or particleinthebox eigenfunctions are commonly used to expand the nuclear wavefunction along the radial degrees of freedom. As to the angular degree(s) of freedom,usually spherical harmonics (Legendre polynomials) basis functions are used. The numerical results of such calculations have been shown to be of exponential convergence ^{[9, 29]}.
The standard fastFourier transform (FFT) technique can be directly applied to propagate the initial wavepacket when the Fourier series are used as basis functions ^{[29]}. The introduction of the FFT ^{[30, 31]} technique led a surge of timedependent wavepacket applications. The FFT technique is easy to implement,even for complex kinetical operators ^{[28]}. At the same time,the obtained physical concept is clear,similar to that in a timedependent classic dynamics.
Several timedependent wavepacket calculations have illustrated that the Fourier functions can also be used as the basis set for the angular variable(s) ^{[28, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41]}. However,using the Fourier basis functions for the angular variable(s),singularities may be encountered. Several groups have investigated this problem and proposed to use a transformed Hamiltonian or to use cos θ instead of θ as the angular variable to avoid the difficulty with singularities ^{[28, 32, 33, 34, 35, 36, 37, 39]}. Moreover,to propagate the wavepacket with the convenient second order split operator method,which is combined with the FFT technique,a transformation of the triatomic Hamiltonian may be required ^{[24, 26, 42]}.
In this work,we point out that once a set of orthogonal basis functions have been chosen,for instance the Fourier basis set,there are limitations on how the Hamiltonians may be transformed. Otherwise the Hamiltonian matrix may become nonhermitian,leading to numerical problems. Tuvi et al. ^{[43]} have investigated the nonhermiticity problem using the Fourier grid Hamiltonian (FGH) discrete variable representation (DVR) method ^{[39]}. They suggested that to improve the numerical convergence,the Hamiltonian operators should be in explicitly symmetric form. If the Fourier basis set is chosen,the Hamiltonian should be kept in an explicitly symmetric form to avoid nonhermiticity resulting from integration errors.
In the following we will explicitly consider timedependent wavepacket calculations using the Fourier basis sets and various triatomic Hamiltonians in bond bondangle and Radau coordinates. The A^{2}A_{1}←X^{2}B_{1} absorption spectrum of the OClO molecule is calculated to illustrate the effect of transforming the Hamiltonian or changing the form of the Fourier basis set to the numerical results. Two types of nonhermiticity will be considered. The first type we will refer to as analytical nonhermitcity. It results from the use of a Hamiltonian and basis set in which the Hamiltonian matrix is nonhermitian even if it can be analytically constructed. The second type we will refer to as numerical nonhermiticity where the nonhermiticity in the Hamiltonian matrix results purely from integration errors.
In this work,we discuss normalization of the wavefunction,review the principles of the FFT technique and discuss how nonhermiticity arises in the Hamiltonian matrix in a Fourier basis set representation. Hamiltonians for triatomic molecules in bond bond angle and Radau coordinates are used for illustration. Numerical illustrations of the issues discussed are presented. This is done by calculating the absorption spectrum of the OClO molecule with a 3D (J=0) timedependent wavepacket model in Eckart and Radau coordinates.
Ⅱ. KINETIC ENERGY OPERATORS,THEIR EVALUATION AND HERMITICITY A. Kinetic energy operators and wavefunction normalizationThe Hamiltonian of the system can be expressed as a sum of the kinetic energy operator $(\hat{T})$ and the potential energy surface (V) in the nonrelativistic limit and within the BornOppenheimer (BO) approximation,
$ \hat{H}=\hat{T}+V $  (1) 
Starting with the Lagrangian form of the kinetic energy in terms of classic velocities,successive transformations lead to
$ T=\displaystyle\frac{1}{2}\mathbf{p}^T\mathbf{G}(\mathbf{q})\mathbf{p} $  (2) 
where p_{i} is the momenta conjugate to the chosen coordinates q_{i} (linearly independent),and $ \mathbf{G}(\mathbf{q}) $ is the wellknown G matrix
$ {G_{jk}} = \sum\limits_i {\frac{1}{{{m_i}}}\frac{{\partial {q_j}}}{{\partial {x_i}}}\frac{{\partial {q_k}}}{{\partial {x_i}}}} $  (3) 
where x_{i} denotes the Cartesian coordinates of atom i with a mass of m_{i}. Then,the Podolsky formalism ^{[44]} yields the quantum mechanical operator
$ \hat T =  \frac{{{\hbar ^2}}}{2}{g^{  1/2}}\sum\limits_{ij} {\frac{\partial }{{\partial {q_i}}}{g^{1/2}}{G_{ij}}\frac{\partial }{{\partial {q_j}}}} $  (4) 
where g=g=G^{1} is the determinant of the metric tensor matrix g ^{[45]}. The operator in Eq.(4) gives the normalization of the wavefunction $\Psi(\mathbf{q})$ as
$ \int \cdots \int {{g^{1/2}}} \prod\limits_i {\rm{d}} {q_i}{\Psi ^*}\Psi =  \int \cdots \int {\rm{d}} \tau {\Psi ^*}\Psi $  (5) 
$ \begin{array}{*{20}{c}} { =  \int \cdots \int {\rm{d}} {x_1}{\rm{d}}{y_1}{\rm{d}}{z_1} \cdots {\Psi ^*}\Psi }\\ { = 1} \end{array} $  (6) 
However this standard normalization of $\Psi$ with the weight factor g^{1/2} (normalized with volume element dτ=dx_{1}dy_{1}dz_{1} in cartesian coordinates ^{[46]}),is not always used. Podolsky ^{[44]} considered a transformed wavefunction normalized with the weight factor equal to 1,that is
$ \int\cdots\int \prod_i\textrm{d}q_i\Psi^*\Psi=1 $  (7) 
with the kinetical operators being
$ \hat T =  \frac{{{\hbar ^2}}}{2}{g^{  1/4}}\sum\limits_{ij} {\frac{\partial }{{\partial {q_i}}}{g^{1/2}}{G_{ij}}\frac{\partial }{{\partial {q_j}}}{g^{  1/4}}} $  (8) 
For triatomic systems,Hamiltonians in for instance bondbond angle,Jacobi and Radau coordinates,with different weight factors have been used. A systematic approach to this has been given by Makarewicz ^{[16]}.
Let us assume that the determinant g(q) can be represented as a product of two independent factors $g({\bf{q}}) = {g_A}({{\bf{q}}_A}){g_B}({{\bf{q}}_B})$ where ${\bf{q}} = ({{\bf{q}}_A},{{\bf{q}}_B})$ and let us transform the wavefunction $\Psi(\mathbf{q})$ to
$ \Phi(\mathbf{q})=g_A^{1/4}\Psi(\mathbf{q}) $  (9) 
which,following from Eq.(5),fulfills the normalization condition:
$ \label{normcc} \int\cdots\int g_B^{1/2}\prod_i\textrm{d}q_i\Phi^*\Phi=1 $  (10) 
For g_{A}=1 and g_{B}=g this equation gives the standard normalization defined by Eq.(5),and for g_{A}=g and g_{B}=1 it gives the normalization with weight factor unity,which Podolsky considered ^{[44]}. Because the factor g_{A}^{1/4} is associated with the transformation of the wavefunction,we will refer to g_{A}^{1/4} as the transformation factor in order to distinguish from the weight factor g_{B}^{1/2}. We note here that,if the weight factor g_{B}^{1/4} is equal to 1 (that is,with transformation factor as g),the operators can always be rearranged in an explicitly symmetric form ^{[8, 44]},as shown in Eq.(8).
Let us categorise triatomic Hamiltonians in bondbond angle coordinates that have appeared in the literature according to the associated transformation (g_{A}^{1/4}) and weight (g_{B}^{1/2}) factors ^{[16]}. The Hamiltonian $\widehat {{H_w}}$ considered by for instance Bardo and Wolfberg ^{[47, 48]} falls in the category of Podolsky's normalization (weight factor g_{B}^{1/4} equal to 1) ^{[16]}. The Hamiltonian can be written ^{[47, 48]},
$ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{{\hat H}_w} = {\rm{ }}  \frac{{{\hbar ^2}}}{2}[\frac{1}{{{\mu _1}}}\frac{{{\partial ^2}}}{{\partial r_1^2}} + \frac{1}{{{\mu _2}}}\frac{{{\partial ^2}}}{{\partial r_2^2}} + (\frac{1}{{2{\mu _1}r_1^2}} + \frac{1}{{2{\mu _2}r_2^2}}  {\rm{ }}}\\ {\frac{{\cos \theta }}{{{\mu _3}{r_1}{r_2}}})\frac{{{\partial ^2}}}{{\partial {\theta ^2}}}]  \frac{{{\hbar ^2}}}{{{\mu _3}}}[\frac{{\sin \theta }}{{{r_1}{r_2}}}\frac{\partial }{{\partial \theta }}  \frac{{\cos \theta }}{2} \cdot } \end{array}}\\ {(\frac{1}{{{r_1}}}\frac{\partial }{{\partial {r_2}}} + \frac{1}{{{r_2}}}\frac{\partial }{{\partial {r_1}}}) + \cos \theta \frac{{{\partial ^2}}}{{\partial {r_1}\partial {r_2}}}  } \end{array}}\\ {\sin \theta \frac{\partial }{{\partial \theta }}(\frac{1}{{{r_1}}}\frac{\partial }{{\partial {r_2}}} + \frac{1}{{{r_2}}}\frac{\partial }{{\partial {r_1}}})] + }\\ {\hat V({r_1},{r_2},\theta ) + \Delta \hat V} \end{array} $  (11) 
and
$ \begin{array}{*{20}{c}} {\Delta \hat V = \frac{{{{\cos }^3}\theta }}{{4{\mu _3}{r_2}{r_2}{{\sin }^2}\theta }}  }\\ {\frac{1}{8}\left( {\frac{1}{{{\mu _1}r_1^2}} + \frac{1}{{{\mu _2}r_2^2}}} \right)\left( {1 + {{\csc }^2}\theta } \right)} \end{array} $  (12) 
where μ_{1},μ_{2},and μ_{3} are the related (reduced) masses. The other variables have their traditional definitions. As we have noted,since this Hamiltonian $\widehat {{H_w}}$ has a corresponding weight factor of unity and transformation factor g,it can be rearranged into an explicitly symmetric form
$ \begin{array}{*{20}{c}} {{{\hat H}_{ws}}{\rm{ }} = {\rm{ }}\frac{{\hat P_1^2}}{{2{\mu _1}}} + \frac{{{{\hat j}^2}}}{{2{\mu _1}r_1^2}} + \frac{{\hat P_2^2}}{{2{\mu _2}}} + \frac{{{{\hat j}^2}}}{{2{\mu _2}r_2^2}} + }\\ {\frac{{{{\hat P}_1}{{\hat P}_2}\cos \theta }}{{{\mu _3}}}  \frac{{{{\hat P}_1}{{\hat P}_\theta }}}{{{\mu _3}{r_2}}}  \frac{{{{\hat P}_2}{{\hat P}_\theta }}}{{{\mu _3}{r_1}}}  }\\ {\frac{{\cos \theta {{\hat j}^2} + {{\hat j}^2}\cos \theta }}{{2{\mu _3}{r_1}{r_2}}} + \hat V\left( {{r_1},{r_2},\theta } \right)} \end{array} $  (13) 
$ {\hat j^2} = {\rm{ }}  {\hbar ^2}{\left( {\frac{1}{{\sin \theta }}} \right)^{  1/2}}\frac{\partial }{{\partial \theta }}\sin \theta \frac{\partial }{{\partial \theta }}{\left( {\frac{1}{{\sin \theta }}} \right)^{  1/2}} $  (14) 
$ {\hat P_i} =  i\hbar \frac{\partial }{{\partial {r_i}}},\qquad i = 1,2 $  (15) 
$ {\hat P_\theta } = {\rm{ }}  i\hbar {\left( {\sin \theta } \right)^{1/2}}\frac{\partial }{{\partial \theta }}{\left( {\sin \theta } \right)^{1/2}} $  (16) 
This Hamiltonian gives the wavefunction normalization condition^{[16, 47, 48]}
$ \int \textrm{d} r_1 \int \textrm{d} r_2 \int \textrm{d}\theta \Psi^*(r_1,r_2,\theta)\Psi(r_1,r_2,\theta)=1 $  (17) 
The Hamiltonians $\hat{H}_w$ and $\hat{H}_{ws}$ are not generally used. However,the Hamiltonian $\hat{H}_c$ in bond bond angle coordinates given by Carter and Handy ^{[49]} is more often used,
$ \begin{array}{*{20}{c}} {{{\hat H}_c} =  \frac{{{\hbar ^2}}}{2}[\frac{1}{{{\mu _1}}}\frac{{{\partial ^2}}}{{\partial r_1^2}} + \frac{1}{{{\mu _2}}}\frac{{{\partial ^2}}}{{\partial r_2^2}} + (\frac{1}{{{\mu _1}r_1^2}} + \frac{1}{{{\mu _2}r_2^2}}  }\\ {\frac{{2\cos \theta }}{{{\mu _3}{r_1}{r_2}}})\frac{1}{{\sin \theta }}\frac{\partial }{{\partial \theta }}\sin \theta \frac{\partial }{{\partial \theta }}] + \frac{{{\hbar ^2}}}{{{\mu _3}}}[  \cos \theta \cdot }\\ {\frac{{{\partial ^2}}}{{\partial {r_1}\partial {r_2}}} + (\frac{1}{{{r_1}}}\frac{\partial }{{\partial {r_2}}} + \frac{1}{{{r_2}}}\frac{\partial }{{\partial {r_1}}}  \frac{1}{{{r_1}{r_2}}}) \cdot }\\ {(\sin \theta \frac{\partial }{{\partial \theta }} + \cos \theta )] + \hat V({r_1},{r_2},\theta )} \end{array} $  (18) 
A rearranged form $\hat{H}_{ca}$ of this Hamiltonian is
$ \begin{array}{*{20}{c}} {{{\hat H}_{ca}} = \frac{{\hat P_1^2}}{{2{\mu _1}}} + \frac{{{{\hat j}^2}}}{{2{\mu _1}r_1^2}} + \frac{{\hat P_2^2}}{{2{\mu _2}}} + \frac{{{{\hat j}^2}}}{{2{\mu _2}r_2^2}} + }\\ {\frac{{{{\hat P}_1}{{\hat P}_2}\cos \theta }}{{{\mu _3}}}  \frac{{{{\hat P}_1}{{\hat P}_\theta }}}{{{\mu _3}{r_2}}}  \frac{{{{\hat P}_2}{{\hat P}_\theta }}}{{{\mu _3}{r_1}}}  }\\ {\frac{{\cos \theta {{\hat j}^2} + {{\hat j}^2}\cos \theta }}{{2{\mu _3}{r_1}{r_2}}} + \hat V\left( {{r_1},{r_2},\theta } \right)} \end{array} $  (19) 
$ {\hat j^2}{\rm{ =  }}{\hbar ^{\rm{2}}}\frac{{\rm{1}}}{{\sin \theta }}\frac{\partial }{{\partial \theta }}\sin \theta \frac{\partial }{{\partial \theta }} $  (20) 
$ {{\rm{\hat P}}_{{\rm{1}},{\rm{2}}}}{\rm{ =  i}}\hbar \frac{\partial }{{\partial {{\rm{r}}_{{\rm{1}},{\rm{2}}}}}} $  (21) 
$ {{\rm{\hat P}}_\theta }{\rm{ =  i}}\hbar \frac{\partial }{{\partial \theta }}\sin \theta $  (22) 
The latter Hamiltonian gives the wavefunction normalization ^{[49]},
$ \label{normc2} \int \textrm{d} r_1 \int \textrm{d} r_2 \int \sin\theta \textrm{d}\theta \Psi^*(r_1,r_2,\theta)\Psi(r_1,r_2,\theta)=1\quad $  (23) 
We note that for the triatomic Hamiltonian $\hat{H}_{ca}$ ($\hat{H}_{c}$) in bondbond angle coordinates,the weight factor g_{B}^{1/2}=sinθ. Therefore,the operators of $\hat{H}_c$ cannot be written into an explicitly symmetric form.
B. Action of the kinetic energy operator using Fourier transformIn a timedependent quantum calculation,the action of the potential operators on the wavefunction is simple multiplication when a spatial representation is used. The evaluation of the action of the kinetic energy operator will be in this representation however involve derivatives of the wavefunction. The derivatives can be efficiently found by the Fourier transform technique,which transforms to momentum space where the kinetic energy operator is local and only multiplication with the wavefunction is required.
The forward Fourier transform $\widehat{FT}(x,k)$ is defined as
$ \Phi(k)=\widehat{FT}(x,k)\Psi(x)=\displaystyle\frac{1}{\sqrt{2\pi}}\int{\Psi(x)\textrm{e}^{ikx}}\textrm{d}x $  (24) 
It transforms the wavefunction from the space representation to the momentum representation.Similarly to Eq.(24),the backward Fourier transform $\widehat{FT}(x,k)$^{1} is given by
$ \begin{array}{*{20}{c}} {\Psi (x) = \widehat {FT}{{(k,x)}^{  1}}\Phi (k)}\\ { = {\rm{ }}\frac{1}{{\sqrt {2\pi } }}\int {\Phi (k){{\rm{e}}^{ikx}}} {\rm{d}}k} \end{array} $  (25) 
It is easy to see that for the derivatives of the wavefunction we can write
$ \frac{{{\rm{d}}\Psi (x)}}{{{\rm{d}}x}} = \frac{1}{{\sqrt {2\pi } }}\int {\Phi (k)ik{{\rm{e}}^{ikx}}{\rm{d}}k} $  (26) 
$ \frac{{{{\rm{d}}^2}\Psi (x)}}{{{{\rm{d}}^2}x}} = \frac{1}{{\sqrt {2\pi } }}\int {\Phi (k){{(ik)}^2}{{\rm{e}}^{ikx}}{\rm{d}}k} $  (27) 
Thus the derivatives are evaluated as simple multiplications in the momentum space. Treating complicated kinetic operators with Fourier transform can be derived straightforwardly. For instance,the action of the operator $ \sqrt{1x^2}\displaystyle\frac{\textrm{d}}{\textrm{d}x}\sqrt{1x^2}$ (this is the angular momentum operator which is used in the later sections,see Eq.(42)) on the wavefunction can be written as
$ \begin{array}{*{20}{c}} {\sqrt {1  {x^2}} \frac{{\rm{d}}}{{{\rm{d}}x}}\sqrt {1  {x^2}} \Psi (x) = \frac{{\sqrt {1  {x^2}} }}{{\sqrt {2\pi } }}\int {\Theta (k)(ik){{\rm{e}}^{ikx}}{\rm{d}}k} }\\ {\Theta (k) = \widehat {FT}(k,x)\phi (x)}\\ { = \frac{1}{{\sqrt {2\pi } }}\int {\phi (x){{\rm{e}}^{ikx}}} {\rm{d}}x} \end{array} $  (28) 
$ \phi (x) = \sqrt {1  {x^2}} \Psi (x) $  (29) 
We note here that when model Hamiltonian is expressed in a suitable form,such as the Hamiltonian $\hat{H}_{ca}$ and $\hat{H}_{ws}$,the action of each kinetic operator on the wavefunction always consists of a combination of two first derivatives. As a result,the factor i (imaginary unit) resulting from a single evaluation of the first derivative using the Fourier transform technique does not appear in the numerical implementation.
C. Hermiticity of the angular momentum operator using the Fourier basis 1. HermiticityPropagating the wavepacket using Fourier transform in a timedependent calculation is principally equivalent to evolving the wavefunction using a Fourier series. In the numerical calculation,the spatial and momentum coordinates have to be discretized as $x_n\rangle$ and $k_m\rangle$ which enables the implementation of the FFT technique. The transformation between the discrete spatial coordinate representation $x_n\rangle$ and the discrete momentum representation $k_m\rangle$ is given by
$ {x_n}\rangle = \frac{1}{{\sqrt {2\pi N} }}\sum\limits_{m = 1}^N {{{\rm{e}}^{i{k_m}{x_n}}}} {k_m}\rangle $  (30) 
$ k_m\rangle=\displaystyle\frac{1}{\sqrt{2\pi N}}\sum_{n=1}^N \textrm{e}^{i k_mx_n}x_n\rangle $  (31) 
Fast Fourier transform algorithms can be employed to affect these transformations and to compute derivative(s) of the wavefunction for realizing the action of the operator on the wavefunction or the wavefunction multiplied by functions of the spatial coordinates. In order to have a good representation of the system,the Hamiltonian matrix should be hermitian in the discrete Fourier space. Tuvi et al. ^{[43]} have shown that writing operators of a onedimensional nuclear Hamiltonian in an explicitly symmetric form guarantees the hermiticity of the resulting Hamiltonian matrix. Their conclusion can be extended directly to multidimensional problems. Therefore the Hamiltonian operators given by the Podolsky normalization condition ^{[44]} are definitely hermitian in a Fourier basis set. We may use the corresponding Hamiltonian form,that is,use Hamiltonian with weight factor g_{B}^{1/4}=1 in our calculation using a Fourier basis set.
Often the weight factor g_{B}^{1/4} of the Hamiltonian being used is not equal to 1 whereby the operators in a Fourier basis set may not be hermitian. For example,the angular momentum operator is often written in a socalled manifestly hermitian form
$ \hat{j^2}=\hbar^2\displaystyle\frac{1}{\sin\theta}\displaystyle\frac{\partial}{\partial\theta}\sin\theta \displaystyle\frac{\partial}{\partial\theta} \label{angular} $  (32) 
but it is easy to see that it is not hermitian in the discrete Fourier basis,$\{ 1/\sqrt{2\pi N}e^{ik_n\theta},k_n=1,N \}$,where k_{n} is a grid point in momentum space. That is,
$ \begin{array}{*{20}{c}} {  \left\langle {{{\rm{e}}^{  i{k_m}\theta }}\left {\frac{1}{{\sin \theta }}\frac{\partial }{{\partial \theta }}\sin \theta \frac{\partial }{{\partial \theta }}} \right{{\rm{e}}^{  i{k_n}\theta }}} \right\rangle \ne \quad \quad }\\ {  {{\left\langle {{{\rm{e}}^{  i{k_n}\theta }}\left {\frac{1}{{\sin \theta }}\frac{\partial }{{\partial \theta }}\sin \theta \frac{\partial }{{\partial \theta }}} \right{{\rm{e}}^{  i{k_m}\theta }}} \right\rangle }^*}} \end{array} $  (33) 
unless n=m,where * indicates its conjugate. As discussed above,this is connected to the fact that the wavefunction $\Psi$ of the angular momentum operator has a weight factor g_{B}^{1/4}=1/sinθ^{1/2} which is not taken into consideration by the given Fourier basis set.
A transformation $\sqrt{\sin\theta}\hat{j}^2{\sin\theta}^{1/2}$ can be made to remove the weight factor sinθ^{1/2}. The transformed operator then becomes
$ (\hat{j^{\prime}})^2=\hbar^2\left(\displaystyle\frac{1}{\sin\theta}\right)^{1/2}\displaystyle\frac{\partial}{\partial\theta}\sin\theta \displaystyle\frac{\partial}{\partial\theta}\left(\displaystyle\frac{1}{\sin\theta}\right)^{1/2} $  (34) 
This means that we write the angular momentum operator using Podolsky's normalization condition ^{[44]} and the weight factor g_{B} becomes unity ^{[16]}. Hence the operator is definitely symmetric. For a triatomic case,this corresponds to transforming the Hamiltonian from $\hat{H}_{ca}$ in Eq.(19) to $\hat{H}_{ws}$ in Eq.(13). When the Fourier basis set $1/\sqrt{2\pi N}\textrm{e}^{ik_n\theta}$ is used in a calculation where the angular momentum operator in Eq.(34) is employed,problems with nonhermiticity do not arise.
The weight factor g_{B}^{1/4}=$1/\sqrt{\sin\theta}$ can also be absorbed by taking cos θ as the variable instead of θ. The operator in Eq.(32) then becomes,
$ \hat{G}=\displaystyle\frac{\partial}{\partial C}\left(1C^2\right) \displaystyle\frac{\partial}{\partial C} \label{px1} $  (35) 
where C=cos θ. The weight factor becomes unity and the normalization condition is
$ \begin{array}{*{20}{c}} {\int_{  1}^1 {{\psi ^*}} (C)\psi (C){\rm{d}}C = \int_{  1}^1 {{\psi ^*}} (\cos \theta )\psi (\cos \theta ){\rm{d}}(\cos \theta )}\\ { = 1} \end{array} $  (36) 
We can then use the discrete Fourier basis set $1/\sqrt{2\pi N}\textrm{e}^{ik_nC}$ in the calculation,whereby the resulting Hamiltonian matrix becomes hermitian. In this way the triatomic Hamiltonian can be obtained in an explicitly symmetric from. From the discussion above,it is clear that to decide whether an operator,particularly the Hamiltonian operator,is hermitian or not,the basis set has to first be decided.
2. Integration errorConsidering the integration error,the Hamiltonian operators used in the numerical calculation should be in an explicitly symmetrical form when the discrete Fourier basis is used ^{[43]}. Sometimes although a Hamiltonian operator is theoretically hermitian in some basis set,the resulting Hamiltonian matrix becomes nonhermitian in the numerical calculations. For instance,the operator in Eq.(35) and
$ \hat{g}=\left(1C^2\right) \displaystyle\frac{\partial^2}{\partial C^2}+2C\displaystyle\frac{\partial}{\partial C} \label{px2} $  (37) 
are equivalent and hermitian in a Fourier basis set when the integration is exact. However because of the integration errors resulting from the discretization,the operator in Eq.(37) is not numerically hermitian in the discrete Fourier basis set. This is because the corresponding Hamiltonian matrix ^{[43]},
$ \begin{array}{*{20}{c}} {{g_{ml}} = \langle {C_m}\hat g{C_l}\rangle }\\ {{\rm{ = }}  \left( {1  C_m^2} \right)\left\langle {{C_m}\left {\frac{{{\partial ^2}}}{{\partial {C^2}}}} \right{C_l}} \right\rangle + }\\ {2{C_m}\left\langle {{C_l}\left {\frac{\partial }{{\partial C}}} \right{C_l}} \right\rangle } \end{array} $  (38) 
and its hermitian conjugate,
$ \begin{array}{*{20}{c}} {g_{lm}^* = {{\langle {C_l}\hat g{C_m}\rangle }^*}}\\ {{\rm{ = }}  \left( {1  C_l^2} \right){{\left\langle {{C_l}\left {\frac{{{\partial ^2}}}{{\partial {C^2}}}} \right{C_m}} \right\rangle }^*} + }\\ {2{C_l}{{\left\langle {{C_l}\left {\frac{\partial }{{\partial C}}} \right{C_m}} \right\rangle }^*}} \end{array} $  (39) 
are only equivalent if there is an odd number of grid points N and m=l ^{[43]},where C_{i}=$i\Delta C$,i=1,$\ldots$,N. Thus ^{[43]},when the model Hamiltonian used in an FFT timedependent calculatition contain operators of the form in Eq.(37),the numerical results may be unreliable.
D. Hermiticity of the J=0 triatomic Hamiltonian in bondbond angle and Radau coordinates in the Fourier basisThe commonly used triatomic Hamiltonian (J=0) in bondbond angle coordinates $(r_{1},r_{2},\theta)$ written in the form of $\hat{H}_{ps}$ with weight factor g_{B} equal to 1 is explicitly symmetric. Therefore the resulting Hamiltonian matrix is hermitian ^{[43]} in the Fourier basis set $\{ 1/\sqrt{2\pi N}\textrm{e}^{i k_l r_{1}},1/\sqrt{2\pi N}\textrm{e}^{i k_m r_{2}},1/\sqrt{2\pi N}\textrm{e}^{i k_n \theta} \}$,below referred to as FBST. However,the rearranged form $\hat{H}_{ca}$ of the Hamiltonian in Eq.(18),with weight factor g_{B}^{1/4}= $1/\sqrt{\sin\theta}$,given by Carter and Handy ^{[49]} does not lead to a hermitian matrix using the FBST basis set and is not suitable for an FFT calculation. As we have discussed,the transformation factor can either be absorbed by a transformation of the Hamiltonian ${\sin\theta}^{1/2}\hat{H}{\sin\theta}^{1/2}$ or by changing the variable θ to cos θ. From the work of Carter and Handy ^{[49]},the Hamiltonian $\hat{H}_{cs}$ in bondbond angle coordinates,taking r_{1},r_{2} and C=cos θ as the variables,can be written as,
$ \begin{array}{*{20}{c}} {{{\hat H}_{cs}} = \frac{{\hat P_1^2}}{{2{\mu _1}}} + \frac{{{{\hat j}^2}}}{{2{\mu _1}r_1^2}} + \frac{{\hat P_2^2}}{{2{\mu _2}}} + \frac{{{{\hat j}^2}}}{{2{\mu _2}r_2^2}} + }\\ {\frac{{{{\hat P}_1}{{\hat P}_2}\cos \theta }}{{{\mu _3}}}  \frac{{{{\hat P}_1}{{\hat P}_\theta }}}{{{\mu _3}{r_2}}}  \frac{{{{\hat P}_2}{{\hat P}_\theta }}}{{{\mu _3}{r_1}}}  }\\ {\frac{{\cos \theta {{\hat j}^2} + {{\hat j}^2}\cos \theta }}{{2{\mu _3}{r_1}{r_2}}} + \hat V\left( {{r_1},{r_2},\theta } \right)} \end{array} $  (40) 
$ {\hat j^2} =  {\hbar ^2}\frac{\partial }{{\partial C}}\left( {1  {C^2}} \right)\frac{\partial }{{\partial C}} $  (41) 
$ {\hat P_\theta } = i\hbar {\left( {1  {C^2}} \right)^{1/2}}\frac{\partial }{{\partial C}}{\left( {1  {C^2}} \right)^{1/2}} $  (42) 
The operators of this Hamiltonian are in explicitly symmetric form and result in a hermitian matrix with the Fourier basis set (FBSC): $\{ 1/\sqrt {2\pi N} {{\rm{e}}^{  i{k_l}{r_1}}},1/\sqrt {2\pi N} {{\rm{e}}^{  i{k_m}{r_2}}},1/\sqrt {2\pi N} {{\rm{e}}^{  i{k_n}C}}\} $.
Similar conclusions hold for the Hamiltonian in Radau coordinates. The usually used Hamiltonian (J=0) in Radau coordinates ($R_1,R_2,\theta$),can be written as ^{[50]},
$ \begin{array}{*{20}{c}} {{{\hat H}_{R\theta a}} =  \frac{{{\hbar ^2}}}{{2{m_1}}}\frac{{{\partial ^2}}}{{\partial R_1^2}}  \frac{{{\hbar ^2}}}{{2{m_2}}}\frac{{{\partial ^2}}}{{\partial R_2^2}}  (\frac{{{\hbar ^2}}}{{2{m_1}R_1^2}} + \quad \quad }\\ {\frac{{{\hbar ^2}}}{{2{m_2}R_2^2}})\frac{1}{{\sin \theta }}\frac{\partial }{{\partial \theta }}\sin \theta \frac{\partial }{{\partial \theta }} + V({R_1},{R_2},\theta ,t)} \end{array} $  (43) 
where m_{1} and m_{2} are the relevant masses. Note that the variable θ here has a different meaning from that in bond bondangle coordinates. This Hamiltonian has a volume element ${\alpha ^{  3}}\sin \theta d{R_1}d{R_2}$dθ where α^{2}=m_{3}/(m_{1}+m_{2}+m_{3}),with weight factor g_{B}^{1/4} equal to sinθ^{1/2}. Its resulted Hamiltonian matrix is not hermitian in the FBST basis set due to the transformation factor sinθ^{1/2}. Similar to the case in bondbond angle coordinates,the Hamiltonian in Eq.(43) can be rewritten as
$ \begin{array}{*{20}{c}} {{{\hat H}_{Rcs}} =  \frac{{{\hbar ^2}}}{{2{m_1}}}\frac{{{\partial ^2}}}{{\partial R_1^2}}  \frac{{{\hbar ^2}}}{{2{m_2}}}\frac{{{\partial ^2}}}{{\partial R_2^2}}  (\frac{{{\hbar ^2}}}{{2{m_1}R_1^2}} + \quad \quad }\\ {\frac{{{\hbar ^2}}}{{2{m_2}R_2^2}})\frac{\partial }{{\partial C}}(1  {C^2})\frac{\partial }{{\partial C}} + V({R_1},{R_2},\theta ,t)} \end{array} $  (44) 
The difficulty arising from the transformation factor $1/\sqrt{\sin\theta}$ has been eliminated when taking $(R_1,R_2,\cos\theta)$ as the variables instead of $(R_1,R_2,\theta)$ and using the FBSC basis set. The Hamiltonian operator is now explicitly symmetric and the resulting Hamiltonian matrix is hermitian.
The Hamiltonian in Radau coordinates can also be transformed ^{[24, 42]}
$ \begin{array}{*{20}{c}} {{{\hat H}_{R\theta s}} =  \frac{{{\hbar ^2}}}{{2{m_1}}}\frac{{{\partial ^2}}}{{\partial R_1^2}}  \frac{{{\hbar ^2}}}{{2{m_2}}}\frac{{{\partial ^2}}}{{\partial R_2^2}}  (\frac{{{\hbar ^2}}}{{2{m_1}R_1^2}} + }\\ {\frac{{{\hbar ^2}}}{{2{m_2}R_2^2}})\left( {\frac{{{\partial ^2}}}{{\partial {\theta ^2}}} + \frac{1}{{4{{\sin }^2}\theta }} + \frac{1}{4}} \right) + }\\ {V({R_1},{R_2},\theta ,t)} \end{array} $  (45) 
with the weight factor being equal to 1. Therefore,the Hamiltonian $\hat{H}_{R\theta s}$ results in a hermitian matrix in the FBST basis set. The Hamiltonian has a volume element d${R_{\rm{1}}}{\rm{d}}{R_{\rm{2}}}$dθ.
We can expect that the convergence of the calculations are good using the Hamiltonian $\hat{H}_{ws}$ in Eq.(13) or $\hat{H}_{R\theta s}$ in Eq.(45) with the FBST basis set and using $\hat{H}_{cs}$ in Eq.(40) or $\hat{H}_{Rcs}$ Eq.(44) with the FBSC basis set. However,calculations using the Hamiltonian $\hat{H}_{ca}$ in Eq.(19) or $\hat{H}_{R\theta a}$ in Eq.(43) with either of the two Fourier basis sets should be only conditionally stable and exhibit numerical difficulties.
Sometimes the Hamiltonian $\hat{H}_{ws}$ in Eq.(13) is used in the expanded form $\hat{H}_{w}$ in Eq.(11),or the Hamiltonian $\hat{H}_{Rcs}$ in Eq.(45) is used in the form,
$ \begin{array}{*{20}{c}} {{{\hat H}_R} =  \frac{{{\hbar ^2}}}{{2{m_1}}}\frac{{{\partial ^2}}}{{\partial R_1^2}}  \frac{{{\hbar ^2}}}{{2{m_2}}}\frac{{{\partial ^2}}}{{\partial R_2^2}}  \left( {\frac{{{\hbar ^2}}}{{2{m_1}R_1^2}} + \frac{{{\hbar ^2}}}{{2{m_2}R_2^2}}} \right) \cdot }\\ {\left[ {(1  {C^2})\frac{{{\partial ^2}}}{{{\partial ^2}C}}  2C\frac{\partial }{{\partial C}}} \right] + V({R_1},{R_2},\varphi ,t)} \end{array} $  (46) 
As we discussed above,the Hamiltonian matrix constructed from the Hamiltonian $\hat{H}_{w}$ in Eq.(11) using the FBST basis set,or from the Hamiltonian $\hat{H}_{R}$ in Eq.(46) using the FBSC basis set,is nonhermitian,due to integration errors,despite the fact that it would be analytically hermitian. We will see later that this kind of nonhermiticity caused purely by integration errors is only conditionally stable and may lead to numerical problems for long time propagations.
Nonhermitian Hamiltonian matrices develop complex eigenvalues which make the propagation unstable. In some cases numerical errors may be seen. Above we have pointed out that to obtain a hermitian Hamiltonian matrix,the forms of the Hamiltonian operators must be carefully constructed and a suitable basis set should be chosen. Therefore,when a derived Hamiltonian is transformed,perhaps for numerical convenience,the basis set used must be reconstructed accordingly.
Wei et al. ^{[8]} in their work pointed out that when the product approximation technique is used to obtain the Hamiltonian matrix in a DVR calculation,the Hamiltonian matrix may exhibit nonhermiticity even though the operator has not been expanded. This results from integration error. Similar problem does not arise in a Fourier basis set calculation which makes it possible also easy to treat with complex kinetical operators with FFT technique,because FFT itself is rigorously antihermitician when the periodical bound condition are fulfilled.
In a timedependent wavepacket calculation,the influence of the nonhermiticity of the Hamiltonian matrix is accumulated which is sensitively indicated by the timedependence of the norm of the wavefunction. In a variational calculation,however,the nonhermiticity of the Hamiltonian matrix may have less impact on the numerical results. When a quantum calculation is done with a nonhermitian Hamiltonian matrix,a series of "ghost" (unphysical) peaks may arise. In some cases,the nonhermitian matrix may change the heights of the physical peaks. The positions of the physical peaks however normally remain correct. The nonhermitian matrix may lead to complex eigenvalues which leads to norm violation. In the next section these issues will be numerically investigated by calculating absorption spectra with a 3D timedependent wavepacket model.
Ⅲ. NUMERICAL EXAMPLESIn this section,we will show absorption spectra of the OClO molecule calculated with the timedependent wavepacket method using Hamiltonians in bondbond angle and Radau coordinates to illustrate the arguments put forward above. Results obtained using the Hamiltonian $\hat{H}_{cs}$ in Eq.(40) and $\hat{H}_{Rcs}$ in Eq.(44) using the FBSC basis set are compared with results calculated using the Hamiltonian $\hat{H}_{ws}$ in Eq.(13) and $\hat{H}_{R\theta s}$ in Eq.(45) using the FBST basis set. These comparisons are used to check the convergence of the calculation using the FFT technique. From the discussion above,we know that these four numerical models should exhibit the best convergence.
The results using the Hamiltonian $\hat{H}_{ca}$ in Eq.(19) and $\hat{H}_{ws}$ in Eq.(13) with the FBST basis set are compared in order to investigate the nonhermiticity problem caused by the weight factor. The calculations using the Hamiltonian $\hat{H}_{w}$ in Eq.(11) with the FBST basis set and $\hat{H}_{R}$ in Eq.(46) with the FBSC basis set are carried out and compared with the results from the best converged models,in order to investigate the nonhermiticity problem caused purely by integration errors when using $\hat{H}_{w}$ and $\hat{H}_{R}$. We find that for the present application even for a short time propagation a nonhermitian Hamiltonian matrix may make the results unreliable.
OClO is a molecule of both experimental and theoretical interest,due to its presumed role in polar stratospheric ozone depletion. Accurate 3D ab initio nearequilibrium potential energy surfaces (PESs) of the X_{2}B_{1} ground state and the excited state,A^{2}A_{2},have been reported ^{[51]}. The PES of the A^{2}A_{1} state has a C_{2v} equilibrium geometry and features strong coupling between the antisymmetric and the symmetric stretch modes. The two surfaces reproduce the experimental absorption spectrum well ^{[25]} and these surfaces are used here. The equilibrium geometries of the both PESs are far from linear with deep bending wells,which enable us to ignore the singularity problem found for linear geometries ^{[32, 33, 34, 35]}.
The absorption spectrum of the OClO molecule is obtained by Fourier transforming the time autocorrelation function of the initial wavefunction. The initial ground vibrational eigenfunction of the X_{2}B_{1} electronic state is obtained by a variational calculation using MorseMorse harmonic wavefunctions ^{[26]}. We note that even though the method used for obtaining the initial wavefunction leads to a nonhermitian Hamiltonian matrix,its influence on the quality of the ground wavefunction is marginal ^{[24]}. The method has been detailed in our previous work ^{[24, 25, 26]}. The absorption spectra below are uniformly broadened with a Lorentzian function of 20 cm^{1} FWHM (full width at half maximum) by damping the autocorrelation functions with an exponential function. In the following calculation,a Chebyshev propagator is applied to evolve the initial wavefunction on the excited PES. The round off error in the Chebyshev polynomial expansion is less than 1×10^{15}. Usually the timestep is chosen to be 0.7 fs and a total of 1024 time steps is enough to obtain a converged spectrum.
A. In Eckart coordinates 1. Hermitian Hamiltonian matrixThe Hamiltonian $\hat{H}_{ws}$ in Eq.(13) using the Fourier basis set FBST results in a Hermitian matrix as does the Hamiltonian $\hat{H}_{cs}$ in Eq.(40) using FBSC since these Hamiltonians are in explicitly symmetric forms ^{[43]}. We can expect that the results obtained with these two numerical models should agree with each other and show good numerical convergence. In the calculations with the Hamiltonian $\hat{H}_{ws}$ we set the grid ranges to [2.3,4.8] in atomic units for r_{1} and r_{2} and [1.3,2.7] in radians for θ. For $\hat{H}_{cs}$ we set the grid ranges [2.3,4.8] in atomic units for r_{1} and r_{2} and and [0.9,0.25] for cos θ. For each calculation,64×64×32 grid points are used (64 grid points along each radial coordinate and 32 along the angular coordinate). Convergence concerning grid ranges and grid spacing has been checked.
A typical part of the absorption spectra of the two numerical models is shown in Fig. 1(a). The corresponding real part of the correlation functions are shown in Fig. 1(b). The results are virtually identical for the two models,indicating good numerical convergel. The norm of the evolved wavepacket is not guaranteed to be conserved in a Chebyshev propagation and the timedependence of the norm gives a measure of convergence. The deviation of the norm from unity using $\hat{H}_{ws}$ in the FBST representation is shown in Fig. 2(a),and that using $\hat{H}_{cs}$ in the FBSC representation is shown in Fig. 2(b). We see that the norm is excellently conserved.
2. Nonhermitian Hamiltonian matrixIn this subsection,we only calculate the OClO absorption spectrum using the Hamiltonian $\hat{H}_{c}$ in Eq.(18) with the FBST basis set parameters used in subsection III (A). As we have discussed,the resulting Hamiltonian matrix is nonhermitian. Although the derivatives of the wavefunction can be evaluated accurately with FFT,the nonhermiticity of the Hamiltonian matrix makes the numerical results unreliable. The calculated absorption spectrum is shown in Fig. 3(b) together with the exact result of Fig. 1. There are a clear deviations. The heights of the peaks have been altered by the nonhermiticity of the Hamiltonian matrix,and a ghost peak has arisen. The timedependent norm of the wavefunction is shown in Fig. 3(a). We see drastical norm violation,which indicates that the nonhermitian matrix develops complex eigenvalues. Thus the numerical model is not stable. We note that the period of the initial oscillation of the timedependent norm approximately agrees with that of the symmetric stretch of the A^{2}A_{2} state.
3. Numerically nonhermitian Hamiltonian matrixFor numerical convenience,sometimes the expanded form in Eq.(11) of the Hamiltonian $\hat{H}_{ws}$ is used in numerical calculation. As Tuvi et al. pointed out ^{[43]},as long as the operators are not explicitly symmetric and the Fourier basis set is applied,the resulting Hamiltonian matrix becomes nonhermitian because of integration errors. The operators
$ \frac{{\cos \theta }}{{{\mu _3}{r_1}{r_2}}}\frac{{{\partial ^2}}}{{\partial {\theta ^2}}},\quad \frac{1}{{{\mu _3}}}\frac{{\sin \theta }}{{{r_1}{r_2}}}\frac{\partial }{{\partial \theta }},\quad \sin \theta \frac{\partial }{{\partial \theta }}\left( {\frac{1}{{{r_1}}}\frac{\partial }{{\partial {r_2}}} + \frac{1}{{{r_2}}}\frac{\partial }{{\partial {r_1}}}} \right) $ 
in Eq.(11) make the Hamiltonian martix numerically nonhermitian. We may expect that this nonhermiticity resulting purely from integration errors should be less drastic than analytic nonhermiticity. However,although the spectrum is satisfactory,which indicates that the nonhermiticity of the Hamiltonian matrix only slightly influences the numerical results,the norm of the wavefunction increases exponentially with propagation time,see Fig. 4. This was not observed using $\hat{H}_{ws}$ with the FBST basis set. The timedependent norm is also sensitive to the values of the Chebyshev parameters,which again indicates the instability of the numerical model. The FBST grid is used here too. The numerical results demonstrate that the nonhermiticity resulting from integration errors only cannot be neglected in some cases,especially for long time propagations.
B. In Radau coordinatesFrom the calculations above,we expect that the numerical models using the Hamiltonian $\hat{H}_{Rcs}$ in Eq.(44) with the FBSC basis set and using the Hamiltonian $\hat{H}_{R\theta s}$ in Eq.(45) with the FBST basis set should show an excellent agreement. The calculated absorption spectra using these two numerical models are shown in Fig. 5(a) and the real part of the correlation function in Fig. 5(b). Results from Fig. 1 are also included and all results agree excellently (to within the thickness of the lines). The norms of the two numerical models are well conserved during the propagation and,for brevity,we don't show the results here.
The calculated absorption spectrum using the Hamiltonian $\hat{H}_{R\theta a}$ in Eq.(43) with the FBST basis set is shown in Fig. 6. The result using Hamiltonian $\hat{H}_{R\theta s}$ in Eq.(45) with the FBST basis set is also shown. The deviations resulting from the nonhermiticity of the matrix caused by the transformation factor are obvious. The norm of the propagated wavefunction shows similar timedependent oscillation to that shown in Fig. 3 using the unsuitable Hamiltonian in Eq.(18) with the FBST basis set and we do not show it here.
Figure 7 illustrates the timedependent norm calculated using the Hamiltonian $\hat{H}_{R}$ in Eq.(46) with the FBSC basis set. The exponential increase of the norm completely destroys the autocorrelation function and its Fourier transform,the absorption spectrum,loses its meaning. We do not show the resulting spectrum but note that the surviving peaks seem to stay at the right energy positions.
In this section,when the Hamiltonian $\hat{H}_{R\theta s}$ in Eq.(45) and $\hat{H}_{R\theta a}$ in Eq.(43) were used. The grid ranges were [1.8,3.9] in atomic units for R_{1} and R_{2} and [1.9,3.1] in radians for θ. When the Hamiltonian $\hat{H}_{Rcs}$ in Eq.(44) is used,the same grid ranges for R_{1} and R_{2} are used but [0.9,0.25] for the variable C,i.e.,for cos θ. For each calculation,64×64×32 grid points are used (64 grid points along each radial coordinate and 32 along the angular coordinate). Convergence concerning grid ranges and grid spacing has been checked.
Ⅳ. DISCUSSIONWe would like to note that the two absorption spectra shown in Fig. 1 are not only virtually identical to each other,they are also virtually identical with those obtained using FFT and Hamiltonians expressed in Jacobi coordinates and in hyperspherical coordinates,where the convergence and the hermiticity of the Hamiltonian matrix has been carefully checked. Further,the lowest vibrational energies obtained using the relaxation method (imaginary time propagation method) ^{[52]} and the Hamiltonians in different coordinates,agree with each other better than 0.001 cm^{1}. This excellent agreement is in surprising contrast to the work of Katz et al. ^{[28]}. This may demonstrate that the FFT is particularly suitable in the case where the PES is well bounded,for the PES of the A^{2}A_{2} state of the OClO molecule exhibits a deep well around the equilibrium geometry ^{[24]}.
Among the numerical models,the calculation with the Hamiltonian in bond bond angle coordinates involves many more kinetic operators which greatly reduces the computation speed. Especially,numerical calculations using the transformed Hamiltonian in Radau coordinates in Eq.(45) are efficient. It only requires three forwardbackward FFTs for each Hamiltonian action,in contrast to the twelve needed for the symmetric Hamiltonian $\hat{H}_{ws}$ or $\hat{H}_{cs}$ in bond bond angle coordinates! Further,the Hamiltonian in Eq.(45) does not mix local and nonlocal operator(s) of the same coordinate,which allows the application of the splitoperator method with the FFT technique to propagate the wavepacket. The numerical model using the Hamiltonian $\hat{H}_{R\theta s}$ in Eq.(45) with the FBST basis set not only gives a simple and efficient way to numerical implementations but also significantly reduces the computation time ^{[24, 26]}.
We emphasize that both the analytic and numerical nonhermiticities of the Hamiltonian matrix discussed here give numerical implementations which are only conditionally stable. Especially in a long time propagation,the wavefunction norm violation may completely destroy the desired physical property of the simulated system. Interestingly,the nonhermiticity resulting from integration errors only always makes the norm of the propagated wavefunction increase exponentially,while the analytic nonhermiticity makes the norm oscillate. This may indicate that more complex and unphysical eigenvalues develop using an unsuitable Hamiltonian,which is reflected in the calculated spectra. We also note that,although the nonhermiticity of the Hamiltonian matrix makes the numerical model unstable and changes the heights of the physical peaks,the positions of the surviving physical peaks seem to be in the correct positions. In a Chebyshev propagation,the timedependent norm of the wavefunction is useful for checking the convergence of numerical model. When a Hamiltonian and a suitable basis set have been correctly constructed,the norm of the timedependent wavepacket should be well conserved.
Ⅴ. SUMMARYIn this work,requirements on the forms of the Hamiltonians to be used with a discrete Fourier basis set is discussed. It is emphasized that for a chosen basis set,the Hamiltonian cannot be arbitrarily transformed. Otherwise,the resulting Hamiltonian matrix becomes nonhermitian,which may lead to numerical errors. Further,it is recommended to use symmetric forms of the operators in numerical calculations as here demonstrated for timedependent wavepacket calculations using the FFT technique. Thus expanded forms of the Hamiltonian should also be avoided since they lead to nonhermitian Hamiltonian matrices due to integration errors. The influence of these two kinds of nonhermiticity of the Hamiltonian matrix,i.e.,resulted from unsuitable Hamiltonian and purely resulted from integration errors,on 3D timedependent wavepacket calculations were numerically investigated in bond bond angle and Radau coordinates. It was found that both of these nonhermiticities lead to the numerical calculation being only conditionally stable. The nonhermiticity problem may be marginal in some cases,perhaps a short time propagation calculation. Although the conclusions drawn here,based on the OClO absorption spectrum calculations using Fourier basis sets,suggests that the numerical errors can be expected to be casedependent,they should be generally applicable to calculations using DVR method.
Ⅵ. ACKNOWLEDGEMENTSThis work was supported by the National Basic Research Program of China (No.2013CB922200),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]  R. Dawes and T. Jr. Carrington, J. Chem. Phys. 122, 134101(2005). 
[2]  B. Poirier and L. C. Light, J. Chem. Phys. 111, 4869(1999). 
[3]  B. Poirier and L. C. Light, J. Chem. Phys. 113, 211(2000). 
[4]  B. Poirier and L. C. Light, J. Chem. Phys. 114, 6562(2001). 
[5]  R. Meyer, J. Chem. Phys. 52, 2053(1970). 
[6]  D. Baye and P. H. Heenen, J. Phys. A 19, 2041(1986). 
[7]  J. C. Light and Z. Bačić, J. Chem. Phys. 87, 4008(1987). 
[8]  H. Wei and T. Carrington Jr., J. Chem. Phys. 101, 1343(1994). 
[9]  J. C. Light and T. Jr. Carrington, Adv. Chem. Phys. 114, 263(2000). 
[10]  D. Yu, S. L. Cong, D. H. Zhang, and Z. G. Sun, Chin. J. Chem. Phys. 26, 755(2013). 
[11]  G. Avila and T. Jr. Carrington, J. Chem. Phys. 135, 064101(2011). 
[12]  X. S. Lin and Z. G. Sun, Chem. Phys. Lett. 621, 35(2015). 
[13]  G. Nyman and H. G. Yu, Rep. Prog. Phys. 63, 1001(2000). 
[14]  X. H. Liu, J. J. Lin, S. Harich, G. C. Schatz, and X. M. Yang, Science 289, 1536(2000). 
[15]  H. S. Lee and J. C. Light, J. Chem. Phys. 120, 5859(2004). 
[16]  J. Makarewicz, J. Phys. B 21, 1803(1988). 
[17]  B. T. Sutcliffe and J. Tennyson, Int. J. Quantum Chem. 91, 183(1991). 
[18]  J. Zúñiga, A. Bastida, and A. Requena, J. Chem. Soc. Faraday Trans. 93, 1681(1997). 
[19]  H. G. Yu and J. T. Muckerman, J. Mol. Spectrosc. 214, 11(2002). 
[20]  R. Schinke, Photodissociation Dynamics Spectroscopy and Fragmentation of Small Polyatomic Molecules, Cambridge University Press (1993). 
[21]  P. Marquetand, A. Materny, N. E. Henriksen, and V. Engel, J. Chem. Phys. 120, 5871(2004). 
[22]  N. E. Henrisken and V. Engel, Int. Rev. Phys. Chem. 20, 93(2001). 
[23]  G. M. Krishnan and S. Mahapatra, J. Chem. Phys. 118, 8715(2003). 
[24]  Z. Sun, N. Lou, and G. Nyman, J. Phys. Chem. A 108, 9226(2004). 
[25]  Z. Sun, N. Lou, and G. Nyman, J. Chem. Phys. 122, 054316(2005). 
[26]  Ğ. Barinovs, N. Marković, and G. Nyman, J. Chem. Phys. 111, 6705(1999). 
[27]  Z. Sun, N. Lou, and G. Nyman, Chem. Phys. 308, 317(2005). 
[28]  G. Katz, K. Yamashita, Y. Zeiri, and R. Kosloff, J. Chem. Phys. 116, 4403(2002). 
[29]  R. Kosloff, J. Phys. Chem. 92, 2087(1988). 
[30]  D. Kosloff and R. Kosloff, J. Comput. Phys. 52, 35(1983). 
[31]  M. D. Feit, J. A. Jr. Fleck, and A. Steiger, J. Comput. Phys. 47, 412(1982). 
[32]  R. N. Dixon, Chem. Phys. Lett. 190, 430(1992). 
[33]  R. N. Dixon, J. Chem. Soc. Faraday Trans. 88, 2575(1992). 
[34]  U. Manthe and H. Köppel, Chem. Phys. Lett. 175, 36(1991). 
[35]  U. Manthe, H. Köppel, and L. S. Cederbaum, J. Chem. Phys. 95, 1708(1991). 
[36]  C. E. Dateo, V. Engel, R. Almeida, and H. Metiu, Comp. Phys. Comm. 63, 435(1991). 
[37]  C. E. Dateo and H. Metiu, J. Chem. Phys. 95, 7392(1991). 
[38]  C. Woywod, Chem. Phys. Lett. 281, 168(1997). 
[39]  G. G. BalintKurti, R. N. Dixon, and C. C. Marston, Int. Rev. Phys. Chem. 11, 317(1992). 
[40]  R. C. Mowrey, Y. Sun, and D. J. Kouri, J. Chem. Phys. 91, 6519(1989). 
[41]  D. Lemoine, J. Chem. Phys. 101, 10526(1994). 
[42]  Ğ. Barinovs, N. Marković, and G. Nyman, J. Phys. Chem. A 105, 7441(2001). 
[43]  I. Tuvi and Y. B. Band, J. Chem. Phys. 107, 9079(1997). 
[44]  B. Podolsky, Phys. Rev. 32, 812(1928). 
[45]  M. L. Boas, Mathematical Methods in the Physical Sciences, New York:Wiley, (1983). 
[46]  P. R. Bunker, Molecular Symmetry and Spectroscopy, Canada:NRC Research Press, (1998). 
[47]  R. D. Bardo and M. Wolfsberg, J. Chem. Phys. 67, 593(1977). 
[48]  G. D. Carney, L. L. Sprandel, and C. W. Kern, Adv. Chem. Phys. 37, 305(1978). 
[49]  S. Carter and N. C. Handy, Mol. Phys. 57, 175(1986). 
[50]  B. R. Johnson and W. P. Reinhardt, J. Chem. Phys. 85, 4538(1986). 
[51]  K. A. Peterson, J. Chem. Phys. 108, 8864(1998). 
[52]  R. Kosloff and H. TalEzer, Chem. Phys. Lett. 127, 223(1986). 