Volume 33 Issue 5
Nov.  2020
Turn off MathJax
Article Contents

Bing-yang Xiao, Jia-bo Xu, Lin-jun Wang. New Energy-Based Decoherence Correction Approaches for Trajectory Surface Hopping†[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 603-612. doi: 10.1063/1674-0068/cjcp2006098
Citation: Bing-yang Xiao, Jia-bo Xu, Lin-jun Wang. New Energy-Based Decoherence Correction Approaches for Trajectory Surface Hopping[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 603-612. doi: 10.1063/1674-0068/cjcp2006098

New Energy-Based Decoherence Correction Approaches for Trajectory Surface Hopping

doi: 10.1063/1674-0068/cjcp2006098
More Information
  • Corresponding author: Lin-jun Wang, E-mail: ljwang@zju.edu.cn
  • Part of the special issue for "the Chinese Chemical Society's 16th National Chemical Dynamics Symposium".
  • Received Date: 2020-06-16
  • Accepted Date: 2020-07-18
  • Publish Date: 2020-10-27
  • Inspired by the branching corrected surface hopping (BCSH) method [J. Xu and L. Wang, J. Chem. Phys. 150 , 164101 (2019)], we present two new decoherence time formulas for trajectory surface hopping. Both the proposed linear and exponential formulas characterize the decoherence time as functions of the energy difference between adiabatic states and correctly capture the decoherence effect due to wave packet reflection as predicted by BCSH. The relevant parameters are trained in a series of 200 diverse models with different initial nuclear momenta, and the exact quantum solutions are utilized as references. As demonstrated in the three standard Tully models, the two new approaches exhibit significantly higher reliability than the widely used counterpart algorithm while holding the appealing efficiency, thus promising for nonadiabatic dynamics simulations of general systems.
  • Part of the special issue for "the Chinese Chemical Society's 16th National Chemical Dynamics Symposium".
  • 加载中
  • [1] S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys. 101, 4657 (1994). doi:  10.1063/1.467455
    [2] K. Song and Q. Shi, J. Chem. Phys. 146, 184108 (2017). doi:  10.1063/1.4982928
    [3] X. Li, D. Hu, Y. Xie, and Z. Lan, J. Chem. Phys. 149, 244104 (2018). doi:  10.1063/1.5048049
    [4] X. Pang, C. Jiang, Y. Qi, L. Yuan, D. Hu, X. Zhang, D. Zhao, D. Wang, Z. Lan, and F. Li, Phys. Chem. Chem. Phys. 20, 25910 (2018). doi:  10.1039/C8CP04762F
    [5] T. Zhang, Y. Fang, X. Song, W. Fang, and G. Cui, J. Phys. Chem. Lett. 11, 2470 (2020). doi:  10.1021/acs.jpclett.0c00294
    [6] X. Gao, H. Geng, Q. Peng, J. Ren, Y. Yi, D. Wang, and Z. Shuai, J. Phys. Chem. C 118, 6631 (2014). doi:  10.1021/jp412782n
    [7] L. Wang, O. V. Prezhdo, and D. Beljonne, Phys. Chem. Chem. Phys. 17, 12395 (2015). doi:  10.1039/C5CP00485C
    [8] J. Spencer, F. Gajdos, and J. Blumberger, J. Chem. Phys. 145, 064102 (2016). doi:  10.1063/1.4960144
    [9] X. Gao, Q. Peng, Y. Niu, D. Wang, and Z. Shuai, Phys. Chem. Chem. Phys. 14, 14207 (2012). doi:  10.1039/c2cp40347a
    [10] J. Huang, L. Du, J. Wang, and Z. Lan, J. Phys. Chem. C 119, 7578 (2015).
    [11] X. Shen, G. Han, D. Fan, Y. Xie, and Y. Yi, J. Phys. Chem. C 119, 11320 (2015).
    [12] A. V. Akimov and O. V. Prezhdo, J. Am. Chem. Soc. 136, 1599 (2014). doi:  10.1021/ja411800n
    [13] L. Wang, Y. Olivier, O. V. Prezhdo, and D. Beljonne, J. Phys. Chem. Lett. 5, 3345 (2014). doi:  10.1021/jz5015955
    [14] R. Crespo-Otero and M. Barbatti, Chem. Rev. 118, 7026 (2018).
    [15] J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
    [16] L. Wang, D. Trivedi, and O. V. Prezhdo, J. Chem. Theory Comput. 10, 3598 (2014). doi:  10.1021/ct5003835
    [17] L. Wang, A. E. Sifain, and O. V. Prezhdo, J. Phys. Chem. Lett. 6, 3827 (2015). doi:  10.1021/acs.jpclett.5b01502
    [18] L. Wang, A. E. Sifain, and O. V. Prezhdo, J. Chem. Phys. 143, 191102 (2015). doi:  10.1063/1.4935971
    [19] J. C. Tully, J. Chem. Phys. 137, 22A301 (2012).
    [20] L. Wang, A. Akimov, and O. V. Prezhdo, J. Phys. Chem. Lett. 7, 2100 (2016). doi:  10.1021/acs.jpclett.6b00710
    [21] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016). doi:  10.1146/annurev-physchem-040215-112245
    [22] W. Feng, L. Xu, X. Li, W. Fang, and Y. Yan, AIP Adv. 4, 077131 (2014). doi:  10.1063/1.4891821
    [23] B. Xie, L. Liu, G. Cui, W. Fang, J. Cao, W. Feng, and X. Li, J. Chem. Phys. 143, 194107 (2015). doi:  10.1063/1.4935800
    [24] C. C. Martens, J. Phys. Chem. Lett. 7, 2610 (2016). doi:  10.1021/acs.jpclett.6b01186
    [25] C. Zhu, A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. 120, 5543 (2004). doi:  10.1063/1.1648306
    [26] C. Zhu, S. Nangia, A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. 121, 7658 (2004). doi:  10.1063/1.1793991
    [27] C. Zhu, A. W. Jasper, and D. G. Truhlar, J. Chem. Theory Comput. 1, 527 (2005). doi:  10.1021/ct050021p
    [28] S. C. Cheng, C. Zhu, K. K. Liang, S. H. Lin, and D. G. Truhlar, J. Chem. Phys. 129, 024112 (2008). doi:  10.1063/1.2948395
    [29] G. Granucci and M. Persico, J. Chem. Phys. 126, 134114 (2007). doi:  10.1063/1.2715585
    [30] Y. Zhang, X. Sun, T. Zhang, X. Liu, and G. Cui, J. Phys. Chem. A 124, 2547 (2020). doi:  10.1021/acs.jpca.0c00791
    [31] H. Uratani and H. Nakai, J. Phys. Chem. Lett. 11, 4448 (2020). doi:  10.1021/acs.jpclett.0c01028
    [32] E. R. Bittner and P. J. Rossky, J. Chem. Phys. 103, 8130 (1995). doi:  10.1063/1.470177
    [33] B. J. Schwartz, E. R. Bittner, O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 104, 5942 (1996). doi:  10.1063/1.471326
    [34] O. V. Prezhdo and P. J. Rossky, J. Chem. Phys. 107, 5863 (1997). doi:  10.1063/1.474312
    [35] J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 024105 (2011). doi:  10.1063/1.3506779
    [36] J. E. Subotnik, J. Phys. Chem. A 115, 12083 (2011). doi:  10.1021/jp206557h
    [37] A. Jain, E. Alguire, and J. E. Subotnik, J. Chem. Theory Comput. 12, 5256 (2016). doi:  10.1021/acs.jctc.6b00673
    [38] N. Shenvi, J. E. Subotnik, and W. Yang, J. Chem. Phys. 134, 144102 (2011). doi:  10.1063/1.3575588
    [39] H. M. Jaeger, S. Fischer, and O. V. Prezhdo, J. Chem. Phys. 137, 22A545 (2012).
    [40] E. J. Heller, J. Chem. Phys. 62, 1544 (1975). doi:  10.1063/1.430620
    [41] J. Xu and L. Wang, J. Chem. Phys. 150, 164101 (2019). doi:  10.1063/1.5090927
    [42] N. Shenvi, J. E. Subotnik, and W. Yang, J. Chem. Phys. 135, 024101 (2011). doi:  10.1063/1.3603447
    [43] N. Shenvi and W. Yang, J. Chem. Phys. 137, 22A528 (2012).
    [44] D. T. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1992). doi:  10.1063/1.462100
    [45] D. E. Manolopoulos, J. Chem. Phys. 117, 9552 (2002).
    [46] T. Gonzalez-Lezana, E. J. Rackham, and D. E. Manolopoulos, J. Chem. Phys. 120, 2247 (2004).
    [47] L. Wang, J. Qiu, X. Bai, and J. Xu, WIREs Comput. Mol. Sci. 10, e1435 (2020).
    [48] G. Granucci, M. Persico, and A. Toniolo, J. Chem. Phys. 114, 10608 (2001).
    [49] S. Fernandez-Alberti, A. E. Roitberg, T. Nelson, and S. Tretiak, J. Chem. Phys. 137, 014512 (2012).
    [50] L. Wang and O. V. Prezhdo, J. Phys. Chem. Lett. 5, 713 (2014).
    [51] X. Bai, J. Qiu, and L. Wang, J. Chem. Phys. 148, 104106 (2018).
    [52] J. Qiu, X. Bai, and L. Wang, J. Phys. Chem. Lett. 9, 4319 (2018).
    [53] J. Qiu, X. Bai, and L. Wang, J. Phys. Chem. Lett. 10, 637 (2019).
    [54] K. Hornik, M. Stinchcombe, and H. White, Neural Netw. 2, 359 (1989).
  • 加载中
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Figures(7)

Article Metrics

Article views(3) PDF downloads(2) Cited by()

Proportional views
Related

New Energy-Based Decoherence Correction Approaches for Trajectory Surface Hopping

doi: 10.1063/1674-0068/cjcp2006098

Abstract: Inspired by the branching corrected surface hopping (BCSH) method [J. Xu and L. Wang, J. Chem. Phys. 150 , 164101 (2019)], we present two new decoherence time formulas for trajectory surface hopping. Both the proposed linear and exponential formulas characterize the decoherence time as functions of the energy difference between adiabatic states and correctly capture the decoherence effect due to wave packet reflection as predicted by BCSH. The relevant parameters are trained in a series of 200 diverse models with different initial nuclear momenta, and the exact quantum solutions are utilized as references. As demonstrated in the three standard Tully models, the two new approaches exhibit significantly higher reliability than the widely used counterpart algorithm while holding the appealing efficiency, thus promising for nonadiabatic dynamics simulations of general systems.

Part of the special issue for "the Chinese Chemical Society's 16th National Chemical Dynamics Symposium".
Bing-yang Xiao, Jia-bo Xu, Lin-jun Wang. New Energy-Based Decoherence Correction Approaches for Trajectory Surface Hopping†[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 603-612. doi: 10.1063/1674-0068/cjcp2006098
Citation: Bing-yang Xiao, Jia-bo Xu, Lin-jun Wang. New Energy-Based Decoherence Correction Approaches for Trajectory Surface Hopping[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 603-612. doi: 10.1063/1674-0068/cjcp2006098
  • In chemistry, physics, biology, and material sciences, there exist many important dynamical processes, e.g., proton transfer [1, 2], photoisomerization [3-5], charge transport [6-8], nonradiative relaxation [9], energy transfer [10], exciton dissociation [11], and singlet fission [12, 13], which all belong to the category of nonadiabatic dynamics. In these circumstances, the traditional Born-Oppenheimer approximation breaks down and the motion of electronic and nuclear degrees of freedom becomes strongly coupled. In the past few decades, mixed quantum-classical dynamics has become a popular approach to simulate nonadiabatic dynamics [14]. Especially, Tully's fewest switches surface hopping (FSSH) algorithm combined with ab initio electronic structure calculations has achieved great success in different research fields [15]. In the standard FSSH and its variants, the electronic wavefunction is propagated quantum mechanically through solving the time-dependent Schrödinger equation (TDSE), the nuclei evolve classically on an active potential energy surface (PES), and surface hops occur stochastically on the basis of hopping probabilities [15-18]. Besides its advantage in formulation and implementation, trajectory surface hopping has generally shown a good balance between efficiency and reliability [19, 20].

    According to quantum mechanics, a quantum system gradually loses coherence when contacting with a classical bath, resulting in the quantum decoherence effect [21]. In quantum measurements, the quantum system is stochastically collapsed onto a pure state [22, 23]. However, this has not been taken into account in the standard FSSH algorithm with independent trajectories because the electronic wavefunction is propagated fully coherently through solving the TDSE [15]. In principle, surface hopping with correlated trajectories can naturally account for the decoherence effect [24]. However, the computational cost for a reliable description of large systems with correlated trajectories is generally prohibitive and additional approximations are required.

    Within the independent trajectory machinery of trajectory surface hopping, decoherence corrections are widely implemented with a predefined decoherence time function related to quantities of individual trajectories [25-39]. Zhu, Truhlar and coworkers proposed a series of decay of mixing methods [25-28] and suggested a general expression for the decoherence time based on the energy difference between adiabatic states and the nuclear kinetic energy along a decoherent direction. This scheme was further simplified by Granucci and Persico to obtain a purely energy-based decoherence time formula [29]. Due to the low computational cost and ease of implementation, this approach has been widely used in the literature for surface hopping simulations of realistic systems [30, 31]. Through adopting the frozen Gaussian ansatz [40] and tracking the exponential decay of overlap between relevant nuclear wave packets (WPs), Rossky, Bittner, Schwartz, and Prezhdo proposed the force-based decoherence time [32-34], which relies on the force difference between different PESs and uses the spatial width of the Gaussian WP as a parameter. Subotnik and coworkers derived a more rigorous expression of the decoherence time from the quantum Liouville equation and proposed the augmented FSSH (AFSSH) algorithm for decoherence correction after making several approximations [35-37]. Shenvi and coworkers proposed the simultaneous-trajectory surface hopping approach, where a trajectory is propagated on each PES simultaneously and the decoherence time is defined according to the decay of off-diagonal components of the electronic density matrix [38]. Prezhdo and coworkers defined the decoherence time according to pure dephasing through the response function for classical fluctuations [39]. In general, available decoherence correction methods require additional parameters or give certain physical limits.

    Recently, we have proposed the branching corrected surface hopping (BCSH) method, which makes decoherence corrections when WP reflection is identified [41]. As benchmarked in a large number of systems, BCSH has shown high accuracy compared with existing decoherence correction approaches. In this study, we propose two new energy-based decoherence time formulas, which can properly describe the decoherence effect due to WP reflection without calculating nuclear forces. Namely, decoherence correction is stronger for smaller nuclear kinetic energy and larger energy difference between states. These formulas achieve more reliable description of nonadiabatic dynamics than the widely used energy-based decoherence correction approach in the literature [29].

  • We consider a general system containing both electronic and nuclear degrees of freedom with coordinates $ {\mathit{\boldsymbol{r}}} $ and $ {\mathit{\boldsymbol{x}}} $, respectively. The total Hamiltonian reads

    where $ T_{\mathit{\boldsymbol{x}}} $ is the kinetic energy operator of the nuclei, and $ {H_0}({\mathit{\boldsymbol{r}}};{\mathit{\boldsymbol{x}}}) $ is the electronic Hamiltonian depending on the nuclear coordinates. At a specified nuclear geometry $ {\mathit{\boldsymbol{x}}} $, the adiabatic states can be obtained through solving the eigenvalue equation

    where $ {E_i}({\mathit{\boldsymbol{x}}}) $ and $ | {{\phi _i}({\mathit{\boldsymbol{r}}};{\mathit{\boldsymbol{x}}})} \rangle $ are the adiabatic potential energies and wavefunctions, respectively. These adiabatic states form a complete orthonormal basis set for the electronic wavefunction and have been widely used in surface hopping simulations.

  • In the traditional FSSH method, a series of independent trajectories are generated to describe nonadiabatic dynamics. For each trajectory, the initial nuclear coordinates $ {\mathit{\boldsymbol{x}}} $ and momenta p are sampled from certain distributions, the electronic wavefunction $ | {\psi ({\mathit{\boldsymbol{r}}})} \rangle $ is set according to the problem under investigation, and an adiabatic state $ a $ is chosen as the active state accordingly. At each time step, the nuclear dynamics is carried out on the active PES by the classical Newton equation,

    and the time evolution of the electronic wavefunction is given by the TDSE,

    In the adiabatic representation, $ | {\psi ({\mathit{\boldsymbol{r}}}, t)} \rangle $ is linearly expanded as

    and thus Eq.(4) can be reformulated to

    where the electronic Hamiltonian in the adiabatic representation is defined as

    and $ {{\mathit{\boldsymbol{d}}}_{ij}}({\mathit{\boldsymbol{x}}}) $ = $ \langle {{\phi _i}\left( {{\mathit{\boldsymbol{r}}};{\mathit{\boldsymbol{x}}}} \right)} |{\nabla _{\mathit{\boldsymbol{x}}}}| {{\phi _j}\left( {{\mathit{\boldsymbol{r}}};{\mathit{\boldsymbol{x}}}} \right)} \rangle $ is the nonadiabatic coupling vector. As a result, the electronic population of the active state $ a $ is given by

    and its time evolution can be expressed as

    where $ {b_{ai}} $ = $ - 2{\mathop{\rm Re}\nolimits} (\rho _{ai}^*\mathit{\boldsymbol{\dot x}} \cdot {{\mathit{\boldsymbol{d}}}_{ai}}) $ is the flux contribution from state $ i $. Within a time step size $ \Delta t $, the hopping probability from the active state to another state $ i $ is given by

    A uniform random number $ \xi $$ \in $$ [0, 1] $ is generated. If $ {{\sum\limits_{i = 1}^{j - 1}} {{g_{ai}}} } $$ < $$ \xi $$ \le $$ {\sum\limits_{i = 1}^j {{g_{ai}}}} $, a surface hop from state $ a $ to $ j $ is expected to happen and the nuclear velocities are adjusted along the direction of $ {{\mathit{\boldsymbol{d}}}_{aj}}\left( {{\mathit{\boldsymbol{x}}}(t)} \right) $ to maintain the total energy. If this cannot be realized, the hopping event is cancelled, and the system remains on the original active state. When the predefined criterion is achieved, the surface hopping simulation stops, and the active states along all trajectories are utilized together to analyze the results of nonadiabatic dynamics.

  • BCSH differs from the standard FSSH in two major aspects: (1) the TDSE is modified to cover the relative phase difference of WPs on different PESs and (2) the electronic wavefunction coefficients are reset according to the judgement of WP reflection [41]. In detail, following the study of Shenvi and coworkers [42, 43], the electronic Hamiltonian in Eq.(7) is replaced by an effective Hamiltonian for phase correction

    where $ {m_\alpha } $ is the mass of the $ \alpha $th nuclear degree of freedom, $ {{\mathit{\boldsymbol{p}}}^\alpha} $ and $ {\mathit{\boldsymbol{p}}}_i^\alpha $ are the corresponding momenta of nuclear WPs on the active and $ i $th PESs under the parallel momenta and energy conservation assumptions. Namely, we assume

    where $ \eta $ is a positive number to be solved. For cases violating Eq.(12) and Eq.(13), $ \eta $ is set to be an infinitesimally small number. Note that Eq.(11) naturally goes back to Eq.(7) when the nuclear kinetic energy is sufficiently high. In a classical manner, the nuclear momentum of the WP on PES $ i $ at time $ t $+$ \Delta t $ can be calculated as

    where $ {\mathit{\boldsymbol{F}}}_i $ is the corresponding adiabatic Hellmann-Feynman force. These WP momenta are used as auxiliary variables for decoherence correction. If the dot products $ {\mathit{\boldsymbol{F}}}_i $($ t $)$ \cdot $$ {\mathit{\boldsymbol{p}}}_i $($ t $) and $ {\mathit{\boldsymbol{F}}}_i $($ t $)$ \cdot $$ {\mathit{\boldsymbol{p}}}_i $($ t $+$ \Delta t $) have different signs, a WP reflection event is supposed to take place on the $ i $th PES during the time interval $ \Delta t $. When the active WP is reflected, we collapse the wavefunction to the active state through resetting the coefficients to $ {c'_{i( \ne a)}} $ = 0 and $ {c'_a} $ = $ {c_a}/| {{c_a}} | $. Otherwise, when the WP on a non-active surface $ j $ is reflected, its population is removed and redistributed to all other states through setting $ {c'_j} $ = 0 and $ {c'_{i( \ne j)}} $ = $ {c_i}/{\left(\sum\limits_{k \ne j} {{{| {{c_k}} |}^2}} \right)^{1/2}} $. If no WP reflection is identified, nothing will happen.

  • In the literature, the energy-based decoherence time has been widely utilized due to its great simplicity and high computational efficiency [30, 31]. This approach was originally proposed by Zhu, Truhlar and co-workers to characterize the decay of mixing [25-28], that is, the tendency of state populations evolves to a statistical mixture of pure states. After the further simplification by Granucci and Persico, the decoherence time between the active state $ a $ and a nonactive state $ i $ is simply expressed as [29]

    where $ E_{ \rm{kin}} $ is the total kinetic energy of the nuclei, and $ A $ is an empirical parameter which is normally set to be 0.1 Hartree. When $ E_{ \rm{kin}} $ is infinitely large, Eq.(15) reduces to $ {\tau _{ai}} $ = $ \hbar /|{E_a} - {E_i}| $, which can be justified by the time-energy uncertainty relation [25]. Apparently, besides the nuclear kinetic energy, the decoherence time in Eq.(15) only relies on the energy difference between adiabatic states, which is easy to obtain through ab inito electronic structure calculations. In this approach, the wavefunction coefficients for all nonactive states are corrected according to $ {c'_{i( \ne a)}} $ = $ {c_i}\exp \left( - \Delta t/{\tau _{ai}}\right) $ at each time step and the amplitude of the active state is set as $ {c'_a} $ = $ {c_a}{\left(1 - {\sum\limits_{i \ne a} {| {{c'_i}} |} ^2}\right)^{1/2}}/| {{c_a}} | $ to conserve the total population.

  • The performance of FSSH, BCSH, and the decoherence-corrected FSSH by Eq.(15) has been systematically benchmarked in a series of diverse scattering models [41]. While BCSH generally reproduces the exact solutions, the widely used energy-based decoherence correction with Eq.(15) does not work well in terms of accuracy. In the picture of BCSH, WP reflection induces trajectory branching and results in strong decoherence [41]. Apparently, WP reflection is easier to take place on surfaces with high potential energies. Thereby, it is reasonable that the decoherence time is inverse proportional to the energy gap between the relevant states as in Eq.(15). However, there exist severe problems regarding the dependence on the nuclear kinetic energy. With the increase of $ E_{ \rm{kin}} $, the decoherence time given by Eq.(15) keeps reducing and approaches the limit of $ {\tau _{ai}} $ = $ \hbar /| {{E_a} - {E_i}} | $, implying that decoherence is still significant even with very high kinetic energies. Unfortunately, this is not correct because WPs on different adiabats propagate very similarly and WP reflection cannot happen in this limit.

    Following the picture of BCSH, we propose two new formulas for the energy-based decoherence time. First, we move $ E_{ \rm{kin}} $ from the denominator in Eq.(15) to the numerator and obtain a general linear function of the nuclear kinetic energy,

    where $ A $ and $ B $ are adjustable parameters. Because the decoherence time for any choice of $ E_{ \rm{kin}} $ should be positive, both $ A $ and $ B $ are nonnegative numbers. Eq.(16) is as simple as Eq.(15), and thus also holds clear advantage in efficiency. Now, the decoherence time increases with the nuclear kinetic energy and reduces with the energy gap between adiabatic states, thus agreeing well with the predictions of BCSH [41]. Considering the fact that WP reflection is mostly determined by the relative magnitude of the kinetic energy and the energy gap between adiabatic states, the decoherence time may be better represented as a function of $ {E_{ \rm{kin}}} $/$ | {{E_a} - {E_i}}| $ directly. A possible exponential formula can be expressed as

    where parameter $ A $ determines the relative strength of decoherence correction and parameter $ B $ characterizes the sensitivity of the decoherence time with $ {E_{ \rm{kin}}} $/$ |{{E_a} - {E_i}}| $. While BCSH needs to calculate time-consuming nuclear forces on each PES to judge WP reflection at each time step, both of our new formulas, i.e., Eqs.(16) and (17), only require the adiabatic state energies. Through tuning the parameters $ A $ and $ B $, we may achieve a better balance between accuracy and efficiency and can potentially give more reliable descriptions of nonadiabatic dynamics even in complex systems.

  • We consider FSSH with the effective Hamiltonian in Eq.(11) for phase-corrected wavefunction propagation and the parameterized decoherence time in Eq.(16) or Eq.(17) for decoherence correction. To find out the optimal parameters in the linear and exponential decoherence time formulas, a systematic parameter scan and benchmark process is utilized. First, we select a set of representative models as the training set and make surface hopping simulations with different $ A $ and $ B $ parameters to obtain a detailed error distribution in the two-dimensional parameter space. Then, we select a certain number of parameter candidates with small errors and carry out surface hopping simulations for a model base with a large number of diverse models to make a reliable assessment of the overall behavior. Finally, the parameter set with the highest accuracy for each formula is taken out and its performance is benchmarked in the three standard Tully models. This approach provides a general and reliable routine to optimize different decoherence time formulas for trajectory surface hopping.

  • In order to systematically examine the performance of the two proposed energy-based decoherence schemes and obtain the optimal parameter sets, a carefully designed model base proposed recently is utilized here [41]. It contains a series of one-dimensional two-level scattering models. As shown in FIG. 1(A), 5 different expressions are utilized to construct the diagonal terms of the corresponding diabatic Hamiltonians,

    Figure 1.  Diabatic potential energy surfaces (blue and red solid lines) to construct the model base with two quantum levels and one classical degree of freedom, and adiabatic potential energy surfaces (red and blue solid lines) and the corresponding nonadiabatic couplings (green dotted lines) of Tully's (B) simple avoided crossing, (C) dual avoided crossing, and (D) extended coupling with reflection models.

    where $ A_1 $ = 0.01, $ B_1 $ = 1.6, $ A_2 $ = 0.02, $ B_2 $ = 0.28, $ E_1 $ = 0, $ E_2 $ = 0.01, and $ E_3 $ = 0.02, and $ x $ is the classical nuclear coordinate. Note that atomic units are adopted throughout this study unless otherwise noted. Besides, we consider a Gaussian function to describe the off-diagonal terms of the diabatic Hamiltonians,

    where $ A_3 $ and $ B_3 $ characterize the strength and width of the diabatic coupling, respectively. Four different $ A_3 $ (i.e., 0.0025, 0.005, 0.02, and 0.04) and two different $ B_3 $ (i.e., 0.4 and 1.5) are considered in the present study. Thereby, there exist 25 and 8 combinations for the diagonal and off-diagonal terms of the diabatic Hamiltonians respectively, and thus a total number of 200 models are constructed in our model base [41]. To simplify the discussions below, we use $ M_{ij}^{{A_3}, {B_3}} $ to express the model using $ H_{11}^{[i]} $, $ H_{22}^{[j]} $ and $ {H_{12}} $ with parameters $ A_3 $ and $ B_3 $ to build the system Hamiltonian. In the nonadiabatic dynamics simulations, we place a Gaussian WP on the negative $ x $ side of the lower adiabatic PES with a positive momentum $ k $,

    where $ x_0 $ and $ \sigma _x $ are the central position and width of the nuclear WP, respectively. For each surface hopping trajectory, the initial coordinate and momentum are sampled by the corresponding Wigner distribution,

    where $ p_0 $ = $ \hbar k $ and $ \sigma_p $ = $ \hbar $/(2$ \sigma_x $). For all models investigated in this study, the spatial width $ \sigma_x $ is set as 10/$ k $, and the nuclear mass is 2000. To reduce the statistical error caused by initial sampling and stochastic surface hops, over 20000 trajectories are used for each surface hopping calculation. The exact quantum dynamics are simulated by the discrete variable representation (DVR) method with transmission-free absorbing potentials [44-46]. For each investigated model in the present study, there are four scattering channels, namely, transmission and reflection on the upper and lower surfaces. Many different choices of $ k $ are utilized for all the models, and their exact solutions are used as references to benchmark the surface hopping studies.

  • In the training process, we choose $ A_3 $ = 0.04 and $ B_3 $ = 1.5 in Eq.(28), which correspond to strong and narrow diabatic couplings. These models generally have large energy gaps between adiabatic states, and thus are associated with strong decoherence effect. As model $ M_{45}^{{A_3}, {B_3}} $ suffers from strong trivial crossings [47], which requires much smaller time step size to get converged results in the adiabatic representation, we move it out and consider the remaining 24 models for training. To track the role of the nuclear kinetic energy, we select 5 representative initial momenta (i.e., $ k $ = 10, 15, 20, 25, and 30). Thereby, we need to carry out 120 calculations for each choice of parameters $ A $ and $ B $. The average error for surface hopping simulations is defined as

    where $ P_{ij}^{{\rm{ SH}}} $ and $ P_{ij}^{{\rm{ \rm{DVR}}}} $ are populations of the $ j $th transmission/reflection channel in the $ i $th calculation by surface hopping and DVR, respectively. The time step size for surface hopping and DVR calculations are set as $ \Delta t $ = 0.05 and 0.02, respectively.

    We first study the linear formula for the decoherence time, i.e., Eq.(16). The error distribution in the two-dimensional parameter space is shown in FIG. 2(A). When $ B $ is small, the decoherence time is dominated by $ A $. In this case, the error of surface hopping is generally large due to the weak kinetic energy dependence of the decoherence time. When either $ A $ or $ B $ becomes large, the error is also evident because the decoherence time is generally too large and cannot properly characterize the decoherence effect. Thereby, the error minimum locates at the region with small $ A $ and intermediate $ B $ to achieve a balanced decoherence strength. We choose several candidate parameter sets, which have the smallest errors. They are indicated by Roman numbers I ($ A $ = 0 and $ B $ = 1000), Ⅱ ($ A $ = 0.5 and $ B $ = 1200), and Ⅲ ($ A $ = 1.0 and $ B $ = 1000) in FIG. 2(A).

    Figure 2.  Error distribution for phase-corrected FSSH with (A) linear (Eq.(16)), (B) exponential (Eq.(17)), and (C) the traditional (Eq.(15)) decoherence time formulas in the training set with 120 chosen calculations. The Roman numbers Ⅰ, Ⅱ, and Ⅲ indicate the parameter sets with the smallest errors in (A), while the Roman numbers IV and V indicate the parameter sets with the smallest errors in (B). The minimum error in (C) is highlighted by the black arrow.

    As shown in FIG. 2(B), the exponential decoherence time formula in Eq.(17) is also investigated. Compared with the linear formula in Eq.(16), the overall errors by the exponential formula are generally close. When $ A $ is large, the decoherence correction is generally too weak due to the exponential increase of the decoherence time. If both $ A $ and $ B $ are small, however, the decoherence correction is too strong, and thus also results in large errors. Thereby, the parameter sets with the highest accuracy are associated with small $ A $ and large $ B $, as highlighted by Roman numbers IV ($ A $ = 25 and $ B $ = 42) and V ($ A $ = 30 and $ B $ = 42) in FIG. 2(B).

    In a similar manner, we can also scan the parameter of the traditional energy-based decoherence time formula in Eq.(15). To prevent negative decoherence time, only non-negative values for parameter $ A $ are considered here. As shown in FIG. 2(C), the error monotonously reduces with increasing $ A $ when $ A $ is small. The minimum error is achieved with $ A $ = 35, which is much larger than the widely adopted value of $ A $ = 0.1 [26, 29]. Besides, this error is very close to that with infinitely large $ A $, which corresponds to the case without any decoherence correction. Thereby, using Eq.(15) for decoherence correction has intrinsic problems even after the parameter is properly adjusted.

  • A systematic benchmark of the candidate decoherence time setups (i.e., Ⅰ, Ⅱ, Ⅲ, Ⅳ, and Ⅴ in FIG. 2) is achieved through calculating the average population error for the whole model base, which includes 200 models and 26 different momenta (i.e., $ k $ = 5, 6, $ \ldots $, 30) per model. Thereby, for each candidate parameter set, we need to perform 5200 individual calculations. In these calculations, a time step size of $ \Delta t $ = 0.05 is used except for models suffering from strong trivial crossings [48-53], where a smaller $ \Delta t $ in the range of 0.0005-0.01 is used instead. Similar to Eq.(31), the average error of the model base for each individual $ k $ is defined as

    where $ P_{ij}^{ \rm{SH}} $ and $ P_{ij}^{ \rm{DVR}} $ are populations of the $ j $th transmission/reflection channel in the $ i $th model by surface hopping and DVR, respectively.

    In FIG. 3(A), the performance of different surface hopping approaches is studied. Compared with the traditional FSSH without any decoherence correction, FSSH with the energy-based decoherence time in Eq.(15) reduces the error slightly for small $ k $ but induces larger error for large $ k $. In comparison, our new decoherence corrections with both Eq.(16) and Eq.(17) systematically reduce the computational error in all cases. However, the performance of the two formulas behaves differently in different regions. Namely, the linear formula works slightly better for $ k $$ < $12, while the exponential formula is superior for $ k $$ > $12. Actually, the error of the exponential formula is almost negligible for large $ k $ because its decoherence time increases quickly with the kinetic energy and decoherence correction is not needed in the high energy limit [41]. Encouragingly, the results of our linear formula are quite close to those of the AFSSH method [41], which is widely recognized as a rigorous decoherence scheme in the literature [35]. In general, the chosen parameter sets show similar results. For the linear formula, $ A $ = 0.5 and $ B $ = 1200 give smaller errors for large $ k $ and thus are chosen as the optimal parameters. For the exponential formula, the errors of parameter sets IV and V are very close and we use $ A $ = 25 and $ B $ = 42 in the following calculations.

    Figure 3.  (A) Average population error as a function of the initial nuclear momentum by FSSH and phase-corrected FSSH with Eq.(15), Eq.(16), and Eq.(17) for the decoherence time. The parameters $ A $ and $ B $ are shown as indicated. The 200 models in the model base with 26 different initial momenta and 4 transmission/reflection channels per model are considered. (B) Distribution of the corresponding population errors for surface hopping calculations with the chosen decoherence correction approaches as indicated.

    In FIG. 3(B), we show the error distribution for the three energy-based decoherence correction approaches with Eqs.(15), (16), and (17). In each of the 5200 individual calculations, we have 4 transmission or reflection channels. Thereby, we obtain a total number of 20800 population errors using the exact quantum results as references. It is evident that our two new decoherence correction approaches exhibit significantly narrower error distribution than the traditional method. For both linear and exponential formulas, less than 0.8% of the errors are larger than 0.1. In contrast, this ratio is about 8% for the traditional formula. Consequently, the new decoherence time formulas are superior regarding both the averaged population error and the error distribution.

  • The two new decoherence time formulas with their optimal parameter sets are further benchmarked with the three standard model systems proposed by Tully [15]. We set $ \Delta t $ as 0.1 and 0.01 for surface hopping and DVR, respectively. The first system is the simple avoided crossing model (see FIG. 1(B)) defined by

    where $ A $ = 0.01, $ B $ = 1.6, $ C $ = 0.005, and $ D $ = 1. In FIG. 4, we show the calculated $ k $-dependent populations of the four transmission/reflection channels with different surface hopping methods. FSSH agrees well with the exact quantum solutions because the decoherence effect in this model is very weak. However, FSSH with the widely used energy-based decoherence time formula in Eq.(15) exhibits evident deviations in the transmission on the lower and upper surfaces for 10$ < $$ k $$ < $22, suggesting that the decoherence correction is too strong. In contrast, the results of both Eq.(16) and Eq.(17) match well with the exact results despite their different formulas and parameters.

    Figure 4.  Transmission populations on (A) the lower and (B) the upper surfaces, reflection populations on (C) the lower and (D) the upper surfaces for the simple avoided crossing model of Tully. Open circles are exact quantum solutions by DVR. The results of FSSH are shown as green solid squares. The phase-corrected FSSH with Eq.(15), Eq.(16) ($ A $ = 0.5 and $ B $ = 1200), and Eq.(17) ($ A $ = 25 and $ B $ = 42) for the energy-based decoherence correction are shown as blue solid triangles, light blue diamonds, and red solid circles, respectively.

    The second system is the dual avoided crossing model (see FIG. 1(C)), which is defined as

    with $ E_0 $ = 0.05, $ A $ = 0.1, $ B $ = 0.28, $ C $ = 0.015, and $ D $ = 0.06. As shown in FIG. 5, FSSH with Eq.(15) shows large errors for $ k $$ > $17 and the oscillations in FSSH and DVR results are smeared out, which means that its decoherence strength is strongly overestimated. For the reflection channel on the lower surface, the results do not improve significantly and the errors are even larger than FSSH for $ k $$ < $10. Our two new decoherence time formulas give a more reasonable relation with the kinetic energy, and thus the results of both methods agree well with the exact quantum solutions. Especially, the reflection probabilities on the lower surface are remarkably improved. The results of the exponential decoherence time formula are almost on top of those of DVR.

    Figure 5.  Transmission populations on (A) the lower and (B) the upper surfaces, reflection populations on (C) the lower and (D) the upper surfaces for the dual avoided crossing model of Tully. The symbols are the same as those in FIG. 4.

    The third system is the extended coupling with reflection model (see FIG. 1(D)), in which the Hamiltonian reads

    where $ A $ = 6$ \times $10$ ^{-4} $, $ B $ = 0.1, and $ C $ = 0.9. The corresponding channel populations are shown in FIG. 6. Although FSSH reproduces the transmission channels, the reflection populations on the lower and upper surfaces exhibit obvious oscillations for small $ k $ due to the absence of decoherence correction. Note that much stronger oscillations are present if the initial sampling of the Wigner distribution is not carried out [15]. The problem is not resolved after the application of the traditional energy-based decoherence correction through Eq.(15), and the description of transmission even becomes worse for larger $ k $. In comparison, the results of our linear and exponential formulas agree very well with the exact results in all channels.

    Figure 6.  Transmission populations on (A) the lower and (B) the upper surfaces, reflection populations on (C) the lower and (D) the upper surfaces for the extended coupling with reflection model of Tully. The symbols are the same as those in FIG. 4.

    Thereby, we have confirmed that the proposed two new decoherence time formulas with the chosen parameter sets do improve significantly the reliability of trajectory surface hopping simulations while holding the high simplicity and efficiency. This can only be achieved when the correct decoherence strength is captured in real time simulations, therefore justifying our modification of the energy-based decoherence time formulas. There are a few points that worth discussing. (i) In principle, the kinetic energy along the NAC direction should be used instead of the total kinetic energy for decoherence time calculations [25-28]. This needs to be considered in real applications with many classical degrees of freedom. However, it is not an issue in the present study because we focus on one-dimensional systems here. (ii) To reduce possible overfitting, the associated parameters in the proposed decoherence time formulas need to be further scanned or trained in more complex systems with more electronic levels and higher dimensions. Better formulas for the energy-based decoherence time may also exist and can be studied with similar ways as proposed in this work. (iii) Since the neural network can generally describe any nonlinear function [54], it could possibly characterize better the complex relation between the decoherence time and related quantities along surface hopping trajectories. Relevant studies are currently under way.

  • In summary, we have proposed linear and exponential formulas of the energy-based decoherence time for trajectory surface hopping. These two new approaches have been designed to reproduce the strong decoherence effect due to WP reflection as predicted by BCSH. We have utilized a general strategy to scan and obtain the optimal parameters for these decoherence time formulas. As demonstrated in hundreds of diverse models and the three standard scattering models of Tully, the new decoherence correction approaches have captured the proper strength of decoherence correction and shown both high efficiency and reliability, and they thus deserve further development for application in surface hopping simulations of complex nonadiabatic dynamics. By comparison, the linear formula exhibits higher accuracy for small nuclear momenta, while the exponential formula behaves better in the large momentum region.

  • This work was supported by the National Natural Science Foundation of China (No.21922305, No.21873080, and No.21703202).

Reference (54)

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return