Chinese Journal of Chemical Physics  2017, Vol. 30 Issue (6): 800-810

#### The article information

Jie Zheng, Yu Xie, Sheng-shi Jiang, Yun-ze Long, Xin Ning, Zheng-gang Lan

Ultrafast Electron Transfer with Symmetrical Quasi-classical Dynamics based on Mapping Hamiltonian and Quantum Dynamics based on ML-MCTDH

Chinese Journal of Chemical Physics, 2017, 30(6): 800-810

http://dx.doi.org/10.1063/1674-0068/30/cjcp1711210

### Article history

Accepted on: December 27, 2017
Ultrafast Electron Transfer with Symmetrical Quasi-classical Dynamics based on Mapping Hamiltonian and Quantum Dynamics based on ML-MCTDH
Jie Zhenga,b, Yu Xiec,d, Sheng-shi Jiangc,d, Yun-ze Longa,b, Xin Ninga, Zheng-gang Lanc,d
Dated: Received on November 14, 2017; Accepted on December 27, 2017
a. Industrial Research Institute of Nonwovens & Technical Textiles, College of Textiles Clothing, Qingdao University, Qingdao 266071, China;
b. College of Physics, Qingdao University, Qingdao 266071, China;
c. Key Laboratory of Biobased Materials, Qingdao Institute of Bioenergy and Bioprocess Technology, Chinese Academy of Sciences, Qingdao 266101, China;
d. University of Chinese Academy of Sciences, Beijing 100049, China
*Author to whom correspondence should be addressed. Yu Xie, E-mail:xieyu@qibebt.ac.cn; Zheng-gang Lan, E-mail:lanzg@qibebt.ac.cn
Abstract: Symmetrical quasi-classical (SQC) method based on mapping Hamiltonian is an efficient approach that is potentially useful to treat the nonadiabatic dynamics of very large systems. We try to evaluate the performance of this method in the ultrafast electron transfer processes involving a few of electronic states and a large number of vibrational modes. The multilayer multiconfigurational time-dependent Hartree (ML-MCTDH) method was used to get the accurate dynamical results for benchmark. Although the population dynamics in the longtime limit show differences in the ML-MCTDH and SQC calculations, the SQC method gives acceptable results.
Key words: Mapping Hamiltonian    Semiclassical dynamics    Nonadiabatic dynamics
Ⅰ. INTRODUCTION

The simulation of the nonadiabatic dynamics of complex systems is very challenging due to the breakdown of the Born-Oppenheimer (BO) approximation and the involvement of a large number of strongly-coupled nuclear/electronic freedoms [1, 2]. Although the theoretical approaches based on quantum dynamics, such as full quantum wavepacket dynamics [1], the multiconfiguration time-dependent Hartree (MCTDH) method [3] and multi-layer MCTDH (ML-MCTDH) method [4-7], can provide the accurate solution of non-BO dynamics, they suffer from the high computational cost. The Gaussian-wavepacket based methods, such as ab initio multiple spawning (AIMS) [8], Gaussian-based MCTDH (G-MCTDH) [9, 10], are promising methods for poly-atomic molecules, while the employment of them in the large systems still requires the significant amount of computational cost. It is thus necessary to develop the semiclassical dynamics for the effective treatment of nonadiabatic dynamics of large systems. Along this research line, several practical theoretical approaches were developed, such as mean-field (Ehrenfest) dynamics and its extensions [1, 2, 11], surface-hopping dynamics [1, 2, 12-14], and so on. After the implementation of these efficient methods in the manner of on-the-fly approaches [15-24], it is possible to simulate the nonadiabatic dynamics of realistic polyatomic systems with the inclusion of all nuclear degrees of freedom. However, we cannot neglect the limitation of these practical methods. For example, the widely-used fewest-switches trajectory surface-hopping method developed by Tully [12] suffers from the improper treatment of the electronic coherence and the appearance of frustrated hops [1, 2, 13, 25, 26].

Alternative approach for the semiclassical approximation of nonadiabatic dynamics was proposed by Meyer and Miller [27, 28], and later on by Stock and Thoss [29]. Analogy of Schwinger transformation, the MMST (Meyer-Miller-Stock-Thoss) model maps $N$ discrete quantum states to $N$ coupled harmonic oscillator. Liu [30] pointed out that such mapping way is not uniquely defined. Previous studies [2, 31] prove that the combination of this mapping Hamiltonian and rigorous semiclassical treatments (such as initial value representation) may provide the accurate description of nonadiabatic dynamics, while the large amount of computational cost limits its wide application of complicated systems. Alternatively, if replacing the quantum operator by their classical counterpart in the mapping Hamiltonian, we obtained a fully classical Hamiltonian. This opens the possibility to run the quasi-classical (QC) dynamics after the proper initial sampling [2]. Overall, the mapping idea received considerable attention and different variants were proposed [32, 33].

In the classical dynamics simulation, the assignment of the final quantum state of the reaction product may be realized by the idea of energy window [34]. For example, if the coordinate and momentum of a harmonic vibrational mode are obtained, it is possible to get its total classical energy. When such energy falls into an energy window belonging to a particular quantum state, we may assign the corresponding quantum state of this vibrational mode. Cotton and Miller [35, 36] suggested that the same window treatment should be considered in the initial sampling step when it is used in the final assignment, resulting in the so-called symmetrical quasi-classical trajectory (SQC) method. Because of its simple implementation and its good performance in a few of model systems [37-40], this mapping-Hamiltonian-based SQC dynamics is assumed to be useful to treat the nonadiabatic dynamics of complicated systems. Due to its efficiency, several studies try to extend the ability of this method for the treatment of different situations, for examples see the work by Tao [41-43].

Through above discussions, it is clearly known that many theoretical methods are available in the study of nonadiabatic dynamics. Particularly, various rigorous semiclassical methods perform extremely well in some model systems, which may give fully correct results consistent with quantum dynamics. However, in reality we need to treat many realistic problems with high complexity. For example, if we wish to treat the nonadiabatic dynamics of extremely large molecular aggregates, a huge number of electronic states and nuclear degrees of freedom are included. Alternatively, if we wish to perform the on-the-fly calculations of polyatomic systems (solvated or biological systems), the computational cost is also very large. The employment of these rigorous semiclassical methods is still beyond the current computational facility, even if these theoretical frameworks work well in model systems. Thus, it is very important to develop some cheaper theoretical methods, which allow us to run the simulation of nonadiabatic dynamics with the balance of computational efficiency and accuracy. Although surface hopping method is such a type of methods and it shows great successes in many applications [1, 2, 12, 13, 25, 26], it suffers from some deficiencies, such as the overestimation of electronic coherence, frustrated hops, and so on. Therefore, it is still essential to "look for" other possible semi-classical or quasiclassical approaches to treat realistic complex systems. Previous studies [36, 39] showed the very good performance of SQC method, even at low temperature. It is reasonable to assume that the SQC may be a possible approach to treat the nonadiabatic dynamics of realistic complex systems. Therefore, extensive benchmark calculations should be done to check its performance on more situations. However, the examination of the accuracy of SQC can only be performed by running calculations in some model systems that can be solved by accurate methods, because it is not possible to get fully numerical correct results for extremely complex system. This is the purpose of this work.

In this work, we try to evaluate the performance of the SQC dynamics with the mapping Hamiltonian in the ultrafast electron transfer processes involving a few of electronic states and a large number of vibrational modes. In our previous work [44], the electron transfer (ET) dynamics from the acene (anthracene) to C$_{60}$ was discussed by using ML-MCTDH based on the diabatic Hamiltonian constructed from electronic structure calculations. In the current work, we wish to take this spin-boson model to test the performance of QC and SQC, because the accurate dynamics results can be obtained by ML-MCTDH. Our purpose here is to check whether QC and SQC can predict the general feature of population dynamics. More importantly, we wish to know whether this highly approximated SQC method can still give the qualitatively and semi-quantitatively reasonable results or not. We wish to know how large the deviation may be generated by the QC and SQC methods. This may provide the valuable information on the further employment of these methods. In addition, because the current model also represents a typical ET model in organic photovoltaics, this work provided the important information on how far we can trust the results of the QC and SQC treatment on similar problems. We believe that this work should help us understand the performance of mapping Hamiltonian and SQC dynamics. This information will be extremely useful in the future study of nonadiabatic dynamics of realistic complicated systems.

Ⅱ. THEORY AND METHODS A. Hamiltonian

The diabatic model in our previous work [44] was taken in the current calculations, which describes the ultrafast acene (Ac)$\rightarrow$fullerene (C$_{60}$) electron transfer (ET) process. The model Hamiltonian includes a donor state (Ac*-C$_{60}$, denoted as $|\psi_d\rangle$) and three acceptor states (Ac$^+$-C$_{60}$$^-, denoted as |\psi_a\rangle). The Hamiltonian in the diabatic representation reads as follows:  \begin{eqnarray} H &=& {T_{\rm{nuc}}} + \left| {{\psi _d}} \right\rangle {V_{dd}}\left\langle {{\psi _d}} \right| + \sum\limits_a {\left| {{\psi _a}} \right\rangle {V_{aa}}\left\langle {{\psi _a}} \right|} + \nonumber\\ &&{}\sum\limits_a {\left( {\left| {{\psi _d}} \right\rangle {V_{da}}\left\langle {{\psi _a}} \right| + \left| {{\psi _a}} \right\rangle {V_{ad}}\left\langle {{\psi _d}} \right|} \right)} \end{eqnarray} (1) where T_{\rm{nuc}} is the kinetic energy. V_{dd} and V_{aa} are the energies of donor and acceptor states, respectively. V_{da}(V_{ad}) are the couplings between donor and acceptor states.  \begin{eqnarray} \begin{array}{l} {V_{dd}} = {E_d} + \displaystyle\frac{1}{2}\sum\limits_i {{\omega _i}Q_i^2} + \sum\limits_i {\kappa _i^{\left( d \right)}{Q_i}} \\ {V_{aa}} = {E_a} + \displaystyle\frac{1}{2}\sum\limits_i {{\omega _i}Q_i^2} + \sum\limits_i {\kappa _i^{\left( a \right)}{Q_i}} \end{array} \end{eqnarray} (2) where E_d and E_a are energies of the donor and acceptor electronic states at the ground-state minimum geometry. Q_i are the dimensionless normal coordinates with associated frequencies \omega_i, and \kappa_i^{(n)} (n=d, a) is the linear vibronic couplings. Here the summation is over all nuclear coordinates. The D-A couplings V_{da}(V_{ad}) take their zero-order couplings. In our previous work [44], all parameters in this four-state model including one donor state (D) and three acceptor states (A1, A2 and A3) were obtained through electronic structure calculations. For illustration, the energies of diabatic states and their couplings between states were given in FIG. 1. For the electron-phonon couplings, we have already discussed all details in previous work [44] and here we only give a short description in Appendix.  FIG. 1 Schematic diagram of the four-state Hamiltonian model used in the current calculations, which includes one donor state and three acceptor states, denoted as D, A1, A2, and A3. The electronic energy of each state and the diabatic couplings between donor state and three acceptor states were given. B. Quasi-classical dynamics method based on mapping Hamiltonian 1. Mapping Hamiltonian The Hamiltonian of a N-state system \hat H is expressed as  \begin{eqnarray} \hat H = \sum\limits_{k, l} {{{\hat h}_{kl}}} \left| {{\phi _k}} \right\rangle \left\langle {{\phi _l}} \right| \end{eqnarray} (3) where \phi_k (or \phi_l) is the k-th (or l-th) electronic state of the system. Because we discuss the system with several nuclear degrees of freedom, the diagonal element of \hat h_kl (k=l) represents the energy (kinetic and potential energies) of an electronic state, which is a function of nuclear coordinate and momentum operators. The off-diagonal element \hat h_kl (k$$\neq$$l) denotes the interstate couplings between the k-th and l-th states. This Hamiltonian is represented in the discrete electronic basis (\phi_k). The main idea of the mapping approach is to find a way to map the current Hamiltonian from the discrete representation to the continuous representation. In the MMST approaches [1, 2, 26-28], a N-state system is mapped to a N coupled harmonic oscillators. The mapping relation is construed based on the annihilation operator \hat a_k and creation operator \hat a_l^+ of the k-th harmonic oscillator, namely  \begin{eqnarray} \begin{array}{l} \left| {{\phi _k}} \right\rangle \left\langle {{\phi _l}} \right| \mapsto \hat a_k^ + {{\hat a}_l}\\ \left| {{\phi _k}} \right\rangle \mapsto \left| {{0_1} \cdot \cdot \cdot {1_k} \cdot \cdot \cdot {0_N}} \right\rangle \end{array} \end{eqnarray} (4) The commutation relation of the annihilation operator and creation operator is given as \left[\hat a_k, \hat a_l^+\right]=\delta_{kl}. By introducing the Cartesian coordinate operator \hat x_k=\left(\hat a_k+\hat a_l^+\right)/\sqrt{2} and momentum operator \hat p_k=i(\hat a_k$$-$$\hat a_l^+)/\sqrt{2} of the harmonic oscillator, the MMST Hamiltonian is obtained, i.e.,  \begin{eqnarray} \hat H = \sum\limits_k {\frac{1}{2}\left( {{{\hat x}_k}^2 + {{\hat p}_k}^2 - 1} \right)} {\hat h_{kk}} + \frac{1}{2}\sum\limits_{k \ne l} {\left( {{{\hat x}_k}{{\hat x}_l} + {{\hat p}_k}{{\hat p}_l}} \right){{\hat h}_{kl}}}\nonumber\\ \end{eqnarray} (5) In this mapping Hamiltonian, the several coupled harmonic oscillators with their Cartesian coordinate operators \hat x_k and momentum operators \hat p_k are used to represent coupled electronic states. For simplicity, when the term of "the coordinate and momentum of electronic part" is used in the below discussion, we refer to the physical quantities defined by the above equation (Eq.(5)), instead of the real coordinates and momenta of electrons. 2. Quasi-classical dynamics If substituting the quantum operators in the MMST Hamiltonian with their corresponding physical variables, a pure classical mapping Hamiltonian is obtained, namely  \begin{eqnarray} {H_{MM}}\left( {x, p, Q, P} \right) &=& \sum\limits_k {\frac{1}{2}\left( {{x_k}^2 + {p_k}^2 - 1} \right)} {H_{kk}}\left( {Q, P} \right) + \nonumber\\ &&{}\frac{1}{2}\sum\limits_{k \ne l} {\left( {{x_k}{x_l} + {p_k}{p_l}} \right){H_{kl}}} \left( {Q, P} \right) \end{eqnarray} (6) where x_k is the classical coordinates and p_k is the classical momentums for the electronic part. Q and P are the coordinates and momentums of nuclei. H_{kk} and H_{kl} (k$$\neq$l), as functions of $Q$ and $P$, are the classical correspondences of the quantum Hamiltonian matrix elements. Starting from this classical mapping Hamiltonian, it is possible to perform the proper initial sampling and run classical dynamics over many trajectories. This provides us a QC trajectory-based approach to treat nonadiabatic dynamics. For each trajectory, the quantum occupation number $n_k$ is represented by $1/2({x_k}^2+{p_k}^2-1)$. After the evolution of many classical trajectories, the electronic population $P^{\rm{pop}}_{{\rm map}, k}$ of the $k$-th state can be obtained by averaging quantum occupation number $n_k$ over all the trajectories

 $\begin{eqnarray} P^{\rm{pop}}_{{\rm map}, k} = {\left\langle {\frac{1}{2}\left( {{x_k}^2 + {p_k}^2 - 1} \right)} \right\rangle _{\rm{traj}}} = {\left\langle {{n_k}} \right\rangle _{\rm{traj}}} \end{eqnarray}$ (7)

where the bracket $\langle\ldots\rangle_{\rm{traj}}$ refers that the average is taken over all trajectories.

The mapping approach naturally provides a way to construct the initial phase-space distribution for the electronic part. The coordinate and momentum are sampled by the action-angle sampling method [2] with fixed electronic action $n_k$=0 (unoccupied) or 1 (occupied).

 $\begin{eqnarray} \begin{array}{l} {x_k} = \sqrt {2{n_k} + 1} \cos \theta \\ {p_k} = \sqrt {2{n_k} + 1} \sin \theta \end{array} \end{eqnarray}$ (8)

where the angle $\theta$ is a random number in the range of [$-\pi$, $\pi$).

For the nuclear part, the current approach does not define a unique way to construct the initial classical distribution function in the phase space. Thus, in principle we may use the action-angle sampling or Wigner sampling [45, 46]. Several previous studiess replying on the mapping approach took the Wigner sampling for the nuclear motion [37, 47, 48] and this way is also widely used in many other semiclassical dynamics [15-25, 49], we thus employed this approach.

Starting from each initial condition, the classical dynamics on the basis of mapping Hamiltonian are calculated. After the dynamics evolution, for a single trajectory the coordinate $x_k$ and momentum $p_k$ for the $k$-th electronic states should be obtained and the corresponding electronic population is given by $n_k$=$1/2\left({x_k}^2+{p_k}^2-1\right)$. Then the electronic population $P^{\rm{pop}}_{{\rm map}, k}$ of the k-th state can be obtained by averaging quantum occupation number $n_k$ of mapping electronic degrees of freedom over all trajectories, namely $P^{\rm{pop}}_{{\rm map}, k}$=$\langle n_k\rangle_{\rm{traj}}$.

3. The symmetrical quasi-classical method

In the quasi-classical dynamics based on the mapping Hamiltonian, it is also possible to obtain the distribution of final quantum states (occupation) by calculating the final values of the action variables accumulated in "bins" (or histograms) centered at 0 or 1. In recent years, Cotton and Miller [35] proposed that such the "bin" idea should be employed for both initial-state sampling and final-state assignment. In this symmetrical quasi-classical (SQC) method, the classical action variables $n_k$ are "quantized" symmetrically (i.e., initially and finally) by defining "window functions" $W_k$,

 $\begin{eqnarray} {W_k}\left( {{n_k}, N} \right) = \frac{1}{{\Delta n}}h\left( {\frac{{\Delta n}}{2} - \left| {{n_k} - N} \right|} \right) \end{eqnarray}$ (9)

where $h(z)$ is the Heaviside function, i.e. $h(z)$=0 when $z$$<0 and h(z)=1 when z$$\geq$0. $\Delta n$ is the window width. The window is centered at the quantum values $N$=1 (occupied) or 0 (unoccupied), i.e., $n_k$ is required to be in the interval [$N$$-$$\Delta n/2$, $N$+$\Delta n/2$]. In this method, $n_k$=1/2(${x_k}^2$+${p_k}^2$$-$$\gamma)$, $\gamma$ is a parameter accounting for the effective zero-point energy [47, 48] and $\gamma$=$\Delta n$. As defined in the last subsection, $x_k$ are the classical coordinates and $p_k$ the classical momentums for the electronic part.

In the SQC approach, the initial sampling of the phase space for the electronic part is performed with the action-angle sampling by setting the action within a window with the width of $\gamma$ and the center at 0 or 1. For the final states, the same binning way is also performed for the assignment of the quantum state of the product. The time-dependent occupations $P^{\rm{occ}}_{{\rm map}, k}$ of electronic states are evaluated by the below average over all trajectories

 $P^{\rm{occ}}_{{\rm map}, k}( t ) = \frac{A}{B}\\ A=\bigg\langle\bigg[W_i( n_i( 0 ), 1 )\prod\limits_{l \ne i}W_l(n_l( 0 ), 0 ) \bigg]\cdot\nonumber\\ \;\;\;\;\;\;\bigg[W_k( n_k( t ), 1)\prod\limits_{l \ne k}W_l(n_l( t ), 0 ) \bigg] \bigg\rangle _{\rm{traj}}\nonumber\\ B=\sum\limits_m^{N_s}\bigg\langle\bigg[W_i( n_i( 0 ), 1 )\prod\limits_{l \ne i}W_l( n_l( 0 ), 0 ) \bigg]\cdot\nonumber\\ \;\;\;\;\;\;\bigg[W_m(n_m( t ), 1 )\prod\limits_{l \ne m}W_l(n_l( t ), 0 ) \bigg] \bigg\rangle _{\rm{traj}}\nonumber$ (10)

in which the average is over all trajectories and the summation of the index $m$ is over all electronic states.

Next, we explain the numerical details of the SQC methods. At time 0, we assume that only one electronic state (state $i$) is occupied and all other states are unoccupied. For the occupied state $i$, $n_i$(0) is in the interval [$1$$-$$\gamma/2$, $1$+$\gamma$/2] and $W_l(n_l(0), 1)$=1/$\gamma$ after the proper initial sampling. For any unoccupied state $l$, $n_l(0)$ is in the interval [0$-$$\gamma/2, 0+\gamma/2] and W_i(n_i(0), 9)=1/\gamma. Thus, \left[{{W_i}\left({{n_i}\left(0 \right), 1} \right)\prod\limits_{l \ne i} {{W_l}\left({{n_l}\left(0 \right), 0} \right)} } \right]=\left({1/\gamma } \right)^{{N_s}}, where N_s is the number of states. At time t, the electronic population n_k(t) of the k-th state can be calculated by {n_k}=1/2\left({{x_k}^2 + {p_k}^2 -\gamma } \right) for a single trajectory. The above relationship can be used again to judge the quantum state by assigning the occupation value according to window functions. When the state k is occupied and the other states are all unoccupied, the product of \left[{{W_k}\left({{n_k}\left(t \right), 1} \right)\prod\limits_{l \ne k} {{W_l}\left({{n_l}\left(t \right), 0} \right)} } \right]={\left({1/\gamma } \right)^{{N_s}}}. In all the other cases, this product equals zero. If two or more states are occupied or no state is occupied at time t, the product of window functions also equals zero. Then the time-dependent occupation P^{\rm{occ}}_{{\rm map}, k} of the k-th electronic state is calculated according to Eq.(10). Following previous discussions [35, 47], the SQC dynamics with \gamma=1 and the SQC dynamics with \gamma=0.732 were considered. For a clear illustration in this work, we always use P^{\rm{pop}}_{{\rm map}, k} and P^{\rm{occ}}_{{\rm map}, k} to represent of the electronic population (Eq.(7)) and state occupation (Eq.(10)) of the k-th state. For the QC dynamics with the mapping Hamiltonian, we label it as "non-window QC dynamics". When the window is used, we always use the SQC approach to represent this situation. The convergence of SQC calculations and the dependence of results on the number of trajectories are discussed in Appendix. C. ML-MCTDH In traditional quantum wave packet dynamics, a group of time-independent basis for each primary coordinate is used to represent the wave function, and the time-dependent coefficients evolve with time. MCTDH employs the time-dependent bases to represent wave function, i.e.  \Psi \left( {{Q_1}, \cdots, {Q_f}, t} \right) =\nonumber\\ \;\;\;\;\;\;\;\;\sum\limits_{{j_1} = 1}^{{n_1}} { \cdots \sum\limits_{{j_f} = 1}^{{n_f}} {{A_{{j_1} \cdots {j_f}}}\left( t \right)} } \prod\limits_{\kappa = 1}^f {\varphi _{{j_\kappa }}^{\left( \kappa \right)}} \left( {{Q_\kappa }, t} \right) (11) where f is the number of degrees of freedom. {{A_{{j_1} \cdots {j_f}}}} is time-dependent coefficient and {\varphi _{{j_\kappa }}^{\left(\kappa \right)}} is the time-dependent basis function namely single particle function (SPF). Q_k refer to the primary coordinates. Using the variational principle, the coefficient and SPFs time evolution of the equations of motion can be obtained. ML-MCTDH is an extension of the standard MCTDH method. In ML-MCTDH, the wave function is represented by a multilayer tree structure (shown in FIG. 2), which is expanded by the multi-dimensional SPFs recursively until the time-independent bases of each primary coordinate are reached. In other words, an SPF in the upper layer is an expansion of SPFs in the lower layer; i.e.,  \begin{eqnarray} &&\varphi _m^{l - 1;{\kappa _1} \cdots {\kappa _{l - 1}}}\left( {Q_{{\kappa _{l - 1}}}^{l - 1;{\kappa _1} \cdots {\kappa _{l - 2}}}, t} \right) =\sum\limits_{{j_1} = 1}^{{n_1}}\cdots \nonumber\\ &&\sum\limits_{{j_{{p_{{\kappa _l}}}}} = 1}^{{n_{{\kappa _l}}}} {A_{m;{j_1} \cdots {j_{{p_{{\kappa _l}}}}}}^{l;{\kappa _1} \cdots {\kappa _{l - 1}}}\left( t \right)}\prod\limits_{{\kappa _l} = 1}^{{p_{{\kappa _l}}}} {\varphi _{{j_{{\kappa _l}}}}^{l;{\kappa _1} \cdots {\kappa _l}}} \left( {Q_{{\kappa _l}}^{l;{\kappa _1} \cdots {\kappa _{l - 1}}}, t} \right)\quad\quad \end{eqnarray} (12)  FIG. 2 Schematic tree structure for the ML-MCTDH wavefunctions. where l denotes the layer depth, and \kappa_1, \ldots, \kappa_{l-1} denote the indices of the logical degrees of freedom of the upper layer down to a lower layer. \varphi_{j_{\kappa+l}}^{l; \kappa_1\ldots\kappa_l} is the j_{\kappa_l}-th SPF for each degree of freedom \kappa_l on the layer l and Q_{{\kappa _l}}^{l; {\kappa _1} \cdots {\kappa _{l -1}}} are collective coordinates of the SPF. Then the equation of motions (EOMs) of ML-MCTDH can be resolved by applying variational principle and the recursive algorithm [4-7]. The construction of the reasonable tree expansion of a wave packet is not a trivial task, which directly affects the accuracy and efficiency of calculation. Because this topic involves deep theoretical discussions beyond the current work, we refer readers to check previous studies [4-7, 50-53]. We tried to run ML-MCTDH calculations to get the accurate dynamical results for benchmark. The ML-MCTDH calculations were performed using the Heidelberg MCTDH package [54]. For convenience, we used P^{\rm{pop}}_{{\rm ml}, k} to represent electronic population of the k-the state in ML-MCTDH. Previous work [44] already provided a proper tree expansion with the reliable description of ultrafast dynamics to 120 fs. Because we wish to know the long-time dynamics behavior of the QC and SQC dynamics, all ML-MCTDH dynamics are re-calculated up to a longer time (150 fs). Because the test calculations show that the population dynamics becomes stable after 120 fs, the same tree expansions employed in previous work were taken in the current study. Ⅲ. RESULTS AND DISCUSSION A. Dynamics of two-state models We first considered a few of two-state models including a donor (D) and an acceptor (A1) state for the evaluation of the performance of the mapping approaches. Various models with different numbers (7, 11, 16, 44 and 246) of vibrational degrees of freedom were considered. The time-dependent electronic populations were calculated using ML-MCTDH and different mapping approaches. Because previous studies recommended the symmetrical window, we do not consider the calculations of state occupation P^{\rm{occ}}_{{\rm map}, k} in the no-window quasi-classical dynamics with mapping Hamiltonian. Only in the SQC dynamics approach, we computed the state occupation using Eq.(10). All results are shown in FIG. 3.  FIG. 3 Time-dependent electronic population (or occupation) of the donor state in various two-state models treated by different dynamical methods. Black solid thin lines denote P^{\rm{pop}}_{{\rm ml}, k}, the electronic populations obtained by ML-MCTDH calculations. Red dashed lines denote P^{\rm{pop}}_{{\rm map}, k}, the electronic populations obtained by quasi-classical dynamics based on the mapping Hamiltonian. Blue solid thick lines denote the electronic occupations P^{\rm{occ}}_{{\rm map}, k} obtained by the SQC dynamics based on mapping Hamiltonian. Labels A, B and C denote that the mapping results were calculated using no-window quasi-classical dynamics, SQC method with \gamma=1 and SQC method with \gamma=0.732, respectively. Labels 1, 2, 3, 4, 5 denote five two-state models including 7, 11, 16, 44 and 246 vibrational modes respectively. In FIG. 3 (A1)-(A5), both ML-MCTDH and no-window quasi-classical dynamics with mapping Hamiltonian provide similar results of the time-dependent electronic population of the donor states, while in all cases the result difference appears in the long-time limit. When more nuclear degrees of freedom are included, the no-window QC dynamics show more deviation on the time-dependent donor-state population with respect to the ML-MCTDH results. Such difference may come from the improper treatment of the detailed balance, because the quasi-classical dynamics mapping essentially share some similarity to the Ehrenfest dynamics [1, 2]. It is well known that the Ehrenfest dynamics does not obey detailed balance in the long-time simulation limit [1, 2, 55]. Thus, the similar problems may also exist in the current approach. In addition, we also noticed that the electronic population P^{\rm{pop}}_{{\rm map}, k} of the donor state drops to less than zero with time being in the two-state models with 11, 16 and 44 modes. Previous work also noticed similar feature [2, 48]. Because the energy of all modes tends to become equal in classical dynamics, the energy flow takes place from the high-frequency modes to the low-frequency modes. In the current system, two electronic states are represented by two coupled harmonic oscillators. The energy of the high-frequency oscillator (corresponding to the higher diabatic electronic state) drops to lower than its zero-point energy (ZPE) due to the energy flow in classical dynamics [47]. Thus, the electronic population of the donor state finally becomes less than zero. Somehow, this ZPE problem is also one of reasons that cause the improper account of detailed balance [47, 55, 56]. The results based on SQC method with \gamma=1 are given in FIG. 3 (B1)-(B5). Roughly speaking, the results of no-window quasi-classical dynamics with the mapping Hamiltonian and the SQC dynamics with \gamma=1 give similar results on the time-dependent electronic populations, except the visible different in the case of the 7-mode model (FIG. 3(B1)). The similar problem of the negative electronic population is observed in the cases of the two-state models with 11, 16 and 44 modes FIG. 3 (B2), FIG. 3 (B3) and (B4)). This negative-population problem is naturally avoided by considering the average electronic occupations in the SQC dynamics. However, the inclusion of the window trick does not strongly improve the consistency between the electronic occupation in the SQC dynamics and the electronic population P^{\rm{pop}}_{{\rm ml}, k} in the ML-MCTDH dynamics (see FIG. 3 (A) and (B)). Stock and coworkers [47, 48] once proposed to lower the zero-point energy to alleviate the above negative population problem. Later on, Cotton and Miller [35] also recommend using the similar idea in the SQC treatment. By checking the SQC results in several test examples, they recommend to set \gamma=0.732. This gives a lower zero-point energy and a narrow energy window. When their approaches were employed, we still noticed the appearance of the negative electronic population P^{\rm{pop}}_{{\rm map}, k} in the long-time dynamics, the consistency between the electronic occupation P^{\rm{occ}}_{{\rm map}, k} in the SQC dynamics and the electronic population P^{\rm{pop}}_{{\rm ml}, k} in the ML-MCTDH dynamic becomes significantly improved. In the model with 246 modes, the deviation still exists although such difference becomes smaller compared with the SQC dynamics with \gamma=1. For reduced models, the shift of potential minimum, \Delta {x_i}=|\kappa _i^{\left(a \right)} -\kappa _i^{\left(d \right)}{\rm{|/}}{\omega _i}, along each mode was used as criterion to characterize the electron-phonon couplings. And \Delta {x_i}$$>$0.6, $\Delta {x_i}$$>0.5, \Delta {x_i}$$>$0.4, and $\Delta {x_i}$$>0.2 correspond to the 7, 11, 16, and 44 modes. Generally speaking, the electronic coherence is gradually dissipated as the number of modes increases. However, as the number of modes increasing, more vibrational modes with frequencies close to the Rabi oscillation of the pure electronic dynamics are incorporated. This may create the complicated electronic dynamics if the number of modes in the model is not large enough. This may explain why in a few cases the short-time electronic coherence is slightly enhanced. When the mode number is large enough, we expect that the coherence should definitely become weaker with the increasing of mode number. Actually, the coherence in the 246-mode model is weaker than that of 44-mode model both in short-time region and long-time region. B. Dynamics of full dimensional model In this section, the SQC dynamics methods with \gamma=1 and \gamma=0.732 were used to treat the full dimensional ET model including 4 states (D, A1, A2, A3) and 246 vibrational modes. All SQC and ML-MCTDH results are shown in FIG. 4.  FIG. 4 Time-dependent diabatic electronic populations and occupations of the full dimensional model including 4 states (D, A1, A2, A3) and 246 vibrational modes. (A) Populations using ML-MCTDH (P^{\rm{pop}}_{{\rm ml}, k}), (B) populations using the SQC dynamics with \gamma=1 (P^{\rm{pop}}_{{\rm map}, k}), (C) occupations using the SQC dynamics with \gamma=1 (P^{\rm{occ}}_{{\rm map}, k}), (D) populations using the SQC dynamics with \gamma=0.732 (P^{\rm{pop}}_{{\rm map}, k}), (E) occupations using the SQC dynamics with \gamma=0.732 (P^{\rm{occ}}_{{\rm map}, k}). As shown in FIG. 4A, the ML-MCTDH population P^{\rm{pop}}_{\rm{ml, donor}} of the donor state decays to 0.5 in very short time (\sim16 fs). After a few weak recurrences, the populations of D reaches a plateau (nearly zero). After 120 fs, the population dynamics becomes stable, all three acceptor states are populated and their ratio is A1:A2:A3=0.63:0.29:0.08. As expected, the lower A state gets more population. In both QC and SQC dynamics, we first discuss the state occupation instead of electronic population. As shown in FIG. 4(C), the occupation P^{\rm{occ}}_{\rm{map, donor}} of the donor state decays to 0.5 within 12 fs in the SQC dynamics with \gamma=1, which is close to the ML-MCTDH results. However, afterward the decay of the donor-state occupation P^{\rm{occ}}_{\rm{map, donor}} becomes much slower and the time-dependent occupation P^{\rm{occ}}_{\rm{map, donor}} was remarkably different from the ML-MCTDH results P^{\rm{pop}}_{\rm{ml, donor}} (FIG. 4(A)). The occupations of all electronic states become stable after 125 fs. The final occupations of D, A1, A2 and A3 states are 0.29, 0.38, 0.23 and 0.10, respectively at 150 fs. The decay of the donor-state occupation does not go to zero, different from ML-MCTDH results. For the acceptor states, the branching ratio of state occupation at the end of dynamics is qualitatively consistent with the ML-MCTDH results, except that the slightly larger difference exists in the occupation of the A1 state. The above difference between SQC and ML-MCTDH results is obviously reduced when \gamma=0.732 was used in the SQC dynamics (FIG. 4(E)). This improvement is consistent with results in the two-state models. The occupation P^{\rm{occ}}_{\rm{map, donor}} of the donor state decays to 0.5 in about 16 fs. Although the decay of P^{\rm{occ}}_{\rm{map, donor}} in the SQC dynamics with \gamma=0.732 also becomes a little slower compared to the ML-MCTDH result, the overall decay feature is similar. At 150 fs, the occupation of the D, A1, A2 and A3 states are 0.20, 0.49, 0.26 and 0.09, respectively. However, due to the lack of detail balance, the electronic occupation P^{\rm{occ}}_{{\rm map}, k} deviates from the electronic population P^{\rm{pop}}_{{\rm map}, k} in the long-time limit. Although the existence of this deficiency, the SQC dynamics with \gamma=0.732 seems to give the consistent results with the ML-MCTDH dynamics, at least at the qualitative and semi-quantitative levels. As shown in FIG. 4 (B) and (D), it is very interesting to notice that the electronic populations P^{\rm{pop}}_{{\rm map}, k} of the SQC dynamics are similar to the ML-MCTDH results, particular in the case of \gamma=0.732. This tendency is different from the situation in the two-state models. However, we should notice that such consistency is only valid for the short-time dynamics. In the long-time SQC dynamics, the electronic population P^{\rm{pop}}_{\rm{map, donor}} of the donor state always drops below zero in both the \gamma=1 (FIG. 5(A)) and \gamma=0.732 (FIG. 5(B)) situations. Thus, as discussed in previous work and the 2-state model in the current work, it seems that the population of the higher electronic state always tends to become negative in the long-time dynamics. As discussed previously [47], this problem is associated with the zero-point energy leaking problem as the results of fully classical simulation. In this sense, it is not suitable to use P^{\rm{pop}}_{{\rm map}, k} to describe the population dynamics in the QC dynamics with mapping Hamiltonian. This also indicates that some additional efforts should be made to improve the performance of QC dynamics or to develop more consistent but cheap approaches.  FIG. 5 Long-time time-dependent diabatic electronic populations and occupations of the full dimensional model including 4 states (D, A1, A2, A3) and 246 vibrational modes. (A) Populations using the SQC dynamics with \gamma=1 (P^{\rm{pop}}_{{\rm map}, k}), (B) populations using the SQC dynamics with \gamma=0.732 (P^{\rm{pop}}_{{\rm map}, k}). FIG. 3 and FIG. 4 indicate that both mapping-Hamiltonian based QC and SQC methods underestimate the oscillation pattern in the population dynamics. It is well known that such oscillation is governed by coherence, thus it means that both approaches underestimate the coherence. Previous studies [37, 39] showed that the SQC might reproduce the coherence effects in some models, while it also underestimates such effects in other cases, such as with strong electronic couplings or at low temperature. At the same time, the situation with overestimation of coherence does not appear at all in both of previous studies and current work. The underline reason may be relevant to the fact that the mapping-Hamiltonian based QC and SQC methods employed the classical dynamics, while the coherence is a quantum effect. Also due to this reason, the underestimation of coherence happens much easier at the low temperature. Further discussion on temperature effect is given in Appendix. Overall, we found that the electronic occupation in the SQC dynamics with \gamma=0.732 can provide the reasonable description on the nonadiabatic dynamics of the current model at least at the qualitative and semi-quantitative levels. Considering its high efficiency, it is worthwhile to employ it in the treatment of other model systems with similar Hamiltonian, or even more general systems. Therefore, more additional efforts should be made in the further development of the QC and semi-classical dynamics method based on mapping Hamiltonian. Ⅳ. CONCLUSION In this work we examined the performance of the quasi-classical dynamics method based on mapping Hamiltonian by checking its results against the accurate ones obtained with the ML-MCTDH method. The electron transfer models in acene/C_{60} system studied in our previous work was used in current work. A few of reduced dimensional two-state models were constructed, which includes different numbers of vibrational modes. The full dimensional model with 246 modes and 4 electronic states was also used. We consider the non-window QC dynamics with the mapping Hamiltonian and the SQC dynamics with the mapping Hamiltonian. In all cases, all quasi-classical dynamics methods can capture the main feature of the nonadiabatic dynamics qualitatively. However, the negative electronic population may appear in the long-time dynamics no matter whether the 'window' trick is applied or not. The electronic occupation obtained in the SQC dynamics is a good alternative choice to measure the electronic dynamics because of its positive feature. Consistent with previous studies, it is also preferable to choose a lower zero-point energy and narrow energy window (\gamma=0.732) to improve the accuracy of the SQC dynamics. Although some difference exists between the results of ML-MCTDH and SQC dynamics, the SQC dynamics with the mapping Hamiltonian give qualitatively reasonable results. Because the SQC dynamics method with the mapping Hamiltonian is quite simple and efficient, thus it shows great potential to treat nonadiabatic dynamics of extremely large systems. Also, due to such reason, additional efforts should be made to make the further improvement of the quasi-classical and semiclassical dynamics based on mapping Hamiltonian. Ⅴ. Acknowledgements This work was supported by the National Natural Science Foundation of China (No.11747170, No.21673266, and No.21503248), the Natural Science Foundation of Shandong Province for Distinguished Young Scholars (JQ201504). The authors also thank Supercomputing Centre, Computer Network Information Center, CAS and the Super Computational Centre of CAS-QIBEBT for providing computational resources. APPENDIX A: Electron-phonon couplings For the reduced models, the shift of potential minimum, \Delta {x_i}=|\kappa _i^{\left(a \right)}$$-$$\kappa _i^{\left(d \right)}|/{\omega _i}, along each mode was used to characterize the electron-phonon couplings. Frequencies and coupling strengths of vibrational modes of the acene/C_{60} system are provided in FIG. 6.  FIG. 6 The shift of potential minima, \Delta x_i, for all normal modes (left column) and along the modes with \omega_i$$>$100 cm$^{-1}$ (right column) for the three acceptor states (denoted as A1, A2, A3). Reproduction with the permission from J. Chem. Phys. 142, 084706 (2015). Copyright 2015 AIP Publishing.
APPENDIX B: Convergence test

Generally speaking, 2000$-$5000 trajectories are needed to obtain converged results in the SQC calculations. We choose the 4-state-246-mode model with $\gamma$=0.732 as an example. Population and occupation results with different number of trajectories are given in FIG. 7. The calculations with 2000 trajectories give reasonably converged results. In all of our mapping calculations, 5000 trajectories were used to make sure to achieve convergence.

 FIG. 7 Time-dependent diabatic electronic populations and occupations of the full dimensional model with different number of trajectories. (A) Populations using the SQC dynamics with $\gamma$=0.732 ($P^{\rm{pop}}_{{\rm map}, k}$), (B) occupations using the SQC dynamics with $\gamma$=0.732 ($P^{\rm{occ}}_{{\rm map}, k}$).
APPENDIX C: Temperature effect

The temperature applied in most calculations of this paper is 0 K. The classical treatment of molecular motion (QC and SQC methods) is more suitable for high temperature and the ML-MCTDH method is more efficient at zero and low temperature.

In the current model system, the ET dynamics is extremely fast ($<$100 fs). This means that very low-frequency modes with slow motion should have minor effects on the population dynamics and the temperature effect should be very weak. This was proven by our previous work [44] using ML-MCTDH method. We also run both QC and SQC calculations at two different temperatures, 0 and 300 K (FIG. 8). For all three methods, the temperature effects are really minor. Certainly, if the dynamics become slower, the temperature effects should become more significant.

 FIG. 8 Time-dependent diabatic electronic populations and occupations of donor state in the full dimensional model at different temperatures of 0 and 300 K. (A) Populations using the SQC dynamics with $\gamma$=0.732 ($P^{\rm{pop}}_{\rm{map, donor}}$), (B) occupations using the SQC dynamics with $\gamma$=0.732 ($P^{\rm{occ}}_{\rm{map, donor}}$).

Because low frequency modes are not relevant to the ultrafast dynamics here, we may think that the current model in fact defines a low-temperature situation considering relevant modes. It is certainly true that the QC and SQC methods may show some deviation in this situation. However, we still wish to run QC and SQC dynamics, because we really wish to know how large of the result deviation may come out by using these highly-approximated methods, compare with the accurate ML-MCTDH results. Particularly, we wish to know whether the results are at least qualitatively/semi-quantitatively reasonable or completely wrong. As discussed in the main manuscript, we found that the electronic occupation in the SQC dynamics with $\gamma$=0.732 can provide the reasonable description on the nonadiabatic dynamics of the current model at least at the qualitative and semi-quantitative levels.

Reference
 [1] W. Domcke, D. R. Yarkony, and H. Köppel, Conical Intersections:Electronic Structure, Dynamics and Spectroscopy. Singapore: World Scientic (2004). [2] G. Stock, and M. Thoss, Adv. Chem. Phys. 131 , 243 (2005). [3] M. H. Beck, A. Jackle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324 , 1 (2000). DOI:10.1016/S0370-1573(99)00047-2 [4] H. B. Wang, and M. Thoss, J. Chem. Phys. 119 , 1289 (2003). DOI:10.1063/1.1580111 [5] U. Manthe, J. Chem. Phys. 128 , 164116 (2008). DOI:10.1063/1.2902982 [6] U. Manthe, J. Chem. Phys. 130 , 054109 (2009). DOI:10.1063/1.3069655 [7] O. Vendrell, and H. D. Meyer, J. Chem. Phys. 134 , 044135 (2011). DOI:10.1063/1.3535541 [8] M. Ben-Nun, and T. J. Martinez, Adv. Chem. Phys. 121 , 439 (2002). [9] I. Burghardt, H. D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 111 , 2927 (1999). DOI:10.1063/1.479574 [10] G. A. Worth, and I. Burghardt, Chem. Phys. Lett. 368 , 502 (2003). DOI:10.1016/S0009-2614(02)01920-6 [11] S. C. Cheng, C. Y. Zhu, K. K. Liang, S. H. Lin, and D. G. Truhlar, J. Chem. Phys. 129 , 024112 (2008). DOI:10.1063/1.2948395 [12] J. C. Tully, J. Chem. Phys. 93 , 1061 (1990). DOI:10.1063/1.459170 [13] G. Granucci, and M. Persico, J. Chem. Phys. 126 , 134114 (2007). DOI:10.1063/1.2715585 [14] M. H. Yang, C. Y. Huo, A. Y. Li, Y. B. Lei, L. Yu, and C. Y. Zhu, Phys. Chem. Chem. Phys. 19 , 12185 (2017). DOI:10.1039/C7CP00102A [15] M. Persico, and G. Granucci, Theor. Chem. Acc. 133 , 1526 (2014). DOI:10.1007/s00214-014-1526-1 [16] L. K. Du, and Z. G. Lan, J. Chem. Theory Comput. 11 , 1360 (2015). DOI:10.1021/ct501106d [17] L. K. Du, and Z. G. Lan, J. Chem. Theory Comput. 11 , 4522 (2015). DOI:10.1021/acs.jctc.5b00654 [18] H. J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, WIREs Comput. Mol. Sci. 2 , 242 (2012). DOI:10.1002/wcms.82 [19] M. Barbatti, G. Granucci, M. Persico, M. Ruckenbauer, M. Vazdar, M. Eckert-Maksic, and H. Lischka, J. Photoch. Photobio. A 190 , 228 (2007). DOI:10.1016/j.jphotochem.2006.12.008 [20] R. Mitric, U. Werner, and V. Bonacic-Koutecky, J. Chem. Phys. 129 , 164118 (2008). [21] E. Fabiano, T. W. Keal, and W. Thiel, Chem. Phys. 349 , 334 (2008). DOI:10.1016/j.chemphys.2008.01.044 [22] T. Nelson, S. Fernandez-Alberti, V. Chernyak, A. E. Roitberg, and S. Tretiak, J. Chem. Phys. 136 , 054108 (2012). DOI:10.1063/1.3680565 [23] N. L. Doltsinis, and D. Marx, Phys. Rev. Lett. 88 , 166402 (2002). DOI:10.1103/PhysRevLett.88.166402 [24] I. Tavernelli, E. Tapavicza, and U. Rothlisberger, J. Mol. Struct. 914 , 22 (2009). DOI:10.1016/j.theochem.2009.04.020 [25] J. E. Subotnik, and N. Shenvi, J. Chem. Phys. 134 , 024105 (2011). DOI:10.1063/1.3506779 [26] U. Müller, and G. Stock, J. Chem. Phys. 107 , 6230 (1997). DOI:10.1063/1.474288 [27] H. D. Meyer, and W. H. Miller, J. Chem. Phys. 70 , 3214 (1979). DOI:10.1063/1.437910 [28] H. D. Meyer, and W. H. Miller, J. Chem. Phys. 71 , 2156 (1979). DOI:10.1063/1.438598 [29] G. Stock, and M. Thoss, Phys. Rev. Lett. 78 , 578 (1997). DOI:10.1103/PhysRevLett.78.578 [30] J. Liu, J. Chem. Phys. 145 , 204105 (2016). DOI:10.1063/1.4967815 [31] H. B. Wang, X. Sun, and W. H. Miller, J. Chem. Phys. 108 , 9726 (1998). DOI:10.1063/1.476447 [32] S. Bonella, and D. F. Coker, J. Chem. Phys. 114 , 7778 (2001). DOI:10.1063/1.1366331 [33] S. Bonella, and D. F. Coker, Chem. Phys. 268 , 189 (2001). DOI:10.1016/S0301-0104(01)00329-9 [34] M. Karplus, and R. D. Sharma, J. Chem. Phys. 43 , 3259 (1965). DOI:10.1063/1.1697301 [35] S. J. Cotton, and W. H. Miller, J. Phys. Chem. A 117 , 7190 (2013). [36] S. J. Cotton, and W. H. Miller, J. Chem. Phys. 139 , 234112 (2013). DOI:10.1063/1.4845235 [37] S. J. Cotton, K. Igumenshchev, and W. H. Miller, J. Chem. Phys. 141 , 084104 (2014). DOI:10.1063/1.4893345 [38] S. J. Cotton, and W. H. Miller, J. Chem. Phys. 145 , 144108 (2016). DOI:10.1063/1.4963914 [39] S. J. Cotton, and W. H. Miller, J. Chem. Theory Comput. 12 , 983 (2016). DOI:10.1021/acs.jctc.5b01178 [40] S. J. Cotton, and W. H. Miller, J. Phys. Chem. A 119 , 12138 (2015). [41] G. Tao, J. Chem. Phys. 144 , 094108 (2016). DOI:10.1063/1.4943006 [42] G. H. Tao, J. Phys. Chem. Lett. 7 , 4335 (2016). DOI:10.1021/acs.jpclett.6b01857 [43] G. Tao, J. Chem. Phys. 147 , 044107 (2017). DOI:10.1063/1.4985898 [44] Y. Xie, J. Zheng, and Z. Lan, J. Chem. Phys. 142 , 084706 (2015). DOI:10.1063/1.4909521 [45] E. Wigner, Phys. Rev. 40 , 0749 (1932). DOI:10.1103/PhysRev.40.749 [46] M. Hillery, R. F. Oconnell, M. O. Scully, and E. P. Wigner, Phys. Rep. 106 , 121 (1984). DOI:10.1016/0370-1573(84)90160-1 [47] G. Stock, and U. Müller, J. Chem. Phys. 111 , 65 (1999). DOI:10.1063/1.479254 [48] U. Müller, and G. Stock, J. Chem. Phys. 111 , 77 (1999). DOI:10.1063/1.479255 [49] B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martinez, Chem. Phys. 347 , 3 (2008). DOI:10.1016/j.chemphys.2008.01.014 [50] M. Schroter, S. D. Ivanov, J. Schulze, S. P. Polyutov, Y. Yan, T. Pullerits, and O. Kuhn, Phys. Rep. 567 , 1 (2015). DOI:10.1016/j.physrep.2014.12.001 [51] Q. Meng, S. Faraji, O. Vendrell, and H. D. Meyer, J. Chem. Phys. 137 , 134302 (2012). DOI:10.1063/1.4755372 [52] J. Zheng, Y. Xie, S. Jiang, and Z. Lan, J. Phys. Chem. C 120 , 1375 (2016). DOI:10.1021/acs.jpcc.5b09921 [53] S. Jiang, J. Zheng, Y. Yi, Y. Xie, F. Yuan, and Z. Lan, J. Phys. Chem. C 121 , 27263 (2017). DOI:10.1021/acs.jpcc.7b08175 [54] O. Vendrell, and H. D. Meyer, The MCTDH Package, Version 8.5 (2011). [55] P. V. Parandekar, and J. C. Tully, J. Chem. Theory Comput. 2 , 229 (2006). DOI:10.1021/ct050213k [56] W. H. Miller, and S. J. Cotton, J. Chem. Phys. 142 , 131103 (2015). DOI:10.1063/1.4916945

a. 青岛大学非织造材料与产业用纺织品创新研究院, 青岛 266071;
b. 青岛大学物理学院, 青岛 266071;
c. 中国科学院青岛生物能源与过程研究所中科院生物基材料重点实验室, 青岛 266101;
d. 中国科学院大学, 北京 100049