Nonadiabatic Quantum Molecular Dynamics in ... - ACS Publications

Jul 25, 2019 - approximation to the multiple cloning method. MCDMS avoids both the ... electronic state basis at every step, either in the adiabatic o...
0 downloads 0 Views 500KB Size
Subscriber access provided by GUILFORD COLLEGE

Spectroscopy and Photochemistry; General Theory

Nonadiabatic Quantum Molecular Dynamics in Dense Manifolds of Electronic States Dmitry A. Fedorov, and Benjamin G. Levine J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.9b01902 • Publication Date (Web): 25 Jul 2019 Downloaded from pubs.acs.org on July 28, 2019

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 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 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.

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 20 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 Quantum Molecular Dynamics in Dense Manifolds of  Electronic States  Dmitry A. Fedorov, Benjamin G. Levine* Department of Chemistry, Michigan State University, East Lansing, Michigan 48824, United States *E-mail: [email protected]

  Abstract  Most nonadiabatic molecular dynamics methods require the determination of a basis of adiabatic or diabatic electronic states at every time step, but in dense manifolds of electronic states such approaches become intractable. A notable exception is Ehrenfest molecular dynamics, which can be implemented without explicit determination of such a basis, but suffers from unphysical behavior when propagation on a mean-field potential energy surface (PES) does not accurately reflect the true dynamics on multiple electronic states. Here we introduce the multiple cloning for dense manifolds of states (MCDMS) method, a systematically improvable approximation to the multiple cloning method. MCDMS avoids both the mean-field PES problem and the need to compute the full electronic spectrum. This is achieved by reformulating multiple cloning to use a subspace of approximate eigenstates constructed from the time-dependent Ehrenfest electronic wave function. By application to model systems, we show that this approach allows a substantial reduction in the size of the required electronic basis without significant loss in accuracy.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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

TOC Graphic 

With forces computed via on-the-fly electronic structure calculations, AIMD enables the exploration of chemical behavior without prior knowledge of the relevant reaction coordinates or mechanisms.1-5 Ab initio nonadiabatic molecular dynamics methods based on trajectory surface hopping6-8 or Gaussian trajectory basis functions4, 9-11 are able to accurately simulate dynamics on multiple electronic states. However, these methods require determination of the electronic state basis at every step, either in the adiabatic or diabatic representation (hereafter “electronic eigenspectrum”). When a very large number of electronic states are potentially populated by the dynamics, explicit calculation of the electronic eigenspectrum at every step renders these methods intractable. There is tremendous interest in problems that fall in this category, however, including semiconductor nanomaterials, plasmonic nanomaterials, molecules in strong laser fields, chemical reactions involving free electrons, molecules on metal surfaces, and radiation damage to molecules and materials. For several of these problems, progress has been made towards the accurate simulation of nonadiabatic molecular dynamics,12-27 but modeling nonadiabatic dynamics in dense manifolds of electronic states remains an area where generally applicable methodological developments are needed.

2 ACS Paragon Plus Environment

Page 2 of 20

Page 3 of 20 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

Methods based on Ehrenfest (mean field) trajectories offer an appealing alternative because they require propagation of a single, time-dependent electronic wave function rather than computation of the full electronic eigenspectrum.14, 28-33 Ehrenfest trajectories provide an accurate description of dynamics in regions of strong coupling, but the mean-field nature of the method results in well-known problems describing dynamics in the regions where full quantum mechanics predicts that the molecular wave function should split into multiple wave packets.34 In such situations, an Ehrenfest trajectory will follow a mean-field potential energy surface (PES), the dynamics on which cannot accurately reflect distinct dynamics occurring on different electronic states. There are numerous approaches to incorporate decoherence into mixed quantum-classical simulations,35-50 but even those based on Ehrenfest trajectories require determination of the electronic eigenspectrum at every step, placing limits on the problems that can be practically solved. In the present work we develop an approach for modeling nonadiabatic dynamics on many electronic states, accurately describing the splitting of the population between states without ever requiring calculation of the full electronic eigenspectrum. Our approach is a variant of the multiple cloning (MC) method.51-55 MC is a clever strategy that allows the accurate splitting of population between states by combining the multiconfigurational Ehrenfest approach56-58 with a basis-set expansion strategy adapted from ab initio multiple spawning.4, 9-10 Here we modify the MC method to eliminate the need to compute the electronic eigenspectrum. The resulting method, multiple cloning for dense manifolds of states (MCDMS), is a systematically improvable approximation to MC dynamics. We will demonstrate the effectiveness of MCDMS by application to several model systems, comparing the results to those of MC simulations performed using the full electronic eigenspectrum.

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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 20

A more detailed description of MC can be found in the literature.51, 55 Here we present only the essential equations with our modifications.

The MC total molecular wave function is

represented as a superposition of vibronic product states, N (t )

  r, R, t    Cn  t  n  r; R, t   n  R; R n  t  , Pn  t  ,  n  t  ,  n  .

(1)

n

Here the nuclear basis functions are represented as  n and the time-dependent Ehrenfest electronic wave functions  n . Nuclear basis functions are frozen multidimensional Gaussians,59 hereafter called trajectory basis functions (TBFs), that can be expressed as products of one-dimensional Gaussians,

 n  R; R n  t  , Pn  t  ,  n  t  ,  n   ei

n

 t t

3N

   R ; R  t  , P  t  ,    ,   1

n

n

n

n

(2)

where 1



n

2  2 n  4     exp   n R  Rn  t   iP n  t  R  Rn  t     









(3)

and Rn and Pn are the phase space center of the nth basis function. The number of nuclear degrees of freedom is N and  n is the time-independent Gaussian width along degree of freedom ρ. The phase space center of each TBF is propagated according to classical equations of motion on Ehrenfest PESs, En  R n  t     n Hˆ el  n

Rn t 

,

(4)

with the force,

Fn  R n  t   

d  n Hˆ el  n dR

Rn t 

.

4 ACS Paragon Plus Environment

(5)

Page 5 of 20 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

Throughout this work the classical equations are integrated using the velocity-Verlet algorithm with a 0.1 au time step. The nuclear mass is that of a proton (1822.0 au). A width of  n = 6.0 is used. The time-dependent quantum amplitudes in eq. (1), C, are obtained by integration of the molecular TDSE,





  i H  iS C SC

(6)

where overlap matrix elements are approximated as Sij   i  j  i  R i (t )   j  R j (t ) 

(7)

and the right-acting time-derivative matrix elements as     Sij   i  j  i  R i (t )   j  R j (t )    i  j  i  R i (t )   j  R j (t )  . t t

(8)

Following past work on MC, off-diagonal potential matrix elements are computed in the bra-ket averaged Taylor expansion (BAT) approximation,51

Vij 



1  i  R i (t )  H el  j  R i (t )    i  R j (t )  H el  j  R j (t )  2



i

j .

(9)

The electronic wave function associated with each TBF is expanded in a basis,

 n  r; R, t    aI  t  I  r; R  .

(10)

I

The electronic time-dependent Schrodinger equation (TDSE) associated with each TBF, ia  H eleca ,

(11)

is integrated in the symplectic split operator approximation,60 with a time step of 0.001 au. Here Helec is the electronic Hamiltonian, a is the vector of electronic wave function amplitudes, and a  a / t . (See Supporting Information for algorithmic details.)

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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 20

The key advance of this work is the development of a strategy for avoiding propagation on unphysical mean-field PESs that does not require explicit knowledge of the electronic eigenspectrum. Each TBF is propagated using the standard Ehrenfest approach until it gets to a region where the mean-field nature of the Ehrenfest trajectory prevents the population from correctly branching between states. In standard MC, this is avoided through adaptive expansion of the set of TBFs (“cloning”). The criterion for triggering a cloning event requires knowledge of the electronic eigenspectrum. To avoid this bottleneck, we require a procedure to approximate the adiabatic eigenstates. To this end, at every nuclear time step, the real and imaginary parts of the electronic

wave

function

from

the

Nstep

previous

electronic

time

steps,

a R ,t , a I ,t , a R ,t t , a I ,t t ,..., a R ,t  N 1t , a I ,t  N 1 t , are saved as columns of a matrix, X. This basis is

analogous to the 2Nstep-order Krylov subspace of the electronic Hamiltonian and a, and therefore should be a good basis in which to approximate the eigenstates that project most strongly onto a. A QR decomposition of X is performed to obtain an orthonormal basis and associated transformation matrix, Q, X  QR

(12)

Matrix Q is used to transform the original electronic wave function a QT a  a k .

(13)

QT H elec Q  H k .

(14)

and Hamiltonian Helec,

A set of 2Nstep approximate eigenstates is determined by diagonalization of the Hamiltonian Hk. If 2Nstep is increased to the total number of electronic states of the system, the approximate eigenvectors converge systematically to the real eigenvectors of the Hamiltonian Helec. In this limit, MCDMS becomes full MC. 6 ACS Paragon Plus Environment

Page 7 of 20 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

For each approximate eigenstate, I, we calculate the absolute difference between the energy of that approximate eigenstate, EIn  R n  t   , and the Ehrenfest energy, En  R n  t   , at each time step. If this difference is larger than a user-defined threshold,  clone , EIn  R n  t    En  R n  t     clone ,

(15)

the cloning procedure is triggered. Note that, though we choose an energy-based criterion for decoherence here, a more sophisticated force-based criterion (as used in past work on MC and many variants of surface hopping35, 49, 61) could be inserted without significant changes to the broader algorithm. Cloning replaces the initial (“parent”) TBF with two TBFs (“clones”) that propagate independently in nuclear phase space, thus allowing the proper splitting of population between electronic states. The first clone TBF is constructed to project onto electronic state I,  'n   n'

aIn n I , aIn

(16)

and the electronic wave function of the second clone TBF is that of the parent wave function after projecting out state I,

 ''n   n''

 J I

1 1 a

n 2 I

aJn Jn .

(17)

Here Jn represents the electronic wave function corresponding to approximate eigenstate J associated with TBF n. To preserve the nuclear norm and the total electronic wave function, the nuclear amplitudes are set according to

Cn'  Cn aIn and 7 ACS Paragon Plus Environment

(18)

The Journal of Physical Chemistry Letters 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

2

Cn''  Cn 1  aIn .

Page 8 of 20

(19)

Additional details of the cloning procedure are presented in Supporting Information. Our ultimate goal is an AIMD scheme, but in this proof-of-concept study the electronic wave function is not explicitly computed. Instead, we will work with model potentials in which diabaticity of the electronic basis is assumed. We choose the model systems to be comprised of a single linear diabatic PES crossing one or more bands containing several parallel linear PESs (see, for example, the colored curves in Figure 1, top panel). The system contains a single nuclear degree of freedom (x). Each model contains a single negatively sloped diabatic basis state, with energy, H11   w1 x .

(20)

Here w1 is a positive number. Hereafter we will refer to this state as diabat 1. The remaining diabatic basis states are positively sloped with slope w2. There are Nupper of these states which are evenly spaced in energy, H ii  w2 x  (i  1) for 1  i  Nupper  1 .

(21)

These comprise what we will call the upper band. In several of the models, a second (lower) band of Nlower states exists, shifted by an additional energy, ε (see for example, the colored curves in Figure 2, top panel), H ii  w2 x  (i  1)   for Nupper  1  i  Nupper  N lower  1 .

(22)

The coupling between diabat 1 and all other states is constant, except where noted below (Model 3), H1i  H i1  k for 1  i . The coupling is zero between all other pairs of states, 8 ACS Paragon Plus Environment

(23)

Page 9 of 20 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

H ij  0 for i  1 and j  1 .

(24)

All models studied in this work have w1=0.25 au, w2=0.025 au, and k=0.005 au. The values of the other parameters are introduced with each model below. In all simulations reported below, population is initiated on the highest adiabatic electronic state (corresponding to diabat 1), with initial position and momentum of -1.0 au and 10.0 au, respectively. All figures in this work contain two sets of lines: colored lines representing the adiabatic state energies as a function of time following the primary (initial) TBF, and black lines reporting the Ehrenfest energy, En  R n  t   , of each TBF in the simulation. The total probability of populating a particular band of states is reported as a percentage. These percentages are computed as the incoherent averages of the projections of the individual TBFs onto the exact electronic eigenstates. The probability of populating the ground electronic state will be referred to as the transmission probability, because ground state population has transmitted diabatically through all intersections in the system. Results labeled MC indicate multiple cloning simulations using the exact electronic eigenspectrum, while simulations labeled MCDMS(M) indicate MCDMS calculations using M approximate eigenstates (M = 2Nsteps). The time domain of all figures is narrowed to only show the crossing region. The first model we study (model 1) contains a single band of 8 states (𝛿 0.01 au, Nupper = 8, and Nlower = 0). Results for three simulations of this system are presented in Figure 1: a traditional Ehrenfest simulation, MC, and MCDMS(2) in the top, middle, and bottom panels, respectively. Very similar transmission probabilities are predicted in these simulations: 78%, 76%, and 77%, respectively. Thus, the MC and MCDMS simulations both faithfully reproduce the Ehrenfest transmission probability. Not shown are MCDMS(4), MCDMS(6), and MCDMS(8) results (see Figures S1-S3); in all of these cases the predicted transmission probability falls between 76% and 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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

79%. Importantly, numerically exact quantum dynamical calculations predict a very similar transmission probability: 78%. A comparison of the final nuclear densities on diabat 1 for model 1 can be found in Figure S4. Agreement is quite good. The densities drawn from MC and MCDMS(N) simulations are nearly identical, regardless of N. Not surprisingly, exact quantum dynamics predicts a wider distribution of positions than the MC and MCDMS simulations, which are based on a frozen-width Gaussian basis. It may be interesting to investigate methods based on a more flexible Gaussian basis for efficient description of wave packet spreading.62 As can be seen in the black lines, after the crossing region the Ehrenfest trajectory propagates on an unphysical mean-field PES. That is, after leaving the crossing region it propagates along a fictional PES corresponding to a linear combination of eigenstates. Both MC and MCDMS(2) simulations exhibit the proper bifurcation of population between adiabatic states, with one TBF following the ground electronic states, and varying numbers of TBFs following states in the band. Note that more TBF states are created in the band when the exact eigenstates are used (8 TBFs) compared to MCDMS(2) (1 TBF). There is a tradeoff here that is inherent to MCDMS. There is a reduction of cost associated with the fact that fewer TBFs are created in the MCDMS(2) simulation, but detailed information about dynamics within the band is not available in this case.

10 ACS Paragon Plus Environment

Page 10 of 20

Page 11 of 20 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

Figure 1. Energies of the exact adiabatic eigenstates (colored lines) and Ehrenfest energies (black lines) for three simulations of Model 1. Simulations are performed via traditional Ehrenfest (top), MC (middle), and MCDMS(2) (bottom). All colored lines correspond to the initial TBF. For MC and MCDMS(2) each black line corresponds to a different TBF.

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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

Figure 2. Energies of the exact adiabatic eigenstates (colored lines) and Ehrenfest energies (black lines) for two simulations of Model 2. Simulations are performed via MC (top) and MCDMS(2) (bottom). All colored lines correspond to the initial TBF. Model 2 is identical to Model 1, except that we introduce a large 0.08 au gap in the middle of the 8-state band, splitting it into two 4-state bands (𝛿 0.01 au, Nupper = 4, Nlower = 4, and ε = 0.08 au). We constructed this model to test how many approximate eigenstates are needed for accurate description of a system where we expect the molecular wavefunction to split into more than two separate wave packets. Similar to model system 1, we vary the number of approximate eigenstates from 8 to 2, but now we look at the population distribution between three possible outcomes: diabat 1, the upper band, and the lower band. The MC and MCDMS(2) cases are shown in Figure 2, top and bottom panels, respectively, while the other MCDMS simulations are

12 ACS Paragon Plus Environment

Page 12 of 20

Page 13 of 20 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

presented in Figure S5-S7. The transmission probabilities predicted for all of these simulations fall within the range 76-82%, while the probabilities of the populations residing in the upper or lower band at the end of the simulations all fall in the ranges 12-14% (upper band) and 6-10% (lower band). It may appear surprising that even with two approximate eigenstates MCDMS is able to describe the splitting of the population between three different outcomes in nearly quantitative agreement to MC simulations with exact eigenstates. This is because interactions with two 4-state bands are separated in time, and in each nonadiabatic event we have significant population transfer only between the initial diabat and a single 4-state band. Thus, at any given time, only two approximate eigenstates are needed to describe the population dynamics. In simulations of systems with many electronic states, it is often true that some pairs of states are not coupled due, for example, to spatial separation or symmetry. Here we illustrate another important advantage of MCDMS, which is that the states that are not populated during the dynamics are automatically excluded from the calculation, eliminating the associated computational expense. Toward this goal we introduce Model 3, a variant of Model 2 in which the coupling between the three lowest states in the lower 4-state band and diabat 1 are set to zero. Figure 3 shows MC and MCDMS(6) results in the top and bottom panels respectively. The branching ratios in these two simulations are identical. This is because the three uncoupled states are never computed in the MCDMS(6) simulation, while the remaining six eigenstates are computed exactly, yielding identical results to full MC.

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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

Figure 3. Energies of the adiabatic eigenstates (colored lines) and Ehrenfest energies (black lines) for two simulations of Model 3. Simulations are performed via MC (top) and MCDMS(6) (bottom). All colored lines correspond to the initial TBF. The exact adiabatic eigenstates are shown for MC, while only the six approximate adiabatic eigenstates are shown for MCDMS(6). Finally, we investigate whether we can still predict accurate branching ratios with a modest number of approximate eigenstates as the density of states is increased. For this purpose, model 4 is defined to have a band of 16 states with a state density four times that of model 1 (𝛿 0.0025 au, Nupper = 16, and Nlower = 0). Results for model 4 computed via MC and MCDMS(2) are shown in Figure 4, top and bottom panels, respectively. The transmission probability remains in nearquantitative agreement, changing from 59% for MC to 61% for MCDMS(2). As in the case of model 1, agreement with the transmission probability predicted by exact quantum dynamical simulations—60%—is excellent. 14 ACS Paragon Plus Environment

Page 14 of 20

Page 15 of 20 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

Figure 4. Energies of the exact adiabatic eigenstates (colored lines) and Ehrenfest energies (black lines) for two simulations of Model 4. Simulations are performed via MC (top) and MCDMS(2) (bottom). All colored lines correspond to the initial TBF. To summarize, we introduced MCDMS, a new method for studying nonadiabatic dynamics in systems with dense manifolds of electronic states. This variant of multiple cloning utilizes an approximation to the electronic eigenstates derived from the time-dependent Ehrenfest electronic wave function to avoid the need to determine the full electronic eigenspectrum. In contrast to traditional Ehrenfest dynamics, long-time propagation on unphysical mean-field PESs is completely avoided, as demonstrated by application to several model problems above. Though

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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 present method is a fully quantum mechanical treatment (based on frozen Gaussian basis functions), the approximate eigenstate strategy demonstrated here can easily be incorporated into mixed quantum-classical treatments. This will be the topic of a future paper. Integration of MCDMS with our GPU-accelerated implementation of TD-CI63 is also underway. This advance will enable the accurate treatment of nonadiabatic molecular dynamics in dense manifolds of electronic states in an AIMD context.

Computational Methods  All calculations were performed with a development version of pySpawn. PySpawn is a modern software implementation of the full multiple spawning method. The procedure used for exact quantum dynamical calculations is described in Supporting Information.

Associated Content  Supporting Information Available: Supplementary figures for additional calculations mentioned

in the text, details of the electronic propagation algorithm, details of the cloning procedure, and details of the exact quantum propagation.

Acknowledgements  We thank Mike Esch, Joe Subotnik, and Sergei Tretiak for useful discussions. This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award # DE-SC0018432. The initial development of the pySpawn software package was funded by the National Science Foundation under grant CHE-1565634.

16 ACS Paragon Plus Environment

Page 16 of 20

Page 17 of 20 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

References  (1) Warshel, A.; Karplus, M. Semiclassical Trajectory Approach to Photoisomerization. Chem. Phys. Lett. 1975, 32, 11-17. (2) Car, R.; Parrinello, M. Unified Approach for Molecular-Dynamics and DensityFunctional Theory. Phys. Rev. Lett. 1985, 55, 2471-2474. (3) Rothlisberger, U.; Klein, M. L. Ab-Inito Molecular-Dynamics Investigation of Singlet C2H2LI2 - Determination of the Ground-State Structure and Observation of LiH Intermediates. J. Am. Chem. Soc. 1995, 117, 42-48. (4) Ben-Nun, M.; Quenneville, J.; Martinez, T. J. Ab Initio Multiple Spawning: Photochemistry from First Principles Quantum Molecular Dynamics. J. Phys. Chem. A 2000, 104, 5161-5175. (5) Worth, G. A.; Robb, M. A.; Burghardt, I. A Novel Algorithm for Non-Adiabatic Direct Dynamics using Variational Gaussian Wavepackets. Faraday Discuss. 2004, 127, 307-323. (6) Tully, J. C., Molecular-Dynamics with Electronic-Transitions. J. Chem. Phys. 1990, 93, 1061-1071. (7) Hammes-Schiffer, S.; Tully, J. C. Proton-Transfer in Solution - Molecular-Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657-4667. (8) Tully, J. C. Mixed Quantum-Classical Dynamics. Faraday Discuss. 1998, 110, 407-419. (9) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martínez, T. J. Implementation of Ab Initio Multiple Spawning in the Molpro Quantum Chemistry Package. Chem. Phys. 2008, 347, 3-16. (10) Curchod, B. F. E.; Martinez, T. J. Ab Initio Nonadiabatic Quantum Molecular Dynamics. Chem. Rev. 2018, 118, 3305-3336. (11) Burghardt, I.; Meyer, H. D.; Cederbaum, L. S. Approaches to the Approximate Treatment of Complex Molecular Systems by the Multiconfiguration Time-Dependent Hartree Method. J. Chem. Phys. 1999, 111, 2927-2939. (12) Rego, L. G. C.; Batista, V. S. Quantum Dynamics Simulations of Interfacial Electron Transfer in Sensitized TiO2 Semiconductors. J. Am. Chem. Soc. 2003, 125, 7989-7997. (13) Craig, C. F.; Duncan, W. R.; Prezhdo, O. V. Trajectory Surface Hopping in the TimeDependent Kohn-Sham Approach for Electron-Nuclear Dynamics. Phys. Rev. Lett. 2005, 95, 163001. (14) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106. (15) Akimov, A. V.; Neukirch, A. J.; Prezhdo, O. V. Theoretical Insights into Photoinduced Charge Transfer and Catalysis at Oxide Interfaces. Chem. Rev. 2013, 113, 4496-4565. (16) Patchkovskii, S.; Schuurman, M. S. Short-Time Dynamics at a Conical Intersection in High-Harmonic Spectroscopy. J. Phys. Chem. A 2014, 118, 12069-12079. (17) Dou, W. J.; Nitzan, A.; Subotnik, J. E. Surface Hopping with a Manifold of Electronic States. II. Application to the Many-Body Anderson-Holstein Model. J. Chem. Phys. 2015, 142, 084110. (18) Kilina, S.; Kilin, D.; Tretiak, S. Light-Driven and Phonon-Assisted Dynamics in Organic and Semiconductor Nanostructures. Chem. Rev. 2015, 115, 5929-5978. (19) Min, S. K.; Agostini, F.; Gross, E. K. U. Coupled-Trajectory Quantum-Classical Approach to Electronic Decoherence in Nonadiabatic Processes. Phys. Rev. Lett. 2015, 115, 073001.

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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

(20) Schleife, A.; Kanai, Y.; Correa, A. A. Accurate Atomistic First-Principles Calculations of Electronic Stopping. Phys. Rev. B 2015, 91, 014306. (21) Askerka, M.; Maurer, R. J.; Batista, V. S.; Tully, J. C. Role of Tensorial Electronic Friction in Energy Transfer at Metal Surfaces. Phys. Rev. Lett. 2016, 116, 217601. (22) Kilina, S. V.; Tamukong, P. K.; Kilin, D. S. Surface Chemistry of Semiconducting Quantum Dots: Theoretical Perspectives. Acc. Chem, Res. 2016, 49, 2127-2135. (23) Mignolet, B.; Curchod, B. F. E.; Martinez, T. J. Communication: XFAIMS-eXternal Field Ab Initio Multiple Spawning for Electron-Nuclear Dynamics Triggered by Short Laser Pulses. J. Chem. Phys. 2016, 145, 191104. (24) Donati, G.; Lingerfelt, D. B.; Aikens, C. M.; Li, X. S. Molecular Vibration Induced Plasmon Decay. J. Phys. Chem. C 2017, 121, 15368-15374. (25) Maurer, R. J.; Jiang, B.; Guo, H.; Tully, J. C. Mode Specific Electronic Friction in Dissociative Chemisorption on Metal Surfaces: H-2 on Ag(111). Phys. Rev. Lett. 2017, 118, 256001. (26) Dou, W. J.; Subotnik, J. E. Perspective: How to Understand Electronic Friction. J. Chem. Phys. 2018, 148, 230901. (27) Schuurman, M. S.; Stolow, A. Dynamics at Conical Intersections. In Annual Review of Physical Chemistry, Vol 69, Johnson, M. A.; Martinez, T. J. Eds. Annual Reviews: Palo Alto, 2018; Vol. 69, pp 427-450. (28) Meyer, H. D.; Miller, W. H. Classical Analog for Electronic Degrees of Freedom in NonAdiabatic Collision Processes. J. Chem. Phys. 1979, 70, 3214-3223. (29) Negele, J. W. The Mean-Field Theory of Nuclear-Structure and Dynamics. Rev. Mod. Phys. 1982, 54, 913-1015. (30) Billing, G. D. On the Use of Ehrenfest's Theorem in Molecular Scattering. Chem. Phys. Lett. 1983, 100, 535-539. (31) Micha, D. A. A Self-Consistent Eikonal Treatment of Electronic-Transitions in Molecular-Collisions. J. Chem. Phys. 1983, 78, 7138-7145. (32) Isborn, C. M.; Li, X.; Tully, J. C. Time-dependent density functional theory Ehrenfest dynamics: collisions between atomic oxygen and graphite clusters. J. Chem. Phys. 2007, 126, 134307. (33) Vacher, M.; Mendive-Tapia, D.; Bearpark, M. J.; Robb, M. A. The Second-Order Ehrenfest Method A Practical CASSCF Approach to Coupled Electron-Nuclear Dynamics. Theor. Chem. Acc. 2014, 133, 1505. (34) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061-1071. (35) Bittner, E. R.; Rossky, P. J. Quantum Decoherence in Mixed Quantum-Classical Systems - Nonadiabatic Processes. J. Chem. Phys. 1995, 103, 8130-8143. (36) Prezhdo, O. V.; Rossky, P. J. Mean-Field Molecular Dynamics with Surface Hopping. J. Chem. Phys. 1997, 107, 825-834. (37) Fang, J. Y.; Hammes-Schiffer, S. Improvement of the Internal Consistency in Trajectory Surface Hopping. J. Phys. Chem. A 1999, 103, 9399-9407. (38) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919-8929. (39) Nielsen, S.; Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Surface Hopping Dynamics. J. Chem. Phys. 2000, 112, 6543-6553.

18 ACS Paragon Plus Environment

Page 18 of 20

Page 19 of 20 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

(40) Zhu, C. Y.; Nangia, S.; Jasper, A. W.; Truhlar, D. G. Coherent Switching with Decay of Mixing: An Improved Treatment of Electronic Coherence for Non-Born-Oppenheimer Trajectories. J. Chem. Phys. 2004, 121, 7658-7670. (41) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-field Dynamics with Stochastic Decoherence (MF-SD): A New Algorithm for Nonadiabatic Mixed Quantum/Classical Molecular-Dynamics Simulations with Nuclear-Induced Decoherence. J. Chem. Phys. 2005, 123, 234106. (42) Jasper, A. W.; Nangia, S.; Zhu, C. Y.; Truhlar, D. G. Non-Born-Oppenheimer Molecular Dynamics. Acc. Chem, Res. 2006, 39, 101-108. (43) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-Induced Surface Hopping. J. Chem. Phys. 2012, 137, 22a545. (44) Subotnik, J. E.; Ouyang, W. J.; Landry, B. R. Can We Derive Tully's Surface-Hopping Algorithm from the Semiclassical Quantum Liouville eEquation? Almost, but only with Decoherence. J. Chem. Phys. 2013, 139, 214107. (45) Cotton, S. J.; Miller, W. H. Symmetrical Windowing for Quantum States in QuasiClassical Trajectory Simulations. J. Phys. Chem. A 2013, 117, 7190-7194. (46) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics: Treatment of Electronic Decoherence. J. Chem. Phys. 2013, 138, 224111. (47) Ouyang, W. J.; Subotnik, J. E. Estimating the Entropy and Quantifying the Impurity of a Swarm of Surface-Hopping Trajectories: A New Perspective on Decoherence. J. Chem. Phys. 2014, 140, 204102. (48) Makri, N. Quantum-Classical Path Integral: A Rigorous Approach to Condensed Phase Dynamics. Int. J. Quantum Chem. 2015, 115, 1209-1214. (49) Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W. J.; Bellonzi, N. Understanding the Surface Hopping View of Electronic Transitions and Decoherence. In Annual Review of Physical Chemistry, Vol 67, Johnson, M. A.; Martinez, T. J., Eds. Annual Reviews: Palo Alto, 2016; Vol. 67, pp 387-417. (50) Min, S. K.; Agostini, F.; Tayernelli, I.; Gross, E. K. U. Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence. J. Phys. Chem. Lett. 2017, 8, 3048-3055. (51) Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V. Ab Initio Multiple Cloning Algorithm for Quantum Nonadiabatic Molecular Dynamics. J. Chem. Phys. 2014, 141, 054110. (52) Freixas, V. M.; Fernandez-Alberti, S.; Makhov, D. V.; Tretiak, S.; Shalashilin, D. An Ab Initio Multiple Cloning Approach for the Simulation of Photoinduced Dynamics in Conjugated Molecules. Phys. Chem. Chem. Phys. 2018, 20, 17762-17772. (53) Green, J. A.; Makhov, D. V.; Cole-Filipiak, N. C.; Symonds, C.; Stavros, V. G.; Shalashilin, D. V. Ultrafast Photodissociation Dynamics of 2-Ethylpyrrole: Adding Insight to Experiment with Ab Initio Multiple Cloning. Phys. Chem. Chem. Phys. 2019, 21, 3832-3841. (54) Makhov, D. V.; Martinez, T. J.; Shalashilin, D. V. Toward Fully Quantum Modelling of Ultrafast Photodissociation Imaging Experiments. Treating Tunnelling in the Ab Initio Multiple Cloning Approach. Faraday Discuss. 2016, 194, 81-94. (55) Makhov, D. V.; Saita, K.; Martinez, T. J.; Shalashilin, D. V. Ab Initio Multiple Cloning Simulations of Pyrrole Photodissociation: TKER Spectra and Velocity Map Imaging. Phys. Chem. Chem. Phys. 2015, 17, 3316-3325.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 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

(56) Saita, K.; Shalashilin, D. V. On-the-Fly Ab Initio Molecular Dynamics with Multiconfigurational Ehrenfest Method. J. Chem. Phys. 2012, 137, 22a506. (57) Shalashilin, D. V. Quantum Mechanics with the Basis Set Guided by Ehrenfest Trajectories: Theory and Application to Spin-Boson Model. J. Chem. Phys. 2009, 130, 244101. (58) Shalashilin, D. V. Nonadiabatic Dynamics with the Help of Multiconfigurational Ehrenfest Method: Improved Theory and Fully Quantum 24D Simulation of Pyrazine. J. Chem. Phys. 2010, 132, 244111. (59) Heller, E. J. Time-Dependent Approach to Semiclassical Dynamics. J. Chem. Phys. 1975, 62, 1544-1555. (60) Blanes, S.; Casas, F.; Murua, A. Symplectic Splitting Operator Methods for the TimeDependent Schrodinger Equation. J. Chem. Phys. 2006, 124, 234105. (61) Jasper, A. W.; Truhlar, D. G. Electronic Decoherence Time for Non-Born-Oppenheimer Trajectories. J. Chem. Phys. 2005, 123, 064103. (62) White, A.; Tretiak, S.; Mozyrsky, D. Coupled Wave-Packets for Non-Adiabatic Molecular Dynamics: a Generalization of Gaussian Wave-Packet Dynamics to Multiple Potential Energy Surfaces. Chem. Sci. 2016, 7, 4905-4911. (63) Peng, W. T.; Fales, B. S.; Levine, B. G. Simulating Electron Dynamics of Complex Molecules with Time-Dependent Complete Active Space Configuration Interaction. J. Chem. Theory Comput. 2018, 14, 4129-4138.

20 ACS Paragon Plus Environment

Page 20 of 20