Quantum Trajectory Mean-Field Method for Nonadiabatic Dynamics in

Subscriber access provided by Nottingham Trent University

Feature Article

Quantum Trajectory Mean-Field Method for Nonadiabatic Dynamics in Photochemistry Lin Shen, Diandong Tang, Bin Bin Xie, and Wei-Hai Fang J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b03480 • Publication Date (Web): 02 Aug 2019 Downloaded from pubs.acs.org on August 10, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Quantum Trajectory Mean-Field Method for Nonadiabatic Dynamics in Photochemistry Lin Shen1, Diandong Tang1, Binbin Xie2 and Wei-Hai Fang1,* 1. Key Laboratory of Theoretical and Computational Photochemistry of Ministry of Education, College of Chemistry, Beijing Normal University, Beijing 100875, P. R. China 2. Hangzhou Institute of Advanced Studies, Zhejiang Normal University, 1108 Gengwen Road, Hangzhou 311231, Zhejiang, P. R. China

ABSTRACT The mixed quantum-classical (MQC) dynamical approaches have been widely used to study nonadiabatic phenomena in photochemistry and photobiology, in which the time evolutions of the electronic and nuclear subsystems are treated based on quantum and classical mechanics, respectively. The key issue is how to deal with coherence and decoherence during the propagation of the two subsystems, which has been the subject of numerous investigations for a few decades. A brief description on Ehrenfest mean field (MEF) and surface-hopping (SH) methods is firstly provided and then different algorithms for treatment of quantum decoherence are reviewed in the present paper. More attentions were paid to quantum trajectory mean-field (QTMF) method under the picture of quantum measurements, which is able to overcome the overcoherence problem. Furthermore, the combined QTMF and SH algorithm is proposed in the present work, which takes advantages of the QTMF and SH methods. The potential to extend the applicability of the QTMF method was briefly discussed, such as the generalization to other type of nonadiabatic transitions, the combination with multiscale computational models, and possible improvements on its accuracy and efficiency by using machine learning techniques.

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

I. INTRODUCTION Understanding the mechanism of photophysical and photochemical processes in multiple electronic states is a challenging task because of the coupling of nuclear and electronic motion in various cases, such as spectroscopy and dynamics,1-3 electron and energy transfer,4-7 scattering at metal surfaces,8-12 and radiationless transitions near conical intersections.13-18 In the ab initio molecular dynamics (MD) simulation, the motion of nuclei is described classically on an adiabatic potential energy surface (PES) that can be obtained by using ab initio quantum chemistry. Despite its great success to study thermochemical reactions, for photochemical processes classical MD simulation is of problematic basis. The adiabatic or Born-Oppenheimer (BO) approximation, which assumes the separation of nuclear and electronic motion, is unable to capture the nonadiabatic dynamics of systems in photochemistry. In the vicinity of conical intersections, for example, electrons no longer adjust to the change of nuclear positions instantaneously. Instead, the change of electronic state is induced by the motion of nuclei, and the nuclei should also evolve in response to quantum transition of electrons. In such a case, the feedback between the nuclear and electronic motions must be considered and the BO approximation fails. The time evolution of a system including nuclei and electrons is governed by the laws of quantum mechanism (QM). A fully QM treatment should be employed in principle.19-24 However, such a simulation on a realistic molecular system is extremely expensive computationally with the increase of degrees of freedom, if not impossible. The global high-dimensional PESs in multiple electronic states are also required prior to propagation, which further hampers the applicability of the fully QM dynamic simulations. As an alternative, the mixed quantum-classical (MQC) MD approach provides a powerful tool to allow a reliable quantum treatment on the degrees of 2


Page 2 of 47

Page 3 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


freedom for electrons and an effective classical treatment on the degrees of freedom for nuclei. During the MQC MD simulations, the electrons are described as quantum particles based on the time-dependent Schrödinger equation (TDSE), while the nuclei are treated as classical particles usually based on the Newton’s second law. In the large family of MQC MD approaches, a variety of theoretical models and practical algorithms have been developed over the last 40 years.25-28 A proper MQC MD method should satisfy four conditions. First, the time evolution on the nuclear and electronic degrees of freedom must be self-consistent. In other words, the nuclear and electronic motion during MD simulation must respond to each other appropriately. Second, the trajectory moving in the region without nonadiabatic effect must evolve as the same as the classical MD evolution. Third, the MQC MD on a molecular system can be implemented on-the-fly. In other words, the construction on the global high-dimensional PES must be unnecessary. Finally, the nonlocal effects of the nuclei, which have been lost because of the classical treatments, can be compensated at least approximately. In this paper, we will summarize two standard MQC MD methods, Ehrenfest mean field (EMF) and surface hopping (SH), in Section II. After a brief review on some advanced approaches incorporating decoherence in Section III, we will introduce the recently developed quantum trajectory mean-field (QTMF) method and propose a variant combined with SH in Sections IV and V, respectively, followed by an outlook in Section VI. II. EHRENFEST MEAN FIELD AND SURFACE HOPPING METHOD The Ehrenfest mean field and surface hopping methods are both well-established to simulate nonadiabatic processes. Taking a factorization or a multi-configuration expansion of the total molecular wavefunctions, the time evolution on nuclear degrees of freedom can be reduced to classical motions approximately, which leads to the 3



Page 4 of 47

equation of motions (EOMs) for nuclei in the EMF or SH, respectively. The EOMs for electrons can be also derived from the TDSE at the classical limit on the nuclei. A short summary is presented here. More details can be seen from Tully’s paper.29 Ehrenfest Mean Field Method The time evolution of a molecule is governed by TDSE as ∂Ψ

(1)

𝑖ℏ ∂𝑡 = 𝐻𝛹 where 𝛹 is the total molecular wavefunction, and 𝐻 is the full Hamiltionian as ℏ2

ℏ2

ℏ2

𝐻 = ― ∑𝛼2𝑀𝛼∇2𝑅𝛼 ― ∑𝛽2𝑚𝛽∇2𝑟𝛽 +𝑉(𝐫,𝐑) = ― ∑𝛼2𝑀𝛼∇2𝑅𝛼 + 𝐻𝑒𝑙(𝐫,𝐑)

(2)

where α denotes nuclei, β denotes electrons, R, r, Mα and mβ are the coordinates of nuclei and electrons, and the masses of nuclei and electrons, respectively, and 𝑉(𝐫,𝐑) includes all interactions between particles. The total molecular wavefunction can be represented as a factorization 𝑖 𝑡 ∫0𝐸𝑒𝑙(𝜏)d𝜏

𝛹(𝐫,𝐑,𝑡) = 𝛷(𝐫,𝑡)𝜒(𝐑,𝑡)𝑒ℏ

(3)

Here 𝛷(𝐫,𝑡) and 𝜒(𝐑,𝑡) are the wavefunctions of electrons and nuclei, respectively, and Eel is the phase factor that is defined as 𝐸𝑒𝑙 = ∫𝛷 ∗ (𝐫,𝑡)𝜒 ∗ (𝐑,𝑡)𝐻𝑒𝑙𝜒(𝐑,𝑡)𝛷(𝐫,𝑡)d𝐫d𝐑

(4)

The phase factor of 𝛷(𝐫,𝑡) and 𝜒(𝐑,𝑡) can be specified as 𝑖ℏ∫𝛷 ∗ (𝐫,𝑡)

∂𝛷(𝐫,𝑡) ∂𝑡

d𝐫 = 𝐸𝑒𝑙

(5)

and 𝑖ℏ∫𝜒 ∗ (𝐑,𝑡)

∂𝜒(𝐑,𝑡) ∂𝑡

d𝐑 = 𝐸

(6)

where E is the expectation value of the full Hamiltonian on the total wavefunction in Eq (1). Applying Eq (3) to (1), we obtain 𝑖ℏ

∂𝛷(𝐫,𝑡) ∂𝑡

ℏ2

= ― ∑𝛽2𝑚𝛽∇2𝑟𝛽𝛷(𝐫,𝑡) + [∫𝜒 ∗ (𝐑,𝑡)𝑉(𝐫,𝐑)𝜒(𝐑,𝑡)d𝐑]𝛷(𝐫,𝑡) 4


(7)

Page 5 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and 𝑖ℏ

∂𝜒(𝐑,𝑡) ∂𝑡

ℏ2

= ― ∑𝛼2𝑀𝛼∇2𝑅𝛼𝜒(𝐑,𝑡) + [∫𝛷 ∗ (𝐫,𝑡)𝐻𝑒𝑙𝛷(𝐫,𝑡)d𝐫]𝜒(𝐑,𝑡)

(8)

Under the picture of MQC, the nuclear degrees of freedom are treated classically. Hereby 𝜒(𝐑,𝑡) in Eq (7) can be replaced by a delta function at 𝐑(𝑡), leading to the EOM of electronic motion as 𝑖ℏ

∂𝛷(𝐫,𝐑(𝑡), 𝑡) ∂𝑡

= 𝐻𝑒𝑙𝛷(𝐫,𝐑(𝑡),𝑡)

(9)

On the other hand, 𝜒(𝐑,𝑡) in Eq (8) can be represented using 𝑖

𝜒(𝐑,𝑡) = 𝐴(𝐑,𝑡)𝑒

ℏ𝑆(𝐑,𝑡)

(10)

Taking the classical limit as ℏ→0, the Hamiltonian-Jacobi equation is then obtained ∂𝑆(𝐑,𝑡) ∂𝑡

1

+ ∑𝛼2𝑀𝛼[∇𝑅𝛼𝑆(𝐑,𝑡)]2 + ∫𝛷 ∗ (𝐫,𝑡)𝐻𝑒𝑙𝛷(𝐫,𝑡)d𝐫 = 0

(11)

Define the classical momentum of nuclei as 𝑝𝛼 = ∇𝑅𝛼𝑆, the above equation is rewritten as 𝑀𝛼𝑅𝛼 =

d𝑝𝛼 d𝑡

= ― ∇𝑅𝛼[∫𝛷 ∗ (𝐫,𝐑(𝑡),𝑡)𝐻𝑒𝑙𝛷(𝐫,𝐑(𝑡),𝑡)d𝐫]

(12)

It indicates that the nuclear motion is governed by the Newton’s second law in the average potential generated by 𝛷(𝐫,𝐑,𝑡). The electronic wavefunction 𝛷(𝐫,𝐑,𝑡) is expanded as a linear combination of basis sets as 𝛷(𝐫,𝐑,𝑡) = ∑𝑗𝑐𝑗(𝑡)𝜑𝑗(𝐫,𝐑)

(13)

where j denote the jth basis function. Without loss of generality, we choose a set of wavefunctions at adiabatic electronic states as basis functions, which satisfies 𝐻𝑒𝑙𝜑𝑗(𝐫,𝐑) = 𝐸𝑗𝜑𝑗(𝐫,𝐑)

(14)

Based on Eqs (9), (13) and (14), the coefficient cj evolves as follows 𝑖

𝑐𝑗 = ― ℏ𝑐𝑗𝐸𝑗 ― ∑𝑘(𝐑 ∙ 𝒅𝑗𝑘)𝑐𝑘 and the EOM of density matrix is 5


(15)


𝑖

𝜌𝑗𝑘 = ― ℏ𝜌𝑗𝑘(𝐸𝑗 ― 𝐸𝑘) ― ∑𝑙𝐑 ∙ (𝒅𝑗𝑙𝜌𝑙𝑘 ― 𝜌𝑗𝑙𝒅𝑙𝑘)

Page 6 of 47

(16)

where djk is the derivative nonadiabatic coupling between two electronic states j and k defined as 𝒅𝛼𝑗𝑘(𝐑) = ∫𝜑𝑗∗ (𝐫,𝐑)∇𝑅𝛼𝜑𝑘(𝐫,𝐑)d𝐫

(17)

We rewrite Eq (12) as 𝑀𝛼𝑅𝛼 = ― ∇𝑅𝛼[∑𝑗𝜌𝑗𝑗𝐸𝑗] = ― ∑𝑗𝜌𝑗𝑗∇𝑅𝛼𝐸𝑗 + ∑𝑗,𝑘𝒅𝛼𝑗𝑘𝜌𝑘𝑗(𝐸𝑗 ― 𝐸𝑘)

(18)

During the MQC MD simulations using EMF, the nuclear motion is governed by Eq (18), and the electronic degrees of freedom propagate according to Eq (16). The nonadiabatic effects are involved intrinsically in the EOMs for nuclei and electrons. Fewest-Switches Surface Hopping Method Unlike the derivation of EMF, the total wavefunction can be represented as a multi-configuration expansion as follows 𝛹(𝐫,𝐑,𝑡) = ∑𝑗𝜑𝑗(𝐫,𝐑)𝜒𝑗(𝐑,𝑡)

(19)

where 𝜑𝑗(𝐫,𝐑) can be selected according to Eq (14), and 𝜒𝑗(𝐑,𝑡) describes the nuclear degrees of freedom when the electrons evolve in state j. Based on the TDSE, 𝑖ℏ

∂𝜒𝑗(𝐑,𝑡) ∂𝑡

ℏ2

ℏ2

= ― ∑𝛼2𝑀𝛼∇2𝑅𝛼𝜒𝑗(𝐑,𝑡) ― ∑𝑘∑𝛼𝑀𝛼𝒅𝑗𝑘 ∙ ∇𝑅𝛼𝜒𝑘(𝐑,𝑡) ℏ2

∑𝑘∑𝛼2𝑀 [∫𝜑𝑗∗ (𝐫,𝐑)∇2𝑅𝛼𝜑𝑘(𝐫,𝐑)d𝐫]𝜒𝑘(𝐑,𝑡) + 𝐸𝑗𝜒𝑗(𝐑,𝑡) 𝛼

― (20)

where djk and Ej are defined above. The nuclear wavefunctions are expressed similarly to Eq (10) as 𝑖

𝜒𝑗(𝐑,𝑡) = 𝐴𝑗(𝐑,𝑡)𝑒ℏ

𝑆𝑗(𝐑,𝑡)

(21)

Therefore, ∂𝑆𝑗(𝐑,𝑡) ∂𝑡

1

+ ∑𝛼2𝑀𝛼[∇𝑅𝛼𝑆𝑗(𝐑,𝑡)]2 + 𝐸𝑗 = 0

and 6


(22)

Page 7 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∂𝐴𝑗(𝐑,𝑡) ∂𝑡

𝑖

+ ∑𝑘[𝐑 ∙ 𝒅𝑗𝑘]𝐴𝑘(𝐑,𝑡)𝑒ℏ

[𝑆𝑘(𝐑,𝑡) ― 𝑆𝑗(𝐑,𝑡)]

+

1

∑𝛼2𝑀 [𝐴𝑗(𝐑,𝑡)∇2𝑅 𝑆𝑗(𝐑,𝑡) + 2∇𝑅𝛼𝐴𝑗(𝐑,𝑡) ∙ ∇𝑅𝛼𝑆𝑗(𝐑,𝑡)] = 0 𝛼 𝛼

(23)

Here the independent trajectory approximation that 𝑝𝛼 = ∇𝑅𝛼𝑆𝑘 regardless of k is employed. Compared with EMF, Eq (22) suggests that nuclei travel classically on one PES in a single adiabatic electronic state rather than the average potential, that is, (24)

𝑀𝛼𝑅𝛼 = ― ∇𝑅𝛼𝐸𝑗

where j denotes the adiabatic electronic state at the current time step. This is the first point of SH. On the other hand, the EOM of the electronic degrees of freedom is not so evident. Curchod et al. introduced a quantum amplitude 𝐴QM 𝑗 (𝐑(𝑡),𝑡), a quantum phase 𝑆QM 𝑗 (𝐑(𝑡),𝑡) and 𝑖 QM 𝑆𝑗 (𝐑(𝑡),𝑡)

ℏ 𝑐𝑗(𝐑(𝑡),𝑡) = 𝐴QM 𝑗 (𝐑(𝑡),𝑡)𝑒

(25)

Replacing Aj and Sj by 𝐴QM and 𝑆QM in Eqs (22) and (23) and ignoring the gradients 𝑗 𝑗 30 of 𝐴QM and 𝑆QM 𝑗 𝑗 , the time evolution of cj reaches the same formula as Eq (15).

Hereby the second point of SH is to share the EOM of the density matrix in EMF as Eq (16). The classical trajectory travels on a single PES and picks up the quantum amplitude of each electronic state, and the nonadiabatic couplings are related to the contribution to quantum interference from other states. It is quite different from the EMF. When the trajectory goes to the weak intersection region, it is reduced to the adiabatic classical MD and rare transitions occur during the simulation with the SH. Using the EMF, however, the nuclear motion keeps on the average PES and thus suffers from the inappropriate combination of adiabatic states. When the trajectory goes around the conical interaction, it just jumps between PESs and still evolves adiabatically between two hops using the SH, while the mean-field approximation in the EMF makes more sense. Both the EMF and SH nonadiabatic simulations may 7



Page 8 of 47

encounter the numerical problem arising from the very large derivative nonadiabatic coupling near a conical intersection, but it is rare for most of trajectories in practice. The final point of SH is an ensemble of independent trajectories. Each trajectory can switch from one PES to another during the time evolution. Switching or so-called surface hopping occurs in infinitesimal time stochastically with a probability extracted from the electronic degrees of freedom. The population of trajectories on state j must be equal to 𝜌𝑗𝑗 as 𝑁𝑗 = 𝜌𝑗𝑗𝑁, where N is the total number of trajectories, and Nj is the occupation number in state j. In the most famous fewest-switches surface hopping (FSSH) method developed by Tully,31-32 the fewest number of nonadiabatic transitions can be obtained with the assumption that 𝑃𝑗(𝑡,𝑡 + d𝑡) =

𝑁𝑗(𝑡) ― 𝑁𝑗(𝑡,𝑡 + d𝑡)

(26)

𝑁𝑗(𝑡)

where 𝑃𝑗(𝑡,𝑡 + d𝑡) is the probability that a trajectory moving in state j at time t switches to any other state between 𝑡 and 𝑡 + d𝑡. We can rewrite Eq (26) as 𝑃𝑗(𝑡,𝑡 + d𝑡) = ―

𝜌𝑗𝑗(𝑡)d𝑡

(27)

𝜌𝑗𝑗(𝑡)

Applying Eq (16) to (27), the transition probability from state j to k during a finite time interval Δ𝑡 can be obtained as

[

𝑃𝑗→𝑘(𝑡,𝑡 + Δ𝑡) = 𝑚𝑎𝑥 0,

2𝑅𝑒(𝑹 ∙ 𝒅𝑗𝑘𝜌𝑗𝑘∗ )d𝑡′

𝑡 + Δ𝑡

∫𝑡

𝜌𝑗𝑗

]

(28)

During the above integration we assume no hop occurs. It makes the population of trajectories on state j unchanged, so does 𝜌𝑗𝑗 according to 𝑁𝑗 = 𝜌𝑗𝑗𝑁. Overall, a swarm of trajectories are employed and propagate independently on a single PES. The nuclear and electronic degrees of freedom of each trajectory progress over time with Eqs (24) and (16), respectively. The trajectory can switch to any other electronic state stochastically at each time step with a transition probability in Eq (28). After surface hopping occurs, the velocities of nuclei are rescaled to maintain the total 8


Page 9 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


energy conservation. In the case of energy forbidden, the hop should be aborted. The SH algorithms such as the FSSH have been applied to study different photochemical and photobiological systems,33-36 varying from small radicals in atmospheric chemical mechanisms to large proteins for diverse bioluminescence phenomena.

III. DECOHERENCE Despite the success and popularity of the EMF and SH MD simulations, the drawbacks on the two methods have been also observed in a variety of situations for decades, especially the overcoherence problem that directly affects the accuracy of nonadiabatic dynamics simulations on photochemical reactions.37 This problem was pointed out in the first paper on FSSH by Tully,31 where the MQC MD simulations fails on the third example so-called extended coupling Hamiltonian with reflection. As shown in Figure 1, a wavepacket entering in the coupling region will spawn into two wavepackets. Because the upper surface is repulsive, the wavepacket on the upper surface will reflect at a low momentum, while that on the lower surface can transmit, resulting in a separation between two wavepackets. It is denoted as the decoherence effect. Schwartz and coworkers revealed the origin of the overcoherence problem based on the comparison between the total density matrix and the electronic density matrix.38 With the classical treatments on nuclei, the total density matrix is reduced to the electronic density matrix, and the contributions of the nuclear degrees of freedom to the decoherence effect are lost. The overestimation on the coherence effects is thus an inevitable result in the EMF and FSSH simulations. It is a general and practical disadvantage that has been observed in applications. For example, the relaxation time for photochemical processes is sometimes underestimated during FSSH simulations without decoherence correction.39-42 In the absence of nuclear wavefunctions in the MQC models, the electrons can 9



Page 10 of 47

be considered as a closed system that obeys the TDSE. The off-diagonal elements of the electronic density matrix must decay to zero to represent decoherence when the nonadiabatic effect can be ignored.43-44 Tully introduced a damping term to the EOM of the electronic density matrix and changed Eq (16) to 𝑖

𝜌𝑗𝑘 = ― ℏ(𝐸𝑗 ― 𝐸𝑘) ― ∑𝑙𝐑 ∙ (𝒅𝑗𝑙𝜌𝑙𝑘 ― 𝜌𝑗𝑙𝒅𝑙𝑘) ―γ(1 ― 𝛿𝑗𝑘)𝜌𝑗𝑘

(29)

where 𝛾 is the empirical decoherence rate.31 In order to determine this parameter, a straightforward idea is to introduce some auxiliary variables or wavepackets to mimic such a quantum nonlocal effect on nuclei. A large number of ad hoc algorithms have been reported over the past years.27-28, 37 One class of these approaches is on the basis of potential energies in electronic states. As proposed by Zhu, Truhlar and coworkers, the decoherence rate is a function of the potential energy difference between two PESs and the classical kinetic energy of nuclei.45-48 It is inconsistent with the intuition that the decoherence rate should depends on the force difference between two PESs rather than energy difference since the separation between two wavepackets moving on lower and upper surfaces is mainly due to different forces. The second class of these approaches is derived based on the overlap between wavepackets. Prezhdo and coworkers obtained the decoherence rate from the overlap evolution of frozen Gaussian wavepackets, which is proportional to the difference in forces.39, 49 Granucci et al. developed a SH method with overlap-based corrections, in which frozen Gaussians are assigned to each electronic state other than the current one on which the classical trajectory travels, and the electronic population can be modified based on the overlap.50 Subotnik et al. proposed the augmented FSSH algorithm using the propagation of auxiliary variables,51 followed by some developments as simultaneous FFSH by Shenvi et al., in which spawned Gaussians on other electronic states are employed during the time evolution to calculate the decoherence rate.52 Akimov et al. 10


Page 11 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reported the coherence penalty functional method.53 Using the estimated decoherence rate, a penalty term is added to Hamiltonian to suppress coherence. The third class of these approaches is some schemes that combine EMF and SH together, such as the MF with SH in which the MF and SH trajectories are propagated simultaneously and compared with each other to capture the coherence loss,54 the MF with stochastic decoherence in which the decoherence rate is employed to stochastically determine whether the decoherence take places at each step of EMF simulations,38 the decay of mixing method in which a mean-field wavefunction can be driven towards a single state accompanied by the surface hopping that occurs over a finite time rather than instantaneously,45-47 and the ab initio multiple cloning method in which a linear combination of configurations in Eq (3) is cloned into two and propagate on a single PES and the mean field, respectively.55 Finally, Gross and coworkers developed a new formalism for exact factorization as well as a time-dependent PES.56-59 Under the picture of exact factorization, a coupled-trajectory method was proposed to introduce quantum forces to the EOMs of electrons and nuclei to address overcoherence.60-62 The quantum force is determined by a swarm of interacting trajectories propagating simultaneously, which means that this model can be reduced to EMF in the case of uncoupled trajectories without quantum force. Beyond the scope of traditional EMF and SH, the methods based on the mixed quantum-classical Liouville equation are more rigorous. Starting from the quantum Liouville equation, a partial Wigner transformation with respect to classical degrees of freedom is performed on the density matrix, and then a set of adiabatic electronic basis is applied, resulting in the evolution equation for the mixed quantum-classical system under the adiabatic representation. Ando, Kapral and Ciccotti reported several variations on quantum-classical Liouville dynamics, all of which are able to capture 11



electronic coherence correctly.63-68 The connection to FSSH has been also discussed, which suggests that the quantum-classical Liouville dynamics can be reduced to FSSH with a large nuclear velocity along the direction of the derivative nonadiabatic coupling and a unique trajectory at any point in phase space at the price of the overcoherence problem in FSSH.37,

69

However, quantum-classical Liouville

dynamics suffers from a huge amount of necessary trajectories and the numerical instability during simulations, even to study simple systems with few degrees of freedom. The multiple spawning method developed by Martínez and coworkers is another attractive computational tool for ab initio MQC MD simulations.70-76 Starting from a multi-configuration expansion of the molecular wavefunction, the time-dependent nuclear wavefunction for each electronic state is represented as a linear combination of frozen Gaussian wavepackets. The EOMs of the coefficients as well as the parameters of Gaussians such as position, momentum and nuclear phase are derived from TDSE. A spawning procedure is employed to modulate the number of Gaussians according to the change on nonadiabatic couplings during simulations, which makes the method practicable to study photochemical systems at an affordable computational cost. The overcoherence problem disappears because many Gaussian wavepackets evolve simultaneously. On the other hand, this feature requires complicated numerical implementation. Some recent algorithms mentioned above, such as the overlap-based decoherence correction by Granucci et al.,50 the simultaneous FFSH by Shenvi et al.52 and the ab initio multiple cloning method by Makhov et al.,55 are also related to the spirit of multiple spawning. Unlike all of the MQC approaches discussed above, a semiclassical description can be applied to both the electronic and nuclear subsystems, for example, in the 12


Page 12 of 47

Page 13 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


semiclassical Meyer-Miller mapping Hamiltonian method77-79. The electronic degrees of freedom are projected onto a series of classical harmonic oscillators. The electronic states are represented using classical coordinates and momenta, and the traditional EMF scheme is employed with few changes. Based on the semiclassical quantization of the initial and final electronic action variables, the deficiency of the EMF method can be overcome using the symmetrical quasi-classical procedure80-82, in which the delta function for quantization is replaced by window functions and the information on electronic state is extracted from classical trajectories. The smoothing effect in the windowing procedure brings the over-coherent EMF dynamics more consistent with quantum results. More benchmarks on realistic photochemical reactions combined with on-the-fly electronic structure calculations are still necessary. For example, how to select an appropriate windowing function with an optimal parameter for different systems and how to deal with the long-time dynamic trajectories outside any window are unclear till now. IV. QUANTUM TRAJECTORY MEAN-FIELD METHOD In recent years we have developed the quantum trajectory mean-field method to address the overcoherence problem in a different way.83 Compared with the EMF and SH where the time evolution of quantum degrees of freedom is derived based on the wavefunction, in this approach we start from the EOM of the density matrix directly because the density matrix representation provides a better view for time-dependent problems. The wavefunction-based MQC MD theory has an obvious flaw from the absence of the classical counterpart of wavefunction. Unlike many existing ad hoc MQC algorithms in which an empirical damping term is introduced to represent the decoherence, we derived the EOM of the density matrix under the picture of quantum measurements. The classical nuclei are viewed as a classical observer that obeys the 13



Page 14 of 47

Newton’s second law, and the nuclear motion is considered as continuous quantum measurements to the electronic motion. The decoherence effect is reproduced with a measurement on the quantum subsystem by the observation of the classical subsystem. The strength of quantum measurement is related to the coupling between electronic and nuclear motion. The feedback of nuclear motion on the time evolution of the degrees of freedom for electrons is thus included intrinsically. The electronic degrees of freedom are expected to collapse to a single adiabatic state in the region without nonadiabatic effect. Therefore, the nonadiabatic effect near the conical intersections and the decoherence effect in the weak coupling region are captured simultaneously. In the QTMF method, we rewrite Eq (16) in a compact form and add two terms to the right hand side. We obtain the EOM of the density matrix as 𝑖 𝜌(𝑡) = ― [𝐻,𝜌(𝑡)] + ℏ

∑Γ

𝑗𝑘𝒟[𝑀𝑗𝑘]𝜌(𝑡)

𝑗≠𝑘

+ ∑𝑗 ≠ 𝑘 𝛾𝐹𝑗𝑘 + 𝛾′𝑗𝑘ℋ[𝑀𝑗𝑘]𝜌(𝑡)𝜉𝑗𝑘(𝑡)

(30)

The first term on the right hand side is equivalent to that of Eq (16) as follows 𝐻𝑗𝑘 = 𝐸𝑗𝛿𝑗𝑘 ―𝑖ℏ𝐑 ∙ 𝒅𝑗𝑘

(31)

The second one is a damping term, which describes the decoherence effect using the Lindblad super-operator 𝒟 as 1

𝒟[𝑀𝑗𝑘]𝜌(𝑡) = 𝑀𝑗𝑘𝜌(𝑡)𝑀𝑗𝑘† ― 2{𝑀𝑗𝑘† 𝑀𝑗𝑘,𝜌(𝑡)}

(32)

𝑀𝑗𝑘 = |𝜑𝑗⟩⟨𝜑𝑗| ― |𝜑𝑘⟩⟨𝜑𝑘|

(33)

where

The third one is a random term that accounts for the back-action effect of the nuclear motion to the electrons using the super-operator ℋ as ℋ[𝑀𝑗𝑘]𝜌(𝑡) = 𝑀𝑗𝑘𝜌(𝑡) +𝜌(𝑡)𝑀𝑗𝑘† ― tr[(𝑀𝑗𝑘 + 𝑀𝑗𝑘† )𝜌(𝑡)]𝜌(𝑡) 14


(34)

Page 15 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and 𝜉𝑗𝑘 denotes white noises. The heuristic consideration towards Eq (30) can be seen from our first paper on QTMF and the reference therein.83-86 A brief description is presented here. Consider two coupled quantum dots occupied by a single electron and measured by a quantum-point-contact. First, two qubit states |1⟩ and |2⟩ as well as a dephasing operator |1⟩⟨1| ― |2⟩⟨2| are defined, and the measurement back-action on the qubit state can be represented using the Lindblad-type master equation. The measurement currents correspond to the qubit electron either in the first or the second dot. Second, an analogy between the quantum measurement and the MQC

dynamics

simulation

is

proposed

with

continuous

quantum

weak

measurements. It makes the output currents stochastic and associated with the information of qubit state. As a result, the Lindblad-type master equation changes to the formula as Eq (30). Finally, the electronic subsystem and the nuclear motion in nonadiabatic dynamics simulations are regarded as the qubit and the measurement, respectively. Similar to the classical currents in the quantum-point-contact measurement, the nuclear forces on different PESs play an essential role to collapse the electronic subsystem to an adiabatic state, leading to the EOM of the density matrix. The damping and random terms are related to the first and the second points on quantum measurements, respectively. The additional terms in Eq (30) extract the information from nuclear motion to identify the electronic degrees of freedom in state j or k and result in a decoherence between j and k. As shown in Figure 1, two wavepackets traveling on two PESs with different forces are possible to diverge, so the decoherence takes place. It suggests the classical force as a key factor on information gain rates. Assume that the classical force is accompanied by stochastic fluctuations and the average satisfies a Gaussian distribution, the force-induced information gain rate denoted as 𝛾𝐹𝑗𝑘 is expressed as a 15



Page 16 of 47

function of classical forces in electronic state j and k in order to identify the two states based on the state distinguishable criterion in the quantum measurement theory, that is, 𝛾𝐹𝑗𝑘

|𝐅𝑗 ― 𝐅𝑘|2

= (|

𝜏 ―1 𝐅𝑗| + |𝐅𝑘|)2 𝑐

(35)

where 𝜏𝑐 is the characteristic time which scales atomic motion compared with the timescale of the electronic motion. Thus, we set 𝜏𝑐―1 = 0.001 in practice. The extra information gain rate in addition to 𝛾𝐹𝑗𝑘 is denoted as 𝛾′𝑗𝑘. For the model systems shown below, this term can be ignored and the QTMF becomes parameter-free. For more complex cases such as the application to high-dimensional molecules or the combination with different MQC MD algorithms, some empirical corrections to the force-induced information gain rate may be beneficial and we can set the parameter 𝛾′𝑗𝑘 as a constant when the nonadiabatic effect is as small as ℏ|(𝐑 ∙ 𝒅𝑗𝑘)/(𝐸𝑗 ― 𝐸𝑘)| < 𝑐

(36)

where c is set as 0.1 in practice. Finally, for the ideal measurements without any information loss we have the overall decoherence rate Γ𝑗𝑘 = 𝛾𝐹𝑗𝑘 + 𝛾′𝑗𝑘. Fortunately, we observed that different values on Γ𝑗𝑘 as well as 𝛾′𝑗𝑘 hardly affect the final results on typical model and realistic systems in our previous work,42,83 making the original QTMF method insensitive to decoherence parameters. The difference between the implementation on EMF and QTMF is the EOMs of the density matrix, in which Eqs (16) and (30) are used, respectively. Because of the random term in Eq (30), the total energy conservation is violated slightly in QTMF, which can be balanced using energy calibration schemes. The QTMF method works excellently on Tully’s three test systems, even for the reflection channels with the extended coupling Hamiltonian in which FSSH fails. 16


Page 17 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We took the photodissociation reaction of the diazirinone (N2CO⟶N2 +CO) as an example to perform MQC MD simulations with QTMF.42 At each MD step, the potential energies, gradients as well as the derivative nonadiabatic couplings were calculated at the highly accurate complete active space self-consistent field (CASSCF) level. The S1 state originates from the excitation of one electron from the oxygen n orbital to the N=N π* orbital. After the excitation to the S1 state at λ > 335 nm, two conical intersections can be achieved. At the lower-energy geometry denoted as CI-1, one of the C-N bonds is broken. However, no C-N bond breaking was observed at the higher-energy geometry denoted as CI-2. The first C-N bond fission is the dominant reaction channel in the S1 state, followed by the internal conversions via CI-1 and the second C-N bond dissociation in the ground state. The time constants for the stepwise C-N bond cleavages were estimated as 73.0 and 137.0 fs, respectively. The QTMF approach makes a good agreement with the ab initio multiple spawning simulations as a reference and surpasses the FSSH algorithm in accuracy. In addition, we replaced 𝛾𝐹 + 𝛾′ in the random term of Eq (30) by the overall decoherence rate Γ and set it as a constant. Interestingly, the influence on the final results is small. We suggested Γ = 0.001 as a reliable empirical parameter in practice. We

also

studied

the

photoisomerization

of

the

acetylacetone

(CH3COCH2COCH3) using QTMF combined with CASSCF calculations.87 The lowest three singlet states were involved in our simulations, that is, the ground state (S0), the nπ* state (S1) and the ππ* state (S2). The S2 state is initially populated upon photoexcitation at ∼265 nm. Three nonadiabatic pathways were observed. Three quarters of trajectories undergo the in-plane proton transfer between two oxygen atoms in the S2 state and relax to the S1 state through a conical intersection S2/S1-1. The rotational isomerization occurs once decaying to the S1 state, followed by an 17



internal conversion to the ground state through a conical intersection S1/S0. In contrast, 17% of trajectories experience the rotational isomerization in the S2 state before reaching another conical intersection S2/S1-2. The subsequent evolution of these trajectories is similar to that in the first pathway. The rest of trajectories move to the three-state intersection S2/S1/S0 via a nonplanar process, stop by in the S1 state, and decay to the S0 state. The time constant for the S2 proton transfer and the S1 lifetime can be well reproduced compared with the experimental data. On the other hand, the nonchelated structure in the S1 state was observed to be too favorable in our simulations, which results in an underestimated ratio of the chelated and non-chelated products. It is possible to further improve the accuracy by performing longer time simulations, if affordable. V. QUANTUM TRAJECTORY MEAN-FIELD METHOD COMBINED WITH SURFACE HOPPING The great popularity of FSSH mainly comes from its simplicity. During the FSSH simulations, nuclei always move on a single PES in an adiabatic electronic state with a probability for nonadiabatic hopping between electronic states. This picture is intuitive in chemistry, but its physical interpretation is controversial and still unclear. Here we propose a combination scheme based on QTMF, in which the decoherence correction is applied straightforwardly in the FSSH simulation. For an individual trajectory, the time evolution on the nuclei and electronic density matrix obey Eqs (24) and (30), respectively, with a transition probability estimated using Eq (28). The results of interest are obtained by averaging over a swarm of independent trajectories after simulations. This method is denoted as QTMF-FSSH. We outline the procedure for one trajectory as follows: (1) Initialize the nuclear coordinates, velocities and electronic density matrix. For 18


Page 18 of 47

Page 19 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a molecular system, the density matrix is set as an adiabatic electronic state based on the initial population after photoexcitation, and the coordinates and the conjugated velocities of nuclei are usually generated randomly from the Wigner distribution of harmonic vibrations in the ground state. (2) Calculate the potential energies and gradients on the current PES as well as the derivative nonadiabatic couplings between electronic states if available. For a molecular system, this is the most time-consuming step on account of the expensive computational cost on ab initio electronic structure calculations (e.g., CASSCF). (3) Propagate the nuclear coordinates and velocities according to Eq (24) with common integration methods such as the velocity Verlet algorithm. (4) Integrate Eq (30) using the fourth-order Runge-Kutta algorithm to evaluate the electronic density matrix. The information gain rates can be chosen based on Eq (35) or empirically as mentioned above. To study a molecular system, we noted that (a) The timescales of nuclear and electronic motion are different. Because of the expensive computational cost on nuclear degrees of freedom in step (2), the multiple time step algorithm is always used such as δ𝑡 = 0.001Δ𝑡, where Δ𝑡 and δ𝑡 are the integration steps in step (3) and (4), respectively. It also requires a linear interpolation procedure to integrate Eq (30) from 𝑡 to 𝑡 + Δ𝑡. (b) If the derivative nonadiabatic coupling is not available, the product 𝐑 ∙ 𝒅𝑗𝑘 in the EOM of the density matrix can be calculated numerically from the overlap of adiabatic electronic wavefunctions between 𝑡 and 𝑡 + Δ𝑡 as follows32

[𝐑 ∙ 𝒅𝑗𝑘](𝑡 + Δ𝑡/2) =

⟨𝜑𝑗(𝑡)│𝜑𝑘(𝑡 + Δ𝑡)⟩ ― ⟨𝜑𝑗(𝑡 + Δ𝑡)│𝜑𝑘(𝑡)⟩ 2Δ𝑡

(37)

Then the interpolation and extrapolation on 𝐑 ∙ 𝒅𝑗𝑘 are performed for the integration from 𝑡 to 𝑡 + Δ𝑡. (5) Calculate the transition probability between electronic states using Eq (28). 19



Applying the multiple time step, the integration should be performed simultaneously to step (4). Given a uniform random number 0 < 𝜉 < 1, the surface hopping to state k 𝑘―1

is attempted if ∑𝑖

𝑘

𝑃𝑗→𝑖 < 𝜉 < ∑𝑖 𝑃𝑗→𝑖, where j denotes the current state.

(6) If no switch occurs, return to step (2). Otherwise, adjust the nuclear velocities to maintain the total energy conservation. If the hop is energetically allowed, accept it and return to step (2). Based on the general semiclassical surface hopping propagator developed by Herman88-89 and the quantum-classical Liouville dynamics mentioned above64, 68, the rigorous way for velocity adjustment is to rescale the momentum along the direction of the derivative nonadiabatic coupling. If the derivative nonadiabatic coupling cannot be obtained on-the-fly, the velocity is often approximately rescaled along its current direction. If the hop is forbidden energetically, reject it and return to step (2). This procedure follows the FSSH algorithm31-32 except step (4), in which the EOM of density matrix changes from Eq (16) to (30). The numerical implementation is as simple as FSSH. Then we employ Tully’s three test systems reported in the original FSSH paper for validation. The exact results using quantum dynamics simulations are extracted from Tully’s paper.31 All parameters below are in atomic units. In this work, all trajectories start at 𝑥 = ― 20 on the lower surface with an incident momentum to the right, propagate with a time step as 0.1, and stop at 𝑥 =  25. The mass was set as 2000. Two different computing strategies are employed to investigate the overall decoherence and information gain rates in QTMF. In the strategy 1, we set the overall decoherence rate Γ𝑗𝑘 = 0.001 and 𝛾𝐹𝑗𝑘 + 𝛾′𝑗𝑘 = Γ𝑗𝑘 in the damping and random terms of Eq (30), while in the strategy 2, 𝛾𝐹𝑗𝑘 is calculated at each step using Eq (35) with 𝜏𝑐―1 = 0.001 and a constant 𝛾′𝑗𝑘 is added according to Eq (36) if the trajectory leaves the intersection regions. For each system, 2000 20


Page 20 of 47

Page 21 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


trajectories were sampled at each incident velocity using each method. The first system is the single-avoided crossing model. The diabatic Hamiltonian is given as 𝑉11(𝑥) =

{

𝐴[1 ― 𝑒 ―𝐵𝑥], 𝑥 > 0 ― 𝐴[1 ― 𝑒𝐵𝑥], 𝑥 < 0

(38a) (38b)

𝑉22(𝑥) = ― 𝑉11(𝑥) 𝑉12(𝑥) = 𝑉21(𝑥) = 𝐶𝑒 ―𝐷𝑥

2

(38c)

where A = 0.01, B = 1.6, C = 0.005 and D = 1.0. The adiabatic PESs on two states as well as the derivative nonadiabatic coupling are shown in Figure 2(a). The results are shown in Figure 3 by using different methods, including FSSH (a), QTMF-FSSH with the strategy 1 (b), and QTMF-FSSH with the strategy 2 with 𝛾′ = 0 (c) or a tunable value as 0.001 (d). All methods work well on this model. The essential resonance phenomenon at 7.7 < k < 8.9 as a reflection on the lower surface can be captured. On the other hand, the quantum tunneling at k < 4.5 cannot be reproduced, which is an intrinsic defect of most of the MQC MD approaches. The second system is the dual-avoided crossing model. The diabatic Hamiltonian is given as (39a)

𝑉11(𝑥) = 0 2

𝑉22(𝑥) = ―𝐴𝑒 ―𝐵𝑥 +𝐸 𝑉12(𝑥) = 𝑉21(𝑥) = 𝐶𝑒 ―𝐷𝑥

(39b) 2

(39c)

where A = 0.10, B = 0.28, C = 0.015, D = 0.06 and E = 0.05. The adiabatic PESs and the derivative nonadiabatic coupling are shown in Figure 2(b), with the results shown in Figure 4 by using FSSH (a), QTMF-FSSH with the strategy 1 (b), and QTMF-FSSH with the strategy 2 with 𝛾′ = 0 (c) or a tunable value as 0.0002 (d). In the low energy regime as ln(𝐸0) < ― 3.0, the FSSH overestimates the probability on 21



Page 22 of 47

the reflection on the lower surface, which can be improved using QTMF-FSSH. The result with the strategy 1 is slightly better than that with the strategy 2. The Stueckelberg oscillations in the high energy regime can be captured using the above methods. However, some errors on the amplitudes appear using the strategy 1, while the strategy 2 is more accurate. The third system is the extended coupling Hamiltonian model. The diabatic Hamiltonian is given as 𝑉11(𝑥) = 𝐴

(40a)

𝑉22(𝑥) = ―𝐴

(40b)

𝑉12(𝑥) = 𝑉21(𝑥) =

{

𝐵[2 ― 𝑒 ―𝐶𝑥], 𝑥 > 0 𝐵𝑒𝐶𝑥, 𝑥 < 0

(40c)

where A = 0.0006, B = 0.1 and C = 0.9. The adiabatic PESs and the derivative nonadiabatic coupling are shown in Figure 2(c), with the results shown in Figure 5 using FSSH (a), QTMF-FSSH with the strategy 1 (b), and QTMF-FSSH with the strategy 2 with 𝛾′ = 0 (c) or a tunable value as 0.10 (d). The poor performance of FSSH on two reflection channels is clearly due to the strong decoherence effect (also see Figure 1). This problem can be addressed excellently combined with QTMF. We examine the influence of the overall decoherence and information gain rates on the population probabilities. First, we repeat the above simulations in the strategy 1 with Γ = 0.0002 and 0.005, respectively. Surprisingly, QTMF-FSSH is sensitive to the value, which is quite different from our previous observations.42, 83 The results with Γ = 0.005 are incorrect seriously for all the three test systems. The simulations using Γ = 0.0002 are similar to FSSH. Therefore, the prediction on the reflection channels in the third system completely fails. Second, we modulate the values of 𝛾′ in the strategy 2. Although the choices on 𝛾′ in the original QTMF algorithm can be 22


Page 23 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


flexible, its influence on the performance of QTMF-FSSH is up and down. For example, the results on the third system can be improved slightly with 𝛾′ = 0.10. In contrast, the value of 𝛾′ > 0.01 hampers the accuracy on the first and second systems. As shown in panel (d) in Figures 3, 4 and 5, the best choice on 𝛾′ is system-dependent in the QTMF-FSSH. Therefore, we suggest taking the strategy 2 with 𝛾′ = 0 and getting 𝛾𝐹 from Eq (35) for applications of QTMF-FSSH in the future. The underlying reasons on the different parameter-dependence between two algorithms also need to be explored. Compared with the original algorithm of FSSH, the use of the EOM of the electronic density matrix in QTMF is able to overcome the overcoherence problem and give satisfactory results on Tully’s three test systems. One difference of SH from EMF is that the derivative nonadiabatic coupling is not necessary in the presence of numerical approximations. The related terms on Eq (16) or (30) can be obtained numerically and much more efficiently.32,

90-92

It makes the QTMF-FSSH attractive

since our final goal is to study complex photochemical and photobiological systems. Applying the QTMF on a realistic molecular system, the analytical calculation on the derivative nonadiabatic coupling required at each time step is very time-consuming. The efficiency can increase significantly in QTMF-FSSH. It also encourages us to employ some advanced electronic structure methods developed very recently93-97 to QTMF MD simulations as long as the potential energies and gradients on different electronic states are available. VI. CONCLUSION AND OUTLOOK In the present article, we have provided an overview of recent progress on the mixed quantum-classical molecular dynamics methods. We focused mainly on the overcoherence problem for Ehrenfest mean field, surface hopping, and quantum 23



trajectory mean-field approaches. The EMF and SH methods have been two standard computational tools to treat nonadiabatic photochemical reactions for several decades. Because of their popularity, people have also recognized some deficiencies on EMF and SH, in particular, the overcoherence problem that originates from the classical treatment on nuclei. It has been fixed up for past years in different ways, for example, the quantum-classical Liouville dynamics, the ab initio multiple spawning simulation and many ad hoc EMF or SH methods in which the decoherence rate is estimated based on the different PESs between electronic states or the overlap between auxiliary nuclear wavepackets. We have developed the QTMF method in recent years. The nuclear degrees of freedom are considered as a classical observer embedding in the mean field generated by the motion of electronic degrees of freedom. The EOM of the electronic density matrix in Eq (16) is replaced by Eq (30) under a different picture of quantum measurements. We have implemented it based on the CASSCF electronic structure calculation. Two applications, the photodissociation of the diazirinone and the photoisomerization of the acetylacetone, demonstrate its success. In this article, we have also developed a new scheme as QTMF-FSSH, in which the EOM of the density matrix in QTMF is introduced to the FSSH algorithm. The overcoherence in Tully’s three test systems can be addressed as excellently as the original QTMF method. The calculation on the derivative nonadiabatic coupling can be replaced by a simple numerical treatment in FSSH, which is an attractive feature compared with the original QTMF. On the other hand, the results on Tully’s three test systems depend on the choice on the information gain rates in the QTMF-FSSH, which is inferior to the original QTMF. We have proposed a simple strategy without further consideration on additional parameters, but special care should be taken case-by-case for applications on realistic molecular systems. 24


Page 24 of 47

Page 25 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We end this article with further improvements and extensions on QTMF. The first one is to generalize the present method to deal with intersystem crossing (ISC) dynamics between electronic states of different multiplicities. The ISC processes are induced by spin-orbit coupling (SOC). On one hand, the diabatic representation, in which the diagonal elements are the common adiabatic electronic wavefunctions and the off-diagonal terms are determined by SOC, can be employed straightforwardly in the existing MQC MD simulations.98 On the other hand, it has been reported that the spin-adiabatic representation is more appropriate to simulate the ISC processes using surface hopping.99-102 However, the adiabatic-to-diabatic transformation should be treated in particular. Both representations are good candidates to describe the spin forbidden transitions under the picture of QTMF, but some nontrivial modifications on the present formulation are necessary. Another extension is the QTMF MD simulation with quantum mechanical and molecular mechanical (QM/MM) approaches for studying the photochemical and photobiological reactions in a complex environment.103-105 In the QM/MM model, a small number of atoms directly related to the nonadiabatic process are selected for QM calculations such as CASSCF, while the contribution of the rest of the system is described with an empirical MM force field. The total potential energy of the whole system in electronic state j in the QM/MM model is written as vdW cov 𝐸 = 𝐸QM + 𝐸ele QM/MM + 𝐸QM/MM + 𝐸QM/MM + 𝐸MM

(41)

The first and the second terms in Eq (41) are treated at the QM level as

⟨ │

𝑞𝑏

│ ⟩

𝑍𝛼𝑞𝑏

𝐸QM + 𝐸ele QM/MM = 𝜑𝑗 𝐻QM + ∑𝑏 ∈ MM|𝐫 ― 𝐑𝑏| 𝜑𝑗 + ∑𝛼 ∈ QM∑𝑏 ∈ MM|𝐑𝛼 ― 𝐑𝑏|

(42)

where α, b and j denote the QM atom, the MM atom and the current electronic state, respectively, Zα is the nuclear charge of QM atom α, qb is the point charge of MM atom b, R denotes the coordinates of the QM (nuclei) and MM atoms, and r denotes 25



the coordinates of the QM electrons. The rest terms in Eq (41) are assumed to be independent of the electronic states and can be obtained using MM force fields. To solve the EOM of the nuclear and electronic degrees of freedom under the QM/MM model, Ej in Eqs (18) and (31) should be calculated using Eq (41), and djk should be calculated using Eq (17) in which the electronic adiabatic wavefunctions 𝜑𝑗 and 𝜑𝑘 are obtained with Eq (42). We also assume that the MM components of djk can be neglected.98 Recently we implemented the QM/MM QTMF MD method to investigate the photo-induced ring-opening process of 2(5H)-thiophenone in the CH3CN solution. The preliminary results are promising. In particular, an ultrafast process along a diabatic pathway was observed, which is quite different from those reported for many other heterocyclic molecules.106-107 The QM/MM QTMF MD method presented here is still far from perfect. First, the polarization of the MM subsystem induced by the QM electrostatic interactions is absent. The polarizable MM force field would be desirable to achieve higher accuracy but at a price of more expensive computational cost.108-109 Second, the influence of the MM subsystem on the damping and random terms on the right hand side of Eq (30) is ignored if the decoherence rate Γ𝑗𝑘 is constant. However, the decoherence induced by nuclear motions should be much more significant in the condensed phase than that in gas phase on account of thousands or millions of classical MM degrees of freedom.38 Further investigations on this issue are needed. Finally, both the original QTMF and QTMF-FSSH methods violate the total energy conservation. In the QM/MM model, the energy discontinuity on the QM subsystem can be balanced by the kinetic energy of the MM environment without artificial energy calibration, but how to perform it during the time evolution in a rigorous way is unclear. Up to now we have not discussed on the multistate electronic structure problem. 26


Page 26 of 47

Page 27 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The MQC MD simulations are performed on-the-fly using the potential energies and gradients in the ground and excited states at the current MD step for propagation. Most of them require the derivative nonadiabatic couplings between electronic states. The calculation on these quantities is quite difficult. Taking account of the accuracy, the multi-configurational self-consistent field (MCSCF) method such as CASSCF is an excellent choice.110 More efficient single-reference method such as time-dependent density functional theory (TDDFT) is frequently employed to study larger systems.111 The electronic structure calculation is the computational bottleneck for all MQC MD simulations on molecular systems. In order to overcome this limitation and achieve high accuracy and efficiency simultaneously, machine learning techniques such as neural network are being more and more attractive in the fields of quantum chemistry and molecular simulation.112-118 It has been shown very recently that kernel ridge regression and deep learning can be used for nonadiabatic dynamics on polyatomic molecules.119-120 On the other hand, some challenging issues are still being under research in our group. For example, the construction on a big database at the CASSCF or higher level is too expensive computationally for many-atom systems. Suffering from insufficient samples at the beginning, one must consider how to estimate the credibility of predictions, how to detect new configurations outside the existing database, and how to adapt the present machine learning model on-the-fly during nonadiabatic dynamics.121-122 Another challenge is an appropriate descriptor on the derivative nonadiabatic coupling because it tends to singularity in the vicinity of conical intersection.123-124 Except the quantitative predictions, another class of machine learning algorithms can be modified and applied to classify the resulting trajectories as well as a large amount of configurations.125-126 More useful information on the mechanism of nonadiabatic photochemical reactions is expected to be extracted 27



effectively. In conclusion, the quantum trajectory mean-field method is a powerful tool for nonadiabatic molecular dynamics simulations in photochemistry. The overcoherence problem on the traditional Ehrenfest mean field and surface hopping methods can be addressed. The combination with the fewest-switches surface hopping provides an alternative and opens up the possibility to incorporate improvements on electronic structure theory and surface hopping algorithm into this method. None of the mixed quantum-classical methods can work well in all cases, which leaves a lot of room in our future work. AUTHOR INFORMATION Corresponding Author *E-mail: [email protected] Notes The authors declare no competing financial interest.

ACKNOWLEDGEMENT We are grateful to financial support from the NSFC (Grant Nos.: 21590801, 21688102, 21520102005, and 21421003) for W.F. and the Fundamental Research Funds for the Central Universities (2018NTST10) for L.S.

REFERENCE 1.

Lin, S. H.; Chang, C. H.; Liang, K. K.; Chang, R.; Shiu, Y. J.; Zhang, J. M.;

Yang, T. S.; Hayashi, M.; Hsu, F. C., Ultrafast dynamics and spectroscopy of bacterial photosynthetic reaction centers. Adv. Chem. Phys. 2002, 121, 1-88. 2.

Niu, Y.; Peng, Q.; Deng, C.; Gao, X.; Shuai, Z., Theory of excited state decays 28


Page 28 of 47

Page 29 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and optical spectra: application to polyatomic molecules. J. Phys. Chem. A 2010, 114, 7817-7831. 3.

Jaeger, H. M.; Hyeon-Deuk, K.; Prezhdo, O. V., Exciton multiplication from first

principles. Acc. Chem. Res. 2013, 46, 1280-1289. 4.

Van Voorhis, T.; Kowalczyk, T.; Kaduk, B.; Wang, L. P.; Cheng, C. L.; Wu, Q.,

The diabatic picture of electron transfer, reaction barriers, and molecular dynamics. Annu. Rev. Phys. Chem. 2010, 61, 149-170. 5.

Blumberger, J., Recent Advances in the Theory and Molecular Simulation of

Biological Electron Transfer Reactions. Chem. Rev. 2015, 115, 11191-11238. 6.

Wu, L.; Cao, X.; Chen, X.; Fang, W.; Dolg, M., Visible-Light Photocatalysis of

C(sp3)-H Fluorination by the Uranyl Ion: Mechanistic Insights. Angew. Chem. Int. Ed. Engl. 2018, 57, 11812-11816. 7.

Rafiq, S.; Scholes, G. D., From Fundamental Theories to Quantum Coherences in

Electron Transfer. J. Am. Chem. Soc. 2019, 141, 708-722. 8.

Head‐Gordon, M.; Tully, J. C., Vibrational relaxation on metal surfaces:

Molecular‐orbital theory and application to CO/Cu(100). J. Chem. Phys. 1992, 96, 3939-3949. 9.

Head‐Gordon, M.; Tully, J. C., Molecular dynamics with electronic frictions. J.

Chem. Phys. 1995, 103, 10137-10145. 10. Shenvi, N.; Roy, S.; Tully, J. C., Dynamical steering and electronic excitation in NO scattering from a gold surface. Science 2009, 326, 829-832. 11. Maurer, R. J.; Jiang, B.; Guo, H.; Tully, J. C., Mode Specific Electronic Friction in Dissociative Chemisorption on Metal Surfaces: H2 on Ag(111). Phys. Rev. Lett. 2017, 118, 256001. 12. Dou, W.; Subotnik, J. E., Perspective: How to understand electronic friction. J. Chem. Phys. 2018, 148, 230901. 13. Yarkony, D. R., Diabolical conical intersections. Rev. Mod. Phys. 1996, 68, 985-1013. 14. Yarkony, D. R., Conical Intersections: The New Conventional Wisdom. J. Phys. Chem. A 2001, 105, 6277-6293. 29



15. Worth, G. A.; Cederbaum, L. S., Beyond Born-Oppenheimer: molecular dynamics through a conical intersection. Annu. Rev. Phys. Chem. 2004, 55, 127-158. 16. Fang, W. H., Ab initio determination of dark structures in radiationless transitions for aromatic carbonyl compounds. Acc. Chem. Res. 2008, 41, 452-457. 17. Yarkony, D. R., Nonadiabatic quantum chemistry--past, present, and future. Chem. Rev. 2012, 112, 481-498. 18. Schuurman, M. S.; Stolow, A., Dynamics at Conical Intersections. Annu. Rev. Phys. Chem. 2018, 69, 427-450. 19. McCullough, E. A.; Wyatt, R. E., Dynamics of the Collinear H+H2 Reaction. I. Probability Density and Flux. J. Chem. Phys. 1971, 54, 3578-3591. 20. Beck, M. H.; Jackle, A.; Worth, G. A.; Meyer, H. D., The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Phys. Rep. 2000, 324, 1-105. 21. Wang, H. B.; Thoss, M., Multilayer formulation of the multiconfiguration time-dependent Hartree theory. J. Chem. Phys. 2003, 119, 1289-1299. 22. Xiao, C.; Xu, X.; Liu, S.; Wang, T.; Dong, W.; Yang, T.; Sun, Z.; Dai, D.; Xu, X.; Zhang, D. H.; Yang, X., Experimental and theoretical differential cross sections for a four-atom reaction: HD + OH --> H2O + D. Science 2011, 333, 440-442. 23. Zhang, D. H.; Guo, H., Recent Advances in Quantum Dynamics of Bimolecular Reactions. Annu. Rev. Phys. Chem. 2016, 67, 135-158. 24. Richings, G. W.; Habershon, S., Direct Quantum Dynamics Using Grid-Based Wave Function Propagation and Machine-Learned Potential Energy Surfaces. J. Chem. Theory Comput. 2017, 13, 4012-4024. 25. Stock, G.; Thoss, M., Classical Description of Nonadiabatic Quantum Dynamics. Adv. Chem. Phys. 2005, 131, 243-375. 26. Yonehara, T.; Hanasaki, K.; Takatsuka, K., Fundamental approaches to nonadiabaticity: toward a chemical theory beyond the Born-Oppenheimer paradigm. Chem. Rev. 2012, 112, 499-542. 27. Gao, L. H.; Xie, B. B.; Fang, W. H., Theories and Applications of Mixed Quantum-Classical Non-adiabatic Dynamics. Chin. J. Chem. Phys. 2018, 31, 12-26. 30


Page 30 of 47

Page 31 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


28. Crespo-Otero, R.; Barbatti, M., Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics. Chem. Rev. 2018, 118, 7026-7068. 29. Tully, J. C., Mixed quantum–classical dynamics. Faraday Discuss. 1998, 110, 407-419. 30. Curchod, B. F.; Rothlisberger, U.; Tavernelli, I., Trajectory-based nonadiabatic dynamics with time-dependent density functional theory. Chemphyschem 2013, 14, 1314-1340. 31. Tully, J. C., Molecular-Dynamics with Electronic-Transitions. J. Chem. Phys. 1990, 93, 1061-1071. 32. Hammes-Schiffer,

S.;

Tully,

J.

C.,

Proton-Transfer

in

Solution

-

Molecular-Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657-4667. 33. Fu, B.; Shepler, B. C.; Bowman, J. M., Three-state trajectory surface hopping studies of the photodissociation dynamics of formaldehyde on ab initio potential energy surfaces. J. Am. Chem. Soc. 2011, 133, 7957-7968. 34. Yue, L.; Lan, Z.; Liu, Y. J., The Theoretical Estimation of the Bioluminescent Efficiency of the Firefly via a Nonadiabatic Molecular Dynamics Simulation. J. Phys. Chem. Lett. 2015, 6, 540-548. 35. Ding, B. W.; Liu, Y. J., Bioluminescence of Firefly Squid via Mechanism of Single Electron-Transfer Oxygenation and Charge-Transfer-Induced Luminescence. J. Am. Chem. Soc. 2017, 139, 1106-1119. 36. Li, Y.; Gong, Q.; Yue, L.; Wang, W.; Liu, F., Photochemistry of the Simplest Criegee Intermediate, CH2OO: Photoisomerization Channel toward Dioxirane Revealed by CASPT2 Calculations and Trajectory Surface-Hopping Dynamics. J. Phys. Chem. Lett. 2018, 9, 978-981. 37. Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N., Understanding the Surface Hopping View of Electronic Transitions and Decoherence. Annu. Rev. Phys. Chem. 2016, 67, 387-417. 38. Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J., Mean-field dynamics with 31



Page 32 of 47

stochastic decoherence (MF-SD): a new algorithm for nonadiabatic mixed quantum/classical

molecular-dynamics

simulations

with

nuclear-induced

decoherence. J. Chem. Phys. 2005, 123, 234106. 39. Schwartz, B. J.; Bittner, E. R.; Prezhdo, O. V.; Rossky, P. J., Quantum decoherence and the isotope effect in condensed phase nonadiabatic molecular dynamics simulations. J. Chem. Phys. 1996, 104, 5942-5955. 40. Larsen, R. E.; Bedard-Hearn, M. J.; Schwartz, B. J., Exploring the role of decoherence in condensed-phase nonadiabatic dynamics: a comparison of different mixed quantum/classical simulation algorithms for the excited hydrated electron. J. Phys. Chem. B 2006, 110, 20055-20066. 41. Jiang, R.; Sibert, E. L., III, Surface hopping simulation of vibrational predissociation of methanol dimer. J. Chem. Phys. 2012, 136, 224104. 42. Xie, B.; Liu, L.; Cui, G.; Fang, W. H.; Cao, J.; Feng, W.; Li, X. Q., Ab initio implementation of quantum trajectory mean-field approach and dynamical simulation of the N2CO photodissociation. J. Chem. Phys. 2015, 143, 194107. 43. Fang, J. Y.; Hammes-Schiffer, S., Improvement of the internal consistency in trajectory surface hopping. J. Phys. Chem. A 1999, 103, 9399-9407. 44. Fang, J. Y.; Hammes-Schiffer, S., Comparison of surface hopping and mean field approaches for model proton transfer reactions. J. Chem. Phys. 1999, 110, 11166-11175. 45. Zhu, C.; Jasper, A. W.; Truhlar, D. G., Non-Born-Oppenheimer trajectories with self-consistent decay of mixing. J. Chem. Phys. 2004, 120, 5543-5557. 46. Zhu, C.; Nangia, S.; Jasper, A. W.; Truhlar, D. G., Coherent switching with decay of mixing: an improved treatment of electronic coherence for non-Born-Oppenheimer trajectories. J. Chem. Phys. 2004, 121, 7658-7670. 47. Zhu, C.; Jasper, A. W.; Truhlar, D. G., Non-Born-Oppenheimer Liouville-von Neumann Dynamics. Evolution of a Subsystem Controlled by Linear and Population-Driven Decay of Mixing with Decoherent and Coherent Switching. J. Chem. Theory Comput. 2005, 1, 527-540. 48. Zhu, C., Restoring electronic coherence/decoherence for a trajectory-based 32


Page 33 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonadiabatic molecular dynamics. Sci. Rep. 2016, 6, 24198. 49. Prezhdo, O. V.; Rossky, P. J., Evaluation of quantum transition rates from quantum-classical molecular dynamics simulations. J. Chem. Phys. 1997, 107, 5863-5878. 50. Granucci, G.; Persico, M.; Zoccante, A., Including quantum decoherence in surface hopping. J. Chem. Phys. 2010, 133, 134111. 51. Subotnik, J. E.; Shenvi, N., A new approach to decoherence and momentum rescaling in the surface hopping algorithm. J. Chem. Phys. 2011, 134, 024105. 52. Shenvi, N.; Subotnik, J. E.; Yang, W., Simultaneous-trajectory surface hopping: a parameter-free algorithm for implementing decoherence in nonadiabatic dynamics. J. Chem. Phys. 2011, 134, 144102. 53. Akimov, A. V.; Long, R.; Prezhdo, O. V., Coherence penalty functional: a simple method for adding decoherence in Ehrenfest dynamics. J. Chem. Phys. 2014, 140, 194107. 54. Prezhdo, O. V.; Rossky, P. J., Mean-field molecular dynamics with surface hopping. J. Chem. Phys. 1997, 107, 825-834. 55. Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V., Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics. J. Chem. Phys. 2014, 141, 054110. 56. Abedi, A.; Maitra, N. T.; Gross, E. K., Exact factorization of the time-dependent electron-nuclear wave function. Phys. Rev. Lett. 2010, 105, 123002. 57. Abedi, A.; Maitra, N. T.; Gross, E. K., Correlated electron-nuclear dynamics: exact factorization of the molecular wavefunction. J. Chem. Phys. 2012, 137, 22A530. 58. Abedi, A.; Agostini, F.; Suzuki, Y.; Gross, E. K., Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface. Phys. Rev. Lett. 2013, 110, 263001. 59. Suzuki, Y.; Abedi, A.; Maitra, N. T.; Yamashita, K.; Gross, E. K. U., Electronic Schrodinger equation with nonclassical nuclei. Phys. Rev. A 2014, 89, 040501. 60. Min, S. K.; Agostini, F.; Gross, E. K., Coupled-Trajectory Quantum-Classical Approach to Electronic Decoherence in Nonadiabatic Processes. Phys. Rev. Lett. 33



Page 34 of 47

2015, 115, 073001. 61. Agostini, F.; Min, S. K.; Abedi, A.; Gross, E. K., Quantum-Classical Nonadiabatic Dynamics: Coupled- vs Independent-Trajectory Methods. J. Chem. Theory Comput. 2016, 12, 2127-2143. 62. Min, S. K.; Agostini, F.; Tavernelli, I.; Gross, E. K. U., Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence. J. Phys. Chem. Lett. 2017, 8, 3048-3055. 63. Kapral, R.; Ciccotti, G., Mixed quantum-classical dynamics. J. Chem. Phys. 1999, 110 , 8919-8929. 64. Nielsen, S.; Kapral, R.; Ciccotti, G., Mixed quantum-classical surface hopping dynamics. J. Chem. Phys. 2000, 112, 6543-6553. 65. Ando, K., Non-adiabatic couplings in Liouville description of mixed quantum-classical dynamics. Chem. Phys. Lett. 2002, 360, 240-242. 66. Ando, K.; Santer, M., Mixed quantum-classical Liouville molecular dynamics without momentum jump. J. Chem. Phys. 2003, 118, 10399-10406. 67. Kapral, R., Progress in the theory of mixed quantum-classical dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129-157. 68. Kapral, R., Surface hopping from the perspective of quantum–classical Liouville dynamics. Chem. Phys. 2016, 481, 77-83. 69. Subotnik, J. E.; Ouyang, W.; Landry, B. R., Can we derive Tully's surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence. J. Chem. Phys. 2013, 139, 214107. 70. Martinez, T. J.; BenNun, M.; Levine, R. D., Multi-electronic-state molecular dynamics: A wave function approach with applications. J. Phys. Chem. 1996, 100, 7884-7895. 71. Ben-Nun, M.; Martinez, T. J., Nonadiabatic molecular dynamics: Validation of the multiple spawning method for a multidimensional problem. J. Chem. Phys. 1998, 108, 7244-7257. 72. Ben-Nun,

M.;

Quenneville,

J.;

Martínez,

T.

J.,

Ab

Initio

Multiple

Spawning: Photochemistry from First Principles Quantum Molecular Dynamics. J. 34


Page 35 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Phys. Chem. A 2000, 104, 5161-5175. 73. Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J., Implementation of ab initio multiple spawning in the MOLPRO quantum chemistry package. Chem. Phys. 2008, 347, 3-16. 74. Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Martinez, T. J., Photodynamics in complex environments: ab initio multiple spawning quantum mechanical/molecular mechanical dynamics. J. Phys. Chem. B 2009, 113, 3280-3291. 75. Yang, S.; Coe, J. D.; Kaduk, B.; Martinez, T. J., An "optimal" spawning algorithm for adaptive basis set expansion in nonadiabatic dynamics. J. Chem. Phys. 2009, 130, 134113. 76. Curchod, B. F. E.; Martinez, T. J., Ab Initio Nonadiabatic Quantum Molecular Dynamics. Chem. Rev. 2018, 118, 3305-3336. 77. Meyera, H. D.; Miller, W. H., A classical analog for electronic degrees of freedom in nonadiabatic collision processes. J. Chem. Phys. 1979, 70, 3214-3223. 78. Meyer, H. D.; Miller, W. H., Analysis and Extension of Some Recently Proposed Classical-Models for Electronic Degrees of Freedom. J. Chem. Phys. 1980, 72, 2272-2281. 79. Cotton, S. J.; Liang, R.; Miller, W. H., On the adiabatic representation of Meyer-Miller electronic-nuclear dynamics. J. Chem. Phys. 2017, 147, 064112. 80. Miller, W. H.; Cotton, S. J., Classical molecular dynamics simulation of electronically non-adiabatic processes. Faraday Discuss. 2016, 195, 9-30. 81. Cotton, S. J.; Miller, W. H., A new symmetrical quasi-classical model for electronically non-adiabatic processes: Application to the case of weak non-adiabatic coupling. J. Chem. Phys. 2016, 145, 144108. 82. Cotton, S. J.; Miller, W. H., A symmetrical quasi-classical windowing model for the molecular dynamics treatment of non-adiabatic processes involving many electronic states. J. Chem. Phys. 2019, 150, 104101. 83. Feng, W.; Xu, L. T.; Li, X. Q.; Fang, W. H.; Yan, Y. J., Nonadiabatic molecular dynamics simulation: An approach based on quantum measurement picture. AIP Adv. 35



2014, 4, 077131. 84. Goan, H. S.; Milburn, G. J., Dynamics of a mesoscopic charge quantum bit under continuous quantum measurement. Phys. Rev. B 2001, 64, 235307. 85. Jin, J. S.; Li, X. Q.; Yan, Y. J., Quantum coherence control of solid-state charge qubit by means of a suboptimal feedback algorithm. Phys. Rev. B 2006, 73, 233302. 86. Wiseman, H. M.; Milburn, G. J., Quantum Measurement and Control. Cambridge University Press: 2010. 87. Xie, B.; Cui, G.; Fang, W. H., Multiple-State Nonadiabatic Dynamics Simulation of Photoisomerization of Acetylacetone with the Direct ab Initio QTMF Approach. J. Chem. Theory Comput. 2017, 13, 2717-2729. 88. Herman, M. F., A Semiclassical Surface Hopping Propagator for Nonadiabatic Problems. J. Chem. Phys. 1995, 103, 8081-8097. 89. Herman, M. F., Toward an accurate and efficient semiclassical surface hopping procedure for nonadiabatic problems. J. Phys. Chem. A 2005, 109, 9196-9205. 90. Fabiano, E.; Keal, T. W.; Thiel, W., Implementation of surface hopping molecular dynamics using semiempirical methods. Chem. Phys. 2008, 349, 334-347. 91. Pittner, J.; Lischka, H.; Barbatti, M., Optimization of mixed quantum-classical dynamics: Time-derivative coupling terms and selected couplings. Chem. Phys. 2009, 356, 147-152. 92. Barbatti, M.; Pittner, J.; Pederzoli, M.; Werner, U.; Mitric, R.; Bonacic-Koutecky, V.; Lischka, H., Non-adiabatic dynamics of pyrrole: Dependence of deactivation mechanisms on the excitation energy. Chem. Phys. 2010, 375, 26-34. 93. Hu, W.; Chan, G. K., Excited-State Geometry Optimization with the Density Matrix Renormalization Group, as Applied to Polyenes. J. Chem. Theory Comput. 2015, 11, 3000-3009. 94. Filatov, M., Spin-restricted ensemble-referenced Kohn-Sham method: basic principles and application to strongly correlated ground and excited states of molecules. Wires Comput Mol Sci 2015, 5, 146-167. 95. Zhang, D.; Peng, D.; Zhang, P.; Yang, W., Analytic gradients, geometry optimization and excited state potential energy surfaces from the particle-particle 36


Page 36 of 47

Page 37 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


random phase approximation. Phys. Chem. Chem. Phys. 2015, 17, 1025-1038. 96. Chattopadhyay, S.; Chaudhuri, R. K.; Mahapatra, U. S.; Ghosh, A.; Ray, S. S., State-specific multireference perturbation theory: development and present status. Wires Comput Mol Sci 2016, 6, 266-291. 97. Lischka, H.; Nachtigallova, D.; Aquino, A. J. A.; Szalay, P. G.; Plasser, F.; Machado, F. B. C.; Barbatti, M., Multireference Approaches for Excited States of Molecules. Chem. Rev. 2018, 118, 7293-7361. 98. Cui, G.; Thiel, W., Generalized trajectory surface-hopping method for internal conversion and intersystem crossing. J. Chem. Phys. 2014, 141, 124101. 99. Richter, M.; Marquetand, P.; Gonzalez-Vazquez, J.; Sola, I.; Gonzalez, L., SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. J. Chem. Theory Comput. 2011, 7, 1253-1258. 100.

Granucci, G.; Persico, M.; Spighi, G., Surface hopping trajectory simulations

with spin-orbit and dynamical couplings. J. Chem. Phys. 2012, 137, 22A501. 101.

Mai, S.; Marquetand, P.; Gonzalez, L., A General Method to Describe

Intersystem Crossing Dynamics in Trajectory Surface Hopping. Int. J. Quant. Chem. 2015, 115, 1215-1231. 102.

Pederzoli, M.; Pittner, J., A new approach to molecular dynamics with

non-adiabatic and spin-orbit effects with applications to QM/MM simulations of thiophene and selenophene. J. Chem. Phys. 2017, 146, 114101. 103.

Brunk, E.; Rothlisberger, U., Mixed Quantum Mechanical/Molecular

Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chem. Rev. 2015, 115, 6217-6263. 104.

Chung, L. W.; Sameera, W. M.; Ramozzi, R.; Page, A. J.; Hatanaka, M.;

Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H. B.; Ding, L.; Morokuma, K., The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678-5796. 105.

Gozem, S.; Luk, H. L.; Schapiro, I.; Olivucci, M., Theory and Simulation of

the Ultrafast Double-Bond Isomerization of Biological Chromophores. Chem. Rev. 2017, 117, 13502-13565. 37



106.

Page 38 of 47

Murdock, D.; Harris, S. J.; Luke, J.; Grubb, M. P.; Orr-Ewing, A. J.; Ashfold,

M. N., Transient UV pump-IR probe investigation of heterocyclic ring-opening dynamics in the solution phase: the role played by nσ* states in the photoinduced reactions of thiophenone and furanone. Phys. Chem. Chem. Phys. 2014, 16, 21271-21279. 107.

Ashfold, M. N. R.; Bain, M.; Hansen, C. S.; Ingle, R. A.; Karsili, T. N. V.;

Marchetti, B.; Murdock, D., Exploring the Dynamics of the Photoinduced Ring-Opening of Heterocyclic Molecules. J. Phys. Chem. Lett. 2017, 8, 3440-3451. 108.

Loco, D.; Cupellini, L., Modeling the absorption lineshape of embedded

systems from molecular dynamics: A tutorial review. Int. J. Quant. Chem. 2019, 119, e25726. 109.

Zuehlsdorff, T. J.; Isborn, C. M., Modeling absorption spectra of molecules

in solution. Int. J. Quant. Chem. 2019, 119, e25719. 110.

Schmidt, M. W.; Gordon, M. S., The construction and interpretation of

MCSCF wavefunctions. Annu. Rev. Phys. Chem. 1998, 49, 233-266. 111.

Casida,

M.

E.;

Huix-Rotllant,

M.,

Progress

in

time-dependent

density-functional theory. Annu. Rev. Phys. Chem. 2012, 63, 287-323. 112.

Morawietz, T.; Singraber, A.; Dellago, C.; Behler, J., How van der Waals

interactions determine the unique properties of water. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 8368-8373. 113.

Shen, L.; Wu, J.; Yang, W., Multiscale Quantum Mechanics/Molecular

Mechanics Simulations with Neural Networks. J. Chem. Theory Comput. 2016, 12, 4934-4946. 114.

Smith, J. S.; Isayev, O.; Roitberg, A. E., ANI-1: an extensible neural network

potential with DFT accuracy at force field computational cost. Chem. Sci. 2017, 8, 3192-3203. 115.

Schneider, E.; Dai, L.; Topper, R. Q.; Drechsel-Grau, C.; Tuckerman, M. E.,

Stochastic Neural Network Approach for Learning High-Dimensional Free Energy Surfaces. Phys. Rev. Lett. 2017, 119, 150601. 116.

Schutt, K. T.; Sauceda, H. E.; Kindermans, P. J.; Tkatchenko, A.; Muller, K. 38


Page 39 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


R., SchNet - A deep learning architecture for molecules and materials. J. Chem. Phys. 2018, 148, 241722. 117.

Shen, L.; Yang, W., Molecular Dynamics Simulations with Quantum

Mechanics/Molecular Mechanics and Adaptive Neural Networks. J. Chem. Theory Comput. 2018, 14, 1442-1455. 118.

Zhang, L.; Han, J.; Wang, H.; Car, R.; E, W., Deep Potential Molecular

Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics. Phys. Rev. Lett. 2018, 120, 143001. 119.

Hu, D.; Xie, Y.; Li, X.; Li, L.; Lan, Z., Inclusion of Machine Learning

Kernel Ridge Regression Potential Energy Surfaces in On-the-Fly Nonadiabatic Molecular Dynamics Simulation. J. Phys. Chem. Lett. 2018, 9, 2725-2732. 120.

Chen, W. K.; Liu, X. Y.; Fang, W. H.; Dral, P. O.; Cui, G., Deep Learning

for Nonadiabatic Excited-State Dynamics. J. Phys. Chem. Lett. 2018, 9, 6702-6708. 121.

Behler, J., First Principles Neural Network Potentials for Reactive

Simulations of Large Molecular and Condensed Systems. Angew. Chem. Int. Ed. Engl. 2017, 56, 12828-12840. 122.

Zhang, P.; Shen, L.; Yang, W., Solvation Free Energy Calculations with

Quantum Mechanics/Molecular Mechanics and Machine Learning Models. J. Phys. Chem. B 2019, 123, 901-908. 123.

Baer, M., Introduction to the theory of electronic non-adiabatic coupling

terms in molecular systems. Phys. Rep. 2002, 358, 75-142. 124.

Dral, P. O.; Barbatti, M.; Thiel, W., Nonadiabatic Excited-State Dynamics

with Machine Learning. J. Phys. Chem. Lett. 2018, 9, 5660-5663. 125.

Li, X.; Xie, Y.; Hu, D.; Lan, Z., Analysis of the Geometrical Evolution in

On-the-Fly Surface-Hopping Nonadiabatic Dynamics with Machine Learning Dimensionality Reduction Approaches: Classical Multidimensional Scaling and Isometric Feature Mapping. J. Chem. Theory Comput. 2017, 13, 4611-4623. 126.

Li, X.; Hu, D.; Xie, Y.; Lan, Z., Analysis of trajectory similarity and

configuration similarity in on-the-fly surface-hopping simulation on multi-channel nonadiabatic photoisomerization dynamics. J. Chem. Phys. 2018, 149, 244104. 39


The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 40


Page 40 of 47

Page 41 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Decoherence on the extended coupling Hamiltonian model (Tully’s third test system). With an incident momentum to the right, the wavepacket propagates on the lower surface and spawns a wavepacket on the upper surface in the intersection region. After travelling to the bifurcation region, the lower wavepacket transmits (blue) while the upper one reflects (red). Finally, two components completely separate away and continue to propagate independently.

41



Figure 2. Adiabatic potential energy surfaces and derivative nonadiabatic couplings for Tully’s three test systems, involving the single-avoided crossing model (a), the dual-avoided crossing model (b) and the extended coupling Hamiltonian model (c). All quantities are in atomic units.

42


Page 42 of 47

Page 43 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Comparisons of simulations using different algorithms (red) with the exact quantum dynamics results (blue) on the first test system. (a) FSSH; (b) QTMF-FSSH using a constant Γ = 0.001; (c) QTMF-FSSH with 𝛾𝐹 as Eq (35) and 𝛾′ = 0; (d) QTMF-FSSH with 𝛾𝐹 as Eq (35) and 𝛾′ = 0.001 according to Eq (36). Different shapes represent the probability on different channels (circle: T1, transmission in the lower surface; triangle: T2, transmission in the upper surface; square: R1, reflection in the lower surface).

43



Figure 4. Comparisons of simulations using different algorithms (red) with the exact quantum dynamics results (blue) on the second system. (a) FSSH; (b) QTMF-FSSH using a constant Γ = 0.001; (c) QTMF-FSSH with 𝛾𝐹 as Eq (35) and 𝛾′ = 0; (d) QTMF-FSSH with 𝛾𝐹 as Eq (35) and 𝛾′ = 0.0002 according to Eq (36). Different shapes represent the probability on different channels (circle: T1, transmission in the lower surface; triangle: T2, transmission in the upper surface; square: R1, reflection in the lower surface).

44


Page 44 of 47

Page 45 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Comparisons of simulations using different algorithms (red) with the exact quantum dynamics results (blue) on the third system. (a) FSSH; (b) QTMF-FSSH using a constant Γ = 0.001; (c) QTMF-FSSH with 𝛾𝐹 as Eq (35) and 𝛾′ = 0; (d) QTMF-FSSH with 𝛾𝐹 as Eq (35) and 𝛾′ = 0.10 according to Eq (36). Different shapes represent the probability on different channels (circle: T1, transmission in the lower surface; triangle: R1, reflection in the lower surface; square: R2, reflection in the upper surface).

45



TOC Graphic

46


Page 46 of 47

Page 47 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Biographies: Lin Shen received his Ph.D. degree in Chemistry under the supervision of Prof. Wei-Hai Fang at Beijing Normal University in 2012. He continued his research as a postdoctoral associate with Prof. Hao Hu at the University of Hong Kong and Prof. Weitao Yang at Duke University. He is an associate professor in the College of Chemistry at Beijing Normal University since 2018. His current research focuses on methodology development for multiscale nonadiabatic dynamics combined with machine learning and enhanced sampling techniques to provide new computational tools in photochemistry and photobiology. Diandong Tang received his B.S. degree in Chemistry from Beijing Normal University in 2017. Currently he is a Ph.D. candidate under the supervision of Prof. Wei-Hai Fang at Beijing Normal University. His research focuses on methodology development on mixed quantum-classical and semiclassical dynamic simulations as well as applications on photophysical and photochemical processes. Binbin Xie received his Ph.D. degree in Chemistry under the supervision of Prof. Wei-Hai Fang at Beijing Normal University in 2017. Currently he is a research assistant at Hangzhou Institute of Advanced Studies at Zhejiang Normal University. His research focuses on understanding diverse photophysical and photochemical processes varying from small molecules to large biological systems by using electronic structure calculations and nonadiabatic dynamic simulations. Wei-Hai Fang earned his Ph.D. degree in theoretical chemistry under the supervision of Professor Ruo-Zhuang Liu at the Beijing Normal University in 1993. He was an Alexander von Humboldt Foundation Scholar and worked with Professor Sigid Peyerimhoff at Bonn University in the year of 1996 - 1998. Then, he joined the faculty of the Department of Chemistry, Beijing Normal University and becomes a full professor in 2000. Currently, his research interests focus mainly on direct ab initio non-adiabatic dynamics simulations of photo-induced chemical and biological processes as well photo-responses of functional materials.

47


Quantum Trajectory Mean-Field Method for Nonadiabatic Dynamics in

Recommend Documents