Simulating Excited State Dynamics in Systems with Multiple Avoided

Subscriber access provided by UNIV OF LETHBRIDGE

Letter

Simulating Excited State Dynamics in Systems with Multiple Avoided Crossings Using Mapping Variable Ring Polymer Molecular Dynamics Jessica Ryan Duke, and Nandini Ananth J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.5b01957 • Publication Date (Web): 01 Oct 2015 Downloaded from http://pubs.acs.org on October 2, 2015

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

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

Page 1 of 16

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

The Journal of Physical Chemistry Letters

Simulating Excited State Dynamics in Systems with Multiple Avoided Crossings Using Mapping Variable Ring Polymer Molecular Dynamics Jessica R. Duke and Nandini Ananth∗ Department of Chemistry and Chemical Biology, Cornell University, Ithaca, New York, 14853, USA E-mail: [email protected]

∗

To whom correspondence should be addressed

1 ACS Paragon Plus Environment


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

Abstract Mapping Variable Ring Polymer Molecular Dynamics (MV-RPMD) is an approximate quantum dynamics method based on imaginary-time path integrals for simulating electronically nonadiabatic photochemical processes. By employing a mapping protocol to transform from a discrete electronic state basis to continuous Cartesian phase-space variables, the method captures electronic state transitions coupled to nuclear motion using only classical MD trajectories. In this work, we extend the applicability of MV-RPMD to simulations of photo-induced excited electronic state dynamics in nonadiabatic systems with multiple avoided crossings. We achieve this by deriving a new electronic state population estimator in the phase space of electronic variables that is exact at equilibrium and numerically accurate in real time. Further, we introduce an efficient constraint protocol to initialize an MV-RPMD simulation to a particular electronic state. We numerically demonstrate the accuracy of this estimator and constraint technique in describing electronic state dynamics from an initial non-equilibrium state in three model systems. We also present simulations of photo-induced excited state dynamics in three model systems that describe photodissociation.

Graphical TOC Entry


Page 2 of 16

Page 3 of 16

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


Developing computational methods capable of simulating photochemical reactions in large, complex molecular systems remains a significant challenge in theoretical chemistry with important applications to several areas of science, including the rational design of materials such as organic photovoltaics for renewable energy harvesting. 1 Scaling limitations of exact quantum simulations with system size have led to the development of many approximate methods, including mixed quantum-classical approaches like surface hopping, 2–6 quantum-classical Liouville dynamics, 7,8 semiclassical (SC) methods, 9–15 and path integral (PI)-based model dynamics. 16–23 One such PI method, Ring Polymer Molecular Dynamics (RPMD), has emerged as an accurate and efficient technique for the direct dynamic simulation of condensed-phase, thermal charge transfer processes. However, the lack of explicit electronic state variables makes RPMD inapplicable to photochemical processes, and recent efforts to extend its applicability 21,22 fail to preserve the exact quantum Boltzmann distribution (QBD). Mapping Variable (MV)-RPMD, 23 derived by one of us, is an analogue of RPMD that retains all the attractive features of the original RPMD formulation 18,20 — favorable scaling in computational effort with system size, a consistent dynamic framework for electronic state transitions and nuclear dynamics, and dynamics that preserve detailed balance— while explicitly incorporating discrete electronic states. In the method’s inaugural paper, 23 we showed that MV-RPMD trajectories correctly capture nonadiabatic dynamics in two-level systems over a wide range of coupling strengths. However, the method lacked a reliable function that reported electronic state populations in time and that allowed us to constrain a system to a particular initial electronic state, necessary features for the simulation of photo-initiated dynamics. In this work, we derive an exact expression for calculating electronic state populations in the MV-RPMD framework, and we use this expression to develop a simple and efficient protocol to initialize a simulation to a particular electronic state. We then numerically demonstrate the accuracy of MV-RPMD trajectories in capturing electronic state dynamics for a series of model nonadiabatic systems with one or more avoided crossings. This work presents the first direct simulations of photo-induced electronic state dynamics using an imaginary-time PI method. First, we briefly review the MV-RPMD formulation for the general case where K electronic



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

Page 4 of 16

states are coupled to an F -dimensional nuclear coordinate. Throughout this work, we use atomic ˆ = units such that ~ = 1. The Hamiltonian for a K-level system is expressed as H PK

n,m=1

ˆTP ˆ P 2M

+

ˆ ˆ ˆ |ψn iVnm (R)hψ m |, where R and P represent nuclear position and momentum operators,

respectively, M is nuclear mass, {ψn } represent diabatic electronic states, and {Vnm (R)} are diabatic potential energy matrix elements. As detailed in Refs. 23 and 24, PI discretization of the quantum canonical partition function (PF) in continuous electronic and nuclear variables is achieved using Meyer-Miller-Stock-Thoss (MMST) mapping 25,26 to introduce Cartesian electronic variables that represent discrete electronic states. Electronic momentum variables are introduced through Wigner-Weyl transformations, 27–31 and nuclear momenta are introduced through normalized Gaussian momentum integrals. The resulting phase-space PI representation for the PF is ˆ

h

Z = Tr e−β H ∝

lim

N →∞

Z

i

(1)

{dRα }

Z

{dPα }

Z

{dxα }

Z

β

{dpα } e− N HN ({Rα },{Pα },{xα },{pα }) sign(Θ) ,

where β = 1/kB T , T is temperature, and we introduce the notation {dξα } = R

R

dξ1 · · · dξN for R

ξ = {R, P, x, p}. The MV-RPMD Hamiltonian in Eq. 1 is

HN =

N X

"

#

N T N Aα + xα xα + pTα pα − ln|Θ|, β β

α=1

(2)

where N is the number of imaginary time slices or “beads,” PTα Pα MN2 T Aα = (Rα − Rα+1 ) (Rα − Rα+1 ) + , 2β 2 2M

Γ=

K X

" N Y

n=1 α=1

#

1 (xα + ipα ) (xα − ipα ) − I M(Rα ) 2 T

(3)

,

(4)

nn

and Θ = Re [Γ]. In Eq. 4, I is the identity matrix, [·]nn refers to a diagonal element of the enclosed


Page 5 of 16

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


K × K matrix, and   β   e− N Vnn (R)

n=m

Mnm (R) = 

(5)

β   − β Vnm (R) e− N Vnn (R)

n 6= m.

N

The variables xα and pα are continuous position and momentum vectors of length K that represent the K electronic states associated with the αth ring polymer bead, and the nuclear positions and momenta of this bead are represented by the F -dimensional vectors Rα and Pα , respectively. Real-time MV-RPMD thermal correlation functions (CFs) are expressed as 23

CAB (t) ≈

hW sign(Θ)A(z0 )B(zt )i , hW sign(Θ)i

(6)

where angular brackets indicate ensemble averages. Initial distributions (t = 0) are generated using path integral Monte Carlo (PIMC) importance sampling with respect to the weighting function β

W = e− N HN , with HN defined in Eq. (2). As emphasized in Ref. 23, the sign term in Eq. (6) originates from the non-positivity of the Boltzmann distribution for K > 2 level systems, and, in all cases tested thus far, numerical convergence with respect to bead number is achieved well before the average value of the sign function approaches zero. In Eq. (6), z0 refers to the set of RP variables, {Rα , Pα , xα , pα }, at time zero, and zt refers to their values at later time t, which are obtained from z0 by integrating classical equations of motion generated by the MV-RPMD Hamiltonian in Eq. (2). In this work, we introduce a population estimator that is exact for equilibrium and, as we will show numerically, reports accurately on population transfer at short times. In order to derive this estimator, we exploit the property that the trace of a product of two operators can be written as the phaseh

i

ˆ = space integral of the product of their Wigner functions, Tr AˆB where the Wigner function is defined as O(x, p) =

R

1 hK

R

R

dx dp A(x, p)B(x, p),

ˆ + ∆x ieipT ∆x ; 31 a detailed d∆x hx − ∆x |O|x 2 2

derivation is provided in the Supporting Information (SI). The resulting expression for calculating



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

Page 6 of 16

the equilibrium average population for electronic state n is D D

E

ˆn = P

W sign(Θ)P¯n hW sign(Θ)i

E

,

(7)

where N 1 X 1 T x −pT p 2 2 K+1 −xα α α ¯ α Pn = [xα ]n + [pα ]n − 2 e , N α=1 2

(8)

which is an average over the Wigner population functions for individual beads and is hereafter referred to as the Wigner estimator. Not surprisingly, the Wigner estimator is similar to population functions used in semiclassical methods that employ ST mapping, 32,33 but here it has been derived exactly in the MV-RPMD framework. We numerically confirm that the Wigner estimator reproduces thermal average electronic state populations exactly; the details of this calculation are presented in the SI. In addition to a reliable population estimator, simulations of many interesting physical processes, especially photo-initiated ones, require a system to be initialized to a particular electronic state. In general, to simulate electronic state dynamics from an initial non-equilibrium state, |ψi i, we evaluate the quantum correlation function h

D

E

ˆ n (t) = P

ˆ ˆ ˆ −iHt Tr |ψi ihψi | eiHt P n e

Tr [|ψi ihψi |]

i

.

(9)

Using the new Wigner estimator, the MV-RPMD approximation to Eq. (9) becomes *

f ({Rα , Pα }) D

E

ˆ n (t) ≈ P

N Q K Q

+

δ(Pk (α)−ak ) P¯n (t)

α=1 k=1

*

f ({Rα , Pα })

N Q K Q

+

.

(10)

δ(Pk (α)−ak )

α=1 k=1

In Eq. (10), the Boltzmann factor from Eq. (6) has been replaced with separate functions, f and δ, that determine the initial nuclear distribution and electronic state populations, respectively. The initial nuclear distribution is generated using PIMC sampling on the function f , which is system-


Page 7 of 16

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


dependent and will be discussed later in the context of specific model systems. The delta functions in population employ ak values 1 or 0 based on the initial occupation of the k th electronic state and constrain the system to a particular electronic state at t = 0. Numerous approaches to implementing this electronic state constraint can be considered; however, we find only one, which is similar in spirit to the focusing approximation in the context of semiclassical simulations, 34 that is stringent enough to ensure the system starts in the desired electronic state: we set each term in the summation in Eq. (8) to 1 if state n is fully occupied and to 0 if it is unoccupied, solving the single-bead equation

T

T

Pn (α) ≡ 2K+1 e−xα xα −pα pα [xα ]2n + [pα ]2n −

    1

1 = an =  2  

(11)

0

for [xα ]2n + [pα ]2n ≡ r2 . For unoccupied states, the solution is r0 ≈ 0.707. For the fully occupied state, there will be two solutions whose values depend on the total number of electronic states in the system. In practice, we use the larger of these two solutions, henceforth called r1 , because the smaller one proves numerically unstable. Thus, to initialize an MV-RPMD simulation such that electronic state n is occupied, electronic variables [xα ]n and [pα ]n for each bead α are sampled on a circle of radius r1 in phase space, while [xα ]m and [pα ]m for unoccupied states (m 6= n) are sampled on circles of radius r0 . MV-RPMD trajectories released from such a constrained distribution move according to the Hamiltonian in Eq. (2). For all models studied here, we use a 4th order Adams-Bashforth-Moulton predictor-corrector integrator. Because of the presence of β throughout Eq. (2), the MV-RPMD Hamiltonian used for dynamics requires a fictitious temperature. Since temperature provides an indirect measure of the energy of a system, we choose this temperature to match the total initial energy, E, of the constrained system using the relationship kB T = E. Bead convergence for all calculations presented here is achieved by converging the equilibrium average energy calculated at this fictitious temperature using PIMC. Finally, real-time electronic state population dynamics are



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

Page 8 of 16

obtained from an ensemble average of the Wigner population estimator, P¯n (t), over MV-RPMD trajectories. We first demonstrate the use of the Wigner estimator and initialization technique using three one-dimensional, two-level models first described in Ref. 3. The functional forms of the potential curves and off-diagonal couplings are provided in the SI and plotted in the left panel of Fig. 1. For each of these models (Models I-III), the exact quantum calculation employs an initial Gaussian 2 /σ 2

wavepacket, ψ(R) = eiP0 R e−(R−R0 )

, for a particle of mass M = 2000 a.u. with a forward

momentum distribution centered at P0 and a nuclear position distribution centered at R0 with width σ =

20 , P0

occupying the lower-energy diabatic state, V11 . 3 To ensure our results are not fortuitously

good, we test three different values of P0 for each model (panels A, B, and C in Fig. 1); the kinetic energies associated with these values are indicated with arrows in the left panel of Fig. 1, and P0 and R0 are reported in Table 1. Table 1: Parameters for Two-State Models (a.u.) Model I

II

III

P0 6.3 8.9 10 10 14.1 16 1.6 2.2 4

R0 -5 -5 -5 -7 -7 -7 -7 -7 -7

N 3 3 3 3 3 3 8 8 6

T (K) 3133 6253 7894 7894 15695 20210 202 382 1263

Trajectories 7200 2400 1920 2400 1920 1920 38,400 21,600 19,200

We ensure that the nuclear phase space distribution for each bead in MV-RPMD matches the corresponding quantum wavepacket probability distribution using the sampling function

f ({Rα , Pα }) =

N Y

2 /σ 2

e(Rα −R0 )

eσ

2 (P

α −P0 )

2

.

(12)

α=1

The system is constrained to electronic state 1 by sampling electronic variables on circles of radius r1 ≈ 1.398 (the solution to Eq. (11) for occupied states when K = 2) and r0 ≈ 0.707. The


Page 9 of 16

temperature used for dynamics is calculated from the initial energy; in these systems, this is simply the kinetic energy,

P02 , 2M

since the lower potential surface is shifted to zero energy. Temperatures

used for the different models are reported in Table1, along with the number of beads and trajectories required to converge our simulations to second decimal place. For comparison, exact quantum dynamics results are obtained using the Discrete Variable Representation (DVR) 35 grid method. For Models I and II, we use a converged grid range of R = − 15 to 15 a.u. and a grid spacing of 0.05 a.u.; Model III requires a grid range of R = − 30 to 30 a.u. and a grid spacing of 0.05 a.u. As shown in Fig. 1, the MV-RPMD results match the exact results well, particularly at short times, where we expect the approximate dynamics to be most accurate. We note that our MV-RPMD results compare favorably with surface hopping and surface hopping with decoherence. 36 0.025

A B

0.010

C

A

B

C

1.0 0.5 0.0

-0.005 A B

0.050

1.0

C

0.000

0.5

-0.050

0.0

0.004

Populations

Potential Energy (a.u.)

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


1.0

A

0.002

0.5 B C

0.000 -10

0.0 -5

0

5

10 0

1000

2000

0

Nuclear Coordinate (a.u.)

1000

2000

0

1000

2000

3000

Time (a.u.)

Figure 1: Diabatic potential energy functions and real-time electronic state populations for Models I, II, and III (top, middle, and bottom rows, respectively). The first column shows potential curves with solid colored lines and coupling functions with black dashed lines (the coupling function for Model III is scaled by 0.02). Kinetic energies associated with the three values of P0 tested for each model are indicated with blue arrows and are labeled to correspond to results in panels A, B, and C. MV-RPMD results in each panel are shown with solid lines, exact quantum results are shown with dashed lines, and lines are colored to match the diabats in the first column. We now turn to a set of model three-state, single-mode systems parametrized to describe photodissocation processes. 33 The physical picture in each case (Models IV-VI) involves an initial, 9 ACS Paragon Plus Environment


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

Page 10 of 16

occupied, excited electronic state generated by optical excitation from a harmonic ground state using a femtosecond laser pulse. Assuming Franck-Condon excitation, the initial non-equilibrium state, |ψi i, is a ground state harmonic oscillator wavefunction centered at R0 with population of 1 in excited electronic state |1i. After photo-excitation, the sub-100 fs dynamics of the electronic state populations are assumed to evolve on a manifold of three coupled excited state surfaces, allowing us to ignore the ground state contribution to these dynamics. The diagonal elements of the diabatic potential energy matrix for the excited state manifold are Morse potentials of the form Vi = Dei (1 − e−βi (R−Rei ) )2 + ci , and the off-diagonal couplings, 2

Vij = Aij e−aij (R−Rij ) , are Gaussian functions centered at the crossing points Rij . Parameters for these functions are reported in Table S1 in the SI, and the potential curves and couplings are plotted in the left panel of Fig. 2. For all three models, the nuclear mass M is 20,000 a.u., and the frequency ω of the harmonic ground state potential is 0.005 a.u. The initial nuclear distribution for the MV-RPMD simulation is generated by PIMC sampling on the ground state Hamiltonian for a harmonic oscillator using the sampling function f ({Rα , Pα }) = β0

0

e N HN , where

HN0

=

N X α=1

"

#

Pα2 MN2 1 2 2 2 + , 2 (Rα − Rα+1 ) + M ω (Rα − R0 ) 0 2M 2 2β

(13)

and β 0 indicates that sampling is performed at 300 K. Electronic variables are, as before, sampled on circles of radius r1 ≈ 1.559 (the solution to Eq. (11) for occupied states when K = 3) and r0 ≈ 0.707. The fictitious temperature used for dynamics is calculated from the initial energy, which for these models is the zero point energy (ZPE) of the ground state plus the potential energy gap between the lowest excited state and the initial occupied excited state at R0 . Parameters for Models IV-VI are reported in Table 2. Exact quantum dynamics results are obtained using the DVR grid method with an initial state |ψi i = e−

Mω (R−R0 )2 2

|1ih1|. We use a converged grid range of

R = 0 to 20 a.u. and a grid spacing of 0.05 a.u. We also compare our results to those obtained using the Linearized SC Initial Value Representation (LSC-IVR), 32,33 a method that also employs


Page 11 of 16

Table 2: Parameters for Three-State Models (a.u.) R0 2.1 3.3 2.9

Model IV V VI

N 3 3 3

T (K) 15288 9605 8843

Trajectories 1200 900 720

a consistent dynamic framework for nuclear and electronic state dynamics and neglects nuclear quantum coherence effects. Our implementation follows Ref. 33 and employs an initial nuclear coherent state of width γ = M ω. Numerical convergence is achieved with 144,000 trajectories for Model IV and 28,800 trajectories for Models V and VI. MV-RPMD, DVR, and LSC-IVR results are shown in the middle and right columns of Fig. 2. 0.06

1.0

0.04 0.5 0.02 0.0 1.0

0.00 0.04

0.5 0.02

Populations

Potential Energy (a.u.)

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


0.0 1.0

0.00 0.04

0.5 0.02 0.00

0.0 2

4

6

8

10

12 0

1000

Nuclear Coordinate (a.u.)

2000

0

1000

2000

3000

Time (a.u.)

Figure 2: Diabatic potential energy curves and real-time electronic state populations for Models IV, V, and VI (top, middle, and bottom rows, respectively). The left column shows potential curves with solid lines and coupling functions with dashed lines; the center of the initial nuclear wavefunction constrained to excited state 1 is indicated by a black arrow. The middle column shows MV-RPMD results, while the right column shows LSC-IVR results, both with solid lines; exact quantum results are shown with dashed lines. Lines are colored to match the diabats in the left column. For all three models, MV-RPMD simulation results match the exact results at short times very well. Model IV involves an avoided crossing at longer times and a likely contribution from nuclear 11 ACS Paragon Plus Environment


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

Page 12 of 16

quantum coherence effects, 33 making it the most challenging for our model dynamics to accurately describe. Models V and VI are better suited to the MV-RPMD approach, with both avoided crossings being accessible at shorter times, as evidenced by the significantly improved agreement with exact quantum results. Fig. 2 also shows that MV-RPMD performs significantly better than LSC-IVR for all three models. In summary, we derive and test a new electronic state population estimator in the MV-RPMD framework, and we introduce an efficient sampling scheme to generate MV-RPMD distributions corresponding to an initial photoexcited electronic state. The numerical accuracy of the MV-RPMD results for all six models makes it clear that the Wigner estimator introduced here is an accurate electronic state population estimator and, in combination with the initialization protocol, can accurately describe photo-induced electronic state dynamics. Future applications of the methodology developed here will include investigating exciton dynamics in organic photovoltaics. We also anticipate that the Wigner estimator will prove useful in developing an MV-RPMD reaction rate theory.

Acknowledgement We sincerely thank Gregory S. Ezra and Thomas F. Miller III for many helpful discussions. N. A. acknowledges funding support through a startup grant from Cornell University, and J. R. D. acknowledges a National Science Foundation Graduate Research Fellowship under Grant No. DGE1144153.

Supporting Information Available Derivation of the Wigner population estimator, details of equilibrium electronic state population calculations using this estimator, potential energy functions for Models I-III, potential energy parameters for Models IV-VI, and MV-RPMD results for N = 1.. This material is available free of charge via the Internet at http://pubs.acs.org/.


Page 13 of 16

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


References (1) Smith, M. B.; Michl, J. Recent Advances in Singlet Fission. Annu. Rev. Phys. Chem. 2013, 64, 361-386. (2) Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H+ with D2 . J. Chem. Phys. 1971, 55, 562-572. (3) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 10611071. (4) Kuntz, P. J. Classical Path Surface-Hopping Dynamics. I. General Theory and Illustrative Trajectories. J. Chem. Phys. 1991, 95, 141-155. (5) Barbatti, M. Nonadiabatic Dynamics with Trajectory Surface Hopping Method. WIREs Comput. Mol. Sci. 2011, 1, 620-633. (6) Shushkov, P.; Li, R.; Tully, J. C. Ring Polymer Molecular Dynamics with Surface Hopping. J. Chem. Phys. 2012, 137, 22A549/1-22A549/13. (7) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919-8929. (8) Kapral, R. Progress in the Theory of Mixed Quantum-Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129-157. (9) Sun, X.; Miller, W. H. Semiclassical Initial Value Representation for Electronically Nonadiabatic Molecular Dynamics. J. Chem. Phys. 1997, 106, 6346-6353. (10) Bonella, S.; Coker, D. F. LAND-Map, A Linearized Approach to Nonadiabatic Dynamics Using the Mapping Formalism. J. Chem. Phys. 2005, 122, 194102/1-194102/13. (11) Ananth, N.; Venkataraman, C.; Miller, W. H. Semiclassical Description of Electronically



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

Page 14 of 16

Nonadiabatic Dynamics via the Initial Value Representation. J. Chem. Phys. 2007, 127, 084114/1-084114/9. (12) Miller, W. H. Electronically Nonadiabatic Dynamics via Semiclassical Initial Value Methods. J. Phys. Chem. A 2009, 113, 1405-1415. (13) Huo, P.; Coker, D. F. Iterative Linearized Density Matrix Propagation for Modeling Coherent Excitation Energy Transfer in Photosynthetic Light Harvesting Systems. J. Chem. Phys. 2010, 133, 184108/1-184108/12. (14) Cotton, S. J.; Miller, W. H. Symmetrical Windowing for Quantum States in Quasi-Classical Trajectory Simulations. J. Phys. Chem. A 2013, 117, 7190-7194. (15) Antipov, S. V.; Ye, Z.; Ananth, N. Dynamically Consistent Method for Mixed QuantumClassical Simulations: A Semiclassical Approach. J. Chem. Phys. 2015, 142, 184102/1184102/9. (16) Voth, G. A. Path-Integral Centroid Methods in Quantum Statistical Mechanics and Dynamics. Adv. Chem. Phys. 1996, 93, 135-218. (17) Jang, S.; Voth, G. A. A Derivation of Centroid Molecular Dynamics and Other Approximate Time Evolution Methods for Path Integral Centroid Variables. J. Chem. Phys. 1999, 111, 2371-2384. (18) Craig, I. R.; Manolopoulos, D. E. Quantum Statistics and Classical Mechanics: Real Time Correlation Functions from Ring Polymer Molecular Dynamics. J. Chem. Phys. 2004, 121, 3368-3373. (19) Huo, P.; Coker, D. F. Semi-Classical Path Integral Non-Adiabatic Dynamics: A Partial Linearized Classical Mapping Hamiltonian Approach. Mol. Phys. 2012, 110, 1035-1052. (20) Habershon, S.; Manolopoulos, D. E.; Markland, T. E.; Miller III, T. F. Ring-Polymer Molecular


Page 15 of 16

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


Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space. Annu. Rev. Phys. Chem. 2013, 64, 387-413. (21) Richardson, J. O.; Thoss, M. Nonadiabatic Ring-Polymer Molecular Dynamics. J. Chem. Phys. 2013, 139, 031102/1-031102/4. (22) Menzeleev, A. R.; Bell, F.; Miller III, T. F. Kinetically Constrained Ring-Polymer Molecular Dynamics for Non-Adiabatic Chemical Reactions. J. Chem. Phys. 2014, 140, 064103/1064103/17. (23) Ananth, N. Mapping Variable Ring Polymer Molecular Dynamics: A Path-Integral Based Method for Nonadiabatic Processes. J. Chem. Phys. 2013, 139, 124102/1-124102/8. (24) Ananth, N.; Miller III, T. F. Exact Quantum Statistics for Electronically Nonadiabatic Systems Using Continuous Path Variables. J. Chem. Phys. 2010, 133, 234103/1-234103/9. (25) 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. (26) Stock, G.; Thoss, M. Semiclassical Description of Nonadiabatic Quantum Dynamics. Phys. Rev. Lett. 1997, 78, 578-581. (27) Weyl, H. Quantum Mechanics and Group Theory. Z. Phys. 1927, 46, 1-46. (28) Wigner. E. On the Quantum Correction For Thermodynamic Equilibrium. Phys. Rev. 1932, 40, 749-759. (29) Wigner. E. On the Penetration of Potential Energy Barriers in Chemical Reactions. Z. Phys. Chem. Abt. B 1932, 19, 203-216. (30) Moyal, J. E. Quantum Mechanics as a Statistical Theory. Math. Proc. Cambridge 1949, 45, 99-124.



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

Page 16 of 16

(31) Case, W. B. Wigner Functions and Weyl Transforms for Pedestrians. Am. J. Phys. 2008, 76, 937-946. (32) Sun, X.; Wang, H.; Miller, W. H. Semiclassical Theory of Electronically Nonadiabatic Dynamics: Results of a Linearized Approximation to the Initial Value Representation. J. Chem. Phys. 1998, 109, 7064-7074. (33) Coronado, E. A.; Xing, J.; Miller, W. H. Ultrafast Non-Adiabatic Dynamics of Systems with Multiple Surface Crossings: A Test of the Meyer-Miller Hamiltonian with Semiclassical Initial Value Representation Methods. Chem. Phys. Lett. 2001, 349, 521-529. (34) Bonella, S.; Coker, D. F. Semiclassical Implementation of the Mapping Hamiltonian Approach for Nonadiabatic Dynamics Using Focused Initial Distribution Sampling. J. Chem. Phys. 2003, 118, 4370-4385. (35) Colbert, D. T.; Miller, W. H. A Novel Discrete Variable Representation for Quantum Mechanical Reactive Scattering via the S-Matrix Kohn Method. J. Chem. Phys. 1992, 96, 1982-1991. (36) Shenvi, N.; Subotnik, J. E.; Yang, W. Simultaneous-Trajectory Surface Hopping: A ParameterFree Algorithm for Implementing Decoherence in Nonadiabatic Dynamics. J. Chem. Phys. 2011, 134, 144102/1-144102/12.


Simulating Excited State Dynamics in Systems with Multiple Avoided

Recommend Documents