Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (4): 608-612

#### The article information

Rui-xue Xu, Xue-cheng Tao, Yao Wang, Yang Liu, Hou-dao Zhang, YiJing Yan

A Hierarchical-Equation-of-Motion Based Semiclassical Approach to Quantum Dissipation

Chinese Journal of Chemical Physics, 2018, 31(4): 608-612

http://dx.doi.org/10.1063/1674-0068/31/cjcp1807172

### Article history

Received on: July 19, 2018
Accepted on: July 22, 2018
A Hierarchical-Equation-of-Motion Based Semiclassical Approach to Quantum Dissipation
Rui-xue Xu, Xue-cheng Tao, Yao Wang, Yang Liu, Hou-dao Zhang, YiJing Yan
Dated: Received on July 19, 2018; Accepted on July 22, 2018
Hefei National Laboratory for Physical Sciences at the Microscale and Department of Chemical Physics and Synergetic Innovation Center of Quantum Information and Quantum Physics and Collaborative Innovation Center of Chemistry for Energy Materials(iChEM), University of Science and Technology of China, Hefei 230026, China
*Author to whom correspondence should be addressed. Rui-xue Xu, E-mail:rxxu@ustc.edu.cn
Abstract: We present a semiclassical (SC) approach for quantum dissipative dynamics, constructed on basis of the hierarchical-equation-of-motion (HEOM) formalism. The dynamical components considered in the developed SC-HEOM are wavepackets' phase-space moments of not only the primary reduced system density operator but also the auxiliary density operators (ADOs) of HEOM. It is a highly numerically efficient method, meanwhile taking into account the high-order effects of system-bath couplings. The SC-HEOM methodology is exemplified in this work on the hierarchical quantum master equation[J. Chem. Phys. 131, 214111 (2009)] and numerically demonstrated on linear spectra of anharmonic oscillators.
Key words: Quantum dissipation    Semiclassical approach    Hierarchical-equation-ofmotion    Anharmonic oscillators    Phase-space moments
Ⅰ. INTRODUCTION

Quantum dissipation plays crucial roles in many fields of modern science. It describes the relaxation dynamics of a centered system to thermal equilibrium embedded in a thermal environment. Assuming Gaussian property for the environment statistics, exact approaches have been developed including the path-integral influence functional formalism [1], its differential equivalences via the hierarchical-equation-of-motion (HEOM) [2-7] or the related stochastic methods [8-12], and a dissipatonequation-of-motion approach as a second-quantization advancement of the HEOM formalism [13-19]. Theoretical developments have also been extended to solventinduced polarization dynamics [14-16] and some classes of non-Gaussian environments [12, 14, 17-19].

The HEOM approach has become nowadays a standard method to simulate non-perturbative and nonMarkovian quantum dissipative dynamics. However, since its dynamical components are a set of density matrices, the HEOM method encounters low numerical efficiency in simulating multi-level systems, even with the state-of-the-art filtering propagator [20] and Padé spectrum decomposition [7, 21, 22] which have largely accelerated its time evolution. For example, the HEOM will be difficult to apply if the primarily interested systems are nuclear motions along certain potentials meanwhile strongly influenced by a thermal environment at a relatively high temperature. This is because it is hard to find a proper representation in advance to describe such a system with a small size of basis set.

On the other hand, occasionally our objective is just the physical observation of dynamic variables such as coordinates or momenta, or other phase-space properties as the nuclear degrees of freedom are concerned. Equations of motion of those interested expectation values can be derived to avoid the time-consuming simulation of density matrices evolutions. In quantum mechanics for general systems, the resulting equations of motion will be coupled equations between different order phase-space moments. This gives at common sense a semiclassical (SC) approach if the order of phase-space moments are truncated at a certain level.

In this work, we apply this strategy to HEOM and construct a SC-HEOM formalism. The dynamical components in SC-HEOM are then different order phasespace moments of not only the primary reduced system density operator but also the auxiliary density operators (ADOs) of HEOM. We take the hierarchical quantum master equation in Ref.[23] as an example to demonstrate the methodology. Numerical illustrations are carried out on linear spectra of anharmonic oscillators with all involved calculations preformed within a minute on a single node normal computer.

Ⅱ. PRELUDE OF HEOM

In this section, we briefly introduce the background of bath statistics and the hierarchical quantum master equation in Ref.[23]. The total Hamiltonian assumes

 $\begin{eqnarray} H_{\rm{T}}(t)=H(t)+QF_{\rm{B}}(t)\ {\rm and}\ H(t)=H_{\rm{S}}-\hat\mu\epsilon(t) \end{eqnarray}$ (2.1)

Here $H(t)$ is the reduced system Hamiltonian with $H_{\rm{S}}$ the time-independent part and $\hat\mu$ the dipole operator through which the external classical field $\epsilon(t)$ acts on the system. $Q$ is a system operator representing a dissipative mode and $F_{\rm{B}}(t)$ is the stochastic Langevin force of bath. In Gaussian statistics, the bath influence on system is completely characterized by the bath correlation

 $\begin{eqnarray} C(t)\equiv\langle F_{\rm{B}}(t)F_{\rm{B}}(0)\rangle_{\rm{B}} \end{eqnarray}$ (2.2)

which is related to the bath spectral density $J(\omega)$ via the fluctuation-dissipation theorem,

 $\begin{eqnarray}\label{FDT} C(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}\!\!\textrm{d}\omega \textrm{e}^{-i\omega t} \frac{J(\omega)}{1-\textrm{e}^{-\beta\omega}} \end{eqnarray}$ (2.3)

Throughout this work we set $\hbar$=1 and $\beta$=$1/(k_\textrm{B}T)$, with $k_\textrm{B}$ the Boltzmann constant and $T$ the temperature.

In this work, the SC-HEOM approach will be exemplified on the hierarchical quantum master equation in Ref.[23], where the Drude bath model is assumed, i.e.

 $\begin{eqnarray} J(\omega)=\frac{2\lambda\gamma\omega}{\omega^2+\gamma^2} \end{eqnarray}$ (2.4)

and the bosonic function is expanded as

 $\begin{eqnarray}\label{exp_3} \frac{1}{1-\textrm{e}^{-\beta\omega}}=\frac{1}{\beta\omega}+\frac{1}{2}+\frac{\beta\omega}{12} +{\cal O}(\beta\omega)^3 \end{eqnarray}$ (2.5)

Using the three lowest expansion terms in Eq.(2.5), the approximated correlation function $C_{\rm{app}}(t)$ is obtained via Eq.(2.3) as [23]

 $\begin{eqnarray}\label{Capp} C_{\rm{app}}(t)=(c_r-ic_i)\textrm{e}^{-\gamma t}+2\Delta\delta(t) \end{eqnarray}$ (2.6)

with

 $\begin{eqnarray} c_r=2\lambda k_\textrm{B}T-\gamma\Delta, \quad c_i=\lambda\gamma, \quad \Delta=\frac{\lambda\gamma}{6k_\textrm{B}T} \end{eqnarray}$ (2.7)

Denote [23]

 $\begin{eqnarray}\label{cripara} \Gamma=\big[6+\sqrt{12+(\beta\gamma)^2}\, \big]/\beta, \quad \kappa=\sqrt{\Gamma/\Delta} \end{eqnarray}$ (2.8)

and $\Omega_{\rm{S}}$ as the characteristic frequency of system. It has been found [23] that the approximation of Eq.(2.6) would lead to accurate simulations on reduced system dynamics when

 $\begin{eqnarray}\label{criterion} \min\{\Gamma/\Omega_{\rm{S}}, \kappa\}>5 \end{eqnarray}$ (2.9)

With Eq.(2.6), the final form of the hierarchical quantum master equation is obtained as [23]

 $\begin{eqnarray}\label{hqme} \dot\rho_n(t)\hspace{-0.12cm}&=&\hspace{-0.12cm}-i[H(t), \rho_n(t)]-n\gamma\rho_n(t)-\Delta[Q, [Q, \rho_n(t)]]-\nonumber\\ &&{}\hspace{-0.12cm}i\sqrt{n/|c_r|}\, \big(c_r[Q, \rho_{n-1}(t)]-ic_i\{Q, \rho_{n-1}(t)\}\big)-\nonumber\\ &&{}\hspace{-0.12cm}i\sqrt{(n+1)|c_r|}\, [Q, \rho_{n+1}(t)] \end{eqnarray}$ (2.10)

Here, $\rho_0$ is just the reduced system density operator $\rho_0$=$\rho$=${\rm tr}_{\rm{B}}\rho_{\rm{tot}}$ and $\{\rho_n;n$=1, 2, 3$\cdots\}$ are a set of well-defined ADOs.

The hierarchical construction of the equations of motion resolves not only the system-bath coupling strength but also the memory time scale of bath correlation. Numerically the hierarchy has to be truncated at a certain tier $L$, which could be large in cases of strong coupling and long-time memory when converged dynamic results are obtained. If a $N$ basis set representation is adopted, each density operator $\rho_n$ will be described as a $N$$\times$$N$-elements matrix. The total scale of Eq.(2.10) in numerical implementation is thus $N$$\times$$N$$\times(L+1). However, Eq.(2.10) is just an exemplification of HEOM for the case of single dissipative mode Q with single exponential term of bath correlation, Eq.(2.6). In general cases, if the number of dissipative modes is M and the number of exponential terms in bath correlation expansion is K, the total number of ADOs will be \displaystyle\frac{(M\times K+L)!}{(M\times K)!L!} [7] where L is the highest tier level. Although the real number of ADOs involved in numerical evolution may be largely decreased by adopting the filtering propagator [20], the total scale of HEOM propagation can still be too huge to be affordable if the basis set dimension N is large. Such systems include, for example, general nuclear motions along anharmonic potentials meanwhile strongly influenced by thermal environments at relatively high temperatures. Since it is hard to determine in advance the characteristic representation to describe such systems being affected by environments, a large size of basis set would have to be adopted in numerical implementations. To overcome these difficulties we develop in the next section a SC-HEOM approach whose dynamic components are phase-space moments, avoiding direct propagations of density operators, thus dramatically increasing the numerical efficiency. Ⅲ. CONSTRUCTION OF THE SC-HEOM We now consider a general nuclear degree of freedom along a general anharmonic potential, under an external field action and embedded in a thermal bath. The system Hamiltonian is  \begin{eqnarray} && H_{\rm{S}}=\frac{p^2}{2M}+V(x) \\ \ \ &&\ \ {\rm with}\ \ V(x)=\sum\limits_{j=1} v_j x^j\nonumber \end{eqnarray} (3.1) The dipolar operator \hat\mu and the dissipative mode Q are assumed as  \begin{eqnarray} \hat\mu(x) = \sum\limits_{j=1} u_j x^j \end{eqnarray} (3.2)  \begin{eqnarray} Q(x) = \sum\limits_{j=1} w_j x^j \end{eqnarray} (3.3) The dynamic components in the SC-HEOM to be constructed are introduced as  \begin{eqnarray}\label{pxn} \langle x^kp^m\rangle_n\equiv{\rm tr}_{\rm{S}}(x^kp^m\rho_n) \end{eqnarray} (3.4) For clarity of the final equation, their complex conjugate counterparts related formulations would not be written explicitly in the following text. For a general dynamic operator \hat A, the time evolution of the associated \langle \hat A\rangle_n can be obtained from Eq.(2.10) as  \begin{eqnarray}\label{Ahqme} \frac{\partial}{\partial t} \langle \hat A\rangle_n \hspace{-0.12cm}&=&\hspace{-0.12cm}-i\langle[\hat A, H_{\rm{S}}]\rangle_n+i\langle[\hat A, \hat\mu]\rangle_n\epsilon(t)-\nonumber\\ & &{}\hspace{-0.12cm}n\gamma \langle \hat A\rangle_n-\Delta\langle[[\hat A, Q], Q]\rangle_n-\nonumber\\ & &{}\hspace{-0.12cm}i\sqrt{n/|c_r|}\, \big(c_r\langle[\hat A, Q]\rangle_{n-1}-ic_i\langle\{\hat A, Q\}\rangle_{n-1}\big)-\nonumber\\ & &{}\hspace{-0.12cm}i\sqrt{(n+1)|c_r|}\, \langle[\hat A, Q]\rangle_{n+1} \end{eqnarray} (3.5) The main task of this work is thus to derive the equations of motion of \{\langle x^kp^m\rangle_n\}. In numerical implementations, the orders m and k need to be truncated at k+m=N_{\rm{trun}}. This represents a type of semiclassical approach. Among those, \{\langle x^kp^m\rangle_0\}, associated with the primary reduced system density operator \rho_0=\rho, gives all interested mechanic information as far as the experimental expectation values are concerned. To derive the SC-HEOM, the following relation would be needed frequently:  \begin{eqnarray}\label{pmxk} [x^k, p^m]=-\sum\limits_{l=1}^{\min(m, k)}(-i)^lC^l_mP^l_k x^{k-l}p^{m-l} \end{eqnarray} (3.6) Here C and P represent the combination and permutation, respectively. Eq.(3.6) can be easily obtained by noticing that p=-i\displaystyle\frac{\partial}{\partial x} and  \begin{eqnarray} \frac{\partial^m}{\partial x^m}[x^kf(x)] =\sum\limits_{l=0}^{\min(m, k)}C^l_mP^l_k x^{k-l}\frac{\partial^{m-l}}{\partial x^{m-l}}f(x) \end{eqnarray} (3.7) Together with [x^kp^m, p^j]=[x^k, p^j]p^m, [x^kp^m, x^j]= -x^k[x^j, p^m], and  \begin{eqnarray} \{x^kp^m, x^j\}=2x^{k+j}p^m-x^k[x^j, p^m] \end{eqnarray} (3.8) the final form of SC-HEOM on \{\langle x^kp^m\rangle_n\} can be obtained from Eq.(3.5) as  \begin{eqnarray}\label{scheom} \frac{\partial}{\partial t}\langle x^kp^m\rangle_n\hspace{-0.12cm} &=&\hspace{-0.12cm} \frac{i}{2M}\sum\limits_{l=1}^{\min(2, k)}(-i)^lC^l_2P^l_k\langle x^{k-l}p^{2+m-l}\rangle_n-\nonumber\\ & &{}\hspace{-0.12cm}n\gamma \langle x^kp^m\rangle_n-i\sum\limits_j[v_j-u_j\epsilon(t)]\times\nonumber\\ & &{}\hspace{-0.12cm} \sum\limits_{l=1}^{\min(m, j)}(-i)^lC^l_mP^l_j\langle x^{k+j-l}p^{m-l}\rangle_n-\nonumber\\ & &{}\hspace{-0.12cm} \Delta\sum\limits_{jj'}w_jw_{j'}\sum\limits_{l=1}^{\min(m, j)}\sum\limits_{l'=1}^{\min(m-l, j')}(-i)^{l+l'}\times\nonumber\\ & &{}\hspace{-0.12cm} C^l_mP^l_jC^{l'}_{m-l}P^{l'}_{j'}\langle x^{k+j+j'-l-l'}p^{m-l-l'}\rangle_n-\nonumber\\ & &{}\hspace{-0.12cm}i\sqrt{(n+1)|c_r|}\, \sum\limits_jw_j\times\nonumber\\ & &{}\hspace{-0.12cm} \sum\limits_{l=1}^{\min(m, j)}(-i)^lC^l_mP^l_j\langle x^{k+j-l}p^{m-l}\rangle_{n+1}-\nonumber\\ & &{}\hspace{-0.12cm} i\sqrt{n/|c_r|}\, (c_r-ic_i)\sum\limits_jw_j\times\nonumber\\ & &{}\hspace{-0.12cm} \sum\limits_{l=1}^{\min(m, j)}(-i)^lC^l_mP^l_j\langle x^{k+j-l}p^{m-l}\rangle_{n-1}-\nonumber\\ & &{}\hspace{-0.12cm}2c_i\sqrt{n/|c_r|}\, \sum\limits_jw_j\langle x^{k+j}p^{m}\rangle_{n-1} \end{eqnarray} (3.9) The total dimension of this equation is (N_{\rm{trun}}+ 1)(N_{\rm{trun}}+2)(L+1)/2 (see comments after Eq.(3.5)). Also note that the complex conjugate counterparts of Eq.(3.9) can be readily obtained. They together with Eq.(3.9) will straightforwardly lead to the equations of motion of \{\langle x^kp^m+p^mx^k\rangle_n, \langle i(x^kp^m$$-$$p^mx^k)\rangle_n\} that are real and of more convenience in numerical implementation. Ⅳ. LINEAR RESPONSE AND NUMERICAL DEMONSTRATION For application, we consider the linear spectrum of an anharmonic oscillator, evaluated via the linear response approach on Eq.(3.9). To the first order of the external field, the linear response relation is written as  \begin{eqnarray} \delta P^{(1)}(t)\hspace{-0.12cm} &=&\hspace{-0.12cm} {\rm tr}_{\rm{S}}{\rm tr}_{\rm{B}}\{\hat\mu[\rho^{(1)}_{\rm{tot}}(t)-\rho_{\rm{tot}}^{\rm eq}]\}\nonumber\\ \hspace{-0.12cm} &=&\hspace{-0.12cm} \langle\hat\mu\rangle^{(1)}_0-\langle\hat\mu\rangle_{\rm eq}\nonumber\\ \hspace{-0.12cm}&\equiv&\hspace{-0.12cm} \int_{-\infty}^t\!\!\textrm{d}\tau\, \chi(t-\tau)\epsilon(\tau)\quad\quad \end{eqnarray} (4.1) Here, \chi(t) is the linear response function and its Fourier transform gives the linear spectrum. To derive the perturbative solution of Eq.(3.9) with respect to the external field, we recast Eq.(3.9) in the matrix-vector form as  \begin{eqnarray}\label{mateom} \dot{\mathbf X}(t)=[{\mathbf\Lambda}+{\mathbf D}\epsilon(t)]{\mathbf X}(t) \end{eqnarray} (4.2) with {\mathbf X}$$\equiv$$\{\langle x^kp^m\rangle_n\} comprised as a vector. {\mathbf\Lambda} and {\mathbf D}, as matrices, correspond to the field-independent and field-dependent counterparts of Eq.(3.9), respectively. Introducing  \begin{eqnarray} {\mathbf X}(t)=\textrm{e}^{{\mathbf\Lambda}t}{\mathbf X}_I(t) \end{eqnarray} (4.3) we have  \begin{eqnarray} \dot{\mathbf X}_I(t)=\textrm{e}^{-{\mathbf\Lambda}t}{\mathbf D}\textrm{e}^{{\mathbf\Lambda}t}{\mathbf X}_I(t)\epsilon(t) \end{eqnarray} (4.4) which leads to  \begin{eqnarray}\label{linearresp} \delta{\mathbf X}^{(1)}(t) = \int_{-\infty}^t\textrm{d}\tau \textrm{e}^{{\mathbf\Lambda}(t-\tau)}{\mathbf D}{\mathbf X}_{\rm eq}\epsilon(\tau)\nonumber\\ \equiv \int_{-\infty}^t\textrm{d}\tau{\chi}(t-\tau)\epsilon(\tau) \end{eqnarray} (4.5) {\mathbf X}_{\rm eq} is the steady state solution of Eq.(3.9) in the absence of external field, satisfying \dot{\mathbf X}_{\rm eq}={\mathbf\Lambda}{\mathbf X}_{\rm eq}=0 with the constraint \langle x^0p^0\rangle_0={\rm tr}_{\rm{S}}(\rho)=1. Eq.(4.5) is thus the linear response solution of Eq.(3.9) or Eq.(4.2) to the first order of the external field \epsilon(t). For numerical demonstration, the potential of the anharmonic oscillator is assumed as  \begin{eqnarray} V(x)=V_\textrm{e}\Big[\Big(\frac{x}{a}\Big)^2-\theta\Big(\frac{x}{a}\Big)^3\Big] \end{eqnarray} (4.6) The parameters are chosen as M=1 amu, a=2 a.u., V_\textrm{e}=2500{\rm{\ cm}}^{-1}, and \theta=0.3 or 0.15 for different anharmonicities. The characteristic frequency is about \Omega_{\rm{S}}=388{\rm{\ cm}}^{-1}. The dissipative mode is chosen as Q=2(x/a) and the dipole is \hat\mu$$\propto$$x. Bath parameters are selected as \gamma=200{\rm{\ cm}}^{-1} and \lambda=150{\rm{\ cm}}^{-1} with the temperature T=300 K. The criterion parameters of Eq.(2.8) are obtained as \Gamma=2000{\rm{\ cm}}^{-1} and \kappa=9. Thus Eq.(2.9) is satisfied. The \delta\langle x\rangle^{(1)}_0-related element of {\chi}$$(t)$ in Eq.(4.5) is evaluated. Its Fourier transform gives the linear spectra of the system, which are demonstrated in FIG. 1 for different $\theta$ and $N_{\rm{trun}}$. The inset is the potentials associated with $\theta$=0.3 (thick) and 0.15 (thin), respectively. Calculations are preformed within a minute using a normal computer with the HEOM tier truncated at $L$=8. We can see that with the anharmonicity increasing, the spectral peaks are red shifted and the effects of higher-order phase-space moments become more important.

 FIG. 1 Linear spectra of anharmonic oscillators evaluated with different phase-space truncation orders: $N_{\rm{trun}}$=6 (solid), $N_{\rm{trun}}$=4 (dotted), $N_{\rm{trun}}$=2 (dashed), and different anharmonicities associated with $\theta$=0.3 (thick) and $\theta$=0.15 (thin), respectively. The inset is the corresponding potentials. See the main text for other numerical details.
Ⅴ. SUMMARY

In this work, we propose a SC-HEOM approach, which is constructed on the basis of the well-established HEOM formalism, with its dynamic components being phase-space moments instead of the density matrices. The objective is to achieve the numerical executability to simulate non-perturbative and non-Markovian quantum dissipative dynamics of anharmonic nuclear systems strongly influenced by a thermal bath at a relatively high temperature. Such systems would require practically a rather large basis set in representation leading to the density-matrix based methods pretty difficult to perform.

For clarity of the SC-HEOM formalism and its construction, the present work is exemplified with the hierarchical quantum master equation of Eq.(2.10) for the case of single dissipative mode with single exponential term bath correlation. The numerical demonstration is made on the linear response of a nuclear oscillator with anharmonicity. More generalized SC-HEOM construction and more complex applications will be carried out in future work. Also the truncation on higher-order phase-space moments need further elaboration considering related physical inherent properties.

Ⅵ. ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China (No.21373191, No.21573202, No.21633006, and No.21703225), and the Fundamental Research Funds for the Central Universities (No.2030020028, No.2060030025, and No.2340000074).

Reference
 [1] R. P. Feynman, and F. L. Vernon Jr, Ann. Phys. 24 , 118 (1963). DOI:10.1016/0003-4916(63)90068-X [2] Y. Tanimura, Phys. Rev. A 41 , 6676 (1990). DOI:10.1103/PhysRevA.41.6676 [3] Y. Tanimura, J. Phys. Soc. Jpn. 75 , 082001 (2006). DOI:10.1143/JPSJ.75.082001 [4] Y. A. Yan, F. Yang, Y. Liu, and J. S. Shao, Chem. Phys. Lett. 395 , 216 (2004). DOI:10.1016/j.cplett.2004.07.036 [5] R. X. Xu, P. Cui, X. Q. Li, Y. Mo, and Y. J. Yan, J. Chem. Phys. 122 , 041103 (2005). DOI:10.1063/1.1850899 [6] J. S. Jin, X. Zheng, and Y. J. Yan, J. Chem. Phys. 128 , 234703 (2008). DOI:10.1063/1.2938087 [7] J. J. Ding, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 136 , 224103 (2012). DOI:10.1063/1.4724193 [8] J. S. Shao, J. Chem. Phys. 120 , 5053 (2004). DOI:10.1063/1.1647528 [9] J. M. Moix, and J. S. Cao, J. Chem. Phys. 139 , 134106 (2013). DOI:10.1063/1.4822043 [10] Y. L. Ke, and Y. Zhao, J. Chem. Phys. 145 , 024101 (2016). DOI:10.1063/1.4955107 [11] C. Y. Hsieh, and J. Cao, J. Chem. Phys. 148 , 014103 (2018). DOI:10.1063/1.5018725 [12] C. Y. Hsieh, and J. Cao, J. Chem. Phys. 148 , 014104 (2018). DOI:10.1063/1.5018726 [13] Y. J. Yan, J. Chem. Phys. 140 , 054105 (2014). DOI:10.1063/1.4863379 [14] H. D. Zhang, R. X. Xu, X. Zheng, and Y. J. Yan, J. Chem. Phys. 142 , 024112 (2015). DOI:10.1063/1.4905494 [15] R. X. Xu, H. D. Zhang, X. Zheng, and Y. J. Yan, Sci. China Chem. 58 , 1816 (2015). DOI:10.1007/s11426-015-5499-2 [16] H. D. Zhang, Q. Qiao, R. X. Xu, and Y. J. Yan, Chem. Phys. 481 , 237 (2016). DOI:10.1016/j.chemphys.2016.07.005 [17] R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, Chin. J. Chem. Phys. 30 , 395 (2017). DOI:10.1063/1674-0068/30/cjcp1706123 [18] R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, J. Chem. Phys. 148 , 114103 (2018). DOI:10.1063/1.4991779 [19] Y. Liu, R. X. Xu, H. D. Zhang, and Y. J. Yan, Chin. J. Chem. Phys. 31 , 245 (2018). DOI:10.1063/1674-0068/31/cjcp1804083 [20] Q. Shi, L. P. Chen, G. J. Nan, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 130 , 084105 (2009). DOI:10.1063/1.3077918 [21] J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 133 , 101106 (2010). DOI:10.1063/1.3484491 [22] J. Hu, M. Luo, F. Jiang, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 134 , 244106 (2011). DOI:10.1063/1.3602466 [23] R. X. Xu, B. L. Tian, J. Xu, Q. Shi, and Y. J. Yan, J. Chem. Phys. 131 , 214111 (2009). DOI:10.1063/1.3268922