Nonadiabatic Dynamics in Atomistic Environments: Harnessing

Subscriber access provided by UNIV OF LETHBRIDGE

Letter

Nonadiabatic Dynamics in Atomistic Environments: Harnessing Quantum-Classical Theory with Generalized Quantum Master Equations William C. Pfalzgraff, Aaron Thomas Kelly, and Thomas E. Markland J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.5b02131 • Publication Date (Web): 12 Nov 2015 Downloaded from http://pubs.acs.org on November 17, 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 23

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

Nonadiabatic Dynamics in Atomistic Environments: Harnessing Quantum-Classical Theory with Generalized Quantum Master Equations William C. Pfalzgraff, Aaron Kelly, and Thomas E. Markland∗ Department of Chemistry, Stanford University, Stanford, California 94305, 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 ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23

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


A theoretical description of processes such as electron and proton coupled electron transfer in enzymes and organic electronics, and exciton transport in light harvesting complexes and photovoltaic materials, requires the development of methods that are able to accurately treat the nonadiabatic quantum dynamics of large condensed phase systems. In many of these problems one is particularly interested in how atomistic details, such as the positioning or interactions of particular chemical groups in the environment, couple to the quantum subsystem and give rise to the resulting dynamics. A commonly employed approach has been to map the fully atomistic problem onto a simpler form for the Hamiltonian, usually by representing the environment as a set of effective harmonic degrees of freedom which are bilinearly coupled to the quantum subsystem. 1,2 Once simplified in this way, the mapped problem can usually be solved exactly. 3 However, by performing this mapping one loses the ability to treat anharmonic environments and nonlinear coupling effects in a self-consistent manner. 4 The main alternative to this approach is to simulate the fully atomistic system using approximate dynamics methods. A number of methods have been previously applied to study the dynamics of quantum subsystems embedded in atomistic environments. Probably the most commonly employed of these is the Fewest Switches Surface Hopping (FSSH) approach, 5,6 which converges with relatively few trajectories and usually approximately satisfies detailed balance, allowing relatively accurate long time populations to be obtained in many cases. However, FSSH is known to suffer from the so-called decoherence problem, and predicts an incorrect scaling of the electron transfer rate as a function of electronic coupling in the Marcus regime. 7 Another promising class of approaches to treat these problems are wavepacket based methods, such as ab initio multiple spawning. 8 While formally exact in principle, obtaining converged results can be prohibitive in problems where the number of possible quantum transitions leads to a large number of spawning events. Recently, another class of approaches has been suggested based on extending the ring polymer molecular dynamics ansatz to formulate ensemble conserving nonadiabatic equations

3



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

of motion. While RPMD itself has been applied to electron dynamics in atomistic systems 9 this direct application has been shown to fail to capture the correct turnover in the electron transfer rate in the Marcus inverted regime. 9 This has led to the recent development of a number of methods in the spirit of RPMD whose dynamics can be shown to conserve 10,11 and approximately conserve 12 the exact quantum Boltzmann distribution as well as those that exactly conserve an approximate distribution. 13 These methods have shown considerable promise in a series of model problems, but how they will perform for atomistic systems currently remains unclear. The hierarchy of approximate quantum-classical (or semiclassical) approaches that can be derived from the Quantum-Classical Liouville Equation (QCLE), 14 the semiclassical initial value representation, 15 or from a direct linearization of the path integral, 16–20 also provide a promising route to performing nonadiabatic quantum dynamics in atomistic environments. In principle, this hierarchy allows one to choose a balance between accuracy and efficiency for a given application by providing a selection of methods that offer an increasing amount of dynamical correlation between the subsystem and environment. However, in practice the higher tiers of this hierarchy require too many trajectories to obtain converged results to be practical for large systems, while the lower tiers suffer from severe detailed balance problems, resulting in incorrect population distributions at long times. As such, currently only fully linearized methods (such as LSC-IVR, 16 LAND-map, 19 and PBME 20 ), and recently introduced partially linearized variants (such as PLDM 21 and FBTS 22 ), are practical for combination with a fully atomistic treatment of the environment. Recent work has shown that one can leverage considerable increases in both the accuracy and efficiency of quantum dynamics approaches by combining quantum-classical approaches with the quantum master equation formalism. 23–28 Using this approach, we have shown that the regime of applicability of highly accurate solutions of the QCLE can be vastly extended, 29 and that lower tier quantum-classical trajectory based methods can be made both more accurate and more efficient when applied to nonequilibrium relaxation processes. 30

4


Page 4 of 23

Page 5 of 23

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


In particular, by combining Ehrenfest mean field theory 31 with the generalized quantum master equation (GQME) formalism, yielding the MF-GQME approach, we demonstrated that one can obtain a quantitative description of the subsystem dynamics throughout a very wide range of parameter space in the spin-boson model, ranging from electron and proton transfer to exciton and energy transfer regimes. Despite its apparent simplicity, there are very few methods that are able to converge accurate results for the spin-boson model for such a wide range of regimes. Although the considerable benefits obtained by combining semiclassical methods with the GQME framework have been extensively demonstrated for the spin-boson problem, it is of paramount importance to assess whether these benefits carry over to atomistic systems that lie beyond linear system-bath coupling and harmonic environments. In this paper, we apply our MF-GQME approach to an atomistic model of charge transfer in water, and demonstrate that it is able to obtain quantitative accuracy over a wide range of charge transfer regimes beyond the confines of harmonic environments and bilinear coupling. These results demonstrate that the accuracy of MF-GQME exhibited in charge transfer regimes of the spin boson problem is equally evident in atomistic environments. This indicates that wide regime of accuracy that MF-GQME exhibited in the spin boson problem will carry over to fully atomistic systems, offering the opportunity to quantitatively describe quantum relaxation processes where the atomistic details of the environment play a vital role. The model, recently introduced by Schwerdtfeger et. al, consists of a donor-acceptor complex that can undergo electronic charge transfer, explicitly solvated by water molecules. 32 The Hamiltonian for the system is given by: ˆ =H ˆs + H ˆb + H ˆ sb . H

(1)

ˆ s = ǫˆ H σz + ∆ˆ σx

(2)

The subsystem Hamiltonian is,

5



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 23

where σ ˆx and σ ˆz are Pauli spin matrices, 2ǫ is the energetic bias between the two states and ∆ is the (constant) electronic coupling between the states. There are two subsystem diabatic states: |1i corresponds to the charge localized on the donor and |2i corresponds to the charge localized on the acceptor. Hence when ǫ > 0, localization of the charge on the donor is higher in energy than doing so on the acceptor. The bath Hamiltonian is, ˆ2 ˆ + Vˆ2 (R) ˆ ˆ b = P + 1 Vˆ1 (R) H 2M 2

(3)

ˆ and R ˆ are the multidimensional bath momentum and position operators respecwhere P ˆ is the potential energy of the bath when the subsystem is localized in tively, and Vˆi (R) diabatic state |ii, obtained from the molecular dynamics force field. 32 The system-bath coupling Hamiltonian is, ˆ − Vˆ1 (R) ˆ σ ˆ sb = − 1 Vˆ2 (R) ˆz . H 2

(4)

Note that due to the dependence of Hsb on the atomistic potential energies of the two bath ˆ and Vˆ2 (R), ˆ the system-bath coupling is not bilinear. states, Vˆ1 (R) The population dynamics of the subsystem reduced density matrix, ρˆs (t), can be obtained by integrating the Nakajima-Zwanzig GQME: 23,24 Z t d ρˆs (t) = −iLs ρˆs (t) + I(t) − dτ K(τ )ˆ ρs (t − τ ), dt 0

(5)

where Ls is the subsystem Liouville operator which generates the free subsystem evolution. The effect of the bath on the subsystem dynamics is encoded in the inhomogeneous term, I(t), and the memory kernel, K. For the initial conditions used here the inhomogeneous term vanishes, (see the Supporting Information) and hence the effect of the bath arises solely from the memory kernel. Although computation of the exact memory kernel would allow one to obtain exact subsystem dynamics, this is intractable for arbitrary atomistic condensed phase systems. How-

6


Page 7 of 23

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


ever, we have recently demonstrated that using an approximate description of the dynamics, Ehrenfest mean field theory, to calculate the kernel allows one to obtain high accuracy results for a wide range of regimes of the spin boson problem. Our MF-GQME method exploits the realization that approximating the memory kernel, which typically decays rapidly in nonadiabatic regimes, and using it to integrate the reduced dynamics often provides a more accurate result than using the same approximation directly. 30 The memory kernel can be constructed from unprojected dynamics by exploiting the Dyson decomposition of the propagator, 25 which allows one to express the kernel in terms of system-bath time correlation functions calculated from dynamics in the full space. The implementation of this method for a general Hamiltonian is described in detail in Refs. ( 25,29,30,33 ), and the expressions for the correlation functions needed to compute the memory kernel for this model of atomistic charge transfer are given in the Supporting Information (SI). In this system the bath degrees of freedom that couple to the charge transfer process are low frequency and hence, unlike our earlier work, 29,30 the bath was sampled from the classical rather than Wigner density using the expressions given in the SI. Since exact solutions for the nonadiabatic quantum dynamics of large fully atomistic systems, such as this one, remain elusive, no exact benchmark to assess the accuracy of our MF-GQME approach exists in most regimes. However, due to its high temperature bath and low to modest electronic coupling, this model can be recognized to be within the regime in which Rips-Jortner (RJ) 34 theory is accurate. This theory thus provides a benchmark to assess the performance of our MF-GQME approach. However, it is important to note that the advantage of the MF-GQME approach is that it can be used to calculate results outside the parameter regimes in which RJ theory is accurate. The RJ expression for the rate of charge transfer for this model is, 4π |∆|2 τL k = 1+ h ¯λ "

#−1

(−2ǫ+λ)2 2π |∆|2 − √ e 4λkB T h ¯ 4πλkB T

7


!

(6)


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

where λ is the reorganization energy, h ¯ is the reduced Planck constant, kB is the Boltzmann constant, T is the temperature, and τL is the solvent longitudinal relaxation time. At low electronic couplings (∆ → 0) the charge transfer rate is limited by the electronic curve crossing between the diabatic states, and hence the first term in this expression approaches unity yielding simply the Marcus rate. 35 However, at higher electronic couplings, where the charge transfer becomes facile, the rate is limited by the solvent relaxation time which is accounted for by the first term. The reorganization energy, λ, and the solvent relaxation time, τL , can be obtained in an atomistic simulation from the fluctuations in the energy gap between the two diabatic states. The reorganization energy was calculated by directly computing the free energy curves as a function of diabatic energy gap in the reactant and product states, 36 and fitting the resulting curves to the standard parabolic form, giving λ = 20.1 kcal/mol. The relaxation time, τL can be obtained by fitting the autocorrelation function of the diabatic energy gap, from which we obtained τL = 0.155 ps. Figure 1 shows the electronic population relaxation in the donor-acceptor charge transfer complex following initialization of the subsystem in the excited state donor diabat, in a regime of large energetic bias and weak electronic coupling where RJ theory is expected to be highly accurate. This system has been previously studied using FSSH, 32 where the population decay (shown in Fig. 1) was found to be over an order of magnitude faster than the RJ rate for this electronic coupling strength. In contrast, the direct Ehrenfest population dynamics in Fig. 1 show a similar relaxation rate to RJ theory, but give long time populations where both states are close to being equally occupied. However, due to the large energetic bias in favor of the acceptor state, the final population of the donor state should be almost zero (e−2ǫ/kB T = 2 × 10−5 ). Our MF-GQME approach, which combines Ehrenfest mean field theory with the GQME framework, is in near perfect agreement with with the RJ result, obtaining both the correct rate and long time population. Figure 2 shows how the charge transfer rate varies as a function of the electronic coupling,

8


Page 8 of 23

Page 9 of 23





Page 10 of 23

Page 11 of 23

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


∆, where the first point corresponds to the population dynamics shown in Fig. 1. MF-GQME again gives results in excellent agreement with the RJ rate across the full range of electronic couplings. In particular, MF-GQME closely tracks the RJ rate at higher electronic couplings, where the system adiabaticity increases and the bath dynamics begin to become rate limiting. This leads to a contribution from the frictional effects of the solvent via the first term in Eq. 6, which corrects the Marcus rate to account for these effects. Direct Ehrenfest theory again yields excellent agreement with the rate, although as in Fig. 1, the relaxation is to incorrect long time populations. This is consistent with previous studies on simple model problems of electron transfer 37 and exciton Hamiltonians, 38 where Ehrenfest rates have been observed to be accurate even when long time populations are given poorly. Hence, treating these processes using direct Ehrenfest theory should be approached with extreme caution if one is interested in where the charge actually ends up. In addition, at the highest electronic couplings Ehrenfest theory is observed to overpredict the RJ rate, failing to capture the full effect of solvent frictional effects as the coupling increases. Fig. 2 also shows the results obtained by Schwerdtfeger et. al 32 using FSSH. 6 The FSSH rate scales linearly with the electronic coupling for this problem, which is consistent with that observed for a two state donor-acceptor complex coupled to a one dimensional Langevin damped harmonic environment, 7 in contrast the the correct quadratic dependence predicted by Marcus theory at low electronic couplings. In this system, these different scaling behaviors lead to rates that are one to two orders of magnitude larger than RJ, consistent with previous observations for FSSH. 39 The most severe divergence in the rates occurs at low electronic coupling, with FSSH also seemingly failing to capture the curvature in the regime where the RJ correction to Marcus theory becomes significant. Decoherence corrections to FSSH might be expected to improve the rates obtained, although simple collapse and Augmented-FSSH have previously been shown to offer little improvement for these regimes of this model. 32 Figure 3 shows the dependence of the rate on the energetic driving force, ǫ, leading to the Marcus turnover in the rate as the driving force is increased. Consistent with what

11




Page 12 of 23

Page 13 of 23




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

number of force evaluations required to converge the property of interest, such as the population dynamics. In MF-GQME, the long-time population dynamics is generated by numerically integrating Eq. 5 using a memory kernel that is truncated at some finite lifetime. Hence if the memory kernel decays rapidly relative to the time required for the full population decay which one wishes to simulate, large speed-ups can be obtained. In practice we have shown MF-GQME to be orders of magnitude more efficient than even direct Ehrenfest for a wide variety of spin boson regimes, with the speed-up increasing as the system becomes more nonadiabatic. Figure 4 shows the charge transfer rate obtained as a function of the memory kernel cut-off time. In all cases, by 75 fs the rate is unchanged as more memory is included in the kernel. However, when the memory kernel is more aggressively truncated, larger errors occur for the higher electronic couplings. This is consistent with our previous observation that the memory kernels decay much faster as the nonadiabaticity of the system is increased. 30 For the lowest electronic coupling the ratio of the memory kernel decay time of ∼ 75 fs, from which the full population decay can be generated by integrating Eq. 5, and the population decay time of 1.5 ns, is thus ∼ 500, 000 fold while for the higher coupling this is reduced to ∼ 400 fold. Figure 5 shows the rate obtained from the integrated dynamics with a memory kernel cutoff off time of 75 fs as the number of trajectories used is increased. For the lowest electronic coupling, even with 500 trajectories the rate is within a factor of 2 of the RJ result and is well converged with around 3000 trajectories (each of length 75 fs). In contrast, for the higher electronic couplings at least 4000 trajectories must be averaged over to reduce the statistical noise in the kernel elements to a small enough level to allow them to yield stable population dynamics using Eq. 5. In addition, as the system becomes more adiabatic (larger ∆) it is clear that more trajectories are required to converge to an accurate rate, with the highest coupling needing at least 20,000 trajectories. With this information in hand, one is able to assess the relative cost, in terms of num-

14


Page 14 of 23

Page 15 of 23




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

ber of force evaluations needed, to generate the full population decay via the traditional direct Ehrenfest approach and our MF-GQME approach. To generate a a single Ehrenfest trajectory for the full decay time at the lowest electronic coupling studied here, ∆ = 0.01 kcal/mol, requires about 12 billion force evaluations (1.5 ns with a 0.25 fs time-step with force evaluations on both electronic surfaces at each time-step), and is typically well converged in 2000 trajectories. In comparison, MF-GQME requires 300 force evaluations per 75 fs trajectory and, as shown in Fig. 5, the memory kernel can be integrated accurately with around 3,000 trajectories and hence only 450,000 force evaluations are needed to obtain a converged rate i.e. significantly less than that for even a single direct Ehrenfest trajectory. Thus, for the lowest electronic coupling, the saving from using MF-GQME is around 2,700 times. For the highest electronic coupling investigated, ∆ = 0.5 kcal/mol, the increased number of trajectories required to produce a converged kernel combined with the larger ratio of the memory kernel time to the population decay time (∼ 0.05), this saving is reduced to ∼ 3 times. In addition, this speed-up is accompanied by the significant improvement in the accuracy of the rate (vs. FSSH) and long time populations (vs. direct Ehrenfest) shown in Fig. 1. The ability to simultaneously achieve more accurate and more efficient results is thus one of the hallmarks of this approach as compared with a direct application of quantum dynamics methods. Finally, we note that the need to perform a large number of very short trajectories in MF-GQME offers considerable computational benefits over a direct Ehrenfest calculation, both in the ability to trivially parallelize over the trajectories, and by reducing the loss due to failure of any single trajectory. In summary, we have shown that MF-GQME provides an accurate and efficient method to treat nonequilibrium charge transfer in anharmonic atomistic environments. The success of this method follows from the well established observation that approximating the fast decaying memory kernel in a master equation is often more accurate and efficient than simulating the dynamics directly, an idea that has been exploited historically in areas such as the development of theories of liquid relaxation and dynamics. 40 This development provides

16


Page 16 of 23

Page 17 of 23

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


new opportunities to obtain quantitative accuracy for nonadiabatic processes that occur in atomistic environments in regimes outside established perturbative limits. Furthermore, this work suggests that other trajectory based methods, such as semiclassics, AIMS or FSSH, could also reap massive benefits in both efficiency and accuracy by combination with the generalized quantum master equation framework, particularly in nonadiabatic regimes.

Simulation methods Simulations of the atomistic model were performed using a custom interface to the LAMMPS molecular dynamics package. 41 The donor-acceptor complex is represented by two spherical ions held fixed at a distance of 0.5 nm, solvated by 433 water molecules. Force field parameters used for the ions and water are provided in the Supporting Information. Periodic boundary conditions were applied using the minimum image convention in a cubic simulation box with sides of length 2.3795 nm. Coulomb interactions were evaluated using particle-particle, particle-mesh Ewald summation 43 with a real-space cutoff of 1.2 nm. The equations of motion were integrated using a 0.25 fs time-step. After 10 ns of equilibration, initial configurations of the isolated bath were generated every 2.0 ps from a 160 ns NVT trajectory at 298K using the stochastic velocity rescaling thermostat 44 with a time constant of 400 fs. These initial configurations were then used to perform Ehrenfest trajectories to obtain the time correlation functions needed to construct the memory kernel. Equation 5 was then numerically integrated using the 2nd order Runge-Kutta method to obtain the subsystem population dynamics. The rate of charge transfer was extracted from an exponential fit to the population dynamics, after an initial transient time. 45

Acknowledgement The authors thank Tim Berkelbach for helpful comments and a thorough reading of this manuscript. This material is based upon work supported by the U.S. Department of Energy,

17



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

Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0014437. T.E.M also acknowledges support from a Cottrell Scholarship from the Research Corporation for Science Advancement, a Terman fellowship, an Alfred P. Sloan Research fellowship and a Hellman Faculty Scholar Fund fellowship. W.C.P acknowledges support from the Melvin and Joan Lane Stanford Graduate Fellowship. A.K acknowledges a postdoctoral fellowship from the Stanford Center for Molecular Analysis and Design.

Supporting Information Available Additional details about the inhomogeneous term of the GQME for this model, expressions for the correlation functions needed to compute the memory kernel, and molecular dynamics simulation parameters.

This material is available free of charge via the Internet at http:

//pubs.acs.org/.

References (1) Rivera, E.; Montemayor, D.; Masia, M.; Coker, D. F. Influence of Site-Dependent Pigment-Protein Interactions on Excitation Energy Transfer in Photosynthetic Light Harvesting. J. Phys. Chem. B. 2013, 117, 5510-5521. (2) Chandrasekaran, S.; Aghtar, M.; Valleau, S.; Aspuru-Guzik, A.; Kleinekath¨ofer, U. Influence of Force Fields and Quantum Chemistry Approach on Spectral Densities of BChl a in Solution and in FMO Proteins. J. Phys. Chem. B. 2015, 119, 9995-10004. (3) Tamura, H.; Burghardt, I.; Tsukuda, M. Exciton Dissociation at Thiophene/Fullerene Interfaces: The Electronic Structures and Quantum Dynamics. J. Phys. Chem. C. 2011, 115, 10205-10210. (4) Gottwald, F.; Ivanov, S.; K¨ uhn, O. Applicability of the Caldeira-Leggett Model to Vibrational Spectroscopy in Solution. J. Phys. Chem. Lett. 2015, 6, 2722-2727.

18


Page 18 of 23

Page 19 of 23

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


(5) 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. (6) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061-1071. (7) Landry, B. R.; Subotnik, J. E. Communication: Standard Surface Hopping Predicts Incorrect Scaling for Marcus Golden-Rule Rate: The Decoherence Problem Cannot Be Ignored. J. Chem. Phys. 2011 135, 191101. (8) Ben-Nun, M.; Mart´ınez, T. J. Nonadiabatic Molecular Dynamics: Validation of the Multiple Spawning Method for a Multidimensional Problem. J. Chem. Phys. 1998 108, 7244-7257. (9) Menzeleev, A.R. Ananth, N. Miller III, T.F. Direct Simulation of Electron Transfer using Ring Polymer Molecular Dynamics: Comparison with Semiclassical Instanton Theory and Exact Quantum Methods . J. Chem. Phys. 2011, 135, 074106. (10) Ananth, N. Mapping Variable Ring Polymer Molecular Dynamics: A Path-Integral Based Method for Nonadiabatic Processes. J. Chem. Phys. 2013, 139, 124102. (11) Duke, J.R. Ananth, N. Simulating Excited State Dynamics in Systems with Multiple Avoided Crossings Using Mapping Variable Ring Polymer Molecular Dynamics. J. Phys. Chem. Lett. 2015, 6, 4219-4223. (12) Richardson, J. O.; Thoss, M. Communication: Nonadiabatic Ring-Polymer Molecular Dynamics. J. Chem. Phys. 2013, 139, 031102 (13) 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, 074106.

19



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

(14) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999 110, 8919-8929. (15) Miller, W. H. Classical S Matrix: Numerical Application to Inelastic Collisions. J. Chem. Phys. 1970, 53, 3578-3587. (16) Sun, X.; Wang, H. B.; 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. (17) Shi, Q.; Geva, E. A Relationship Between Semiclassical and Centroid Correlation Functions. J. Chem. Phys. 2003, 118, 8173-8184. (18) Poulsen, J. A.; Nyman, G.; Rossky, P. J. Practical Evaluation of Condensed Phase Quantum Correlation Functions: A Feynman-Kleinert Variational Linearized Path Integral Method. J. Chem. Phys. 2003, 119, 12179-12193. (19) Bonella, S.; Coker D. F. LAND-map, a Linearized Approach to Nonadiabatic Dynamics Using the Mapping Formalism. J. Chem. Phys. 2005, 122, 194102. (20) Kim; H.; Nassimi, A.; Kapral., R. Quantum-Classical Liouville Dynamics in the Mapping Basis. J. Chem. Phys. 2008, 129, 084102. (21) Huo, P.; Coker, D. F. Communication: Partial Linearized Density Matrix Dynamics for Dissipative, Non-Adiabatic Quantum Evolution. J. Chem. Phys. 2011, 135, 201101. (22) Hsieh C.-Y.; Kapral., R. Nonadiabatic Dynamics in Open Quantum-Classical Systems: Forward-Backward Trajectory Solution. J. Chem. Phys. 2012, 137, 22A507. (23) Nakajima, S. On Quantum Theory of Transport Phenomena: Steady Diffusion. Prog. Theor. Phys. 1958, 20, 948-959. (24) Zwanzig, R. Ensemble Method in the Theory of Irreversibility. J. Chem. Phys. 1960, 33, 1338-1341. 20


Page 20 of 23

Page 21 of 23

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


(25) Shi, Q.; Geva, E. A New Approach to Calculating the Memory Kernel of the Generalized Quantum Master Equation for an Arbitrary System-Bath Coupling. J. Chem. Phys. 2003, 119, 12063-12076. (26) Shi, Q.; Geva,E. A Semiclassical Generalized Quantum Master Equation for an Arbitrary System-Bath Coupling. J. Chem. Phys. 2004, 120, 10647-10658. (27) Berkelbach, T. C.; Reichman, D. R.; Markland, T. E. Reduced Density Matrix Hybrid Approach: An Efficient and Accurate Method for Adiabatic and Non-Adiabatic Quantum Dynamics.. J. Chem. Phys. 2012, 136, 034113. (28) Berkelbach, T. C.; Markland, T. E.; Reichman, D. R. Reduced Density Matrix Hybrid Approach: Application to Electronic Energy Transfer. J. Chem. Phys. 2012, 136, 084104. (29) Kelly, A.; Markland, T. E. Efficient and Accurate Surface Hopping for Long Time Nonadiabatic Quantum Dynamics. J. Chem. Phys. 2013, 139, 014014. (30) Kelly, A.; Brackbill, N. J.; Markland, T. E. Accurate Nonadiabatic Quantum Dynamics on the Cheap: Making the Most of Mean Field Theory with Master Equations. J. Chem. Phys. 2015, 142, 094110. (31) McLachlan, A. D. A Variational Solution of the Time-Dependent Schrodinger Equation. Mol. Phys. 1964, 8, 39-44. (32) Schwerdtfeger, C. A.; Soudackov, A. V.; Hammes-Schiffer, S. Nonadiabatic Dynamics of Electron Transfer in Solution: Explicit and Implicit Solvent Treatments that Include Multiple Relaxation Time Scales. J. Chem. Phys. 2012, 140, 034113. (33) Zhang, M.-L.; Ka, B. J.; Geva, E. Nonequilibrium Quantum Dynamics in the Condensed Phase via the Generalized Quantum Master Equation. J. Chem. Phys. 2006, 125, 044106. 21



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

(34) Rips, I.; Jorter, J. Dynamic Solvent Effects on OuterSphere Electron Transfer. J. Chem. Phys. 1987, 87, 2090-2104. (35) Marcus, R.A. Chemical and Electrochemical Electron-Transfer Theory. Annu. Rev. Phys. Chem. 1964, 15, 155-196. (36) Blumberger, J.; Sprik, M. Redox free-energies from vertical energy gaps: ab initio MD implementation, in Lecture Notes in Physics. Ferrario, M.; Ciccotti, G.; Binder, K, eds., Springer: Berlin Heidelberg 2006; vol. 704; pp. 481-506 (37) Wang, H. B. ; Thoss, M. ; Miller, W. H. Systematic Convergence in the Dynamical Hybrid Approach for Complex Systems: A Numerically Exact Methodology. J. Chem. Phys. 2001, 115, 2979-2990. (38) Ishizaki, A.; Fleming, G. R. On the Interpretation of Quantum Coherent Beats Observed in Two-Dimensional Electronic Spectra of Photosynthetic Light Harvesting Complexes. J. Phys. Chem. B. 2011, 115, 62276233 . (39) Jaeger, H. M.; Green, J. R.; Prezhdo, O. V. Decoherence Reduces Thermal Energy Loss in Graphene Quantum Dots. Appl. Phys. Lett. 2013, 103, 073111. (40) Zwanzig, R. Nonequilibrium Statistical Mechanics Oxford University Press: New York; 2001. (41) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 1995, 117, 1-19. (42) Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D; Impey, R.W.; Klein, M.L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. (43) Hockney, R.W.; Eastwood, J.W Computer Simulation Using Particles, Adam Hilger: New York 1989 22


Page 22 of 23

Page 23 of 23

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


(44) Bussi, G.; Donadio, D.; Parinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (45) M¨ uhlbacher, L.; Egger, R. Electron Transfer Rates for Asymmetric Reactions. Chem. Phys. 2004, 296, 193-199.

23


Nonadiabatic Dynamics in Atomistic Environments: Harnessing

Recommend Documents