Chinese Journal of Chemical Physics  2019, Vol. 32 Issue (3): 333-342

#### The article information

Hai-mei Shi, Guang-hai Guo, Zhi-gang Sun

Numerical Convergence of the Sinc Discrete Variable Representation for Solving Molecular Vibrational States with a Conical Intersection in Adiabatic Representation

Chinese Journal of Chemical Physics, 2019, 32(3): 333-342

http://dx.doi.org/10.1063/1674-0068/cjcp1812275

### Article history

Accepted on: January 14, 2019
Numerical Convergence of the Sinc Discrete Variable Representation for Solving Molecular Vibrational States with a Conical Intersection in Adiabatic Representation
Hai-mei Shia,b , Guang-hai Guoa , Zhi-gang Sunb
Dated: Received on December 6, 2018; Accepted on January 14, 2019
a. School of Mathematics and Physics, Qingdao University of Science and Technology, Qingdao 266061, China;
b. State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical Computational Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian 116023, China
Abstract: Within the Born-Oppenheimer (BO) approximation, nuclear motions of a molecule are often envisioned to occur on an adiabatic potential energy surface (PES). However, this single PES picture should be reconsidered if a conical intersection (CI) is present, although the energy is well below the CI. The presence of the CI results in two additional terms in the nuclear Hamiltonian in the adiabatic presentation, i.e., the diagonal BO correction (DBOC) and the geometric phase (GP), which are divergent at the CI. At the same time, there are cusps in the adiabatic PESs. Thus usually it is regarded that there is numerical difficulty in a quantum dynamics calculation for treating CI in the adiabatic representation. A popular numerical method in nuclear quantum dynamics calculations is the Sinc discrete variable representation (DVR) method. We examine the numerical accuracy of the Sinc DVR method for solving the Schrödinger equation of a two dimensional model of two electronic states with a CI in both the adiabatic and diabatic representation. The results suggest that the Sinc DVR method is capable of giving reliable results in the adiabatic representation with usual density of the grid points, without special treatment of the divergence of the DBOC and the GP. The numerical uncertainty is not worse than that after the introduction of an arbitrary vector potential for accounting the GP, whose accurate form usually is not easy to obtain.
Key words: Discrete variable representation    Conical intersection    Adiabatic and diabatic representaton    Geometric phase
Ⅰ. INTRODUCTION

Non-adiabatic couplings are ubiquitous in polyatomic molecular excited state dynamics, such as in the photodissociation, photochemsitry, and isomerization processes [1-5]. With developments of the computational methods and computer resources, more and more attention has been payed to non-adiabatic molecular dynamics involving more than one electronic states [1, 6-13]. Conical intersection (CI) is one of the important non-adiabatic interactions, where two or more potential energy surfaces (PESs) are degenerate and the non-adiabatic couplings between these states cannot be removed simply in adiabatic representation. In the vicinity of the CIs, the Born-Oppenheimer (BO) approximation breaks down and the couplings between electronic and nuclear motion become important.

The CIs of multi-electronic states play a role as "funnels" which convert rapidly the extra electronic energy into nuclear motions [14, 15]. In the adiabatic picture, the presence of a CI results in two additional operators into the Hamiltonian: one is the geometric phase (GP) and the other is the diagonal BO correction (DBOC) [16-18]. The introduction of the GP explains a sign change of adiabatic electronic wave functions when it is conveyed along a closed path of nuclear configurations encircling the CI seam, which makes the electronic wave function double valued [19, 20]. Since the total wave function must be single valued, the double valued electronic wave function renders double valuedness of the nuclear wave function as well. In the adiabatic representation, the GP is associated with nuclear motion, which has important consequences in spectroscopy and dynamics for molecules influenced by the CIs, even in situations when the nuclear wave function is localized far from the region of the CI [3, 18, 20-49]. The introduced DBOC term is singular and repulsive around the CI, which accounts for the little populations of the adiabatic states in the region around the CIs.

Recently, it was found that during the photodissociation process the GP leads to an additional phase accumulation for fragments of the nuclear wave packet that transports around the CI to opposite sides [18, 28, 31, 50, 51]. This results in destructive interference, resulting in either spontaneously localizing the wave function or slowing down the nuclear dynamics compared with that in the case where the GP is not considered, which alters the lifetime of the dissociated state much [28, 31, 50, 51]. And, the GP shifts the spectrum of a bound state system by altering the pattern of nodes in the nuclear wave function [16, 22, 40, 43-45]. An understanding of the GP effects in a chemical reaction has developed only recently, owing mainly to a series of theoretical investigations and experiments on the hydrogen-exchange H+H$_2$$\toH_2+H reaction and its isotopologues [26, 27, 34, 35, 46, 48, 49, 52-54]. Very recently, using the high-resolution velocity map imaging (VMI) apparatus, in combination with a newly developed state-to-state quantum wavepacket reactive scattering theory using only the reactant Jacobi coordinate including the GP, the GP effect in the H+HD\toH_2+D reaction at a collision energy of 2.77 eV was clearly identified [54]. It was found that the experimentally measured fast angular oscillations around forward scattering of some specific product ro-vibrational states only agree with the theoretical predictions with inclusion of the GP, but is out-of-phase with the theoretical predictions without inclusion of the GP. The topological argument was developed by Althorpe and coworkers [48, 49], they explained how the GP effect was clearly observed in their work [50]. For understanding the role of CIs in a molecular dynamics, both in diabatic and adiabatic representation, numerical quantum dynamics investigations are essential. A popular numerical method for solving Schrödinger equation is the discrete variable representation (DVR) method. The DVR method based on classical polynomial or triangular functions using Gauss quadrature exhibits exponential convergence when the PESs are smooth, since essentially it is a spectral or pesudo-spectral method. However, such a spectral or pesudo-spectral method converges slowly when there are singularities or cusps, as those arise in a molecule with CIs in adiabatic representation. Thus it is a question that if the DVR method is capable of giving accurate results when it is simply applied with usual coordinates. In the work by Guo and his co-workers, cylindrical coordinates are applied to deal with the singularities and the cusps arising in the molecular Hamiltonian with a CI in the adiabatic presentation, and the DVR method was not applied [18]. In this work, the accuracy of the Sinc DVR method for dealing with the CI in normal coordinates both in adiabatic and diabatic representation is investigated, with and without inclusion of the GP and DBOC terms. It is found that the Sinc DVR method for solving the Schrödinger equation with the CI in diabatic representation gives exponential convergence, as expected. Although the Sinc DVR method cannot give numerical results with exponential convergence for solving the Schrödinger equation with the CI in adiabatic representation, especially for those states of energy around the CI, it does give the results with convergence of high order finite difference method. With the usual grid density of the Sinc DVR method, it is also found that the errors of the calculated vibrational energies in the adiabatic representation, which includes the GP and DBOC operators, are less than that introduced by the function forms of the vector potential (or the "mixing" angle for switching the diabatic and adiabatic representation) accounting for the GP. Since accurate function form of the mixing angle or the vector potential for a polyatomic molecule usually is unknown, in a practical calculation in adiabatic representation for accounting for the GP, an arbitrary function, which is simply capable of making the nuclear wavefunction doubleness when transported along the CI one circle, is applied. We conclude that the Sinc DVR method is accurate enough to investigate the molecular quantum dynamics with CIs in the adiabatic representation. Ⅱ. THEORY The well-known two-dimensional (2D) Jahn-Teller model [18, 19, 31, 55] was adopted in our present work, where the CI between two interacting electronic states is simply a point. The CI model of 2D represents a special case, where adiabatic non-adiabatic couplings can be completely removed by the adiabatic-to-diabatic transformation. The diabatic Hamiltonian assumes the following form (\hbar=1):  $$\hat{H}^\textrm{d}=\hat{H}_N+\hat{H}_e=\hat{T}_N \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}+ \begin{pmatrix} V_{11} & V_{12}\\ V_{21} & V_{22} \end{pmatrix} \label{dia}$$ (1) where  \hat{T}_N=-\frac{1}{2}\frac{\partial^2}{\partial x^2}-\frac{1}{2}\frac{\partial^2}{\partial y^2} which is the kinetic energy operator while the potential energy operator is a 2\times2 matrix. The corresponding adiabatic Hamiltonian can be written as  $$\hat{H}^{\textrm{a}}=\mathbf{U}^{\dagger}\hat{H}^\textrm{d}\mathbf{U}=\hat{T}_N \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}+\mathbf{U}^{\dagger}\cdot[\hat{T}_N, \mathbf{U}]+ \begin{pmatrix} V_{-} & 0\\ 0 & V_{+} \end{pmatrix} \label{ad}$$ (2) which is obtained by a unitary adiabatic-to-diabatic transformation  \begin{matrix} \mathbf{U}= \begin{pmatrix} \cos\alpha & \sin\alpha\\ -\sin\alpha & \cos\alpha \end{pmatrix} \end{matrix} where \alpha is the mixing angle between two diabatic states,  \alpha=\frac{1}{2}\arctan\frac{2V_{12}}{V_{22}-V_{11}} and \dagger denotes Hermitian conjugation. In the adiabatic representation, V_{\pm} is expressed as  V_{\pm}=\frac{1}{2}(V_{11}+V_{22})\pm \frac{1}{2}\sqrt{(V_{11}-V_{22})^2+4V{_{12}}^2} which is diagonal. The non-adiabatic couplings \mathbf{U}^{\dagger}\cdot[\hat{T}_N, \mathbf{U}] can be written as  \begin{eqnarray} \mathbf{U}^{\dagger}\cdot[\hat{T}_N, \mathbf{U}]=\begin{matrix} \begin{pmatrix} \hat{\tau}_{11} & i\hat{\tau}_{12}\\ -i\hat{\tau}_{21} & \hat{\tau}_{22} \end{pmatrix} \end{matrix}\nonumber \end{eqnarray} with the DBOC as  \hat{\tau}_{11}=\hat{\tau}_{22}=\frac{1}{2}\nabla\alpha\cdot\nabla\alpha and the derivative couplings as  \hat{\tau}_{12}=\hat{\tau}_{21}=\frac{1}{2}\left[(-i\nabla)^{\dagger}\cdot\nabla\alpha+\nabla\alpha\cdot(i\nabla)\right] In the presence of CI, although the adiabatic and diabatic representation are related by a unitary transformation, the diabatic Hamiltonian in Eq.(1) and the adiabatic Hamiltonian in Eq.(3) required different boundary conditions, i.e., single-valued vs. double-valued boundary conditions. To impose the double-valued boundary condition of the adiabatic wave function, the vector potential method proposed by Mead and Truhlar may be adopted [20]. In the method, a position-dependent phase factor \textrm{e}^{il\eta(x, y)} (l is an integer) is introduced into the adiabatic wave function, which makes the total wave function single valued and leads to the following Hamiltonian with a vector potential:  \begin{array}{l} \hat H_{{\rm{GP}}}^{\rm{a}} = {{\rm{e}}^{ - il\eta }}{{\hat H}^{\rm{a}}}{{\rm{e}}^{il\eta }} = {{\hat T}_N}\left( {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{{\hat \zeta }_{11}}}&{i{{\hat \zeta }_{12}}}\\ { - i{{\hat \zeta }_{21}}}&{{{\hat \zeta }_{22}}} \end{array}} \right) + \\ \left( {\begin{array}{*{20}{c}} {{V_ - }}&0\\ 0&{{V_ + }} \end{array}} \right) \end{array} (3) where  \begin{eqnarray} \hat{\zeta}_{11} = \hat{\zeta}_{22}\hspace{-0.05cm}=\hspace{-0.05cm}\hat{\tau}_{11}\hspace{-0.05cm}+\hspace{-0.05cm}\frac{l^2}{2}\nabla\eta\cdot\nabla\eta\hspace{-0.05cm}+\hspace{-0.05cm}\frac{l}{2}[(-i\nabla)^{\dagger}\cdot\nabla\eta\hspace{-0.05cm}+\hspace{-0.05cm}\nabla\eta\cdot (i\nabla)] \nonumber \\ \hat{\zeta}_{12} = \hat{\zeta}_{21}\hspace{-0.05cm}=\hspace{-0.05cm}\hat{\tau}_{12}+\frac{l}{2}(\nabla\alpha\cdot\nabla\eta+\nabla\eta\cdot\nabla\alpha)\nonumber \end{eqnarray} The form of function \eta(x, y) may be arbitrary, but only enforces the double-valued boundary condition of the adiabatic wave function. In practice, usually the exact solution of \alpha (x, y) is unknown, since diabatization of the potential energy surfaces is not easy to do. A unified form of the function, which is arbitrary but capable of making the nuclear wavefunction doubleness when transported along the CI one circle for \alpha and \eta, may be applied. We note that only for l as odd numbers, the GP effects are correctly included. Sometimes the adiabatic excited state can be well ignored, then often only the adiabatic ground state need to be considered. The Hamiltonian operators for the adiabatic ground state can be written as  \begin{eqnarray} \hat{H}^{(\textrm{a})}_1 = \hat{T}+V_-+\zeta_{11}\nonumber\\ = \hat{T}+V_-+\hat{\tau}_{11}+\frac{l^2}{2}\nabla\eta\cdot\nabla\eta+\nonumber\\ \frac{l}{2}[(-i\nabla)^{\dagger}\cdot\nabla\eta+\nabla\eta\cdot (i\nabla)] \label{adia1} \end{eqnarray} (4) where the GP effects are included for l as odd numbers. Quite often the GP effects are completely ignored in practical calculations, then the corresponding Hamiltonian operators can be written as  $$\hat{H}^{(\textrm{a})}_1=\hat{T}+V_-+\hat{\tau}_{11} \label{adia1a}$$ (5) In the following numerical calculations, the results calculated with the above Hamiltonian, i.e., two-state diabatic, one-state and two-state adiabatic with/without the GP and the DBOC operators, are compared with each other to illustrate their roles. In practical calculations, Hamiltonian operators in the adiabatic presentation involves divergent derivatives of \eta or \alpha and the non-differentiable cusps in the PESs (V_{\pm}), thus usually they are regarded to be difficult to deal with using modern DVR method in a simply way. On the other hand, since the exact form of the position-dependent phase factor \textrm{e}^{il\eta(x, y)} and the form of its resulted vector potential is unknown, is the numerical uncertainty caused by the usual DVR method in the adiabatic representation larger than that caused by the arbitrary forms of the vector potential? What is the exact convergence behaviour of the DVR method for dealing with the CI in the adiabatic representation including the GP effect? In this work, the 2D Jahn-Teller model defined in Ref.[18] is taken as the numerical example and the Sinc DVR method is used to answer these questions. The Sinc DVR is based on a uniform grid with spacing h, x_j=x_0+(j-1)h, and basis functions are chosen to be the Sinc interpolation function S_{j, h}(x) as [56-58]  \begin{eqnarray} S_{j, h}(x) = \frac{{\rm Sinc}(x-jh)}{h}\nonumber\\ = \begin{cases} \displaystyle\frac{\sin[\pi(x-jh)/h]}{\pi(x-jh)/h} & \quad \text{for } x\ne x_j\\ 1 & \quad \text{for } x = x_j \end{cases}\nonumber \end{eqnarray} The grid interval between x_0 and x_1 is (N-1)h, where N is the number of the grid points. The matrix form of first and second derivative operators on this uniform grid are given by  \frac{\partial}{\partial x} = \begin{cases} \displaystyle\frac{(-1)^{i-j}}{(i-j)h} & \quad \text{for } i\ne j\\ 0 & \quad \text{for } i = j \end{cases} and  \frac{\partial^2}{\partial x^2} = \begin{cases} -\displaystyle\frac{2(-1)^{i-j}}{(i-j)^2h^2} & \quad \text{for } i\ne j\\ -\displaystyle\frac{\pi^2}{3h^2} & \quad \text{for } i = j \end{cases} The Sinc functions are a set of basis functions orthogonal with unit weight function. With the first and second derivative operators given above, the Schrödinger equation can be written in matrix form directly and solved. Ⅲ. RESULTS AND DISCUSSION We first examine the well-studied 2D Jahn-Teller model [18, 19, 31, 55]:  V_{11} = \frac{\omega{_1}^2}{2}\left(x+\frac{a}{2}\right)^2+\frac{\omega{_2}^2}{2}y^2 (6)  V_{22} = \frac{\omega{_1}^2}{2}\left(x-\frac{a}{2}\right)^2+\frac{\omega{_2}^2}{2}y^2-\Delta (7)  V_{12} = cy (8) with \omega_1=1, \omega_2=1, a=6, \Delta=0.01, c=3. The PESs in adiabatic representation are presented in FIG. 1. It is seen that the CI is quite low in energy, thus it would influence the low lying vibrational states significantly. This model will be denoted as model-Ⅰ in the following. The grid ranges for x and y in the calculations are set as [-12.0, 12.0].  FIG. 1 The adiabatic potential energy surfaces of model-Ⅰ with parameters (in atomic units) as \omega_1=1, \omega_2=1, a=6, \Delta=0.01, c=3. The lowest two eigenfunctions and their corresponding eigenenergies are presented in FIG. 2. The plots in the left column are calculated with Eq.(1), the plots in the middle column are calculated with Eq.(5), whereas the plots in the right column are calculated with Eq.(4). In the calculations, \alpha (x, y) and \eta(x, y) are chosen as the same exact form, which can be written out explicitly  \begin{eqnarray} \alpha=\eta = \frac{1}{2}\arctan\frac{2V_{12}}{V_{22}-V_{11}}\nonumber\\ = -\frac{1}{2}\arctan\left(\frac{2cy}{\Delta+a\omega_1x}\right) \end{eqnarray} (9)  FIG. 2 Eigenfunctions for n=0 and 1 of model-Ⅰ for (a, d) diabaitic PES, (b, e) one-state adiabatic PES without the GP and DBOC, and (c, f) one state adiabatic PES with the GP and DBOC. The corresponding eigenvalues are noted in the panels. The wave functions in panel (a) and (d) agree well with those in panel (c) and (f), but with one more node than that in panel (b) and one less node than that in panel (e). Comparing the plots in the middle and right columns of FIG. 2, we see that the patterns of the eigenfunctions have been altered much with inclusion of the GP effects. There is an extra node for the ground state, but one node is missed for the first excited state. Only with inclusion of the GP effects, the calculated wave functions agree with those calculated in the "exact" diabatic representation. The convergence of typical vibrational states as a function of number of the grid points (n=n_x=n_y) are displayed in FIG. 3 for different sets of calculations: (ⅰ) diabatic two-state (Eq.(1)), (ⅱ) adiabatic two-state with even l and odd l (Eq.(3)), (ⅲ) adiabatic one-state with GP with even l and odd l /without GP and with/without DBOC. The error in the plots is defined as  {\rm{Error}} = {\log _{10}}\left( {\frac{{|{E_n} - E_n^0|}}{{E_n^0}}} \right) (10)  FIG. 3 Numerical convergence of eigenvalues E_n with n=0, 10, and 20 for model-Ⅰ PES with respective to the number of grid points n_x(=n_y), for (a) diabatic, (b) one-state adiabatic without GP and DBOC, (c) one-state adiabatic without GP but with DBOC, (d) one-state adiabatic with GP but without DBOC, (e) one-state adiabatic with GP and DBOC, and (f) two-state adiabatic calculations. The error is defined as log_{10}$$\displaystyle{\left(\frac{| E_n-E^0_n|}{E^0_n}\right)}$. The reference values $E^0_n$ are calculated with $n_x$=$n_y$=100. The diabatic calculations exhibit exponential convergence and all of the other calculations exhibit slower convergence.

where $E_n^0$ is calculated using $n_x$=$n_y$=100. $E_i^0$ are not expected to be the exactly converged values, however, they can be taken as good references to examine the numerical convergence. The plots in FIG. 3(a) tell us that the calculations in diabatic representation using the Sinc DVR exhibit exponential convergence. However, the results in other panels converge in a much slower way, especially for the ground state. For the one-state adiabatic without the GP and the DBOC in FIG. 3(b) and one-state adiabatic without the GP but with DBOC, the numerical convergences for state $n$=20 are particularly slow.

For understanding their different convergence behaviours, the eigenfunctions of $n$=10 and 20 for one-state adiabatic without the GP and the DBOC and one-state adiabatic with GP and DBOC are presented in FIG. 4. It is observed that there is considerable amplitude for the wave function around the CI for $n$=20 of one-state adiabatic without the GP and the DBOC. Therefore, this wave function feels the cusp and divergence of the CI strongly and converges slowly. However, the wave functions for $n$=20 avoid the region around the CI for one-state adiabatic with GP and DBOC. The wave functions for $n$=10 avoid the region around the CI for one-state adiabatic both with and without GP and DBOC. Thus these three states converge faster. Anyway, from the numerical convergence shown in FIG. 3 it needs to point out that inclusion of the GP and the DBOC in the adiabatic calculations does not introduce more of the numerical difficulty for the Sinc DVR method, as compared with the calculations without the DBOC and GP.

 FIG. 4 Eigenfunctions with $n$=10 and 20 for (a, b) one-state adiabatic without GP and DBOC, (c, d) for one state-adiabatic with GP and DBOC. There are considerable amplitudes around the region of the conical intersection for wave functions in panel (a) and (b). However, there is negligible amplitude of the wave function in the other two panels (c, d) around the conical intersection region, especially for that in panel (d).

For obtaining clearer impression of the numerical convergence of all these kinds of calculations, numerical results with many digits calculated using $n_x$=$n_y$=40 (black) and $n_x$=$n_y$=100 (blue and with states denoted by Roman numbers) are listed in Table Ⅰ. The digits which are different in these two sets of calculations are in bold font. For most of the results, the convergence is better than the fifth decimal number after the decimal point.

Table Ⅰ Eigenvalues (in atomic units) for diabatic and adiabatic models. The parameters (in atomic units) used in model-Ⅰ PESs: $\omega_1$=1, $\omega_2$=1, $a$=6, $\Delta$=0.01, $c$=3.

Since the CI of model Ⅰ is of low energy, most of the low lying states are influenced by the CI strongly. For a model with CI of much higher energy, one may expect that the numerical convergence for the Sinc DVR would be much better. In the following calculations, we will examine the numerical performance of the Sinc DVR for the model with $\omega_1$=1, $\omega_2$=1, $a$=10, $\Delta$=0.01, $c$=5. The PESs in adiabatic representation are displayed in FIG. 5. It is seen that the CI now is of much higher energy. We will denote this model as model Ⅱ. The grid ranges for $x$ and $y$ in the calculations are set as [-12.0, 12.0] again.

 FIG. 5 The adiabatic potential energy surfaces of model-Ⅱ with parameters (in atomic units) as $\omega_1$=1, $\omega_2$=1, $a$=10, $\Delta$=0.01, $c$=5.

Similar to FIG. 2, two typical eigenfunctions and their corresponding eigenenergies are presented in FIG. 6. Similarly, the plots in the left column are calculated with Eq.(1), the plots in the middle column are calculated with Eq.(5), whereas the plots in the right column are calculated with Eq.(4). In the calculations, $\alpha (x, y)$ and $\eta(x, y)$ are chosen as the same form, which can be written out explicitly as Eq.(6). From FIG. 6 we see again that, with inclusion of the GP effects, the eigenfunctions have been altered. There is an extra node for the ground state, but one node is missed for the first excited state, which is the same as those in FIG. 2. However, now there is an additional node along the radial direction taken CI as the centre of a circle.

 FIG. 6 Eigenfunctions for $n$=14 and 15 of model-Ⅱ for (a, d) diabaitic PES, (b, e) for $n$=15 and 16 for one-state adiabatic PES without the GP and DBOC, and (c, f) for $n$=14 and 15 one state adiabatic PES with the GP and DBOC. The corresponding eigenvalues are noted in the panels. The wave functions in panel (a) and (d) agree well with those in panel (c) and (f), but quite different from those in panel (b) and (e), which are similar to the plots in FIG. 2.

Similar to FIG. 3, the convergence of typical vibrational states as a function of number of the grid points ($n$=$n_x$=$n_y$) is displayed in FIG. 7 for different sets of calculations: (ⅰ) diabatic two-state (Eq.(1)), (ⅱ) adiabatic two-state with even $n$ and odd $n$ (Eq.3), (ⅲ) adiabatic one-state with even $l$ and odd $l$ /without the GP and with/without the DBOC. The error in the plots is again defined as Eq.(10), where $E_n^0$ is also calculated using $n_x$=$n_y$=100. The plots in FIG. 7(a) tell us that the calculations in diabatic representation using the Sinc DVR exhibit exponential convergence. Interestingly, the results in other panels now converge in a similar way. This must be due to the high energy of the CI, which is inaccessible to the states studied in the plots, and thus all of them exhibit fast numerical convergence as a function of number of the grid points. Thus, in the case where the CI is of high energy, the Sinc DVR method indeed works well in the adiabatic representation, regardless including the GP and/or the DBOC or not. Anyway, we should note that only with inclusion of the GP, the eigenenergies of the adiabatic models agree with those of the diabatic two state calculations, as we see the wave functions and their energies in FIG. 6 and the results shown in Table Ⅱ.

 FIG. 7 Numerical convergence of eigenvalues $E_n$ with $n$=0, 10 and 20 for model-Ⅱ PES with respective to the number of grid points $n_x$(=$n_y$), for (a) diabatic, (b) one-state adiabatic without GP and DBOC, (c) one-state adiabatic without GP but with DBOC, (d) one-state adiabatic with GP but without DBOC, (e) one-state adiabatic with GP and DBOC, (f) two-state adiabatic calculations. The error is defined as log$_{10}\left( {\frac{{|{E_n} - E_n^0|}}{{E_n^0}}} \right)$. The reference values $E^0_n$ are calculated with $n_x$=$n_y$=100. All of the calculations, irrespective of the divergence of GP and DBOC in the adiabatic calculations, exhibit exponential convergence.
Table Ⅱ Eigenvalues (in atomic units) for two-state and one-state adiabatic models with different vector potentials, where the GP and DBOC are both included. The parameters (in atomic units) used in model-Ⅰ PESs: $\omega_1$=1, $\omega_2$=1, $a$=6, $\Delta$=0.01, $c$=3 and model-Ⅱ PESs: $\omega_1$=1, $\omega_2$=1, $a$=10, $\Delta$=0.01, $c$=5.

Quite often, in a practical calculation, only the ground PES in adiabatic representation is known. Thus, $\alpha (x, y)$ and $\eta(x, y)$, whose exact expressions are required for the GP and the DBOC calculations, usually are unknown. In this case, there is no unique definition of the angles, $\alpha (x, y)$ and $\eta(x, y)$, which can be any function of the coordinates, for which $\alpha (x, y)$$\to$$\alpha(x, y)$+2$\pi$ and $\eta(x, y)$$\to$$\eta(x, y)$+2$\pi$, describe a complete path around the CI. For the 2D model we studied here, for convenience we may choose $\alpha'$=$\eta'$=$\displaystyle{\frac{1}{2}\arctan\frac{y}{x}}$ to examine the effect of the function forms of the mixing angle to the adiabatic calculations for including the GP and the DBOC terms.

Typical eigenenergies of model Ⅰ and Ⅱ in two-state and one-state adiabatic representation for the GP and the DBOC, with $\alpha'$=$\eta'$=$\displaystyle{\frac{1}{2}\arctan\frac{y}{x}}$ and $\alpha$=$\eta$=$-\displaystyle{\frac{1}{2}\arctan\frac{2cy}{\Delta+a\omega_1x}}$, are listed in Table Ⅱ. The numbers in black colour are calculated with $n_x$=$n_y$=40 and those in blue with $n_x$=$n_y$=100. For the numbers in blue, the digits after the decimal point, which are the same for $\alpha'$ and $\alpha$, are in bold font. As we compare the results calculated with $n_x$=$n_y$=40 and $n_x$=$n_y$=100, it is seen that the uncertainty introduced by the $\alpha'$, which is different from exact $\alpha$, is larger than the numerical error introduced by the Sinc DVR method. Therefore, we may conclude that the Sinc DVR method works well for a molecular system with CIs in a usual coordinate. Correspondingly, it is not surprising to see that the Sinc DVR method could give accurate state-to-state results for H+H$_2$ and H+O$_2$ reaction and its isotopologues in adiabatic representation including the GP and the DBOC terms [34, 35, 54, 59, 60].

Finally, we want to emphasise that the GP influences the states of a molecular quantum system much, regardless of the energy of the CI. However, different energies of the CI, do result in different effects on the numerical convergence of the Sinc DVR method. The Sinc DVR method converges fast for the states far from the CI. For the states of energies around the CI, the best choice is to set up the "exact" diabatic representation to describe. Anyway, in the case where the diabatic representation is unavailable, the numerical uncertainty introduced by the Sinc DVR is not larger than that introduced by the arbitrary function forms of the vector potential accounting for the GP effects.

Ⅳ. CONCLUSION

In this work, we examine the numerical convergence of the Sinc DVR method for calculating vibrational eigenstates of 2D Jahn-Teller model both in diabatic and adiabatic representation. In the adiabatic representation, models both with inclusion of the GP and DBOC Hamiltonian operators and without inclusion of them were considered. At the same time, the effects of the different function forms of the mixing angle for transforming the diabatic and adiabatic representation, and the effect of the resulted vector potential for including the GP, on the eigenstates in the adiabatic representation are examined. It was found that the Sinc DVR method works well for 2D Jahn-Teller model not only in diabatic representation, but also in the adiabatic representation with inclusion of the GP and the DBOC, even there is cusp and divergent terms in the Hamiltonian operators. At the same time, different function forms of the mixing angle for transforming the diabatic and adiabatic representation, which also gives the vector potential accounting the GP and the DBOC operator, do have small influence on the eigenstate of the system, which is even larger than the numerical error introduced by the Sinc DVR method with the modest number of grid points.

Ⅴ. ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No.21733006 and No.21825303) and NSFC Center for Chemical Dynamics (No.21688102), the Strategic Priority Research Program of Chinese Academy of Sciences (No.XDB17000000), and the Chinese Academy of Sciences, and the Key Research Program of the Chinese Academy of Sciences.

Reference
 [1] W. Domcke, D. R. Yarkony, and H. Köppel, Conical Intersections: Theory, Computation and Experiment, Singapore: World Scientific, (2011). [2] M. Barbatti, and H. Lischka, J. Am. Chem. Soc. 130, 6831(2008). DOI:10.1021/ja800589p [3] M. N. R. Ashfold, B. Cronin, A. L. Devine, R. N. Dixon, and M. G. D. Nix, Science 312, 1637(2006). DOI:10.1126/science.1125436 [4] W. K. Chen, X. Y. Liu, W. H. Fang, P. O. Dral, and G. L. Cui, J. Phys. Chem. Lett. 9, 6702(2018). DOI:10.1021/acs.jpclett.8b03026 [5] X. Y. Liu, G. L. Cui, and W. H. Fang, Theo. Chem. Acc. 136, 8(2017). DOI:10.1007/s00214-016-2029-z [6] J. W. Zwanziger, and E. R. Grant, J. Chem. Phys. 87, 2954(1987). DOI:10.1063/1.453083 [7] H. Koizumi, and S. Sugano, J. Chem. Phys. 102, 4472(1995). DOI:10.1063/1.469495 [8] W. Domcke, and G. Stock, Adv. Chem. Phys. 100, 1(1997). [9] G. A. Worth, and L. S. Cederbaum, Ann. Rev. Phys. Chem. 55, 127(2004). DOI:10.1146/annurev.physchem.55.091602.094335 [10] K. Kleinermanns, N. Nachtigallová, and M.S. de Vries, Int. Rev. Phys. Chem. 32, 1(2013). DOI:10.1080/0144235X.2012.729331 [11] B. G. Levine, and T. J. Martínez, Ann. Rev. Phys. Chem. 57, 613(2007). [12] D. R. Yarkony, Rev. Mod. Phys. 68, 985(1996). DOI:10.1103/RevModPhys.68.985 [13] A. V. Akimov, and O. V. Prezhdo, J. Chem. Theo. Comput. 9, 4959(2013). DOI:10.1021/ct400641n [14] B. Balzer, S. Hahn, and G. Stock, Chem. Phys. Lett. 379, 351(2003). DOI:10.1016/j.cplett.2003.08.052 [15] S. Hahn, and G. Stock, J. Phys. Chem. B 104, 1146(2000). DOI:10.1021/jp992939g [16] A. F. Izmaylov, J. R. Li, and L. Joubert-Doriol, J. Chem. Theo. Comput. 12, 5278(2016). DOI:10.1021/acs.jctc.6b00760 [17] S. Mahapatra, H. Koppel, and L. S. Cederbaum, J. Phys. Chem. A 105, 2321(2001). DOI:10.1021/jp003784b [18] C. J. Xie, D. R. Yarkony, and H. Guo, Phys. Rev. A 95, 022104(2017). DOI:10.1103/PhysRevA.95.022104 [19] H.C. Longuet-Higgins, U. Opik, M. H. L. Pryce, and R. A. Sack, Proc. R. Soc. London Ser. A 244, 1(1958). DOI:10.1098/rspa.1958.0022 [20] C. A. Mead, and D. G. Truhlar, J. Chem. Phys. 70, 2284(1979). DOI:10.1063/1.437734 [21] V. J. Barclay, C. E. Dateo, I. P. Hamilton, B. Kendrick, R. T. Pack, and D. W. Schwenke, J. Chem. Phys. 103, 3864(1995). DOI:10.1063/1.470039 [22] C. A. Mead, Rev. Mod. Phys. 64, 51(1992). DOI:10.1103/RevModPhys.64.51 [23] S. C. Althorpe, J. Chem. Phys. 124, 084105(2006). DOI:10.1063/1.2161220 [24] Z. G. Lan, W. Domcke, V. Vallet, A. L. Sobolewski, and S. Mahapatra, J. Chem. Phys. 122, 224315(2005). DOI:10.1063/1.1906218 [25] C. A. Mead, and D. G. Truhlar, J. Chem. Phys. 77, 6090(1982). DOI:10.1063/1.443853 [26] B. K. Kendrick, J. Hazra, and N. Balakrishnan, Phys. Rev. Lett. 115, 153201(2015). DOI:10.1103/PhysRevLett.115.153201 [27] B. K. Kendrick, J. Hazra, and N. Balakrishnan, Nat. Comm. 6, 7918(2015). DOI:10.1038/ncomms8918 [28] F. Bouakline, Chem. Phys. 442, 31(2014). DOI:10.1016/j.chemphys.2014.02.010 [29] I. G. Ryabinkin, L. Joubert-Doriol, and A. F. Izmaylov, J. Chem. Phys. 140, 214116(2014). DOI:10.1063/1.4881147 [30] L. Joubert-Doriol, I. G. Ryabinkin, and A. F. Izmaylov, J. Chem. Phys. 139, 234103(2014). [31] I. G. Ryabinkin, and A. F. Izmaylov, Phys. Rev. Lett. 111, 220406(2013). DOI:10.1103/PhysRevLett.111.220406 [32] R. N. Dixon, T. A. A. Oliver, and M. N. R. Ashfold, J. Chem. Phys. 134, 194303(2011). DOI:10.1063/1.3585609 [33] M. G. D. Nix, A. L. Devine, R. N. Dixon, and M. N. R. Ashfold, Chem. Phys. Lett. 463, 305(2008). DOI:10.1016/j.cplett.2008.08.085 [34] F. Bouakline, S. C. Althorpe, and D. Pelaéz Ruiz, J. Chem. Phys. 128, 124322(2008). DOI:10.1063/1.2897920 [35] J.C. Juanes-Marcos, and S. C. Althorpe, J. Chem. Phys. 122, 204324(2005). DOI:10.1063/1.1924411 [36] D. Babikov, B. K. Kendrick, P. Zhang, and K. Morokuma, J. Chem. Phys. 122, 044315(2005). DOI:10.1063/1.1824905 [37] Z. Xu, M. Baer, and A. J. C. Varandas, J. Chem. Phys. 112, 2746(2000). DOI:10.1063/1.480848 [38] B. K. Kendrick, J. Chem. Phys. 112, 5679(2000). DOI:10.1063/1.481143 [39] A. J. C. Varandas, and H. G. Yu, J. Chem. Soc. Faraday Trans. 93, 819(1997). DOI:10.1039/a605777b [40] B. K. Kendrick, Phys. Rev. Lett. 79, 2431(1997). DOI:10.1103/PhysRevLett.79.2431 [41] B. K. Kendrick, and R. T. Pack, J. Chem. Phys. 104, 7475(1996). DOI:10.1063/1.471460 [42] R. Baer, D. M. Charutz, R. Kosloff, and M. Baer, J. Chem. Phys. 105, 9141(1996). DOI:10.1063/1.472748 [43] B. E. Applegate, T. A. Barckholtz, and T. A. Miller, Chem. Soc. Rev. 32, 38(2003). DOI:10.1039/A910269H [44] M. S. Child, Adv. Chem. Phys. 124, 1(2002). [45] H. von Busch, Dev Vas, H. A. Eckel, S. Kasahara, J. Wang, W. Demtröder, P. Sebald, and W. Meyer, Phys. Rev. Lett. 81, 4584(1998). DOI:10.1103/PhysRevLett.81.4584 [46] F. Bouakline, S. C. Althorpe, P. Larregaray, and L. Bonnet, Mol. Phys. 108, 969(2010). DOI:10.1080/00268971003610218 [47] B. J. Rao, R. Padmanaban, and S. Mahapatra, Chem. Phys. 333, 135(2007). DOI:10.1016/j.chemphys.2007.01.012 [48] J.C. Juanes-Marcos, S. C. Althorpe, and E. Wrede, J. Chem. Phys. 126, 044317(2007). DOI:10.1063/1.2430708 [49] J.C. Juanes-Marcos, S. C. Althorpe, and E. Wrede, Science 309, 1227(2005). DOI:10.1126/science.1114890 [50] J. Schön, and H. Köppel, J. Chem. Phys. 103, 9292(1995). DOI:10.1063/1.469988 [51] C. J. Xie, J. Y. Ma, X. L. Zhu, D. R. Yarkony, D. Q. Xie, and H. Guo, J. Am. Chem. Soc. 138, 7828(2016). DOI:10.1021/jacs.6b03288 [52] B. K. Kendrick, J. Phys. Chem. A 107, 6739(2003). [53] B. K. Kendrick, J. Hazra, and N. Balakrishnan, New J. Phys. 18, 123020(2016). DOI:10.1088/1367-2630/aa4fd2 [54] D. F. Yuan, Y. F. Guan, W. T. Chen, H. L. Zhao, S. R. Yu, C. Luo, Y. X. Tan, T. Xie, X. A. Wang, Z. G. Sun, D. H. Zhang, and X. M. Yang, Science 362, 1289(2018). DOI:10.1126/science.aav1356 [55] H. Köppel, W. Domcke, and L. S. Cederbaum, Adv. Chem. Phys. 57, 59(1984). [56] D. T. Colbert, and W. H. Miller, J. Chem. Phys. 96, 1982(1992). DOI:10.1063/1.462100 [57] P. Amore, J. Phys. A: Math. Gen. 39, L349(2006). DOI:10.1088/0305-4470/39/22/L01 [58] J. Lo, and B. D. Shizgal, J. Chem. Phys. 125, 194108(2006). DOI:10.1063/1.2378622 [59] B. Kendrick, and R. T. Pack, J. Chem. Phys. 104, 7502(1996). DOI:10.1063/1.471461 [60] B. Kendrick, and R. T. Pack, J. Chem. Phys. 106, 3519(1997). DOI:10.1063/1.473449

a. 青岛科技大学数理学院，青岛 266061;
b. 中国科学院大连化学物理研究所，分子反应动力学国家重点实验室和理论计算中心，大连 116023