Mapping State Space to Quasiclassical Trajectory Dynamics in

Article pubs.acs.org/JPCA

Mapping State Space to Quasiclassical Trajectory Dynamics in Coherence-Controlled Nonadiabatic Simulations for Condensed Phase Problems Guohua Tao* and Na Shen School of Advanced Materials, Peking University Shenzhen Graduate School, Shenzhen, China 518055 Shenzhen Key Laboratory of New Energy Materials by Design, Peking University, Shenzhen, China 518055 S Supporting Information *

ABSTRACT: Recently a coherence controlled (CC) approach to nonadiabatic dynamics was proposed by one of the authors based on the mapping between the decomposed classical state space and different types of nuclear dynamics. Here we elaborate the state-space decomposition scheme and the corresponding state-space-todynamics mapping of the CC approach in a general high-dimensional framework. In the CC formalism, dynamical properties such as the full electronic matrix can be evaluated by means of the ensemble of trajectories in the active state space, which consists of single-state domains and coherence domains. The feasibility of the state space decomposition and related mappings and the performance of the CC approach are demonstrated by the implementation to benchmark problems of nonadiabatic molecular dynamics in condensed phase including the spin-boson model and the excitation energy transfer problem in photosynthesis. The results obtained from the CC approach are in reasonably good agreement with exact or benchmark calculations, and it is also shown that the CC approach satisfies the detailed balance approximately and is capable of efficiently describing condensed phase nonadiabatic molecular dynamics at reasonable accuracy.

I. INTRODUCTION Nonadiabatic dynamics could play a vital role in many fundamental molecular processes in condensed phase,1−4 such as charge transfer and energy transfer in chemical reactions, energy materials or photosynthesis. Despite extensive research having been done, the accurate and efficient description on the coupled electronic-nuclear motions and the complicated many-body interactions in the presence of the environment in condensed phase nonadiabatic dynamics is still challenging. So far the most popular theoretical approaches for nonadiabatic dynamics simulations are perhaps the mixed quantum-classical (MQC) methods, such as Ehrenfest5 and surface hopping,6 which are essentially based on classical mechanics and therefore have practical advantages of treating complicated molecular systems over the full quantum approaches. Recent development7−14 of the MQC methods has made significant progresses in improving the accuracy and expanding the scope of their applications. However, the major problem with the MQC methods is the improper treatment of quantum coherence,2 due to the inconsistent mixed quantum-classical treatment. For example, Ehrenfest dynamics features a mean field model for nuclear motions in which nuclei always propagate on an averaged potential energy surface of multiple states. Therefore, the classical treatment of nuclear dynamics results in the overestimation of electronic coherence7 and unphysical asymptotic limits. The popular surface hopping treatment involves an essentially Born−Oppenheimer (BO) trajectory switched © XXXX American Chemical Society

among multiple states and the inconsistence in treating electronic-nuclear dynamics results in an effective overestimation of nuclear coherence.15 One way to consistently treat electronic and nuclear degrees of freedom (DOFs) at equal footing is the Meyer and Miller (MM) model,16 which was proved to be exact at the quantum mechanical level.17,18 Practically the MM model could be implemented at the semiclassical or quasiclassical levels, depending on the approximations to be made. Although the MM nuclear dynamics is still Ehrenfest, electronic coherence is well described either by the semiclassical treatment via the semicalssical initial value representation (IVR),19 or by the windowing technology in the symmetrical quasiclassical (SQC) approach.20,21 The MM formalism, in particular when combined with the SQC approach, has been demonstrated to be extremely efficient and reasonably accurate in describing nonadiabatic dynamics in both simple models and complex molecular systems.21−26 As in the conventional Ehrenfest dynamics, the mean field description on nuclear motions in the original MM model or the recent SQC approach may suffer from the numerical instability due to the instant negative electronic population variables. To better describe nuclear dynamics, we proposed a coherence Received: October 31, 2016 Revised: February 8, 2017 Published: February 8, 2017 A

DOI: 10.1021/acs.jpca.6b10936 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A controlled (CC) scheme27 to consistently switch nuclear dynamics to be BO-like in the single state regime or to be Ehrenfest in the coherence regime. The CC approach treats electronic DOFs the same way as the SQC approach, but creates a mapping between the decomposed classical state space with different types of nuclear dynamics, resulting in a simple picture of quasi-classical treatment on coupled electronic-nuclear motions. How to incorporate the BO-like single state dynamics and the coherent Ehrenfest dynamics into a consistent framework is a long-standing problem. Fruitful results were achieved from early efforts made in solving the problem. For example, Prezhdo and Rossky28 suggested a switching criterion based on the perturbative measurement to the deviation of Ehrenfest dynamics from the active single state dynamics. Truhlar and co-workers7,29 introduced a decoherence term into electronic dynamics to account for the decay of the system to pure state associated with a parametric decay time. However, the implementation of these sophisticated ideas to general complicated problems might be difficult. For example, the convergence of the perturbative switching criteria may become problematic if the single state dynamics are significantly different from each other, whereas explicitly adding decoherence terms to nonadiabatic dynamics involving many states may require a set of decay times with many possibilities of combination. The CC approach27 suggests a switching criterion solely based on the time-dependent population variables from the selfconsistent propagation of coupled electronic-nuclear trajectories with no perturbation or adding terms to the true dynamics at the quasiclassical level. Therefore, the generalization of the CC algorithm to complex molecular systems is straightforward and practical. In previous work,27 the CC approach was applied to several benchmark systems including two or three states but only one nuclear DOF. Here we demonstrate that the state space decomposition and the mappings between electronic density matrix, nuclear dynamics and the decomposed state space are applicable in high dimensional systems, and the CC approach could retain its simplicity in treating nonadiabatic dynamics in large systems with many states coupled with thermal bath. The switching criterion for high dimensional systems is provided, the evaluation of the full electronic matrix is presented, and the validity of the detailed balance is discussed. This paper is organized as the follows: Section II provides the theoretical framework of the CC approach, followed by a description on models and computational details in Section III, Section IV shows the results and discussions, and Section V concludes.

limit. Vkk and Vkl are the diagonal and off-diagonal Hamiltonian matrix elements, and F is the number of states. The electronic dynamics can be cast into a compact form by introducing a classical electronic variable matrix B ≡ [Bkl], i.e. F

̇ = iℏBmk

F

F

with 1 [(xkxl + pk pl ) + i(xkpl − pk xl)] 2

Bkl =

(2b)

We denote the classical electronic variable of population, coherence, and flux as follows: hk = Re[Bkk ], ckl = Re[Bkl ], fkl = Im[Bkl ]

(2c)

The nuclear dynamics is Hamilton dynamics, i.e., dQ k dt

=

∂H ∂Pk

(3a)

dPk ∂H =− dt ∂Q k

(3b)

Equation 2 and eq 3 therefore constitute classical equations of motion (EOM) of coupled electronic-nuclear dynamics. The advantage of the MM dynamics is that both electronic and nuclear DOFs are treated consistently on an equal footing, because the electronic variables are no longer discrete numbers but continuous variables. The price to pay for this consistency is that the instant population of a particular state, say nkk, is in general a real number and could become negative for unoccupied (or insignificant) states. Furthermore, nuclear dynamics in eq 3 is Ehrenfest, which is governed by a mean-field force including contributions from all states. Consequently, the MM trajectory may be overcorrected by the force contributed from the insignificant states with negative population, and the algorithm could become unstable. Considerable improvement on the MM model may be achieved when electronic population is requantized using a symmetrical windowing technology proposed in the SQC approach.21 By introducing a histogram window function for both initial and final electronic variables, i.e., Wk(hk , N ) =

⎞ 1 ⎛⎜ Δn h − |hk − Nk|⎟ ⎝ ⎠ 2 Δn

(4a)

with Nk (0 or 1) the electronic quantum number, Δn = 2γ the window width, and h(t) the Heaviside function, i.e. 0 (t < 0) h(t ) = , the electronic state associated with the 1 (t ≥ 0) system can be determined by using the windowing function for the full system, i.e.,

F

∑ nkkVkk(Q) + ∑ ∑ nklVkl(Q) k=1

(2a)

l=1

II. THEORY We start from the coupled electronic-nuclear Hamiltonian in the standard MM form:16,17 P2 + /(x , p, Q, P) = 2μ

∑ [Bml Vlk − VmlBlk ]

{

k=1 k≠1

(1)

with μ the nuclear reduced mass, (x, p) the classical electronic coordinates and momenta, and (Q, P) the nuclear 1 coordinates and momenta. Here nkk ≡ 2 (xk2 + pk2 ) − γ and

F

W ({hk }, {Nk}) =

1 (x x 2 k l

nkl ≡ + pk pl ) are defined as the population of state k and the coherence between state k and state l, respectively, and the parameter γ represents for the classical measure of zero point energy, which is 0.5 in the full quantum-classical correspondence

∏ Wk(hk , Nk) k=1

(4b)

with the set of electronic population variables {hk} = (h1, h2, ...hk, ...hF) and electronic quantum numbers {Nk} = (N1, N2, ...Nk, ...NF). For example, if the system is in state k for a particular B


Article

The Journal of Physical Chemistry A

Figure 1. State space decomposition for a three-state model. Colored cubes correspond to the single-state domains for state 1 (D11, red), state 2 (D22, golden), state 3 (D33, blue), and coherence domains (other colors), which are mapped to Born−Oppenheimer-like and Ehrenfest dynamics (see the text), respectively. (a) Standard scheme; (b) unprojected scheme, i.e., all coherence domains forms a single union.

the quantization of nuclear dynamics according to eq 5, the mean-field force includes only the contributions within the subspace. The switching criterion provided by eq 6 assures the nuclear dynamics outside the coherence domain to be correctly described in a single state flavor, which greatly improves the trajectory stability in the nonadiabatic dynamics simulations (see below). It would be more illustrative to introduce the CC approach in terms of a mapping between the classical decomposed state space and different types of nuclear dynamics, as we did in the previous work.27 The MM model eq 1 establishes a mapping between discrete quantum states and continuous classical electronic variables, the latter of which forms a continuous classical state space with the basis consisting of the set of classical population variables, i.e., (h1, h2, ...hk, ...hF). Quantum mechanically the state k corresponds to a single point of (0, 0, ...hk = 1, ...0), while classically the electronic variables may not be exactly the integer 0 or 1, and a histogram bin or window function is used to define a regime in the state space, which is associated with a particular state. Such trajectory histograms in the state space could produce the distribution of electronic states, which can be used to evaluate the electronic population and related properties. Figure 1 illustrates the state space of a three-state model, in which each single state is associated with a cubic domain defined by the window function eq 4 and the window width to be the quantum value γ = 0.5. The state space is therefore decomposed into single-state domains (red, yellow, and blue for three single states in a three-state model), and coherence domains (other colors), which is defined in a similar way as the single-state domain but with a window width γc (here we adopt γc = γ = 0.5). The union of the single-state domains and coherence domains are called the active space, and the complementary space is called the void space (not colored). We suggested that such a picture of the state space decomposition would be illuminating for understanding nonadiabatic dynamics. For example, it appears reasonable that the symmetrical windowing technology used in the SQC approach is superior to the earlier asymmetric binning treatment,21 because each state point in the state space should be on equal footing and associated with a single-state domain with an identical size so that the microscopic reversibility holds. Moreover, the decomposition scheme in Figure 1 not only provides a different mapping between the full electronic density matrix and the classical state space, but offers a straightforward and consistent way to map the classical state space to different

MM trajectory, i.e., {Nk} = (0, 0, ...Nk = 1, ...0), we should expect that Wk(hk, Nk) = 1 and Wm(hm, Nm) = 1 with Nk = 1, and Nm = 0 for all m ≠ k. The ensemble average of trajectories results in the desired time-dependent distribution of electronic states. The application of symmetrical windowing functions in the SQC approach is equivalent to quantize the time dependent electronic population variables for each MM trajectory when performing related calculations. The SQC treatment recovers the detailed balance and considerably improves the evaluation of electronic dynamics. However, the electronic population variables themselves and the EOM bear no change at all, therefore the nuclear dynamics is still Ehrenfest. The recently proposed CC scheme27 takes one step forward by quantizing the nuclear dynamics eq 3 according to the timedependent electronic population variables {hk}, i.e., ⎧ 1, occupied state nkk = W ({hk }, {Nk}) = ⎨ ⎩ 0, unoccupied state ⎪

⎪

(5)

If the system is on the occupied state k = occ, then nkk → 1, and nll → 0 for all l ≠ k. The CC nuclear dynamics therefore reads: Ṗ = −

∂V (Q) ∂/ = − kk − ∂Q ∂Q

F

∑

ckl

l ≠ k , k = occ

∂Vkl(Q) ∂Q

(6a)

and k=1

Ṗ = −

∂V (Q) ∂/ = − ∑ hk kk − ∂Q ∂Q F

k=1 l≠k

∑ ∑ ckl F

F

∂Vkl(Q) ∂Q (6b)

when the system is in the normal coherent superposition of multiple states, i.e., nkk = 0 for all k. If the system enters the super coherent regime, i.e., ∃ dim({nkk}) > 1, nkk = 1, (or nkk = 1 for more than one state), the nuclear EOM follows: F

Ṗ = −

∂V (Q) ∂/ = − ∑ nkk hk kk − ∂Q ∂Q k=1

F

F

∑ ∑ nkknllckl k=1 l≠k

∂Vkl(Q) ∂Q (6c)

Note that here nkk is 1 or 0 (see eq 5), by contrast, hk and ckl are continuous variables. Equation 6a is essentially Born− Oppenheimer dynamics modified by interstate interactions, whereas eq 6b is Ehrenfest. Equation 6c is a projected Ehrenfest dynamics that means if the system is projected into a subspace in C


Article


the potential term n1V1. The nuclear dynamics is dominated by the BO term, i.e., the first term on the right-hand side of eq 6a. Interactions with other states are considered but in an equally distributed way, for example c12V12 and c21V21 are assigned to state 1 and 2, respectively (see eq 6a). If the trajectory enters the cube centered at (1,1,0), the nuclear dynamics is governed by contribution from the state 1 and state 2 and their coherence (see eq 6c with n11 = n22 = 1, n33 = 0), and we map this cube to the potential term n12V12. For the cube centered at (1,1,1), all three states contribute, and we map nuclear dynamics to standard Ehrenfest (eq 6b), unlike eq 6c, which projects Ehrenfest dynamics to the corresponding subspace. There is no potential term in eq 1 to be mapped to this cube since the third order term n111V111 is neglected. One special cube centered at (0,0,0) is connected to every single state and constitutes a full coherence domain in a sense that for a subspace consisting of two states, every dual off-diagonal electronic matrix elements or potential terms (say ρ21 or n21V21) may be mapped into this domain. Therefore, we map the corresponding nuclear dynamics to standard Ehrenfest eq 6b. Alternatively, we can map the entire coherence domain (including all coherence cubes) to Ehrenfest dynamics (also to all off-diagonal electronic matrix elements), resulting in an unprojected CC scheme (Figure 1b). In general, the one-to-one mapping between the decomposed state space and nuclear dynamics does not exist or is more complicated and thus more difficult to be implemented than the proposed CC scheme. The mapping between the state space and nuclear dynamics is not always well-defined if the single-state domains and the coherence domain overlap or are separated, except for some limiting cases. For example, when the coherence window width γc → ∞, the coherence domain overlaps with single-state domains, and nuclear dynamics is Ehrenfest everywhere, exactly the case in the standard SQC approach.21 In the opposite limit, i.e., γc → 0, there is no coherence domain. Nuclear dynamics is always Born−Oppenheimer-like, and the CC approach becomes to the augmented image (AI) version of the multistate trajectory (MST) approach (unpublished results). Once the state space is decomposed, the individual offdiagonal element can be evaluated in the following way:

types of nuclear dynamics. The motivation of the state decomposition scheme proposed in the CC approach27 is to establish a one-to-one mapping between the state space and the nuclear dynamics, i.e., ⎛ BO‐like, eq.6a ⎞ ⎛ noncoherence domain ⎞ ⎟⎟ ⎜ ⎟ ↔ ⎜⎜ ⎝ coherence domain ⎠ ⎝ Ehrenfest, eq.6b ⎠

(7)

Here the noncoherence domain is the union of the single-state domain and the void space. The mapping between the diagonal electronic density matrix elements (population) and the single-state domains is obvious one-to-one. On the other hand, the mapping between the offdiagonal elements and coherence domains is not so obvious. Indeed for a system of F states, the number of full electronic matrix elements is F2, while the number of (high dimensional) cubes in the active state space is 2F. Apparently associating each off-diagonal element of electronic matrix to a particular cube does not constitute a one-to-one mapping in general. This inconsistency may be solved if we consider a generalized MM Hamiltonian including many-state interactions, i.e., /(x , P, Q, P) = F

+

F

P2 + 2μ

F

F

F

∑ nkVk(Q) + ∑ ∑ nklVkl(Q) k=1

k=1 l≠k

F

∑∑ ∑

nmklVmkl(Q) + ... + nk1k 2...kF Vk1k 2...kF(Q)

k=1 l≠k m≠l≠k

(8)

Here each term in the potential represents for the (coherent) interactions involving the corresponding number of states, say the first order potential term nkVk may be associated with the contribution from the single state k, and the second order term nklVkl corresponds to the contribution from two states k and l, etc., exactly the same way in the binomial expansion. So the total number of potential terms in eq 8 is 2F − 1, and the decomposed active state space can be precisely mapped into the general form of MM Hamiltonian (eq 8), with the special cube centered at the origin to be mapped into the full coherence domain (see below). Practically high order terms are neglected, and eq 8 reduces to the conventional MM Hamiltonian eq 1. The potential terms corresponding to single state and two-state coherence represent the diagonal (F elements) and off-diagonal electronic 1 matrix elements [ 2 F(F − 1) elements considering no distinction between the symmetrical terms, say ρ12 ↔ ρ21], respectively. Note that the number of terms in eq 1 are indeed F(F + 1) 1 CF1 + CF2 = F + 2 F(F − 1) = 2 , where C1F, and C2F are the coefficients of the binomial expansion. The state-space-to-dynamics mapping for high dimensional systems is basically the same as eq 7, except that the coherence domain needs further classifications. Aside from the single-state domains which map to the first order potential terms in eq 8, coherence domains correspond to high order potential terms, which indicates that which states are involved according to the correspondence between the state space point at the center of cubes (h1, h2, ...hk, ...hF) and the coefficients of the potential terms nk1k2...kF. Let us take the three-state model in Figure 1a as an example to illustrate the CC scheme. Assuming that the MM Hamiltonian can be truncated at the second order terms and higher orders terms are negligible, we have the standard MM model. If the trajectory enters the cube centered at (1,0,0), the contribution from the state 1 dominates, and we map this cube to

cij

Re[ρij ] =

W ({hk }, {Nk}, t = 0)Woff ({hk }, {Nk}, t ) |c | α ij

F ∑k = 1 ⟨W ({hk },

{Nk}, t = 0)W ({hk }, {Nk}, t )⟩ (9a) fij

W ({hk }, {Nk}, t = 0)Woff ({hk }, {Nk}, t ) |f | α Im[ρij ] =

ij

F ∑k = 1 ⟨W ({hk },

{Nk}, t = 0)W ({hk }, {Nk}, t )⟩ (9b)

Equations 9a and 9b are the real and imaginary parts of the off-diagonal matrix element ρij. Here W({hk}, {Nk}, t = 0) is the window function evaluated at time t = 0 corresponding to the sampling of initial configurations, and Woff({hk}, {Nk}, t) is the window function of the off-diagonal regions (Doff) evaluated at time t, i.e., ⎧1, if {hk } ∈ Doff Woff ({hk }, {Nk}, t ) = ⎨ ⎩ 0, otherwise ⎪

⎪

D

(9c)


Article


directly from the trajectories in single-state domains based on the classical electronic variable matrix (eq 2). Note that here the magnitudes of the coherence and flux terms (or action variables) are quantized according to eq 9e, exactly the same strategy to quantize electronic population as in the standard SQC approach.21 Scheme D in Figure 2d extends the off-diagonal region into the diagonal regions therefore including contributions from the trajectories not only when they are in the singlestate regions but also when they move outside. The mapping schemes in Figure 2 connect the state space to the electronic density matrix, which has nothing to do with dynamics but only for the evaluation of electronic variables and related properties. However, in some cases, the CC approach,27 which proposes a mapping between the decomposed state space and nuclear dynamics may coincide with a particular state-spaceto-density-matrix mapping (equivalent to Scheme A but with a simpler symmetry). However, any other mapping schemes like those in Figure 2 may be combined with the CC approach or other dynamics methods (such as SQC) to obtain dynamical information on the electronic density matrix. Consistently the SQC may be taken as a special case of the CC treatment, in which the whole state space is mapped to Ehrenfest dynamics. It is worth noting that the normalization factor I in eqs 9a and 9b is the total number of trajectories detected by all singlestate windows, i.e.,

The weighting parameter ⎧ ⎛1 ⎞⎛ 1 ⎞ ⎪ + γ ⎟ , if {hk } ∉ Dll ⎪ ⎜⎝ + γ ⎟⎜ ⎠⎝ 2 ⎠ α=⎨ 2 ⎪ ⎪ γ(1 + γ ) , if {h } ∈ D ⎩ k ll

(9d)

which is determined by the quantization of the action variables, i.e., |cij| = |fij | =

quantization

(ni + γ )(nj + γ ) ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ α

(9e)

Finally the ensemble average (denoted by ⟨...⟩) of the relevant quantities is normalized by the total number of trajectories associated with all possible final states. Besides the CC scheme, there are other possibilities to map the off-diagonal electronic density matrix elements to state space. Figure 2 exhibits four representative mapping schemes

F

I=

∑ ⟨W ({hk}, {Nk}, t = 0)W ({hk}, {Nk}, t )⟩ k=1

(10a)

Alternatively, the window function associated with the offdiagonal elements may be used to obtain the normalization factor, i.e., F

I=

∑ ⟨W ({hk}, {Nk}, t = 0)Woff ({hk}, {Nk}, t )⟩ k=1

(10b)

This choice is definitely reasonable in the mathematical sense. However, physically trajectory outside the single-state domains defined by the corresponding window functions is not real. Here the single-state window functions serve as detector, and the number of trajectories captured by the window function represents for the corresponding state probability. The sum of all trajectories detected by all single-state windows resembles the partition function in statistical mechanics, and the total probability of finding a trajectory associated with any particular state should be unit. Therefore, there is no room for trajectory outside the single-state domains to enter the normalization process, and in principle any dynamical quantities could be obtained from the partition function eq 10a generated by the trajectories associated with single-state domains, i.e., from the active space trajectory ensemble. The generalization of the CC approach to complex systems of multiple states is straightforward. The window functions of single state are high-dimensional cubes centered at the corresponding state point in the state space, e.g. (0, 0, ...hk = 1, ...0) for the state k. We take the window width 2γ = 1 from the quantum number. The active space is a bigger cube centered 1 1 1 1 at 2 , 2 , ...hk = 2 , ... 2 , with a side length to be 4γ. The

Figure 2. Schemes for mapping the full electronic density matrix in the state space for a two-state model. Colored cubes correspond to regions for state 1 (D11, blue), state 2 (D22, golden), and the region for the off-diagonal elements (green). The regions with forward slashes or backslashes correspond to different values of α in eq 9d.

(schemes A−D) for the evaluation of the off-diagonal matrix elements, taking a two-state system as an example. Scheme A in Figure 2a is a mapping with no overlap between regions corresponding to diagonal and off-diagonal matrix elements, which allows for a further one-to-one mapping between the decomposed state space and nuclear dynamics. In fact, scheme A is physically equivalent to the standard CC scheme because trajectories never enter the blank triangle regions in the whole active space in Figure 2a, required by the conservation of total population. Figure 2b shows a mapping scheme B for the full electronic density matrix essentially the same as the one proposed by Miller and co-workers,23 except that here we use a different window width. Scheme C in Figure 2c shows a straightforward way to calculate the off-diagonal matrix elements

(

)

complementary space to the single-state domain in the active space is the coherence domain. The complementary space of the super cube of the active space in the whole space is the void space, and for completeness we map the void space to BO-like dynamics (eq 7). Note that trajectories in the void space do not contribute E


Article

The Journal of Physical Chemistry A to the dynamical computations. Equivalently here we define an ensemble of trajectories in the active state space, i.e., the active space trajectory ensemble. For nonadiabatic dynamics of complex molecular systems such as condensed phase problems, there may have many DOFs or atoms associated with each single state. Obviously the CC approach to control nuclear dynamics based on the state space decomposition does not depend on the nuclear DOFs of each state, and the algorithm would become no more complicated than its two opposite limits, namely, the standard SQC or MSTAI approaches. Another important issue about the condensed phase nonadiabatic dynamics is if the detailed balance is maintained. Because it has been shown22 that the SQC approach satisfies the detailed balance, we may expect that the CC approach holds the detailed balance condition at least approximately. Besides the switch in the CC scheme, nuclear dynamics is either BO-like or Ehrenfest (eqs 6a, 6b and 6c), all of which along with electronic dynamics are Hamilton dynamics. Therefore, the total energy of the system is at least piecewise conserved. The switch may introduce the change in the total energy, which could lead to the breakdown of the detailed balance. However, if the switch occurs outside the strong coherence region, i.e., the interstate couplings Vkl → 0, the change in the total energy is mainly the difference between the average potential weighted by the fractional population nkk, and the single state potential. Moreover, on the ensemble average level, the number of trajectories associated with the single state is equal to the fractional population, so the total energy is approximately conserved upon the switch (see below).

case II:

case III:

case IV:

j

c j2 mjωj

α = 0.09,

β Δ = 5,

ωc = 2.5 Δ

ε = 1,

α = 0.1,

β Δ = 5,

ωc = 2.5 Δ

ε = 1,

α = 0.1,

β Δ = 0.25,

ωc =1 Δ (13d)

2⎤ ⎡ 2 ⎛ Pki cki ⎞ ⎥ 1 2 ⎢ ⎟ Hii = εi + ∑ + mkiωki ⎜Q ki − ⎢ 2 mkiωki2 ⎠ ⎥⎦ ⎝ ki = 1 ⎣ 2mki Nb

F

+

Nb

⎡ P2 kj

∑∑⎢ j ≠ i kj = 1

⎢⎣ 2mkj

+

⎤ 1 mkjωkj2Q kj2 ⎥ , ⎥⎦ 2

Hij = Δij = const. (14)

Here εi, Δij are the electronic energy of the state i and the interstate coupling between states i and j. (Qki, Pki) is the phase point of the kth bath mode associated with the state i, and mki, ωki, cki are the corresponding mass, frequency, and the coupling constant. The bath frequency and the coupling constant are determined again by the discretization of the continuous spectral density of a Debye bath, i.e., ηωω J(ω) = 2 c 2 ω + ωc (15) Here η = 2λ is the system-bath coupling strength, λ is the reorganization energy, and ωc is the characteristic frequency of the bath. Nb = 200 is the number of bath modes associated with each state. The exact quantum results of the spin boson model are obtained from ref 21. The hierarchical equations of motion (HEOM) approach for the reduced density matrix was supposed to provide an accurate description on the time-dependent population and related energy transfer dynamics in its applications. The HEOM results in this work were obtained by using the PHI package34 unless specified otherwise. The CC approach and its limiting cases, namely, the MST-AI and SQC approaches, produce converged results by averaging 96 000 trajectories.

(11)

IV. RESULTS AND DISCUSSION Figure 3 displays the results of the population difference between the two states in the spin boson model, i.e., D(t) = P1(t) − P2(t), obtained by using the CC approach, the SQC, the MST-AI, and the linearized semiclassical initial value representation (LSC-IVR) methods, in comparison with the exact quantum calculations. The overall agreement between the CC and the exact quantum results for all four cases in eq 13 is excellent, especially at high temperature for both the symmetric unbiased (Figure 3a) and asymmetric biased spin boson problem (Figure 3d). The performance of the CC approach is still fairly good even at quite a low temperature, i.e., βΔ = 5, although the

δ(ω − ωj) (12a) 33

via a discretization scheme on the Ohmic bath, i.e., π J(ω) = αωe−ω / ωc 2

ε = 0,

The model photosynthetic system was used by Ishizaki and Fleming to study the excitation energy transfer in the pigment/ protein complex of photosynthesis. The Hamiltonian matrix elements are given by30,31

Here ε and Δ are the electronic energy and the interstate coupling. V0(Q) and V1(Q) are the zero and first-order nuclear potentials, N N 1 given by V0(Q) = ∑ j =b 1 2 mjωj2Q j2 , and V1(Q) = ∑ j =b 1 cjQ j . Here Nb = 100 is the number of bath modes, mj and Qj are the mass and coordinate of the mode j, and ωj and cj are the corresponding frequency and coupling constant determined by the spectra density

∑

ωc = 2.5 Δ

(13c)

H22(Q) = V0(Q) − V1(Q) − ε ;

π 2

β Δ = 0.1,

(13b)

H11(Q) = V0(Q) + V1(Q) + ε ;

J(ω) =

α = 0.09,

(13a)

III. MODELS AND COMPUTATIONAL DETAILS The CC approach was tested on the simple two-state or threestate models in previous work.27 Here we extend its implementation to general condensed phase nonadiabatic dynamics, taking the well-studied spin-boson model21 and the excitation energy transfer model for photosynthesis30−32 as benchmark problems. The considered spin boson model consists of two electronic states coupled with harmonic thermal bath, and the Hamiltonian matrix elements are given by21

H12(Q) = H12(Q) = Δ

ε = 0,

case I:

(12b)

We consider the same four settings of the spin boson model representing dynamics in different regimes, namely, F


Article


Figure 3. Spin boson problems of a two-state model. Population difference D(t) = P1(t) − P2(t) is plotted for four cases in eq 13. Results obtained from the coherence controlled nonadiabatic dynamics approach (CC, black solid lines) are compared with quantum calculations (QM, circles), MST augmented image (AI, dotted lines), SQC approaches (dot-dashed lines), and linearized semiclassical initial value representation (LSC, red solid lines). (a) Case I, (b) case II, (c) case III, and (d) case IV.

magnitude of coherence in the unbiased case (Figure 3b) and the asymptotic limit in the biased problem (Figure 3c) are not reproduced quantitatively well. The SQC approach shows the best overall performance, and in particular it predicts the correct asymptotic limit, i.e., the thermal equilibrium of the relaxation in the asymmetric case III (Figure 3c). However, the CC approach is at least as good as (or even slightly better than) the SQC method in the other three cases (Figure 3a,b,d). To check whether the CC approach satisfies the detailed balance, we examine the long-time behavior of the population difference for the asymmetric spin boson problem, namely case III and IV. Figure 4 indicates that the CC approach predicts the stable asymptotic value of the population difference nearly identical to the MST-AI results, which is very close to the SQC prediction. At high temperature (Figure 4b), the difference between the SQC and the CC results becomes even smaller than the low temperature case (Figure 4a). The excitation energy transfer model is another benchmark system to test the performance of the CC approach. Figure 5 shows the time-dependent population of state 1 in a two-state model for excitation energy transfer with a fast thermal bath, calculated by using the CC approach, the HEOM, MST-AI, and SQC methods. The agreement of the CC results and those from other methods are fairly good. The CC predictions are nearly identical to the MST-AI results for all cases of different reorganization energies. Both of them are almost the same as the SQC result but predict a slightly stronger coherence than HEOM for small reorganization energy (Figure 5a). While the SQC prediction agrees better with HEOM in the intermediate systembath coupling regimes (Figure 5b,c). In the strong system-bath coupling regime (Figure 5d), the CC approach shows the best agreement with the HEOM. Figure 6 presents the results for the same problem as in Figure 5, but for a slower thermal bath, i.e., τc = 500 fs. In this case the electronic system has more time to relax before nuclear

Figure 4. Long-time behavior of the population difference D(t) for the two-state spin boson problems. The results obtained from the coherence controlled nonadiabatic dynamics approach (CC, solid lines) are compared with the MST augmented image (AI, dotted lines) and SQC approaches (dot-dashed lines). Results from the exact quantum calculations are also plotted (QM, circles). (a) Case III and (b) case IV.

motions reach to a new equilibrium. So the differences between the treatments of nuclear dynamics in the CC, SQC, and MST-AI approaches diminish, and all these three approaches produce almost the same results for all four cases in Figure 6a−6d, which again show a little stronger coherence than the HEOM predictions. We attribute this discrepancy to the implicit treatment on nuclear dynamics in the HEOM approach. However, this explanation warrants further investigation. G


Article


Figure 6. Same as Figure 5, but for τc = 500 fs.

Figure 5. Time-dependent population of state 1 for the two-state model photosynthetic system. The results obtained from the coherence controlled nonadiabatic dynamics approach (CC, dot-dashed lines) are compared with the hierarchical equation of motion (HEOM, black solid lines, from ref 30), the MST augmented image (AI, dotted lines) and SQC approaches (red solid lines). The characteristic time of the thermal bath τc = 100 fs, and T = 300 K. The difference in the site energies ε1 − ε2 = 100 cm−1, and J12 = 100 cm−1 is interstate electronic coupling. (a) λ = J12/50; (b) λ = J12/5; (c) λ = J12; and (d) λ = 5J12.

(state populations), but also the off-diagonal matrix elements, because the phase information is stored in the electronic variables Bkl in eq 2, and the only difference between the two methods are the treatment on nuclear dynamics. However, the quality of the prediction depends on how we define the coherence region corresponding to the off-diagonal elements. Figure 8 compares the results of the off-diagonal matrix element ρ21 for the two-state model photosynthetic system in Figure 5b, computed by using the HEOM method and the CC and SQC approaches according to four different mapping schemes for the off-diagonal electronic matrix elements shown in Figure 2. The CC approach produces nearly identical results for the imaginary part of ρ21 with the SQC calculations and a little more positive real part of ρ21 than the latter. Both methods show a fairly good overall agreement with the HEOM results. The CC results obtained from the standard electronic mapping scheme (Figure 8a, equivalent to scheme A) and the SQC derived scheme proposed by Miller and co-workers23 (scheme B, Figure 8b) predict similar short time behavior of both real and imaginary part of ρ21 with the HEOM calculations. Good agreement is also seen for the long time imaginary part of ρ21, while the long time real part of ρ21 obtained by the CC approach is small in magnitude. The CC computations according to the single-state scheme (scheme C, Figure 8c) and the mixed singlestate/coherence scheme (scheme D, Figure 8d) show excellent agreement with the HEOM results on the long-time behaviors of both real and imaginary parts of ρ21, whereas both treatments predict a larger magnitude of ρ21 at short times. The best results (in comparison with the HEOM calculations) for this two state problem seem to be the ones obtained by the standard SQC approach, and the CC approach with the single-state scheme.

We also applied the CC approach to the seven-state model for excitation energy transfer in photosynthesis, and compared the results of time-dependent site populations in Figure 7 with those obtained by using the HEOM, the SQC and MST-AI methods. All theoretical methods predict the relatively longlived coherence (up to several hundred femto-seconds, fs), and the energy transfer from the initial state 1 to other states along a particular pathway. The CC approach resembles the MST-AI treatment, and is in better agreement with the HEOM method than the SQC approach on the computation of the site populations. Specifically in comparison with the CC approach, HEOM predicts a little stronger coherence and less amount of transfer from the initial state to weak coupling states such as states 5, 6, and 7. The SQC approach shows relatively preferential transfer to the strong coupling states, namely states 3 and 4. Since the energy transfer in this multistate system involve different transfer pathways, the quantum coherence between pathways24 could make the overall dynamics quite complicated. Therefore, it is not straightforward to explain why HEOM predicts weaker coherence for the two-state problem but stronger coherence for the seven-state problems than other three methods. It would be interesting to see how these methods agree with realistic experiments if such a comparison were allowable. Similar to the SQC approach,23 the CC method can be used to compute not only the diagonal electronic matrix elements H


Article


Figure 7. Time-dependent site populations of the 7-state model photosynthetic system. The results obtained from (a) the hierarchical equation of motion (HEOM, from ref 31), (b) the coherence controlled nonadiabatic dynamics approach (CC); (c) the SQC approach; and (d) the MST augmented image (AI). The characteristic time of the thermal bath τc = 50 fs, λ = 35 cm−1, and T = 300 K. The site energies and interstate electronic couplings are from ref 31.

Figure 8. Off-diagonal electronic matrix elements in the two-site excitation energy transfer model in Figure 5b. The results are computed by using the HEOM, (solid lines) CC (dashed lines)and SQC (dotted lines) approaches for nuclear dynamics, and different schemes in Figure 2 for the mapping between the electronic matrix elements and the state space. (a) CC scheme (equivalent to scheme A in Figure 2a); (b) the standard SQC scheme (scheme B, mc); (c) single state scheme (scheme C, ss); (d) the mixed single-state/coherence scheme (scheme D, sc). Real part of ρ21: Re; imaginary part of ρ21: Im.

subsystem of the 7-state photosynthetic model system in Figure 7, which includes the first three states only, with all parameters the same. We consider four electronic mapping

To further evaluate the performance of the electronic mapping scheme, we consider the high dimensional problem by taking a three state model as an example. This model system is the I


Article


Figure 9. Off-diagonal electronic matrix elements in the three-state excitation energy transfer model. The results are computed by using the HEOM (solid lines), CC (dashed lines), and SQC (dotted lines) approaches for nuclear dynamics, and different schemes for the mapping between the electronic matrix elements and the state space. (a) CC scheme, (b) scheme B (mc), (c) single state scheme (ss), and (d) modified scheme B (mc2). Real part of ρ12: Re; imaginary part of ρ12: Im.

Figure 10. Same as Figure 9 except that the results are for ρ13.

single-state cubes in the active space; (3) the standard SQC scheme, for which the off-diagonal window function corresponds to the cubes centered at the midpoint of the related two singlestate cubes; for example, the mapping domain for ρ12 is the

schemes: (1) the CC scheme, as shown in Figure 1a, for which the off-diagonal window function in eq 10b corresponds to the off-diagonal cubes in the active space; (2) the single state scheme, for which the off-diagonal window function corresponds to the J


Article


Figure 11. Same as Figure 9 except that the results are for ρ23.

Figure 12. Off-diagonal electronic matrix elements in the 7-state excitation energy transfer model. The results are computed by using the HEOM, (solid lines) and CC approaches for nuclear dynamics, and different schemes for the mapping between the electronic matrix elements and the state space: single state scheme (dashed lines, ss), and CC scheme (dotted lines). (a) ρ12; (b) ρ13; (c) ρ14; (d) ρ23. Real part: Re; imaginary part: Im.

Figure 9 displays the results of ρ12 obtained by using the HEOM, CC, and SQC approaches for the four different electronic mapping schemes. The CC scheme (Figure 9a) reproduces most of the HEOM results with substantial underestimation on the long time magnitude of the real part of

cube centered at (1/2, 1/2, 0) with an side length of 2γ; and (4) the modified standard SQC scheme, for which the off-diagonal window function corresponds to the domain that extends the 2D off-diagonal region (Figure 2b) to the full 3D space. K


Article

The Journal of Physical Chemistry A ρ12, and moderate underestimation on the imaginary part of ρ12 in magnitude. The standard SQC scheme (Figure 9b) does not perform as well as the CC scheme, so is the modified version (Figure 9d), which is slightly better than the standard scheme. The best results are provided by using the CC approach with the single state scheme (Figure 9c). The results of ρ13 are shown in Figure 10. Again the CC scheme (Figure 10a) and the single state scheme (Figure 10c) perform better than other treatments (Figure 10b,d). Figure 11 shows the results of ρ23. The CC scheme (Figure 11a) produces the results in fairly good agreement with the HEOM calculations with overestimation on the magnitude of real part of ρ23 and that of the short time imaginary part. By contrast, the SQC related schemes (Figure 11b,d) underestimate the magnitude of real part of ρ23 and that of the short time imaginary part. The single state scheme shows the best overall agreement with the HEOM results. Figure 12 displays the representative off-diagonal electronic matrix elements of the 7-state model system (see also Figure 7 and Figure S1), calculated by using the CC approach for nuclear dynamics with the single-state and the CC electronic mapping schemes. The results show that the single state scheme predicts the time-dependent electronic coherence in reasonably good agreement with the HEOM method, in support of our previous argument that dynamical information should be able to achieved from the active state space trajectory ensemble. Moreover, the CC electronic mapping scheme performs at a very similar level in accuracy to the single state scheme, indicating that the mapping between the decomposed state space and the electronic density matrix could be useful. This finding is not surprising since in our framework the trajectory travels in the single state domain or in the coherence domain belongs to the same active state space trajectory ensemble. Our idea is different from the trajectory ensemble treatment suggested by Martens,35 in which each density matrix element corresponds to a different trajectory ensemble. The high-dimensional state-space decomposition scheme (see for example Figure 1) can be easily implemented to condensed phase systems regardless of the number of nuclear DOFs involved on any particular state. The unprojected scheme shown in Figure 1b is even simpler than the standard scheme in that the former does not project the full Ehrenfest dynamics into any subspace, therefore all coherence terms are retained in a united coherence domain. Figure 13 compares the results of the site

populations in the seven-state model obtained from the CC approach using the standard and the unprojected schemes. It turns out that the effect of different decomposition schemes is small at least for this case. However, it is suggested to always perform the check on this effect for other problems if the unprojected scheme is used. One advantage of the CC approach is its excellent efficiency and numerically stability. It is well-known that the mean field dynamics in the MM model-based nonadiabatic simulations may cause the divergence of trajectory. For example, considering the three state photoreaction model system in our previous work,27 about 72% trajectories in the standard SQC calculation diverge during the time widowing we considered (see Figure S2). By contrast, the diverged trajectory is only 10−4 in the CC simulation, whereas no trajectory diverges in the AI treatment, as expected. Here we obtained well-converged results by using 96 000 trajectories for the seven-state model system consisting of 1400 bath DOFs. The CC approach and the closely related MST-AI are apparently suitable for simulating nonadiabatic molecular dynamics in condensed phases featuring a favorable scaling factor with the system size, as demonstrated in this work and previous studies on complex molecular systems.22,24−26 In principle, nonadiabatic dynamics should not depend on the chosen representation, either adiabatic or diabatic, which is normally taken as a criterion for evaluating the performance of a simulation method. Although our state space decomposition and the coherence controlled approach are proposed based on the MM model in the diabatic representation, the transformation from the diabatic to adiabatic representation is straightforward with the interstate coupling terms in eq 1 replaced by the nonadiabatic coupling in the adiabatic representation, and the following mappings remain the same. Figure S3 shows the state transmission probabilities for single (Figure S3a) and dual (Figure S3b) avoided crossing problems calculated by the CC approach in both diabatic and adiabatic representations, which exhibit about the same accuracy in comparison with the exact quantum DVR results. Energy conservation is an important issue in dynamics simulations. Figure 14 illustrates the energy fluctuation of a representative CC trajectory for the two-state model system described in Figure 5b. Upon the transition between the single

Figure 14. Energy deviation of a representative trajectory for the twosite excitation energy transfer model in Figure 5b simulated by the CC approach. Plotted are the nuclear-electronic energy of the system when the trajectory in the single state domain (SSD, red), or in the coherence domain (CD, blue), and the total energy of the system (black), which also include the pure electronic energy.

Figure 13. Effect of using different decomposition schemes for coherence domains on site populations of the 7-state model. The CC results obtained from the standard scheme (ex. Figure 1a, colored lines) and the unprojected scheme (ex. Figure 1b, dotted lines). L


Article


in comparison with the HEOM method, which we give our preference in future applications. The applicability of the single state scheme justifies the fact that dynamical information could be obtained from the active state space trajectory ensemble while the success of the CC scheme demonstrates the feasibility of the state space decomposition and the state-space-toelectronic-matrix mapping. The consistent treatment of Born− Oppenheimer-like dynamics in single-state domains and (projected) Ehrenfest dynamics in coherence domains in the CC approach ensures the algorithm to be extremely efficient and numerically stable, as demonstrated by the performance on the benchmark problems with up to seven electronic states and 1400 nuclear DOFs. Therefore, the CC approach will provide a practical way to simulating nonadiabatic molecular dynamics in condensed phase systems.

state domain and the coherence domain, nuclear dynamics is switched between BO-like and Ehrenfest types. The corresponding nuclear-electronic energy appears to be piecewise smooth. However, once taking pure electronic energy into account, the total energy of the system is pretty smooth, and its fluctuation seems to agree with the normal dynamics simulation estimate. For example, the standard deviation of the total energy of the trajectory in Figure 14 is only 1.6 × 10−4 for the time step dt =10 au.

V. CONCLUSIONS The recently proposed coherence controlled (CC) approach to nonadiabatic dynamics was elaborated and implemented to condensed phase problems. The state-space decomposition scheme in general high-dimensional systems was illustrated, and the mapping between the state space and the full electronic matrix and that between the decomposed state space and different treatments of nuclear dynamics were specified by means of a generalized MM Hamiltonian including many-state interactions. Dynamical properties such as the full electronic matrix may be obtained from the active space trajectory ensemble, and the CC dynamics obeys the detailed balance approximately. The CC approach was applied to the benchmark spin boson system and showed a reasonably good agreement with the exact quantum treatment. The agreement between the CC approach and the related SQC and MST-AI methods was also good, and the long-time simulation confirmed that the CC approach satisfied the detailed balance approximately, at least for the current application. The CC approach was also implemented to the excitation energy transfer problems. The agreement between the CC results for a variety of model systems with those obtained using the HEOM, SQC, and MST-AI methods is remarkable. Different electronic mapping schemes for calculating the off-diagonal electronic matrix elements were evaluated, and the results show that both the single state scheme and the CC scheme predict electronic coherence in reasonable accuracy

■ ■

APPENDIX A schematic flowchart of the CC dynamics is provided in Figure 15. ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b10936. Some other off-diagonal electronic matrix elements in the 7-state model system, the population of diverged trajectories in a three-state system, and the results of single and dual avoided crossing models in both diabatic and adiabatic representations (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Guohua Tao: 0000-0002-5374-2513 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS G.T. thanks Professor Akihito Ishizaki for providing HEOM data from refs 30 and 31. We acknowledge the support from the National Science Foundation of China (Grant No. 51471005, and 21673012) and Peking University Shenzhen Graduate School.

■

REFERENCES

(1) Miller, W. H. Electronically Nonadiabatic Dynamics via Semiclassical Initial Value Methods. J. Phys. Chem. A 2009, 113, 1405−1415. (2) Tully, J. C. Perspective: Nonadiabatic Dynamics Theory. J. Chem. Phys. 2012, 137, 22A301. (3) 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. (4) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011−2015. J. Phys. Chem. Lett. 2016, 7, 2100−2112. (5) Ehrenfest, P. Bemerkung Ü ber Die Angenäherte Gültigkeit Der Klassischen Mechanik Innerhalb Der Quantenmechanik. Eur. Phys. J. A 1927, 45, 455−457. (6) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (7) Jasper, A. W.; Nangia, S.; Zhu, C.; Truhlar, D. G. Non-BornOppenheimer Molecular Dynamics. Acc. Chem. Res. 2006, 39, 101−108.

Figure 15. Flowchart of the CC dynamics. M


Article


Transfer: Reduced Hierarchy Equation approach. J. Chem. Phys. 2009, 130, 234111. (31) Ishizaki, A.; Fleming, G. R. Theoretical Examination of Quantum Coherence in a Photosynthetic System at Physiological Temperature. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 17255−17260. (32) Tao, G.; Miller, W. H. Semiclassical Description of Electronic Excitation Population Transfer in a Model Photosynthetic System. J. Phys. Chem. Lett. 2010, 1, 891−894. (33) Wang, H.; Thoss, M.; Sorge, K.; Gelabert, R.; Gimenez, X.; Miller, W. H. Semiclassical Description of Quantum Coherence Effects and their Quenching: A Forward−backward Initial Value Representation Study. J. Chem. Phys. 2001, 114, 2562−2571. (34) Strumpfer, J.; Schulten, K. Open Quantum Dynamics Calculations with the Hierarchy Equations of Motion on Parallel Computers. J. Chem. Theory Comput. 2012, 8, 2808. (35) Martens, C. C. Communication: Fully Coherent Quantum State Hopping. J. Chem. Phys. 2015, 143, 141101.

(8) Wu, Y.; Herman, M. F. Nonadiabatic Surface Hopping HermanKluk Semiclassical Initial Value Representation Method Revisited: Applications to Tully’s Three Model Systems. J. Chem. Phys. 2005, 123, 144106. (9) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-induced Surface Hopping. J. Chem. Phys. 2012, 137, 22A545. (10) Subotnik, J. E.; Shenvi, N. A New Approach to Decoherence and Momentum Rescaling in the Surface Hopping Algorithm. J. Chem. Phys. 2011, 134, 024105. (11) 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. (12) 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. (13) Subotnik, J. E. Fewest-Switches Surface Hopping and Decoherence in Multiple Dimensions. J. Phys. Chem. A 2011, 115, 12083−12096. (14) Shushkov, P.; Li, R.; Tully, J. C. Ring Polymer Molecular Dynamics with Surface Hopping. J. Chem. Phys. 2012, 137, 22A549. (15) Prezhdo, O. V.; Rossky, P. J. Evaluation of Quantum Transition Rates from Quantum-Classical Molecular Dynamics Simulations. J. Chem. Phys. 1997, 107, 5863−5878. (16) Meyer, H. D.; Miller, W. H. A Classical Analog for Electronic Degrees of Freedom in Nonadiabatic Collision Processes. J. Chem. Phys. 1979, 70, 3214−3223. (17) Stock, G.; Thoss, M. Semiclassical Description of Nonadiabatic Quantum Dynamics. Phys. Rev. Lett. 1997, 78, 578−581. (18) Thoss, M.; Stock, G. Mapping Approach to the Semiclassical Description of Nonadiabatic Quantum Dynamics. Phys. Rev. A: At., Mol., Opt. Phys. 1999, 59, 64−79. (19) Miller, W. H. The Semiclassical Initial Value Representation: A Potentially Practical Way for Adding Quantum Effects to Classical Molecular Dynamics Simulations. J. Phys. Chem. A 2001, 105, 2942− 2955. (20) Cotton, S. J.; Miller, W. H. Symmetrical Windowing for Quantum States in Quasi-Classical Trajectory Simulations. J. Phys. Chem. A 2013, 117, 7190−7194. (21) Cotton, S. J.; Miller, W. H. Symmetrical Windowing for Quantum States in Quasi-Classical Trajectory Simulations: Application to Electronically Non-Adiabatic Processes. J. Chem. Phys. 2013, 139, 234112. (22) Miller, W. H.; Cotton, S. J. Communication: Note on Detailed Balance in Symmetrical Quasi-Classical Models for Electronically NonAdiabatic Dynamics. J. Chem. Phys. 2015, 142, 131103. (23) Miller, W. H.; Cotton, S. J. Communication: Wigner Functions in Action-angle Variables, Bohr-Sommerfeld Quantization, the Heisenberg Correspondence Principle, and a Symmetrical Quasiclassical Approach to the Full Electronic Density Matrix. J. Chem. Phys. 2016, 145, 081102. (24) Tao, G. Electronically Non-adiabatic Dynamics in Singlet Fission: a Quasi-Classical Trajectory Simulation. J. Phys. Chem. C 2014, 118, 17299−17305. (25) Tao, G. Bath Effect in Singlet Fission Dynamics. J. Phys. Chem. C 2014, 118, 27258−27264. (26) Tao, G. Understanding Electronically Non-Adiabatic Relaxation Dynamics in Singlet Fission. J. Chem. Theory Comput. 2015, 11, 28−36. (27) Tao, G. Coherence-Controlled Nonadiabatic Dynamics via StateSpace Decomposition: A Consistent Way To Incorporate Ehrenfest and Born-Oppenheimer-Like Treatments of Nuclear Motion. J. Phys. Chem. Lett. 2016, 7, 4335−4339. (28) Prezhdo, O. V.; Rossky, P. J. Evaluation of Quantum Transition Rates from Quantum-classical Molecular Dynamics Simulations. J. Chem. Phys. 1997, 107, 825−834. (29) Hack, M. D.; Truhlar, D. G. A Natural Decay of Mixing Algorithm for Non-Born-Oppenheimer Trajectories. J. Chem. Phys. 2001, 114, 9305−9314. (30) Ishizaki, A.; Fleming, G. R. Unified Treatment of Quantum Coherent and Incoherent Hopping Dynamics in Electronic Energy N


Mapping State Space to Quasiclassical Trajectory Dynamics in

Recommend Documents