Chinese Journal of Chemical Physics  2016, Vol. 29 Issue (1): 112-122

#### The article information

De-quanYu, HeHe, GunnarNyman, Zhi-gangSun
Hermiticity of Hamiltonian Matrix using the Fourier Basis Sets in Bond-Bond-Angle and Radau Coordinates
Chinese Journal of Chemical Physics, 2016, 29(1): 112-122
http://dx.doi.org/10.1063/1674-0068/29/cjcp1507141

### Article history

Accepted on July 13, 2015
Hermiticity of Hamiltonian Matrix using the Fourier Basis Sets in Bond-Bond-Angle and Radau Coordinates
De-quanYua, HeHea, b, GunnarNymanc, Zhi-gangSuna, d
Dated: Received on July 4, 2015; Accepted on July 13, 2015
a. State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical and Computa-tional Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Science, Dalian 116023, China;
b. School of Physics and Electronic Technology, Liaoning Normal University, Dalian 116029, China;
c. Department of Chemistry, Physical Chemistry, Göteborg University, SE-412 96 Göteborg, Sweden;
d. Center for Advanced Chemical Physics and 2011 Frontier Center for Quantum Science and Technol-ogy, University of Science and Technology of China, Hefei 230026, China
Part of the special issue for "the Chinese Chemical Society's 14th National Chemical Dynamics Symposium".
Author: Zhi-gangSun, E-mail: zsun@dicp.ac.cn
Abstract: In quantum calculations a transformed Hamiltonian is often used to avoid singularities in a certain basis set or to reduce computation time. We demonstrate for the Fourier basis set that the Hamiltonian can not be arbitrarily transformed. Otherwise, the Hamiltonian matrix becomes non-hermitian, which may lead to numerical problems. Methods for correctly constructing the Hamiltonian operators are discussed. Specific examples involving the Fourier basis functions for a triatomic molecular Hamiltonian (J=0) in bond-bond angle and Radau coordinates are presented. For illustration, absorption spectra are calculated for the OClO molecule using the time-dependent wavepacket method. Numerical results indicate that the non-hermiticity of the Hamiltonian matrix may also result from integration errors. The conclusion drawn here is generally useful for quantum calculation using basis expansion method using quadrature scheme.
Key words: Discrete variable representation    Hermiticity    Time-dependent wavepacket method    Absorption spectra
Ⅰ. INTRODUCTION

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 ro-vibrational spectra [13, 14, 15, 16, 17, 18, 19]. However,exact quantum calculations beyond triatomic and tetra-atomic molecules are still a formidable task. Improving the computational methods will be of continued interest in the theoretical field.

Quantum dynamics calculations solving the time-dependent molecular Schrödinger equation have proven to be invaluable for understanding dynamical process like photodissociation and femtosecond real-time 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 bond-bond angle coordinates [25, 28]. Fourier series or particle-in-the-box 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 fast-Fourier 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 time-dependent 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 time-dependent classic dynamics.

Several time-dependent 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 non-hermitian,leading to numerical problems. Tuvi et al. [43] have investigated the non-hermiticity 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 non-hermiticity resulting from integration errors.

In the following we will explicitly consider time-dependent wavepacket calculations using the Fourier basis sets and various triatomic Hamiltonians in bond bond-angle and Radau coordinates. The A2A1X2B1 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 non-hermiticity will be considered. The first type we will refer to as analytical non-hermitcity. It results from the use of a Hamiltonian and basis set in which the Hamiltonian matrix is non-hermitian even if it can be analytically constructed. The second type we will refer to as numerical non-hermiticity where the non-hermiticity 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 non-hermiticity 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) time-dependent wavepacket model in Eckart and Radau coordinates.

Ⅱ. KINETIC ENERGY OPERATORS,THEIR EVALUATION AND HERMITICITY A. Kinetic energy operators and wavefunction normalization

The 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 Born-Oppenheimer (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 pi is the momenta conjugate to the chosen coordinates qi (linearly independent),and $\mathbf{G}(\mathbf{q})$ is the well-known 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 xi denotes the Cartesian coordinates of atom i with a mass of mi. 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 g1/2 (normalized with volume element dτ=dx1dy1dz1 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 bond-bond 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 gA=1 and gB=g this equation gives the standard normalization defined by Eq.(5),and for gA=g and gB=1 it gives the normalization with weight factor unity,which Podolsky considered [44]. Because the factor gA1/4 is associated with the transformation of the wavefunction,we will refer to gA1/4 as the transformation factor in order to distinguish from the weight factor gB1/2. We note here that,if the weight factor gB1/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 bond-bond angle coordinates that have appeared in the literature according to the associated transformation (gA1/4) and weight (gB1/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 gB1/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 bond-bond angle coordinates,the weight factor gB1/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 transform

In a time-dependent 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{1-x^2}\displaystyle\frac{\textrm{d}}{\textrm{d}x}\sqrt{1-x^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. Hermiticity

Propagating the wavepacket using Fourier transform in a time-dependent 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 one-dimensional nuclear Hamiltonian in an explicitly symmetric form guarantees the hermiticity of the resulting Hamiltonian matrix. Their conclusion can be extended directly to multi-dimensional 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 gB1/4=1 in our calculation using a Fourier basis set.

Often the weight factor gB1/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 so-called 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 kn 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 gB1/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 gB 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 non-hermiticity do not arise.

The weight factor gB1/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(1-C^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 error

Considering 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 non-hermitian in the numerical calculations. For instance,the operator in Eq.(35) and

 $\hat{g}=-\left(1-C^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 Ci=$i\Delta C$,i=1,$\ldots$,N. Thus [43],when the model Hamiltonian used in an FFT time-dependent calculatition contain operators of the form in Eq.(37),the numerical results may be unreliable.

D. Hermiticity of the J=0 triatomic Hamiltonian in bond-bond angle and Radau coordinates in the Fourier basis

The commonly used triatomic Hamiltonian (J=0) in bond-bond angle coordinates $(r_{1},r_{2},\theta)$ written in the form of $\hat{H}_{ps}$ with weight factor gB 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 gB1/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 bond-bond angle coordinates,taking r1,r2 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 m1 and m2 are the relevant masses. Note that the variable θ here has a different meaning from that in bond bond-angle coordinates. This Hamiltonian has a volume element ${\alpha ^{ - 3}}\sin \theta d{R_1}d{R_2}$dθ where α2=m3/(m1+m2+m3),with weight factor gB1/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 bond-bond 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 non-hermitian,due to integration errors,despite the fact that it would be analytically hermitian. We will see later that this kind of non-hermiticity caused purely by integration errors is only conditionally stable and may lead to numerical problems for long time propagations.

Non-hermitian 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 non-hermiticity 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 anti-hermitician when the periodical bound condition are fulfilled.

In a time-dependent wavepacket calculation,the influence of the non-hermiticity of the Hamiltonian matrix is accumulated which is sensitively indicated by the time-dependence of the norm of the wavefunction. In a variational calculation,however,the non-hermiticity of the Hamiltonian matrix may have less impact on the numerical results. When a quantum calculation is done with a non-hermitian Hamiltonian matrix,a series of "ghost" (unphysical) peaks may arise. In some cases,the non-hermitian matrix may change the heights of the physical peaks. The positions of the physical peaks however normally remain correct. The non-hermitian 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 time-dependent wavepacket model.

Ⅲ. NUMERICAL EXAMPLES

In this section,we will show absorption spectra of the OClO molecule calculated with the time-dependent wavepacket method using Hamiltonians in bond-bond 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 non-hermiticity 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 non-hermiticity 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 non-hermitian 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 near-equilibrium potential energy surfaces (PESs) of the X2B1 ground state and the excited state,A2A2,have been reported [51]. The PES of the A2A1 state has a C2v equilibrium geometry and features strong coupling between the anti-symmetric 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 X2B1 electronic state is obtained by a variational calculation using Morse-Morse harmonic wavefunctions [26]. We note that even though the method used for obtaining the initial wavefunction leads to a non-hermitian 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 time-step 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 matrix

The 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 r1 and r2 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 r1 and r2 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 time-dependence 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.

 FIG. 1 (a) An expanded portion of the calculated absorption spectra and (b) the real part of the correlation functions $\langle\phi(0)|\phi(t)\rangle$ using $\hat{H}_{ws}$ in Eq.(13) with FBST (solid lines) and $\hat{H}_{cs}$ in Eq.(40) with FBSC (dotted lines).
 FIG. 2 Time-dependent norm deviation from its initial value 1 of the wavefunction using $\hat{H}_{ws}$ in Eq.(13) with FBST (a) and $\hat{H}_{cs}$ in Eq.(40) with FBSC (b).
2. Non-hermitian Hamiltonian matrix

In 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 non-hermitian. Although the derivatives of the wavefunction can be evaluated accurately with FFT,the non-hermiticity 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 non-hermiticity of the Hamiltonian matrix,and a ghost peak has arisen. The time-dependent norm of the wavefunction is shown in Fig. 3(a). We see drastical norm violation,which indicates that the non-hermitian matrix develops complex eigenvalues. Thus the numerical model is not stable. We note that the period of the initial oscillation of the time-dependent norm approximately agrees with that of the symmetric stretch of the A2A2 state.

 FIG. 3 (a)The time-dependent norm of the wavefunction resulting from using the non-hermitian matrix of the Hamiltonian $\hat{H}_c$ in Eq.(18) with FBST. (b)The solid line is the expanded portion of the calculated absorption spectra using the Hamiltonian $\hat{H}_c$ in Eq.(18) with FBST and the dotted line is the exact result shown in Fig. 1.
3. Numerically non-hermitian Hamiltonian matrix

For 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 non-hermitian 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 non-hermitian. We may expect that this non-hermiticity resulting purely from integration errors should be less drastic than analytic non-hermiticity. However,although the spectrum is satisfactory,which indicates that the non-hermiticity 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 time-dependent 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 non-hermiticity resulting from integration errors only cannot be neglected in some cases,especially for long time propagations.

 FIG. 4 (a) Time-dependent norm of the wavefunction which results from using $\hat{H}_w$ in Eq.(11) such that purely to integration errors occur. (b) The solid line is the expanded portion of the calculated absorption spectrum using the Hamiltonian $\hat{H}_{w}$ in Eq.(11) with FBST. The dotted line is the exact results shown in Fig. 1.

From 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.

 FIG. 5 (a) The solid line is the expanded portion of the calculated absorption spectra using the Hamiltonian $\hat{H}_{R\theta s}$ in Eq.(45) with FBST; the dotted line is from the model using the Hamiltonian $\hat{H}_{Rcs}$ in Eq.(44) with FBSC; and the dashed line corresponds with the numerical results shown in Fig. 1. (b) The real part of the auto-correlation function $\langle\phi(0)|\phi(t)\rangle$. The solid and dotted lines are the results using the Hamiltonians $\hat{H}_{R\theta s}$ and $\hat{H}_{Rcs}$ with FBST and FBSC,respectively. The dashed line corresponds with the numerical results shown in Fig. 1.

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 non-hermiticity of the matrix caused by the transformation factor are obvious. The norm of the propagated wavefunction shows similar time-dependent 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.

 FIG. 6 The solid line is the expanded portion of the calculated absorption spectra using Hamiltonian $\hat{H}_{R\theta a}$ in Eq.(43) with FBST and the dotted line is the exact result shown in Fig. 5.

Figure 7 illustrates the time-dependent 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 auto-correlation 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.

 FIG. 7 The time-dependent norm of the wavefunction obtained from evolving the non-hermitian matrix rsulting from the expanded Hamiltonian $\hat{H}_{R}$ in Eq.(46) with FBSC.

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 R1 and R2 and [1.9,3.1] in radians for θ. When the Hamiltonian $\hat{H}_{Rcs}$ in Eq.(44) is used,the same grid ranges for R1 and R2 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.

Ⅳ. DISCUSSION

We 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 A2A2 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 forward-backward 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 non-local operator(s) of the same coordinate,which allows the application of the split-operator 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 non-hermiticities 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 non-hermiticity resulting from integration errors only always makes the norm of the propagated wavefunction increase exponentially,while the analytic non-hermiticity 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 non-hermiticity 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 time-dependent 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 time-dependent wavepacket should be well conserved.

Ⅴ. SUMMARY

In 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 non-hermitian,which may lead to numerical errors. Further,it is recommended to use symmetric forms of the operators in numerical calculations as here demonstrated for time-dependent wavepacket calculations using the FFT technique. Thus expanded forms of the Hamiltonian should also be avoided since they lead to non-hermitian Hamiltonian matrices due to integration errors. The influence of these two kinds of non-hermiticity of the Hamiltonian matrix,i.e.,resulted from unsuitable Hamiltonian and purely resulted from integration errors,on 3D time-dependent wavepacket calculations were numerically investigated in bond bond angle and Radau coordinates. It was found that both of these non-hermiticities lead to the numerical calculation being only conditionally stable. The non-hermiticity 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 case-dependent,they should be generally applicable to calculations using DVR method.

Ⅵ. ACKNOWLEDGEMENTS

This 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.

References
 [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. Balint-Kurti, 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. Tal-Ezer, Chem. Phys. Lett. 127, 223(1986).