Letter pubs.acs.org/JPCL
Surface Hopping by Consensus Craig C. Martens* Department of Chemistry, University of California, Irvine, California 92697-2025, United States ABSTRACT: We present a new stochastic surface hopping method for modeling molecular dynamics with electronic transitions. The approach, consensus surface hopping (CSH), is a numerical framework for solving the semiclassical limit Liouville equation describing nuclear dynamics on coupled electronic surfaces using ensembles of trajectories. In contrast to existing techniques based on propagating independent classical trajectories that undergo stochastic hops between the electronic states, the present method determines the probabilities of transition of each trajectory collectively with input from the entire ensemble. The full coherent dynamics of the coupled system arise naturally at the ensemble level and ad hoc corrections, such as momentum rescaling to impose strict trajectory energy conservation and artificial decoherence to avoid the overcoherence of the quantum states associated with independent trajectories, are avoided.
T
justifying the FSSH approach in the context of the foundation provided by this rigorous semiclassical theory.5,18 In this Letter, we describe a new approach to surface hopping, which we call consensus surface hopping (CSH). Rather than adopting the FSSH method as the starting point for corrections and improvements, we build the approach anew on the rigorous foundation of the semiclassical Liouville formalism with the goal of solving the underlying Liouville equation as accurately as possible in the context of an ensemble of stochastic trajectories. The key new element of the method is that trajectories do not evolve as independent “lone wolves”, but collectively represent the evolution of the underlying density matrix elements of the Liouville equation. Both populations and coherences are accurately represented by the evolving ensemble. The stochastic hopping of the individual trajectories between quantum states is dictated by the properties of the instantaneous density matrix supported by the evolving trajectories, and thus involves a consensus among the members of the entire ensemble. We avoid the constraints imposed on the FSSH algorithm such as energy conservation of individual trajectories, and do not introduce empirical corrections such as artificial decoherence to the method. Rather, we allow the system populations and coherence to emerge naturally from the collective nonclassical nature of the trajectory ensemble as it represents the underlying quantum system. We begin by briefly reviewing the semiclassical Liouville approach to nonadiabatic dynamics.9−13 The quantum mechanical Liouville equation for the density operator ρ̂(t) is given by19
rajectory surface hopping is a popular and effective method for simulating molecular dynamics with quantum electronic transitions.1−8 It is a mixed quantum-classical approach that combines a trajectory-based description of the nuclear motion with stochastic transitions between coupled electronic states. Most current activity in developing and applying the surface hopping formalism builds on the seminal contribution by Tully,1 where the fewest switches surface hopping (FSSH) method was introduced and tested on several simple benchmark systems. The FSSH method is a simple, reasonable, practicalbut approximateapproach to treating the full coupled electronic-nuclear quantum dynamics of molecules undergoing nonadiabatic transitions in a mixed quantum-classical framework. From its introduction, it was recognized that the FSSH method has shortcomings and limitations. In particular, the approach relies on an ensemble of independent trajectories, each of which carries with it a pure state quantum density matrix to determine the probabilities of hopping. The phenomenon of decoherence does not arise naturally in this framework, and so FSSH is well-known to be “over-coherent”. In addition, the independence of the trajectories comprising the ensemble in FSSH suggests a strict imposition of energy conservation at the individual trajectory level. While ensuring a correct asymptotic energy budget, this feature of FSSH leads to inaccuracies such as “frustrated hops”, where strict energy conservation precludes channels that are open in a full quantum mechanical treatment. Following the introduction of FSSH, alternative and more rigorous methods for treating molecular dynamics with electronic transitions have been developed and implemented. An example is the mixed quantum-classical method based on a semiclassical limit of the multistate quantum Liouville equation, which have been pursued by several groups.9−17 More recently, considerable attention has been focused on deriving and © 2016 American Chemical Society
Received: May 31, 2016 Accepted: June 26, 2016 Published: June 26, 2016 2610
DOI: 10.1021/acs.jpclett.6b01186 J. Phys. Chem. Lett. 2016, 7, 2610−2615
Letter
The Journal of Physical Chemistry Letters
diagonal surfaces while at the same time collectively representing the off-diagonal quantum coherence. The populations of states 1 and 2 are together represented by a single ensemble of N trajectories, each of which is characterized by a point in phase space Γj(t) = (qj(t),pj(t)) and a binary integer σj(t), which can take on the values 1 or 0, indicating whether the trajectory is associated with quantum state 1 or 2, respectively. The phase space densitites are then given in terms of the trajectory ensemble by
d ρ ̂( t ) = [Ĥ , ρ (̂ t )] (1) dt ̂ ̂ where H is the Hamiltonian of the system and [H, ρ̂] denotes the commutator of the operators Ĥ and ρ̂. For dynamics on a single potential surface, the classical limit of eq 1 is the wellknown classical Liouville equation of nonequilibrium statistical mechanics,20 iℏ
∂ρ = {H , ρ} ∂t
(2)
ρ11(Γ, t ) =
where ρ(q,p,t) and H(q,p,t) are now functions of the classical phase space variables (q,p) and time t, and {H,ρ} is the Poisson bracket of H and ρ. This result can be derived systematically from eq 1 by performing a Wigner-Moyal expansion21,22 of the quantum mechanical Liouville equation. To lowest order in ℏ, this involves replacing commutators by Poisson brackets: [Â , B̂ ] → iℏ{A,B} + O(ℏ2). The semiclassical limit of eq 2 can be generalized to two coupled quantum states (e.g., diabatic electronic surfaces) with Hamiltonian and density matrix ⎛ Ĥ ⎛ ρ ̂ (t ) ρ ̂ (t ) ⎞ V̂ ⎞ 11 12 ⎟ρ ̂(t ) = ⎜ 11 ⎟ Ĥ = ⎜⎜ ⎜ ρ ̂ (t ) ρ ̂ (t )⎟ ⎟ ̂ ̂ ⎠ ⎝ 21 ⎝ V H22 ⎠ 22
∂t ∂ρ22 ∂t
= {H11 , ρ11} + {V , α} −
2V β ℏ
= {H22 , ρ22 } + {V , α} +
2V β ℏ
∑ σj(t )δ(Γ − Γj(t ))
1 N
∑ (1 − σj(t ))δ(Γ − Γj(t ))
ρ22 (Γ, t ) =
(10)
where hq and hp are the widths in position and momentum, respectively. The widths are determined using the methodology of density estimation, a technique for reconstructing underlying probability distributions from sampled data.29 For an ensemble of N trajectories sampled from a distribution with widths wi (i = 1, ..., d) in d dimensions, the density estimation methodology suggests the smoothing widths hi = wi [4/N(d+2)]1/(d+4). For the one degree of freedom problems treated below the phase space dimension d = 2, giving hq = wq (1/N)1/6 and hp = wp (1/ N)1/6, where we take (wq,wp) as the initial widths of the ensemble in phase space. We represent the coherence ρ12(Γ,t) in terms of the same ensemble as the populations; this is appropriate, as the coherence occupies the same region in phase space as the union of the state 1 and 2 populations. Unlike the populations, however, the coherence is a complex quantity and thus the coefficients of the trajectories representing the coherence are complex numbers. The populations and the real and imaginary parts of the coherence are expressed in terms of the smoothed trajectory ensemble as
(5)
∂β V = {H0 , β} − ωα + (ρ11 − ρ22 ) ∂t ℏ
(7)
(9)
j=1
⎛ (q − q )2 (p − pj )2 ⎞ 1 j ⎟ exp⎜ − − ⎜ 2πhqhp 2hq2 2hp2 ⎟⎠ ⎝
g (Γ − Γj) =
(4)
(6)
N
The δ functions represent the discrete points of the trajectories in phase space, while the coefficients σj denote which state the trajectory is currently occupying. In the numerical implementation involving finite trajetory ensembles, the δ functions are broadened using phase space Gaussians of the form
(3)
∂α 1 = {H0 , α} + ωβ + {V , ρ11 + ρ22 } ∂t 2
(8)
j=1
and
With the replacement of the quantum mechanical operators by the corresponding classical phase space functions this becomes a set of coupled classical-like Liouville equations:9−13,23−28 ∂ρ11
N
1 N
Here, we have written the coherence ρ12(Γ, t) = α(Γ, t) + iβ(Γ, t) in terms of its real and imaginary parts and have defined the average Hamiltonian H0 = (H11+H22)/2 and the transition frequency ω = (H11−H22)/ℏ. All higher order terms in ℏ have been discarded, leading to a classical-limit formalism that retains only the most important quantum corrections. Equations 4−7 highlight the role played by the coherence ρ12 in mediating the transfer of population between the quantum states: the coherence, together with the coupling V, acts as a population source or sink term. The evolving population densities, in turn, feed back on this process by acting as a source term for generating quantum coherence. We have presented the semiclassical Liouville formalism in a representation where the coupling between the electronic states appears as an off-diagonal potential energy term V, corresponding to a diabatic representation. The method can be equally well be implemented in the adiabatic representation, where coupling appears through off-diagonal kinetic energy coupling; 11,16,17 we illustrate both diabatic and adiabatic implementations below. The CSH method solves these coupled semiclassical Liouville equations using an ensemble of N trajectories. These trajectories undergo stochastic hops between the
ρ11(Γ, t ) =
ρ22 (Γ, t ) =
α (Γ , t ) =
β (Γ , t ) = 2611
1 N 1 N
N
∑ σj(t )g(Γ − Γj(t )) (11)
j=1 N
∑ (1 − σj(t ))g(Γ − Γj(t )) (12)
j=1 N
1 N
∑ αj(t )g(Γ − Γj(t ))
1 N
∑ βj(t )g(Γ − Γj(t ))
(13)
j=1 N j=1
(14) DOI: 10.1021/acs.jpclett.6b01186 J. Phys. Chem. Lett. 2016, 7, 2610−2615
Letter
The Journal of Physical Chemistry Letters In our implementation of CSH, the coefficients σj(t) are stochastic binary integers, while αj(t) and βj(t) (j = 1,2,...,N) are continuous real numbers. It is not necessary to make such a distinction: in a recent paper, we described an alternative approach to the general problem of quantum state hopping in condensed phases that represented both populations and coherence in terms of separate stochastic processes entangled with the nonstationary response of a nonequilibrium environment.30 The equations of motion for the trajectories Γj(t) and state parameters (σj(t), αj(t), βj(t)) (j = 1,2,..., N) are determined by subtituting eqs 11−14 into the semiclassical Liouville eqs 4−7. In the uncoupled (V = 0) case, the evolution of the populations reduces to purely classical Liouvillian dynamics corresponding to trajectory ensemble evolution under the appropriate electronic surface Hamiltonian. In the presence of coupling, two types of nonclassical terms appear. The first are Poisson bracket terms {V,α}, which occur symmetrically in the equations for both ρ11 and ρ22. These interactions modify the shape of the evolving distributions, but do not change the total state populations. Conservation of population under these terms results from the fact that the trace of a Poisson bracket vanishes for functions that vanish at the boundary of phase space. It is possible to include these nonclassical terms in the dynamics of the trajectories at the cost of computational complexity, and formally these terms are important for the energy budget of the coherent quantum-classical evolution.31 In the present work, we introduce an approximation and ignore these terms. The second type of couplings are the sink and source terms ±2Vβ/ℏ, which are responsible for the population transfer between states. These are included in the CSH implementation. The evolution of the populations is detemined by integrating numerically the ordinary differential equations for the trajectories and the coefficients. Each time step of duration Δt is divided into two parts. First, the coefficients are updated and then the phase space trajectories are propagated forward in time. The trajectories Γj(t) in the population ensemble each evolve under the Hamiltonian correponding to their currently occupied state according to the instantaneous value of σj(t): Hj(t ) = σj(t )H11 + [1 − σj(t )]H22
1 N
N
1 N
N
∑ σkΔσkg(Γ − Γk) = − ∑ k=1
k=1
⟨ρ11⟩j =
⟨β⟩j =
N
k=1
(Γj − Γk(t ))Δt
N
k=1
ℏ
k=1
(19)
1 N
N
∑ βk(t )g(Γj − Γk(t )) k=1
(20)
1 2V (Γj) ⟨β⟩j Δt ⟨ρ11⟩j ℏ
(21)
For trajectories currently evolving on surface 2, the corresponding result is Δσj(t ) = −
1 2V (Γj) ⟨β⟩j Δt ⟨ρ22 ⟩j ℏ
(22)
We emphasize the essential difference between this approach and the FSSH method. Here, the ensemble collectively determines the stochastic hopping probabilities of each of its members. The local densities ⟨ρ11⟩j and ⟨ρ2⟩j and the coherence ⟨β⟩j at point Γj depend on the entire ensemble of evolving trajectories Γk(t) (k = 1,2,...,N). Quantum transitions are thus determined by a consensus among the members of the ensemble representing the full entangled electronic-nuclear quantum state, rather than by the individual “lone wolf” trajectories of FSSH. It should also be emphasized that no momentum rescaling is performed when transitions occur. Individual trajectories representing a quantum system have no requirement to separately conserve energy,31 and they do not do so in this method. We therefore do not rescale the momenta of trajectories when they undergo hops between surfaces. Strict energy conservation at the individual trajectory level is invoked on physical grounds in surface hopping methods based on FSSH. However, momentum rescaling violates the phase space locality of the underlying coupled Liouville equations. This locality is apparent in eqs 11−12: the terms responsible for probability transition appear symmetrically in the equations, so that every element of population that is induced to leave surface 1 by the term −2V(q) β(q,p)/ℏ appears on surface 2 as +2V(q) β(q,p)/ ℏ at the same point (q,p) in phase space. Momentum rescaling to impose energy conservation induces a spurious shift of the probability in phase space. Conservation of energy is only important at the full state level, and emerges formally from the semiclassical equations of motion without momentum rescaling.31 These equations form the basis of a stochastic hopping algorithm. In particular, a random number ξ between 0 and 1 is generated at each time step and compared with the appropriate value of |Δσj| corresponding to the occupied state, which is interpreted as the hopping probability. The value of σj(t) is changed or kept at its current value depending on the outcome.
(15)
2V (qj)
∑ σkg(Γj − Γk)
Δσj(t ) = −
with a similar expression for surface 2. We can evaluate the left and right sides of these expressions at each of the trajectory points of interest, Γj, yielding coupled equations for the change in coefficients. This gives, for the surface 1 coefficients, N
N
1 N
The equation for updating the coefficients of trajectories currently on surface 1 becomes
2V (q) βk (t )g (Γ − Γk(t ))Δt ℏ
∑ σkΔσkg(Γj − Γk) = − 1 ∑
(18)
This approximation gives results in very close agreement to those obtained by solving the linear algebra problem directly, at much less numerical cost. Similarly, we define the value of the coherence at point j as
(16)
1 N
k=1
where ⟨ρ11⟩j, the local density at point Γj on surface 1, is given by
for j = 1,2···,N. To derive a probabilistic algorithm for updating the coefficients, we consider the subsets of trajectories on surface 1 and surface 2 separately. For surface 1, we have 1 N
N
∑ σkΔσkg(Γj − Γk) ≃ ⟨ρ11⟩j Δσj
βk (t )g (17)
for j = 1,2,...,N. In general, this presents a linear algebra problem for determination of the Δσj. We can simplify its solution by making the following approximation: 2612
DOI: 10.1021/acs.jpclett.6b01186 J. Phys. Chem. Lett. 2016, 7, 2610−2615
Letter
The Journal of Physical Chemistry Letters
Figure 1. Time dependence of the state 1 population ρ11(t) and the real and imaginary parts of the coherence ρ12(t) for Tully’s single crossing model. The CSH results (blue dashed lines) are compared with exact wavepacket calculations (black solid lines). (a) ℏk = 10. (b) ℏk = 30.
= ∫ ψi(q,t) ψj*(q,t) dq. Nearly quantitative agreement is observed for these cases, with the CSH and quantum time series often falling on top of each other. It is particularly noteworthy that the coherences are very accurately modeled. The collective nature of the full ensemble faithfully tracks the coherence of the entangled nuclear-electronic dynamics. This is even true for the ℏk = 30 case, where significant coherence persists in the system after the population transfer is complete until the end of the simulation. (We note that FSSH methods that incorporate artificial dephasing would destroy this important feature of the electronic-nuclear dynamics.) In Figure 2 the final state population as a function of initial momentum ℏk is presented. The agreement is excellent except
The electronic coherence evolves in parallel withand coupled tothe dynamics of the population densities. A similar analysis that includes the approximate neglect of terms that leave the trace of ρ12 unchanged yield expressions that describe the evolution of the coefficients over the time step Δt: Δαj = ω(Γj)βjΔt
(23)
⎡ ⎤ 1 Δβj = ⎢ −ω(Γj)αj + V (Γj)(2σj − 1)⎥Δt ⎣ ⎦ ℏ
(24)
We emphasize that no artificial decoherence is added to the evolving system. Coherence and its decay is treated accurately through the collective nature of the method as highlighted by eq 20 and demonstrated below. (We note that adding pure dephasing to a coupled system induces population and energy relaxation, and so in principle should be avoided for systems that should satisfy energy conservation.) Together with integration of Hamilton’s equations for the Γj(t), eqs 21−24 constitute the numerical implementation of CSH. We apply the CSH method to two model systems originally proposed by Tully as benchmark problems and adopted universally by the surface hopping community as test cases.1 The first is Tully’s single crossing system, which is a onedimensional model corresponding to surfaces undergoing a single crossing in the diabatic representation. The system has been treated in many previous publications. We adopt the potentials and parameters from Tully’s original paper1 and treat ensembles of 2000 trajectories sampled randomly from an initial minimum uncertainty Gaussian distribution with spatial width wq = 1.0 and with mean initial position qo = −6. A range of ensembles with varying initial momenta p = ℏk are generated and propagated using the CSH method. The surface hopping results are compared with corresponding quantum wave packet calculations performed in the diabatic representation using the method of Kosloff.32 In Figure 1 we show the initial state population and coherence dynamics computed for Tully’s single crossing model by the CSH method and compare with exact wave packet calculations. Results are calculated and presented in the diabatic represention. In Figure 1a we show result for an initial momentum ℏk = 10 while panel b shows the results for ℏk = 30. The CSH quantities shown are the traces of the corresponding distribution functions, so that ρ11(t) = (1/N) ∑j σj(t), Re ρ12(t) = (1/N)∑j αj(t), and Im ρ12(t) = (1/ N)∑jβj(t). The corresponding exact quantum results are determined from the evolving numerical wave packets: ρij(t)
Figure 2. Final state 1 population ρ22,final versus initial wave vector k for Tully’s single crossing model. Comparison of CSH results (blue dots) with exact wave packet calcuations (solid black line). Results obtained and presented in the diabatic electronic representation.
near threshold, where the CSH algorithm underestimates the quantum population transfer. We also tested the CSH method on Tully’s dual crossing model.1 The diabatic potentials and parameters were taken from Tully’s original paper. For this system, we implement the CSH methodology in the adiabatic electronic representation. The diagonal adiabatic potentials are found by computing the eigenvalues of the diabatic potential matrix as a function of coordinate. The off-diagonal coupling responsible for inducing population transfer is now given by the purely imaginary kinetic term − iℏd(q)p/m, where the nonadiabatic coupling element d(q) = −ϕ′(q)/2, and the mixing angle ϕ(q) is given in terms of the diabatic coupling V(q) and potentials E1(q) and E2(q) by tan ϕ(q) = 2V(q)/[E1(q) − E2(q)].11,16,17,31 As in the single crossing case, we include only terms with nonzero phase space trace, as described above. Writing now the adiabatic coherence as ρ+− = α + iβ we have 2613
DOI: 10.1021/acs.jpclett.6b01186 J. Phys. Chem. Lett. 2016, 7, 2610−2615
Letter
The Journal of Physical Chemistry Letters
Figure 3. Time dependence of the lower adiabatic state population ρ−−(t) and the real and imaginary parts of the coherence ρ+−(t) for Tully’s dual crossing model. The CSH results (blue dashed lines) are compared with exact wavepacket calculations (black solid lines). (a) ℏk = 30. (b) ℏk = 40.
Δσj(t ) = −
pj 2 d(qj) ⟨α⟩j Δt (upper surface) ⟨ρ++ ⟩j m
(25)
Δσj(t ) = −
pj 2 d(qk ) ⟨α⟩j Δt (lower surface) ⟨ρ−− ⟩j m
(26)
pj ⎤ ⎡ Δαj = ⎢ω(Γj)βj + d(qj) (2σj − 1)⎥Δt ⎦ ⎣ m
(27)
Δβj = −ω(Γj)αjΔt
(28)
results are nearly quantitative, while deviations become more pronounced as the energy is lowered. In the results shown here, we have used ensembles with N = 2000 trajectories, which are the same size as in Tully’s original paper.1 For ensembles as small as N = 300, good but noisier results are obtained. The numerical expense of CSH is higher than FSSH, as the local values of the phase space densities are represented by the ensemble of trajectories rather than resorting to the independent trajectory approximation of conventional surface hopping. Evaluating these local quantities requires extra effort, and large enough ensembles must be employed to represent the underlying states accurately. We considered one-dimensional systems in the present study, and found that we could perform all computations interactively on a notebook computer. The scaling of the method with dimensionality will be investigated in future studies. In this letter, we have presented a new approach to modeling molecular dynamics with electronic transitions using ensembles of trajectories. The method, consensus surface hopping (CSH), is based on representing the full coupled time-dependent electronic-nuclear density matrix using an ensemble of trajectories that evolve classically and undergo stochastic transitions between the electronic states. The ensemble collectively models the quantum population transfer as well as the evolving electronic coherence. In sharp contrast with the independent trajectory foundation of the fewest-switches surface hopping (FSSH) approach introduced by Tully and adopted by many groups, the CSH method determines individual hopping probabilities holistically from the entire ensemble. Our approach does not presuppose the FSSH method as a framework on which to build corrections, such as addition of artificial decoherence or the imposition of strict energy conservation at the individual trajectory level. Rather, we derive anew a stochastic approach that is built firmly on the rigorous foundation of the full electronic-nuclear dynamics of the density matrix in the semiclassical limit. The resulting CSH algorithm honors the underlying properties of the semiclassical Liouville dynamics, such as phase space locality of population transfer and definition of quantum coherence at the ensemblerather than individual trajectorylevel. The method models the dynamics of population transfer and the generation, evolution, and decay of quantum coherence using a “by consensus” algorithm, in which the trajectory ensemble collaboratively represents the underlying phase space functions that correspond in the semiclassical limit to quantum populations and coherences, and which in turn guide the system evolution through their collective effect on the evolution of each individual member of the ensemble of trajectories.
The equations for the binary integers σj are solved using the stochastic hopping algorithm, while the coherence coefficients are integrated in time, as above. The classical trajectories are propagated on the appropriate adiabatic potentials between hops. The initial condition is a wave packet on the upper state with the same width as for Tully’s single hopping problem and an initial position of qo = −10. In Figure 3 we show the lower state population and adiabatic representation coherence as a function of time for two initial wave packets. In Figure 3a) results for ℏk = 30 are shown, while Figure 3 b) gives results for ℏk = 40. The results are compared with wave packet calculations that have been transformed to the adiabatic basis. We again see excellent agreement between the CSH and exact results. It is noteworthy that even subtle time-dependent features of both the populations and coherences are faithfully captured by the CSH method. In Figure 4 we show the final lower state population as a function of initial momentum ℏk. At high momentum, the
Figure 4. Final lower state population ρ−−,final vs initial wave vector k for Tully’s dual crossing model. Comparison of CSH results (blue dots) with exact wave packet calcuations (solid black line). Results obtained and presented in the adiabatic electronic representation. 2614
DOI: 10.1021/acs.jpclett.6b01186 J. Phys. Chem. Lett. 2016, 7, 2610−2615
Letter
The Journal of Physical Chemistry Letters
(9) Martens, C. C.; Fang, J. Y. Semiclassical-Limit Molecular Dynamics on Multiple Electronic Surfaces. J. Chem. Phys. 1997, 106, 4918−4930. (10) Donoso, A.; Martens, C. C. Simulation of Coherent Nonadiabatic Dynamics Using Classical Trajectories. J. Phys. Chem. A 1998, 102, 4291−4300. (11) Donoso, A.; Martens, C. C. Semiclassical Multi-State Liouville Dynamics in the Adiabatic Representation. J. Chem. Phys. 2000, 112, 3980−3989. (12) Donoso, A.; Kohen, D.; Martens, C. C. Simulation of Nonadiabatic Wavepacket Interferometry Using Classical Trajectories. J. Chem. Phys. 2000, 112, 7345−7354. (13) Donoso, A.; Martens, C. C. Classical Trajectory-Based Approaches to Solving the Quantum Liouville Equation. Int. J. Quantum Chem. 2002, 90, 1348−1360. (14) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919. (15) Hanna, G.; Kapral, R. Nonadiabatic Dynamics of Condensed Phase Rate Processes. Acc. Chem. Res. 2006, 39, 21−27. (16) Ando, K. Non-Adiabatic Couplings in Liouville Description of Mixed Quantum-Classical Dynamics. Chem. Phys. Lett. 2002, 360, 240−242. (17) Ando, K.; Santer, M. Mixed Quantum-Classical Liouville Molecular Dynamics Without Momentum Jump. J. Chem. Phys. 2003, 118, 10399−10406. (18) Nielsen, S.; Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Surface Hopping Dynamics. J. Chem. Phys. 2000, 112, 6543. (19) Cohen-Tannoudji, C.; Diu, B.; Laloe, F. Quantum Mechanics; Wiley: New York, 1977. (20) McQuarrie, D. A. Statistical Mechanics; HarperCollins: New York, 1976. (21) Takahashi, K. Distribution Functions in Classical and Quantum Mechanics. Prog. Theor. Phys. Suppl. 1989, 98, 109−156. (22) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford University Press: Oxford, 1995. (23) Riga, J. M.; Martens, C. C. Simulation of Environmental Effects on Coherent Quantum Dynamics in Many-Body Systems. J. Chem. Phys. 2004, 120, 6863−6873. (24) Riga, J. M.; Fredj, E.; Martens, C. C. Quantum Vibrational StateDependent Potentials for Classical Many-Body Simulations. J. Chem. Phys. 2005, 122, 174107. (25) Riga, J. M.; Martens, C. C. Environmental Decoherence of Many-Body Quantum Systems: Semiclassical Theory and Simulation. Chem. Phys. 2006, 322, 108−117. (26) Riga, J.; Fredj, E.; Martens, C. Simulation of Vibrational Dephasing of I2 in solid Kr Using the Semiclassical Liouville Method. J. Chem. Phys. 2006, 124, 064506. (27) Hogan, P. A.; Fredj, E.; Martens, C. C. Simulation of Vibrational Dephasing in Liquid Water Using the Semiclassical Liouville Method. Chem. Phys. Lett. 2011, 510, 208−211. (28) Roman, E.; Martens, C. C. Semiclassical Liouville Method for the Simulation of Electronic Transitions: Single Ensemble Formulation. J. Chem. Phys. 2004, 121, 11572. (29) Silverman, B. W. Density Estimation for Statistics and Data Analysis; Chapman and Hall: London, 1986. (30) Martens, C. C. Communication: Fully Coherent Quantum State Hopping. J. Chem. Phys. 2015, 143, 141101. (31) Martens, C. C. Nonadiabatic Dynamics in the Semiclassical Liouville Representation: Locality, Transformation Theory, and the Energy Budget. Chem. Phys. Submitted for publication. (32) Kosloff, R. Propagation Methods for Quantum Molecular Dynamics. Annu. Rev. Phys. Chem. 1994, 45, 145. (33) Donoso, A.; Martens, C. C. Quantum Tunneling Using Entangled Classical Trajectories. Phys. Rev. Lett. 2001, 87, 223202. (34) Donoso, A.; Zheng, Y.; Martens, C. C. Simulation of Quantum Processes Using Entangled Trajectory Molecular Dynamics. J. Chem. Phys. 2003, 119, 5010.
Electronic coherence and the energy budget of the system emerge naturally in CSH at the ensemble level, and ad hoc corrections such as artificial decoherence and implicitly nonlocal momentum rescaling of transferred population densities are avoided. A single classical trajectory has no real quantum mechanical meaning, except in its relation with other trajectories.23,25,33,34 Concepts such as uncertainty, nonlocality, and entanglement are manifestations of this fact in other contexts. In the present method, the coherence of a quantum system emerges as the interrelationships between trajectories in the ensemble. In the CSH approach, problems such as the overcoherence of FSSH surface hopping do not arise. The CSH method was applied to benchmark systems introduced by Tully1 and employed extensively as test problems for nonadiabatic simulation methodology. Excellent agreement was obtained for both time-dependent populations and the evolving electronic coherence of the full coupled system. For nonadiabatic dynamics near threshold, where quantum effects become more important, the agreement was not as quantitative. In the present implementation, additional nonclassical terms that do not induce transitions but modify the dynamics between transitions were neglected for the sake of computational simplicity. The potential for these additional effects to give improved accuracy will be explored in future publications.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS
■
REFERENCES
We are grateful to Shaul Mukamel, Vladimir Mandelshtam, and Greg Ezra for helpful conversations. Stimulating discussions at the UCI Liquid Theory Lunch (LTL) and the Telluride Science Research Center (TSRC) are also acknowledged. This work is an outgrowth of research supported by the National Science Foundation under CHE-0614005 and CHE-0802913.
(1) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061. (2) Tully, J. C. Perspective: Nonadiabatic Dynamics Theory. J. Chem. Phys. 2012, 137, 22A301. (3) Subotnik, J. E.; Shenvi, N. A New Approach to Decoherence and Momentum Rescaling in the Surface Hopping Algorithm. J. Chem. Phys. 2011, 134, 024105. (4) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-Induced Surface Hopping. J. Chem. Phys. 2012, 137, 22A545. (5) 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. (6) Curchod, B. F. E.; Tavernelli, I. On Trajectory-Based Nonadiabatic Dynamics: Bohmian Dynamics versus Trajectory Surface Hopping. J. Chem. Phys. 2013, 138, 184112. (7) Chen, H.-T.; Reichman, D. R. On the Accuracy of Surface Hopping Dynamics in Condensed Phase Non-Adiabatic Problems. J. Chem. Phys. 2016, 144, 094104−10. (8) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011−2015. J. Phys. Chem. Lett. 2016, 7, 2100−2112. 2615
DOI: 10.1021/acs.jpclett.6b01186 J. Phys. Chem. Lett. 2016, 7, 2610−2615
