The article information
 Jie Zheng, Yu Xie, Shengshi Jiang, Yunze Long, Xin Ning, Zhenggang Lan
 郑杰, 谢宇, 江圣世, 龙云泽, 宁新, 兰峥岗
 Ultrafast Electron Transfer with Symmetrical Quasiclassical Dynamics based on Mapping Hamiltonian and Quantum Dynamics based on MLMCTDH^{†}
 使用基于MAPPING哈密顿的对称准经典动力学方法和基于MLMCTDH的量子动力学方法研究超快电子转移
 Chinese Journal of Chemical Physics, 2017, 30(6): 800810
 化学物理学报, 2017, 30(6): 800810
 http://dx.doi.org/10.1063/16740068/30/cjcp1711210

Article history
 Received on: November 14, 2017
 Accepted on: December 27, 2017
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
The simulation of the nonadiabatic dynamics of complex systems is very challenging due to the breakdown of the BornOppenheimer (BO) approximation and the involvement of a large number of stronglycoupled nuclear/electronic freedoms [1, 2]. Although the theoretical approaches based on quantum dynamics, such as full quantum wavepacket dynamics [1], the multiconfiguration timedependent Hartree (MCTDH) method [3] and multilayer MCTDH (MLMCTDH) method [47], can provide the accurate solution of nonBO dynamics, they suffer from the high computational cost. The Gaussianwavepacket based methods, such as ab initio multiple spawning (AIMS) [8], Gaussianbased MCTDH (GMCTDH) [9, 10], are promising methods for polyatomic 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 meanfield (Ehrenfest) dynamics and its extensions [1, 2, 11], surfacehopping dynamics [1, 2, 1214], and so on. After the implementation of these efficient methods in the manner of onthefly approaches [1524], 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 widelyused fewestswitches trajectory surfacehopping 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 (MeyerMillerStockThoss) model maps
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 socalled symmetrical quasiclassical trajectory (SQC) method. Because of its simple implementation and its good performance in a few of model systems [3740], this mappingHamiltonianbased 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 [4143].
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 onthefly 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 semiclassical 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
The diabatic model in our previous work [44] was taken in the current calculations, which describes the ultrafast acene (Ac)
$ \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
$ \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
In our previous work [44], all parameters in this fourstate 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 electronphonon couplings, we have already discussed all details in previous work [44] and here we only give a short description in Appendix.
B. Quasiclassical dynamics method based on mapping Hamiltonian 1. Mapping HamiltonianThe Hamiltonian of a
$ \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
In the MMST approaches [1, 2, 2628], a
$ \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
$ \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
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
$ \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
The mapping approach naturally provides a way to construct the initial phasespace distribution for the electronic part. The coordinate and momentum are sampled by the actionangle sampling method [2] with fixed electronic action
$ \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
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 actionangle 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 [1525, 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
In the quasiclassical 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 initialstate sampling and finalstate assignment. In this symmetrical quasiclassical (SQC) method, the classical action variables
$ \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
In the SQC approach, the initial sampling of the phase space for the electronic part is performed with the actionangle sampling by setting the action within a window with the width of
$ 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
Next, we explain the numerical details of the SQC methods. At time 0, we assume that only one electronic state (state
Following previous discussions [35, 47], the SQC dynamics with
In traditional quantum wave packet dynamics, a group of timeindependent basis for each primary coordinate is used to represent the wave function, and the timedependent coefficients evolve with time. MCTDH employs the timedependent 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
MLMCTDH is an extension of the standard MCTDH method. In MLMCTDH, the wave function is represented by a multilayer tree structure (shown in FIG. 2), which is expanded by the multidimensional SPFs recursively until the timeindependent 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) 
where
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 [47, 5053].
We tried to run MLMCTDH calculations to get the accurate dynamical results for benchmark. The MLMCTDH calculations were performed using the Heidelberg MCTDH package [54]. For convenience, we used
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 longtime dynamics behavior of the QC and SQC dynamics, all MLMCTDH dynamics are recalculated 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 twostate modelsWe first considered a few of twostate 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 timedependent electronic populations were calculated using MLMCTDH and different mapping approaches. Because previous studies recommended the symmetrical window, we do not consider the calculations of state occupation
In FIG. 3 (A1)(A5), both MLMCTDH and nowindow quasiclassical dynamics with mapping Hamiltonian provide similar results of the timedependent electronic population of the donor states, while in all cases the result difference appears in the longtime limit. When more nuclear degrees of freedom are included, the nowindow QC dynamics show more deviation on the timedependent donorstate population with respect to the MLMCTDH results. Such difference may come from the improper treatment of the detailed balance, because the quasiclassical 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 longtime 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
The results based on SQC method with
Stock and coworkers [47, 48] once proposed to lower the zeropoint 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
For reduced models, the shift of potential minimum,
In this section, the SQC dynamics methods with
As shown in FIG. 4A, the MLMCTDH 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
The above difference between SQC and MLMCTDH results is obviously reduced when
As shown in FIG. 4 (B) and (D), it is very interesting to notice that the electronic populations
FIG. 3 and FIG. 4 indicate that both mappingHamiltonian 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 mappingHamiltonian 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
In this work we examined the performance of the quasiclassical dynamics method based on mapping Hamiltonian by checking its results against the accurate ones obtained with the MLMCTDH method. The electron transfer models in acene/C
In all cases, all quasiclassical dynamics methods can capture the main feature of the nonadiabatic dynamics qualitatively. However, the negative electronic population may appear in the longtime 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 zeropoint energy and narrow energy window (
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 CASQIBEBT for providing computational resources.
APPENDIX A: Electronphonon couplingsFor the reduced models, the shift of potential minimum,
Generally speaking, 2000
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 MLMCTDH method is more efficient at zero and low temperature.
In the current model system, the ET dynamics is extremely fast (
Because low frequency modes are not relevant to the ultrafast dynamics here, we may think that the current model in fact defines a lowtemperature 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 highlyapproximated methods, compare with the accurate MLMCTDH results. Particularly, we wish to know whether the results are at least qualitatively/semiquantitatively reasonable or completely wrong. As discussed in the main manuscript, we found that the electronic occupation in the SQC dynamics with
[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/S03701573(99)000472 
[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. BenNun, 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/S00092614(02)019206 
[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/s0021401415261 
[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. EckertMaksic, 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. BonacicKoutecky, 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. FernandezAlberti, 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/S03010104(01)003299 
[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/03701573(84)901601 
[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 
b. 青岛大学物理学院, 青岛 266071;
c. 中国科学院青岛生物能源与过程研究所中科院生物基材料重点实验室, 青岛 266101;
d. 中国科学院大学, 北京 100049