Surface Hopping Dynamics beyond Nonadiabatic Couplings for

Letter Cite This: J. Phys. Chem. Lett. 2018, 9, 1097−1104

pubs.acs.org/JPCL

Surface Hopping Dynamics beyond Nonadiabatic Couplings for Quantum Coherence Jong-Kwon Ha, In Seong Lee, and Seung Kyu Min* Department of Chemistry, School of Natural Science, Ulsan National Institute of Science and Technology (UNIST), 50 UNIST-gil, Ulsan 44919, Republic of Korea S Supporting Information *

ABSTRACT: Description of correct electron−nuclear couplings is crucial in modeling of nonadiabatic dynamics. Within traditional semiclassical or mixed quantum−classical dynamics, the coupling between quantum electronic states and classical nuclear trajectories is governed by nonadiabatic coupling vectors coupled to classical nuclear momenta. This enables us to develop a very powerful nonadiabatic dynamics algorithm, namely, surface hopping dynamics, which can describe the splitting of nuclear wave packets and detailed balance. Despite its efficiency and practicality, it suffers from the lack of quantum decoherence due to incorrect accounts for the electron−nuclear coupling. Here we present a new surface hopping algorithm based on the exact electron−nuclear correlation from the exact factorization of molecular wave functions. This algorithm demands comparable computational costs to existing surface hopping methods. Numerical simulations with two-state models and a multidimensional multistate realistic molecule show that the electron−nuclear coupling beyond the nonadiabatic coupling terms can describe the quantum coherence properly.

D

classical trajectory independently, and its stochastic behavior provides statistically meaningful results from a large number of trajectories. However, this independent trajectory raises an important problem, the so-called “lack of quantum decoherence”. Interactions between a quantum system and a large reservoir make the system recover the purity of the quantum state, i.e., the off-diagonal element of the reduced density matrix should disappear eventually. Recently, thanks to the exact factorization of a full molecular wave function, we can obtain exact equations for coupled electronic and nuclear dynamics.45−47 Here the correct electron−nuclear coupling includes not only the traditional coupling between electronic nonadiabatic coupling vectors (NACVs) and classical nuclear momenta but also the effect of quantum nuclei, i.e., nonlocal nuclear density. A coupledtrajectory approach based on this electron−nuclear coupling has been developed, and it has been tested for several model systems as well as realistic molecules with an ab initio quantum chemistry program.48−50 As a result, the quantum (de)coherence can be captured correctly with the coupled-trajectory approach based on exact factorization. However, apparently this coupled-trajectory approach requires more computational cost compared to an independent-trajectory approach. In this sense, it would be very intriguing to develop an algorithm of surface hopping dynamics with a correct electron−nuclear coupling term derived from exact factorization. In this Letter, we show

ynamics involving important quantum (quasi-)particles coupled to large (classical) enviornments plays a key role in many intriguing phenomena. For example, exciton energy transfer in photosynthesis,1−3 exciton transport in photovoltaic materials,4−7 electron transfer in various chemical systems,8−12 and polaron transport in thermoelectric materials13 are under high investigation in various fields due to their potential applications. For the design of highly efficient photovoltaic/ photosynthetic/thermoelectric materials, theoretical understanding of the coupling between quantum and classical variables is very demanding. The most accurate way is to deal with the entire system quantum mechanically based on the time-dependent Schrödinger equation or the quantum Liouville equation. However, this is only feasible for a system with a few degrees of freedom. For a large-scale system such as an extended molecular system or solids, one of the most promising methods is the so-called mixed quantum−classical (MQC) approach where we treat fast particles quantum mechanically while we consider slowly moving particles with classical trajectories. In trajectory-based approaches, the correct capture of the coupling between quantum and classical particles is crucial. So far, many trajectory-based MQC methods have been developed such as quantum−classical Liouville equations,14−21 the linearized path-integral approaches,22−24 nonadiabatic Bohmian dynamics,25−27 as well as numorous versions of Ehrenfest dynamics28−31 and surface hopping dynamics.32−44 Probably Tully’s fewest switch surface hopping (FSSH) algorithm and its derivatives are by far the most powerful schemes due to their simplicity and practicality of algorithm. Trajectory-based surface hopping dynamics allows us to run a © XXXX American Chemical Society

Received: January 8, 2018 Accepted: February 13, 2018 Published: February 13, 2018 1097

DOI: 10.1021/acs.jpclett.8b00060 J. Phys. Chem. Lett. 2018, 9, 1097−1104

Letter

The Journal of Physical Chemistry Letters

rithm is fully compatible with modern electronic structure methods to obtain adiabatic solutions. However, one major missing ingredient in FSSH is description of quantum decoherence. When system variables are coupled to bath variables, persistent couplings between them make the system state recover purity of the state. If the trajectory (I) propagates on the adiabatic PES EL, then the quantum decoherence implies the collapse of the electronic wave function on the adiabatic state L (i.e., C(I) = 1). L Unfortunately, it does not usually happen in FSSH dynamics because ρ(I) LL does not change after nonadiabatic coupling regions, i.e., ρ̇(I) LL = 0. The correct capture of (de)coherence in FSSH dynamics is crucial to understand a situation where the coherence becomes important as in exciton energy transfer of light-harvesting systems.1,53 For many decades, people have tried to modify the FSSH algorithm to describe the behavior of coherence correctly.36,37,40,41,54,55 One simple way is ad hoc implementation of decoherence correction based on analysis of the decoherence time,41 which requires an adjustment of electronic coefficients during propagation, or one can add a dissipative Hamiltonian at some level for the electronic time evolution from the quantum Liouville equation37,40 or Meyer− Miller−Stock−Thoss formalism.54 The decoherece correction can be calculated either from a semiclassical propagation of Gaussian wave packets to include effects of quantum nuclei or from fluctuations in nuclear degrees of freedom. The algorithms37,54 focus on the off-diagonal elements of the reduced density matrix, and thus, they may result in the lack of internal consistency in representing state occupations while accurate quantum coherences can be obtained from averages. Moreover, we want to mention that one can derive a nonHermitian Hamiltonian, especially a skew-Hermitian Hamiltonian for the dissipative part, from the quantum Liouville equation based on the Born−Huang expansion.40 Even though this is a true nature of coupled electron−nuclear dynamics, the algorithm again requires an additional normalization because the dissipative terms for each state do not compenstate for each other. In the following, we show that our new FSSH scheme based on the direct electronic equation of motion from the exact factorization descibes the quantum coherence while satisfying the internal consistency without any renormalization of the reduced density matrix. According to the exact factorization of a full molecular wave function, we can obtain the exact equations for coupled electronic and nuclear dynamics.45−47 The exact factorization framework enables us to derive coupled-trajectory MQC (CTMQC) dynamics, whereas the coupling between nuclear trajectories is inherited by the spatial shape of a nuclear density.48,49 In CT-MQC dynamics, the electronic motion coupled to a classical trajectory R(I )(t ) is governed by

an efficient algorithm for surface hopping dynamics and its applications for several chemical systems including a realistic molecule to describe quantum decoherence based on the exact factorization. In the surface hopping algorithm, the nuclear motion is treated by a set of classical trajectories, while the electronic wave function evolves by the time-dependent Schrödinger equation for electrons. The nuclear trajectory propagates on a given adiabatic potential energy surface (PES), but it can switch to another PES at any time with a given transition probability or “hopping” probability. In particular, the electronic and nuclear degrees of freedom evolve according to (I ) iℏΦ̇ (r, t ) = ĤBO(r, R(I )(t ))Φ(I )(r, t ) and M ν R̈ ν( I ) (t) = −∇ ν E L( I ) (t), respectively, where r = {r1, r2, ...} and R = {R1, R 2, ...} represent electronic and nuclear degrees of freedom, Φ and Ĥ BO are the electronic wave function and Born−Oppenheimer (BO) Hamiltonian, and Mν is the nuclear mass for the νth nucleus. Here R(I )(t ) is the Ith nuclear trajectory, which moves on the adiabatic PES E(I) L , and L denotes the current or “force” state at time t. Note that any quantity with the superscript (I) means the quantity defined at the (I)th nuclear trajectory, i.e., Φ(I )(r, t ) = Φ(r, t ; R(I )(t )). By expanding the time-dependent electronic wave function Φ on the basis of the adiabatic wave function ϕK at a given nuclear configuration R(I )(t ), i.e., Φ(I )(r, t ) = ∑K CK(I )(t )ϕK(I )(r) with (I) (I) (I) Ĥ (I) BOϕK = EK ϕK , the electronic equation of motion can be transformed equivalently into i (I ) CK̇ (t ) = − EK(I )(t )CK(I )(t ) − ℏ

∑ ∑ d(KJI)ν·Ṙ (νI)(t )C(J I)(t ) J

ν

(1)

where CK’s are the expansion coefficients of the time-dependent electronic wave function and dKJν = ⟨ϕK|∇νϕJ⟩ and Ṙ ν are the NACV and the nuclear velocity for the νth nuclei, respectively. The last term in eq 1 is a nonadiabatic coupling term that describes electron−nuclear couplings from a classical-trajectory approximation. The key ingredient in the FSSH algorithm is to determine the hopping probability from one adiabatic PES, EL, to another one, EK, by considering the time variation of the electronc state population ρ(I) LL(t) (diagonal element of the (I) (I) adiabatic electronic density matrix ρ(I) LK = CL *CK ). According to the original FSSH stochastic algorithm, the transition probability P(I) L→K within time interval [t,t + Δt] can be represented as PL(I→) K =

(I ) I) ̇ (I ) ⎤ 2ℜ⎡⎣ρLK (t ) ∑ν d(LK ν·R ν (t )⎦ (I ) (t ) ρLL

Δt .

We usually perform the FSSH dynamics with a large set of nuclear trajectories and obtain averaged quantities over all trajectories. Due to its stochastic behavior to determine the “force” state, the FSSH dynamics can describe spatial separation or branching of nuclear wave packets in excited-state dynamics successfully. From the practical point of view, the FSSH dynamics has several benefits: (1) The FSSH trajectories are independent of each other, and thus, they can be handled efficiently with parallel computing. (2) In addition, instead of computing NACVs, which requires huge computational costs, we can calculate nonadiabatic coupling matrix elements (NACMEs) σLK ≡ ∑ν dLKν · Ṙ ν(t) exploiting the relation d σIJ(t ) = ⟨ΦI (t )| dt ΦJ (t )⟩.51,52 (3) Moreover, the FSSH algo-

i (I ) CK̇ (t ) = − EK(I )(t )CK(I )(t ) − ℏ +

∑∑ J

ν

1 ∇ν |χ | M ν |χ |

∑ ∑ d(KJI)ν(t )·Ṙ (νI)(t )C(J I)(t ) J

ν

· {f (JIν)(t ) R(I )(t )

− f (KIν)(t )}|C(J I )(t )|2 CK(I )(t )

(2)

where χ (R , t ) is a time-dependent nuclear wave function and fIν is related to a phase of the Ith state, which comes from the spatial derivative of CI(t). The detailed derivation of eq 2 can be found elsewhere.49 Compared to eq 1, we have additional 1098


Letter


Figure 1. (a) BOPESs with nonadiabatic coupling (shifted by −0.3) between them. (b) pI for the Ith state from DISH-XF and FSSH algorithms and ⟨ρII⟩ from quantum dynamics (QD). (c) ⟨|ρ12|2⟩ for DISH-XF, FSSH, and QD. (d) ρII for a chosen single trajectory in DISH-XF and FSSH dynamics. DISH-XF n represents a trajectory on the nth state at the end. We present only one trajectory for FSSH because ρII’s are all similar in FSSH. (e,f) Position (R) and potential energy of nuclei at 38.4 fs from FSSH and DISH-XF, respectively. The distribution of nuclear trajectories is depicted as a histogram.

except for the electronic evolution, which requires additional computations for the decoherence correction. First, to calculate the quantum momentum, DISH-XF exploits fictitious trajectories coupled to each independent (real) trajectory. We generate an auxiliary trajectory R K (t ′) as R K (t ′) = R L(t ′) on each BO state that is not equal to the running state (K ≠ L) when the BO population ρKK becomes nonzero at time t′. We assume that the auxiliary trajectory follows uniform velocity motion during the time interval [t′,t′ + Δt], and we obtain the velocity Ṙ K (t ′) by uniform rescaling of ̇ t ′) to satisfy the energy conservation law as R(

electron−nuclear coupling terms, the so-called decoherence term, depending on the spatial shape of nuclear distribution, i.e., quantum momentum ∇ν |χ (R , t )| /|χ (R , t )|. In the previous coupled-trajectory approach, we calculate the quantum momentum from a set of classical nuclear trajectories by giving a Gaussian-shaped nuclear density to each nuclear configuration. Thus, all nuclear trajectories are coupled via this quantum momentum. Due to its coupled nature, unlike FSSH dynamics, we should run dynamics with multiple nuclear trajectories simultaneously to construct a nuclear density at each time. As a result, the coupled-trajectory algorithm can describe the correct behavior of quantum coherence successfully with a smaller number of nuclear trajectories compared to FSSH dynamics.48−50 Now we provide a practical FSSH algorithm exploiting the decoherence correction from exact factorization, what we call decoherence-induced surface hopping based on exact factorization (DISH-XF), which requires as low of computational costs as the conventional FSSH algorithm. One should note that this DISH-XF algorithm is not a variation of another decoherenceinduced surface hopping (DISH) method developed by Prezhdo and co-workers.36 In the FSSH algorithm, we solve electronic problems with the BO Hamiltonian at a given nuclear configuration R to obtain BO electronic energies EI as well as NACMEs σIJ among electronic states. In addition, the nuclear force is determined solely by the gradient of the PES for a certain force state L, i.e., −∇νEL. In our DISH-XF approach, we follow the same procedure as the FSSH algorithm

∑ ν

1 Mν Ṙ (KIν)2(t ) = 2

∑ ν

1 Mν Ṙ (νI )2(t ) + EL(I )(t ) − EK(I )(t ) 2 (3)

If the kinetic energy term for the Kth state becomes negative, the velocity Ṙ (KI )(t ′) is set to zero. Finally, if ρ(I) KK becomes (numerically) zero, we destroy the auxiliary trajectory on the Kth state. In this evolution scheme, the averaged trajectory over the auxiliary trajectories mimics the Ehrenfest trajectory. During the time evolution of auxiliary trajectories, we can calculate the quantum momentum ∇ν|χ|/|χ| at the R(I )(t ) by constructing Gaussian nuclear density centered at each R (KI )(t ) with a uniform variance σ as 2 (I ) 2 2 χK (R (t ))| = NK ∏ν exp(−|R ν − R Kν| /2σ ) where the am1099


Letter

The Journal of Physical Chemistry Letters plitude NK satisfies |χK (R , t )|2 /|χ (R , t )|2 = ρKK . Then, the

Instead, DISH-XF shows clear decoherence as quantum dynamics. We can understand effects of the decoherence term from the time evolution of diagonal density matrix elements for a single nuclear trajectory (Figure 1d). From eqs 2 and 4, we obtain

quantum momentum at R(I )(t ) yields ∇ν |χ | |χ |

=− R(I )(t )

1 [R (νI )(t ) − ⟨R (νI )(t )⟩] 2σ 2

(4)

ρII(̇ I ) (t ) =

(I) (I) where ⟨R(I) ν (t)⟩ = ∑KρKK(t)RKν(t). Next, the phase term fIν(t) is derived as a time integration of the BO force according to the previous works,48,49 i.e., fIν(t) = −∫ t ∇νϵIBO(t′) dt′. Alternatively, the phase term can be approximated to a time integration of momentum change at the Kth state because −∇νϵBO I (t′)Δt is the momentum change of the νth nucleus at the Ith state for an infinitesimal time change Δt. Thus, the phase term at time t becomes

∑∑ J

ν

(f (JIν)(t )

1 [R (νI )(t ) − ⟨R (νI )(t )⟩]· M νσ 2 − f (IIν)(t ))ρJJ(I )(t )ρII(I )(t )

(6)

in the absence of nonadiabatic couplings. Hence, the decoherence term can describe the change of ρII(t) after coupling regions. The behavior of the decoherence term can be explained by the energy conservation law because motions of auxiliary trajectories are governed by eq 3. In this model situation, for example, if the force state is equal to 1 after the AC (solid lines in Figure 1d), R(t) is larger than ⟨R(t)⟩ because the velocity of the auxiliary trajectory at state 2 is smaller to that of the real trajectory at state 1. Menawhile, from eq 5, the phase term at time t can be reduced to the difference between momentum at time t and t0, where t0 is the time that the auxiliary trajectory at state 2 is created. After the AC, the energy of state 2 increases and the energy of state 1 decreases. Then, the velocity of the trajectory at state 2 becomes smaller, while the velocity of the trajectory at state 1 becomes larger (i.e., f1 > 0, f 2 < 0, and f 2 − f1 < 0), which yields ρ̇11 > 0 and ρ̇22 < 0. For the case with the running state as 2 after the AC, ρ̇22 > 0 and ρ̇11 < 0 because R(t) < ⟨R(t)⟩ and f1 − f 2 > 0 for the same reason. Thus, in the DISH-XF algorithm, a population change of the force state becomes positive, and it becomes zero as the system is completely decohered. In addition, it is interesting to analyze a decoherence time (τJK) at this stage where the decoherence time is a decay time constant of absolute value of off-diagonal densty matrix elements (ρJK) contributed from the decoherence term as37,40

f (KIν)(t ) = f (KIν)(t − Δt ) + Mν[Ṙ (KIν)(t ) − Ṙ (KIν)(t − Δt )] (5)

where Ṙ (I) Kν(t) is obtained from eq 3. When the auxiliary trajectory is created on the Kth state at t′, the phase term fKν(t′) is set to zero initially. With this approximation, we can obtain phase terms at time t without calculating gradients of BO energy surfaces. Summarizing our DISH-XF algorithm, the step-by-step procedure is the following: (i) With a given nuclear configuration R(I )(t ) at a time t, we calculate a set of BO (I) potential energies E(I) I (t) as well as NACMEs σIJ (t). (ii) We evolve the time-dependent BO coefficients C(I) (t) (or density I matrix elements ρIJ(I)(t)) according to eq 2, while the decoherence term can be obtained by eqs 4 and 5. (iii) Nuclear dynamics is governed by the Newtonian equation. (iv) We determine the stochastic transition probability in the same way as the original FSSH algorithm. This approach requires the same amount of quantum mechanical calculation as FSSH with a few additional computations for the decoherence term. Moreover, one can easily prove that the norm, ∑J ρ(I) JJ (t), is always conserved in the DISH-XF approach. As the first application of the DISH-XF algorithm, we prepare a one-dimensional charge transfer model system that was proposed by Shin and Metiu.56 In this model Hamiltonian, we consider two BO potential energy surfaces (BOPESs) that have an avoided crossing (AC) between them as in Figure 1a. Here we prepare the initial nuclear trajectories at near R = −6 au on the first excited state, and after some time, the nuclear wave packet splits into different running states passing the AC with a partial electronic transition. Then, we compare two types of BO populations, pI and ⟨ρII⟩, as explained in the Computational details. pI is obtained from the numbers of trajectories running on the Ith state, while ⟨ρII⟩ is the averaged BO population ρ(I) II over all trajectories. In addition, we calculate 2 ⟨|ρIJ|2⟩, which is the averaged |ρ(I) IJ | over all trajectories to indicate the decoherence. ⟨|ρIJ|2⟩ becomes finite if the system is cohered, and ⟨|ρIJ|2⟩ becomes zero if it is decohered. While both FSSH and DISH-XF describe the correct nuclear BO population pI after a single nonadiabatic coupling region (Figure 1b), the behaviors of electronic decoherence ⟨|ρ12|2⟩ are completely different (Figure 1c). As a reference, quantum dynamics simulation shows that the electronic states are coherent (i.e., ρ12 ≠ 0) as the nuclear wave packet is in the AC region while the electronic decoherence appears (i.e., ρ12 = 0) after the nuclear wave packet splits in nuclear configurational space. The FSSH algorithm shows completely wrong results because we have finite ρ11 and ρ22 even after the AC region.

d ln|ρLK | 1 =− τLK dt 1 =∑ (R ν − ⟨R ν⟩) ·(2f Jν − f Lν − f Kν)ρJJ 2 2 M νσ J ,ν

(7)

Equation 7 implies that the difference in nuclear configurations and phases for each state promotes the decoherence. For the two-state model, the decoherence time becomes 1/τ12 = (R − ⟨R⟩)( f1 − f 2)(ρ11 − ρ22)/2Mσ2. Therefore, when two trajectories become far apart while propagating on different PESs, |τ12| decreases gradually. When all trajectories evolve almost on top of each other as in the high-temperature limit, we recover the Ehrenfest regime (τ12 = ∞) because R = ⟨R⟩ and f1 = f 2. Also, we want to point out the R, ⟨R⟩, and f L are determined not by just instant quantities like energies or forces at a given time but by the entire history of trajectories. Figure 1e,f shows the nuclear positions with potential energies (R,ϵ) at 38.4 fs from FSSH and DISH-XF, respectively. The colors represent ρ11 for each trajectory, i.e., the red is for ρ11 = 1 (ρ22 = 0), the purple is for ρ11 = 0 (ρ22 = 1), and the orange is for finite ρ11 and ρ22. The trajectories from both FSSH and DISH-XF show the same behavior (at least after one AC region). However, in FSSH, the electronic coherence remains after the AC, and we can expect that eventually this will also cause the wrong nuclear dynamics when the trajectories pass ACs multiple times. 1100


Letter


Figure 2. (a) BOPESs and nonadiabatic coupling between them. (b) pI for the Ith state from the DISH-XF and FSSH algorithm. (c) ⟨ρII⟩ for the Ith state from the DISH-XF and FSSH algorithm. QD results are shown for the reference in (b) and (c). (d) ρII from DISH-XF and FSSH for a single trajectory that passes the avoided crossing five times and hops to the ground state after the last avoided crossing.

base, which is regarded as a reference molecule for excited-state molecular dynamics.43,44,59 We consider three BO states, S0, S1, and S2, as well as NACVs among them based on MRCI/SA-4CASSCF(6,4)/6-31G* calculation with the Columbus program package.60 We conduct our DISH-XF dynamics with 80 trajectories on S2 at 300 K after thermal equilibrium (see the Supporting Information). Typical trajectories, with a constant temperature of 300 K, meet the S2/S1 conical intersection elongating the C−N bond and hop to the S1 state (t < 10 fs). After passing the first conical intersection, a pyramidalization at the N side starts to occur and pass the S1/S0 coupling region rotating around the C−N bond (20−60 fs), and finally, they are relaxed to the S0 ground state with isomerized geometry within 100 fs. In Figure 3a,b, we show BO populations for DISH-XF and FSSH. Both DISH-XF and FSSH algorithms can show similar BO population exchanges based on pI (Figure 3a). Trajectories meet a conical intersection between S2 and S1 at around 3 fs and meet another nonadiabatic coupling region between S1 and S0 from 20 to 60 fs. The first electronic transition occurs very rapidly, while the second electronic transition occurs gradually. As in the previous NaI case, ⟨ρII⟩ from FSSH and DISH-XF show completely different behavior (Figure 3b), and DISH-XF shows the internal consistency (pI ≈ ⟨ρII⟩). In addition, we calculate the coherence ⟨|ρIJ|2⟩ for both FSSH and DISH-XF in Figure 3c. We have two distinctive peaks (or bands) at 3 fs for ⟨|ρ12|2⟩ and 20−60 fs for ⟨|ρ01|2⟩ in DISH-XF. The first narrow peak appears because most trajectories pass through the conical intersection between S1 and S2 almost simultaneously ,while the second broad band appears because trajectories encounter another coupling region between S0 and S1 at different times.

Next, we test a situation with multiple ACs to model the femtosecond pump−probe expriment for a gaseous sodium iodide (NaI) molecule.57,58 By a pump pulse, a nuclear wave packet goes up to the excited state, and it meets an AC region between dissociative (neutral) and attractive (ionic) potential energy curves (Figure 2a). If a nuclear wave packet keeps moving on the upper state, it encounters the AC multiple times. In this case, the same behavior of pI and ⟨ρII⟩, known as the internal consistency, is evidence of the correct capture of decoherence.35 Figure 2b,c shows BO populations for DISHXF, FSSH, and quantum dynamics. Both DISH-XF and FSSH results show nice agreement with the result from quantum dynamics after a single AC (∼0.5 ps) as a previous model system. However, only DISH-XF can describe the correct long time behavior after multiple ACs for both pI and ⟨ρII⟩. The reason is clear if we analyze ρII for a single trajectory that passes the AC multiple times. Figure 2d shows ρII for a single trajectory that passes the AC five times and goes to the dissociative channel finally. Whenever the trajectory passes the AC, ρII(t) in FSSH gives an error accumulation due to the lack of decoherece, while DISH-XF recovers the decohered state, i.e., ρ11 = 0 and ρ22 = 1, after some time and shows the dramatic behavior that shows ρ11 = 1 and ρ22 = 0 after the final AC (t > 2 ps). The wrong decription of coherence indeed leads to the wrong nuclear dynamics. This result shows very clear limitation of FSSH and advantage of DISH-XF in terms of the description of decoherence, which comes from the correlation between nuclear motion and the electronic states. Finally we test a realistic molecule with multiple BO states in a multidimensional nuclear configurational space. Here we choose a molecule CH2NH+2 , the smallest protonated Schiff 1101


Letter


Figure 3. (a) BO populations pI from DISH-XF and FSSH. (b) BO populations ⟨ρII⟩ from DISH-XF and FSSH. (c) Electronic coherence ⟨|ρIJ|2⟩. (d) Diagonal electronic density matrix elements ρII for a single trajectory. (e,f) Potential energy profile in time for DISH-XF and FSSH, respectively. Black represents the potential energy of a molecule, while red, green, and blue lines represent the ground, first, and second excited states, respectively. We use 80 nuclear trajectoreis sampled from the ground-state molecular dynamics simulation at 300 K.

However, FSSH shows the overcohered behavior. ⟨|ρ12|2⟩ remains finite after the first peak at 3 fs, and the finite BO populations cause the coherence ⟨|ρ02|2⟩ even though there is no significant coupling between S0 and S2 states. Figure 3d shows the change of diagonal element ρII of a single trajectory. Again, while DISH-XF shows complete decoherences, FSSH shows finite ρII’s during the dynamics even though the nuclear dynamics (Figure 3a) for both FSSH and DISH-XF are similar to each other because FSSH chooses the correct running state (Figure 3e,f). However, we can expect, as in the previous sodium iodide case and Figure 3c, that results may be different for a long time simulation after multiple passages of nonadiabatic coupling regions. As a summary, we have proposed a new algorithm for surface hopping dynamics to account for decoherence in nonadiabatic phenomena. This new algorithm has two major aspects. First, in theoretical aspects, the coupling between electrons and nuclei was derived from the exact factorization of a molecular wave function that is the exact solution of a time-dependent Schrödinger equation of a molecule. The electron−nuclear coupling contains not only nonadiabatic couplings between electronic states and classical nuclear momenta but also decoherence terms from quantum nuclear momenta. This decoherence term depends on the current position of a given nuclear trajectory R(I )(t ) as well as a history or the time integration of a phase term f(I) Jν (t). As explained in eq 6, the decoherence term plays a role to give an additional change of electronic populations only if the system is cohered, and it becomes zero naturally as the system is decohered. As

implemented on top of surface hopping dynamics, the population of a running state approaches 1. This feature leads to internal consisitency with fractions of trajectories, pI, and average occupations, ⟨ρII⟩, which a surface hopping algorithm should meet. Moreover, the norm of electronic populations, a trace of a reduced density matrix, is always conserved automatically. Hence, an additional renormalization procedure is unnecessary. Second, in computational aspects, DISH-XF requires almost the same computational cost as the original surface hopping dynamics. The FSSH approach requires the computation of many-body electronic states ϕ K(I) with corresponding electronic energies E(I) K for a given nuclear (I) position R(I )(t ) as well as NACVs, d(I) IJν(t), or NACMEs, σIJ . According to eqs 3−5, the decoherence term in DISH-XF can be calculated approximately from the same ingredients exploiting the concept of auxiliary trajectories. Thus, DISHXF can use the benefits of independent-trajectory approaches such as high parallel efficiency for larger molecular systems. Finally, we emphasize that the implementation of the current approach is straightforward on top of the FSSH dynamics.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.8b00060. Computational details (PDF) 1102


Letter


■

(17) Kapral, R. Progress in the Theory of Mixed Quantum-Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129−157. (18) Grunwald, R.; Kim, H.; Kapral, R. Surface-Hopping Dynamics and Decoherence with Quantum Equilibrium Structure. J. Chem. Phys. 2008, 128, 164110. (19) Jang, S. Nonadiabatic Quantum Liouville and Master Equations in the Adiabatic Basis. J. Chem. Phys. 2012, 137, 22A536. (20) Hsieh, C.-Y.; Kapral, R. Analysis of the Forward-Backward Trajectory Solution for the Mixed Quantum-Classical Liouville Equation. J. Chem. Phys. 2013, 138, 134110. (21) Kelly, A.; Markland, T. E. Efficient and Accurate Surface Hopping for Long Time Nonadiabatic Quantum Dynamics. J. Chem. Phys. 2013, 139, 014104. (22) Dunkel, E. R.; Bonella, S.; Coker, D. F. Iterative Linearized Approach to Nonadiabatic Dynamics. J. Chem. Phys. 2008, 129, 114106. (23) Bonella, S.; Coker, D. F. LAND-map, a Linearized Approach to Nonadiabatic Dynamics Using the Mapping Formalism. J. Chem. Phys. 2005, 122, 194102. (24) Huo, P.; Coker, D. F. Consistent Schemes for Non-Adiabatic Dynamics Derived from Partial Linearized Density Matrix Propagation. J. Chem. Phys. 2012, 137, 22A535. (25) Curchod, B. F. E.; Tavernelli, I.; Rothlisberger, U. TrajectoryBased Solution of the Nonadiabatic Quantum Dynamics Equations: an On-the-Fly Approach for Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2011, 13, 3231−3236. (26) Curchod, B. F. E.; Tavernelli, I. On Trajectory-Based Nonadiabatic Dynamics: Bohmian Dynamics Versus Trajectory Surface Hopping. J. Chem. Phys. 2013, 138, 184112. (27) Prezhdo, O. V.; Brooksby, C. Quantum Backreaction Through the bohmian Particle. Phys. Rev. Lett. 2001, 86, 3215−3219. (28) Prezhdo, O. V. Mean Field Approximation for the Stochastic Schrödinger Equation. J. Chem. Phys. 1999, 111, 8366−8377. (29) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-Field Dynamics with 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. (30) Alonso, J. L.; Clemente-Gallardo, J.; Cuchí, J. C.; Echenique, P.; Falceto, F. Ehrenfest Dynamics is Purity Non-Preserving: A Necessary Ingredient for Decoherence. J. Chem. Phys. 2012, 137, 054106. (31) 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. (32) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (33) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011−2015. J. Phys. Chem. Lett. 2016, 7, 2100− 2112. (34) Jasper, A. W.; Nangia, S.; Zhu, C.; Truhlar, D. G. Non-BornOppenheimer Molecular Dynamics. Acc. Chem. Res. 2006, 39, 101− 108. (35) Granucci, G.; Persico, M. Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping. J. Chem. Phys. 2007, 126, 134114. (36) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-Induced Surface Hopping. J. Chem. Phys. 2012, 137, 22A545. (37) Subotnik, J. E.; Shenvi, N. A New Approach to Decoherence and Momentum Rescaling in the Surface Hopping Algorithm. J. Chem. Phys. 2011, 134, 024105. (38) 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. (39) 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.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Seung Kyu Min: 0000-0001-5636-3407 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (2016R1C1B2015103).

■

REFERENCES

(1) Scholes, G. D.; et al. Using Coherence to Enhance Function in Chemical and Biophysical Systems. Nature 2017, 543, 647−656. (2) Romero, E.; Novoderezhkin, V. I.; van Grondelle, R. Quantum Design of Photosynthesis for Bio-Inspired Solar-Energy Conversion. Nature 2017, 543, 355−365. (3) Kaucikas, M.; Maghlaoui, K.; Barber, J.; Renger, T.; Van Thor, J. J. Ultrafast Infrared Observation of Exciton Equilibration from Oriented Single Crystals of Photosystem II. Nat. Commun. 2016, 7, 13977. (4) Blancon, J. C.; et al. Extremely Efficient Internal Exciton Dissociation Through Edge States in Layered 2D Perovskites. Science 2017, 355, 1288−1292. (5) Nah, S.; Spokoyny, B.; Stoumpos, C.; Soe, C.; Kanatzidis, M.; Harel, E. Spatially Segregated Free-Carrier and Exciton Populations in Individual Lead Halide Perovskite Grains. Nat. Photonics 2017, 11, 285−288. (6) Huang, J.; Yuan, Y.; Shao, Y.; Yan, Y. Understanding the Physical Properties of Hybrid Perovskites for Photovoltaic Applications. Nat. Rev. Mater. 2017, 2, 17042. (7) Chistyakov, A. A.; Zvaigzne, M. A.; Nikitenko, V. R.; Tameev, A. R.; Martynov, I. L.; Prezhdo, O. V. Optoelectronic Properties of Semiconductor Quantum Dot Solids for Photovoltaic Applications. J. Phys. Chem. Lett. 2017, 8, 4129−4139. (8) Chábera, P.; et al. A Low-Spin Fe (III) Complex with 100-ps Ligand-to-Metal Charge Transfer Photoluminescence. Nature 2017, 543, 695−699. (9) Xia, Z.; He, C.; Wang, X.; Duan, C. Modifying Electron Transfer Between Photoredox and Organocatalytic Units via Framework Interpenetration for β-Carbonyl Functionalization. Nat. Commun. 2017, 8, 361. (10) Fukushima, T.; Gupta, S.; Rad, B.; Cornejo, J. A.; Petzold, C. J.; Chan, L. J. G.; Mizrahi, R. A.; Ralston, C. Y.; Ajo-Franklin, C. M. The Molecular Basis for Binding of an Electron Transfer Protein to a Metal Oxide Surface. J. Am. Chem. Soc. 2017, 139, 12647−12654. (11) Lou, Z.; Fujitsuka, M.; Majima, T. Two-Dimensional AuNanoprism/Reduced Graphene Oxide/Pt-Nanoframe as Plasmonic Photocatalysts with Multiplasmon Modes Boosting Hot Electron Transfer for Hydrogen Generation. J. Phys. Chem. Lett. 2017, 8, 844− 849. (12) Zhang, M.; Wang, L.; Shu, S.; Sancar, A.; Zhong, D. Bifurcating Electron-Transfer Pathways in DNA Photolyases Determine the Repair Quantum Yield. Science 2016, 354, 209−213. (13) Craven, G. T.; Nitzan, A. Electrothermal Transistor Effect and Cyclic Electronic Currents in Multithermal Charge Transfer Networks. Phys. Rev. Lett. 2017, 118, 207201. (14) Kapral, R. Surface Hopping from the Perspective of QuantumClassical Liouville Dynamics. Chem. Phys. 2016, 481, 77−83. (15) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919−8929. (16) Nielsen, S.; Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Surface Hopping Dynamics. J. Chem. Phys. 2000, 112, 6543−6553. 1103


Letter

The Journal of Physical Chemistry Letters (40) Gao, X.; Thiel, W. Non-Hermitian Surface Hopping. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2017, 95, 013308. (41) Granucci, G.; Persico, M.; Zoccante, A. Including Quantum Decoherence in Surface Hopping. J. Chem. Phys. 2010, 133, 134111. (42) Craig, C. F.; Duncan, W. R.; Prezhdo, O. V. Trajectory Surface Hopping in the Time-Dependent Kohn-Sham Approach for ElectronNuclear Dynamics. Phys. Rev. Lett. 2005, 95, 163001. (43) Tapavicza, E.; Tavernelli, I.; Rothlisberger, U. Trajectory Surface Hopping Within Linear Response Time-Dependent Density-Functional Theory. Phys. Rev. Lett. 2007, 98, 023001. (44) Doltsinis, N. L.; Marx, D. Nonadiabatic Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2002, 88, 166402. (45) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. Phys. Rev. Lett. 2010, 105, 123002. (46) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Correlated ElectronNuclear Dynamics: Exact Factorization of the Molecular WaveFunction. J. Chem. Phys. 2012, 137, 22A530. (47) Gidopoulos, N. I.; Gross, E. K. U. Electronic Non-Adiabatic States: Towards a Density Functional Theory Beyond the BornOppenheimer Approximation. Philos. Trans. R. Soc., A 2014, 372, 20130059. (48) Min, S. K.; Agostini, F.; Gross, E. K. U. Coupled-Trajectory Quantum-Classical Approach to Electronic Decoherence in Nonadiabatic Processes. Phys. Rev. Lett. 2015, 115, 073001. (49) Agostini, F.; Min, S. K.; Abedi, A.; Gross, E. K. U. QuantumClassical Nonadiabatic Dynamics: Coupled- vs Independent-Trajectory Methods. J. Chem. Theory Comput. 2016, 12, 2127−2143. (50) 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. (51) Hammes-Schiffer, S.; Tully, J. C. Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657−4667. (52) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods. J. Phys. Chem. Lett. 2015, 6, 4200−4203. (53) Ramanan, C.; Ferretti, M.; van Roon, H.; Novoderezhkin, V. I.; van Grondelle, R. Evidence for Coherent Mixing of Excited and Charge-Transfer States in the Major Plant Light-Harvesting Antenna, LHCII. Phys. Chem. Chem. Phys. 2017, 19, 22877−22886. (54) 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. (55) Shenvi, N.; Subotnik, J. E.; Yang, W. Phase-Corrected Surface Hopping: Correcting the Phase Evolution of the Electronic Wavefunction. J. Chem. Phys. 2011, 135, 024101. (56) Shin, S.; Metiu, H. Nonadiabatic Effects on the Charge Transfer Rate Constant: A Numerical Study of a Simple Model System. J. Chem. Phys. 1995, 102, 9285−9295. (57) Rose, T. S.; Rosker, M. J.; Zewail, A. H. Femtosecond RealTime Probing of Reactions. IV. The Reactions of Alkali Halides. J. Chem. Phys. 1989, 91, 7415−7436. (58) Engel, V.; Metiu, H. A Quantum Mechanical Study of Predissociation Dynamics of NaI Excited by a Femtosecond Laser Pulse. J. Chem. Phys. 1989, 90, 6116−6128. (59) Barbatti, M.; Aquino, A. J. A.; Lischka, H. Ultrafast Two-Step Process in the Non-Adiabatic Relaxation of the CH2NH+2 Molecule. Mol. Phys. 2006, 104, 1053−1060. (60) Lischka, H. et al. COLUMBUS, an Ab Initio Electronic Structure Program, release 7.0; 2015.

1104


Surface Hopping Dynamics beyond Nonadiabatic Couplings for

Recommend Documents