An Improved Ehrenfest Approach to Model Correlated Electron

Spectroscopy and is published by the Photochemistry; General American Chemical Society. 1155 Sixteenth Theory Street N.W., Washington, DC 20036by Iowa State Subscriber access provided Published by American University | Library Chemical Society. Copyright © American Chemical Society. However, no copyright

An Improved Ehrenfest isApproach published by the American Chemical Society. 1155 Sixteenth to Model Correlated Street N.W., Washington, DC 20036by Iowa State Subscriber access provided Published by American University | Library Chemical Society. Copyright © American Chemical Society. However, no copyright

Electron-Nuclear Dynamics is published by the American Chemical

Roman Baskov, Alexander J Society. 1155 Sixteenth Street N.W., Washington, White, and Dmitry Mozyrsky DC 20036

Subscriber access provided by Iowa State Published by American University | Library Chemical Society. Copyright © American Chemical Society. However, no copyright

J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.8b03061 is published by the• American Publication Date (Web): 09Chemical Jan 2019

Society. 1155 Sixteenth

DownloadedStreet fromN.W., http:// Washington, DC 20036by 9, pubs.acs.org onprovided January 2019 Subscriber access Iowa State

Published by American University | Library Chemical Society. Copyright © American Chemical Society. However, no copyright

Just Accepted

published by the “Just Accepted” is manuscripts have been pe American Chemical online prior to technical editing, formatting Society. 1155 Sixteenth Washington, Society provides Street “JustN.W., Accepted” as a serv DC 20036by Iowa State access provided ofSubscriber scientific material as soon as Published by American possible

University | Library

Chemical Society. Copyright © American Chemical Society. However, no copyright

full in PDF format accompanied by an HT peer reviewed, but should not be conside is published(DOI®). by the Digital Object Identifier “Just Acc American Chemical the “Just Accepted” Web site may not inc Society. 1155 Sixteenth N.W., Washington, a manuscript is Street technically edited and for DC 20036by Iowa State Subscriber access provided site and published as by anAmerican ASAP article. Published University | Library

Chemical Society. Copyright © American Chemical Society. However, no copyright

to the manuscript text and/or graphics w ethical guidelines that apply to the journ is published by the consequences arising from the use of inf American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036by Iowa State Subscriber access provided Published by American University | Library Chemical Society. Copyright © American Chemical Society. However, no copyright

0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

a)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.0 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.01.0 0.80.8 0.60.6 0.40.4 0.20.2 0.00.0 -0.2-0.2 -0.4-0.4

0.8 0.6 0.4

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

+ Page The 1 of Journal 30 of Physical Chemistry Letters + + 0

+

& '(-%0.2. 0 +'+ , 0.0 -0.2 -0.4 -0.60 0 x (a.u.) x (a.u.) -0.8 -1.0

-0.6-0.6 -0.8-0.8 -1.0-1.0

0

& '(-% . 0 +'+ ,

+

0

!" ⋅ $0 % + " ⋅ 2 − 2% Δ50 0 !" ⋅ $% + " ⋅ 2 −0 2% Δ5 x (a.u.)

x (a.u.)

(a.u.) x x(a.u.)

x (a.u.)

1 2 3 4 5 6 7 8

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.41.0 -0.60.8 -0.80.6 0.4 -1.00.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

1.0 0.8

& '(-%

b)

0

0.6 1.0 1.0 0.4 1.0 0.8 0.8 0.2 0.8 0.6 . / +'+ 0.6 ,0.0 0.6 0.4 -0.2 0.4 -0.4 0.4 0.2 0.0 -0.6 0.2 0.2 -0.2 -0.8 -0.4 0.0 -1.0 0.0 -0.6 -0.2 -0.2 -0.8 1.0 1.0 -1.0 -0.4 -0.4 0.8 0.8 -0.6 -0.6 1.0 0.6 0.6 0.8 -0.8 -0.8 1.0 0.4 0.4 0.6 0.8 -1.0 0.2 0.4 0.2 -1.0

x (a.u.) 0

x (a.u.)

0.6

x (a.u.)

!"x (a.u.) ⋅ $% & '()* 0

& 0

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

0.2 0.4 0.0 0.2 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

00

x x(a.u.) (a.u.)

x (a.u.)

1.4 1.2 1.0 0.8 0.6 0.4 0.2 1.0 0.8 0.0 0.6 -0.2 0.4 -0.4 0.2 -0.6 0.0 -0.2 -0.8 -0.4 -1.0

0

0.2 0.0 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

1.0 0.8

+ +

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

&

0.00.2

&

(a.u.) x (a.u.) -0.2 x -0.20.0 -0.20.0 -0.2 -0.2 -0.4 -0.4 -0.4 -0.4 -0.4 -0.6 -0.6 -0.6 -0.6 -0.6 -0.8 -0.8 -0.8 -0.8 -0.8 -1.0 -1.0 -1.0 -1.0 0 -1.0 00 0 0 0 0 (a.u.) x (a.u.) xx0(a.u.) x (a.u.) (a.u.) xx(a.u.) 0.00.2

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8

0.2 0.0 -0.2 0.0 -0.2 -0.2 -0.4 -0.4 -0.4 -0.6 -0.6-0.6 -0.8 -0.8-0.8 -1.0 -1.0

0.0 x (a.u.) -0.2 -0.2 -0.4-0.4 -0.6-0.6 -0.8-0.8 -1.0 0 -1.0 0 x (a.u.)0 x (a.u.)

0

00

0 (a.u.) xx(a.u.)

x (a.u.) x (a.u.)

& '(-% . / +'+ ,

%

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

'()* +'+ ,

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

& '(-% . / +'+ ,

x (a.u.)

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

-1.0

%

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 1.0 -0.6 0.8 -0.8 0.6 -1.0

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

& '(-% . / +'+ ,

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 0 x (a.u.) -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 0.4 -0.2 0.2 -0.4 0.0 -0.6 -0.2 -0.8 -0.4 -1.0 -0.6 -0.8 0

+ +

-1.0 x (a.u.)

0 x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

|8′⟩ 0

0

x (a.u.)

time 0

x (a.u.)

1.0 0.8 1.0 0.6 0.8 0.4 0.6 0.2 0.4 0.0 0.2 -0.2 0.0 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0 0

'()* +'+ ,

!" ⋅ $ & !" ⋅ $ & Plus Environment ACS Paragon

+'+ ,

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.) 0

x (a.u.)

1.4 1.2 1.0 0.8 0.6 0.4 1.0 0.2 0.8 1.0 0.8 0.0 0.6 0.6 -0.2 0.4 -0.4 0.4 0.2 -0.6 0.0 0.2 -0.2 -0.8 -1.0 0.0 -0.4

0.6 0.6 1.01.0 0.6 0.41.0 0.41.0 0.4 0.80.8 0.20.8 0.20.8 0.2 0.6 '(-% . /0.60.6 +'+ , 0.00.6 0.00.4 '(-% . / +'+0.0 ,0 -0.2 -0.20.4 x (a.u.) -0.2 0.40.4 0.2 -0.4 -0.40.2 -0.4 0.0 0.0 -0.6 0.2 0.2 -0.6 -0.2 -0.6 -0.2 -0.8 -0.8 -0.4 -0.8 -0.4 0.00.0 0 -1.0 -1.0 -0.6 -0.6 0 -0.6 -0.6 -1.0 0 -0.8 0 -0.2 -0.2x (a.u.) -0.2 -0.8 x (a.u.) -0.8 -0.8 -1.0 -1.0 x (a.u.) x (a.u.) -1.0 -1.0 1.0 1.0 1.0 -0.4 -0.4 -0.41.0 0 0 0 0 0.8 x 0.8 (a.u.) 0.8 x-0.6 (a.u.) x (a.u.) x (a.u.) -0.6 -0.60.8 1.0 1.0 1.0 0.6 1.0 0.6 0.6 1.0 0.6 0.8 0.81.0 0.8 0.8 -0.8 -0.8 -0.8 1.0 0.8 1.0 0.4 0.61.0 0.4 0.6 0.4 0.60.8 0.6 0.4 0.8 0.8 0.8 -1.0 -1.0 +'+ -1.00.2 0.6 +'+0.2 0.6 0.4 0.2 0.40.6 0.4 0.2 ,0.40.6 '(-%. 0.6 '(-%. 0.4 0 0 , 0 0 0 0.20.4 0.2 0.0 0.4 0.4 0.2 0.0 0.20.4 0.0 0.2

&

|8 ⟩

!" ⋅ $% + " ⋅ 2 − 2% Δ5

x (a.u.) 1.0 0.8

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

00

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

|8 ⟩

x (a.u.) x (a.u.)

0

x (a.u.)

x (a.u.) 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6

|8′⟩

E2

A12/15

E The Journal of Physical Chemistry Letters 2

90

90

E1

%

60

60

30

a) 20

2

Difference from exact (%)

1

A12

Refl. 2

E2

Trans. 1

20

Refl. 1

E1

b)

30

k (a.u.)

Page 2 of 30

Trans. 2

40

%

30

10

Trans. 1

E1

% Trans. 1

20

30

c) 40

k (a.u.)

%

15

Trans. 1

30

45

k (a.u.)

%

90

90

V33

V11

60 V12

0 Momentum (a.u)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

A12/100

10

15

20

25

30

30 0

15

0

10

V13

V22

60 E2

-0.04

30

E1

-0.08 -0.12

d)

5 0 10

12

14

16

k (a.u.)

18

e) 4 0 ACS Paragon0 4Plus -4 Environment

f)

-4

20

30

k (a.u.)

0 1000

2000

time (a.u.)

3000

-0.6

-1.0 0

-1.0

0.6 -0.8 0.4

0 0.4 0.2 x (a.u.)

x (a.u.)0.2 -1.0

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0 0 e Journal Page 3 of Physical 30 Chemistry Lett x (a.u.) x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0

1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6

x

1.0 0.8 0.6 0.4 0 0.2 0.0 (a.u.) -0.2 -0.4 -0.6 -0.8 -1.0

1 !" ⋅ $% & '()* +'+, 2 3 !" ⋅ $ % x (a.u.) 4 ⋅ . − .% Δ1 -1.0 +" 0 0 5 x (a.u.) x (a.u.) ACS 6 Paragon Plus Environment + 7 8 x

0

1.0 0.4 0.8 0.6 0.2 0 0.4 0.0 (a.u.)0.2 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0

x

0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

x (a.u.)

1.0 0.8 0.6 0 0.4 0.2 (a.u.) 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

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

An Improved Ehrenfest Approach to Model Correlated Electron-Nuclear Dynamics Roman Baskov,† Alexander White,∗,‡ and Dmitry Mozyrsky∗,‡ †Institute of Physics of the National Academy of Sciences of Ukraine, pr. Nauky 46, Kyiv-28, MSP 03028, Ukraine ‡Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA E-mail: [email protected]; [email protected]

1

ACS Paragon Plus Environment

Page 4 of 30

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

Abstract Mixed quantum-classical mechanical descriptions are critical to modeling coupled electron-nuclear dynamics, i.e. non-adiabatic molecular dynamics, relevant to photochemical and photophysical processes. We introduce an efficient description of such dynamics in terms of an effective Hamiltonian that not only properly captures electronnuclear correlation effects but also helps developing an efficient computational method. In particular we introduce a coupled Gaussian wavepacket parameterization of the nuclear wavefunction, which generalizes the Ehrenfest approach to account for electronnuclei correlations. We test this new approach, Ehrenfest-Plus, on a suite of model problems that probe electron-nuclear correlation in non-adiabatic transitions. The high accuracy of our approach, combined with mixed-quantum classical efficiency, opens a path for improved simulation of non-adiabatic molecular dynamics in realistic molecular systems.

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

Graphical TOC Entry

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0

0

x (a.u.) 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.8 0.4 0.6 0.2 0.0 0.4 0.2 -0.2 -0.4 0.0 -0.6 -0.2 -0.8 -0.4 -1.0 -0.6 -0.8 -1.0 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 1.0 0.4 0.8 0.6 0.2 0 0.4 0.0 x (a.u.)0.2 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0 1.0

0

x (a.u.)

+

0 1.0 0.8 0.6 0.4 0 0.2 0.0 (a.u.) -0.2 -0.4 -0.6 -0.8 -1.0

x (a.u.)

x

0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0 0.4 0.2 (a.u.) 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

0

x (a.u.)

0

0

0

x (a.u.)

x (a.u.)

x (a.u.)

2

x

x

x

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.4 1.2 1.0 0.8 0.6 0 0.4 0.2 (a.u.) 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

!" ⋅ $% & '()* +'+,

!" ⋅ $% x (a.u.) +" ⋅ . − .% Δ1 0

x

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.)

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 1.0 0.8 -0.6 0.6 0 0.4 -0.8 x (a.u.)0.2 -1.0


1.0 0.8 0.6 0.4 0 0.2 1.0 (a.u.) 0.0 0.8 -0.2 -0.4 0.6 -0.6 0.4 -0.8 0.2 -1.0

0.0 -0.2 -0.4 -0.6 -0.8 0 -1.0

x (a.u.)

1.0 0.8 0.6 00.4 (a.u.) 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 0

x (a.u.)

x

1.0 0.8 0.6 00.4 0.2 (a.u.) 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 0

0

x (a.u.)

x (a.u.)


Ab-initio methods play an important role in the simulation of electronic and coupled electronic-nuclear processes in atomic, molecular, soft and condensed matter systems. 1–3 Light absorption and emission, 4,5 charge and energy transfer, 6–9 photoisomerization and photochemistry, 4,10 non-radiative relaxation, 11–13 swift-ion stopping, 14–16 etc. are all influenced or controlled by the interaction of electronic excitation and nuclear motion. These non-adiabatic processes are particularly difficult to model, due to the interaction of the quantum electrons and the nuclei. For the above processes, the two cannot be separated as in the traditional Born-Oppenheimer approximation. A fully quantum mechanical treatment of all the nuclei is computationally prohibitive for most systems. 17–19 Mixed quantum-classical mechanics, which attempt to treat the nuclei semi-classically are therefore highly desirable. One of the most successful first-principles mixed quantum-classical algorithms is the mean-field Ehrenfest approach, 20–23 in which the nuclei are subject to a classical force determined by the instantaneous average density of the electrons. However, as a mean-field approach, it cannot properly describe correlated electron-nuclear process. Ad-hoc methods have been developed to treat these correlation effects, while maintaining computational efficiency. 24 Such methods, however, are unpredictable in terms of their accuracy and no consensus exists on how to treat various situations. However, great efforts have been applied to correct the approaches. 25–44 First-principles based, semi-classical approximations have long been developed, but are often prohibitively expensive for high dimension systems like polyatomic molecules or require inaccurate approximations. 45–63 The development of abinitio non-adiabatic molecular dynamics (NAMD) approaches that closely tie to accessible mixed-quantum classical techniques can provide insight into these ad-hoc algorithms, and potentially provide more reliable, but still numerically tractable, simulations. The ab initio multiple spawning and cloning approachs, 64–68 where the exact time dependent Schrödinger equation (TDSE) is solved in the nonorthogonal basis of fixed width Gaussians, has been highly successful in this regard. 2,69–71 While there are many approaches to generate these basis functions, 36,64,66,72–74 they are ad-hoc, and possibly problem dependent. However, the

3


Page 6 of 30

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


full spawning limit, where the basis is complete by definition, may be prohibitively expensive for many systems of interest. 75 Methods in mixed quantum-classical dynamics have been recently reviewed in Refs. 2,29,36,49,76 . Here we develop an approach that is designed for on-the-fly dynamics where electronic structure is calculated for a specific nuclear configuration at a given time, rather than as a continuous function of nuclear configurations. First we provide a general path integral framework using local adiabatic eigenstates. From this general framework the mean-field Ehrenfest approach can be derived, see Supporting Information. Our extended approach, Ehrenfest-Plus, accounts for fluctuations caused by electron-nuclear correlations. We test this new approach on a set of standard models, which probe electron-nuclear correlation in non-adiabatic quantum scattering. However, we first begin by briefly reviewing the conventional approach to non-adiabatic processes in molecular systems.

Non-adiabatic Hamiltonians A traditional description of NAMD is based on the matrix vector-potential picture, where the transitions between electronic states with energies En (x) parametrically dependent on 3N -component coordinate vector x for the nuclear positions are described in terms of an effective “velocity gauge” Hamiltonian, 77,78

ˆ vg (x) = H

3N ∑

1 ˆ [ˆ pµ − iAˆµ (x)]2 + E(x), 2M µ µ=1

(1)

where the non-adiabatic coupling vector (NACV), Ann′ ,µ (x) = ⟨n(x)|∂µ n′ (x)⟩, and potential energy surface (PES) scalar, Enn′ (x) = En (x)δnn′ , potentials are matrices in the subspace spanned by the adiabatic electronic states |n(x)⟩ (eigenstates of the electronic Hamiltonian, ˆ e ). pˆµ is the conventional momenta operator acting in the x-space, pˆµ = −i∂µ . En (x) is H commonly referred to as the potential energy surface (PES) of state n. One can arrive at Eq. (1) by expanding the full molecular state |Ψ(x)⟩ as a superposition

4



Page 8 of 30

of these adiabatic electronic eigenstates |n(x)⟩,

Ψn (x, t) =

∑

ψn (x, t) |n(x)⟩ ,

(2)

n

by substituting Eq. (2) into the time-dependent Schrödinger equation. Then, after standard manipulations, 77 one arrives at the Scrödinger equation for ψn ’s with Hamiltonian given by Eq. (1). However, for more than a few nuclear degrees of freedom, calculation of the PES and NACV matrices is numerically prohibitive. Thus, for most systems, on-the-fly ab-initio MD methods, where classical or nearly classical trajectories guide the calculation of the PES and NACV matrices, are desirable. In these MD simulations the basis states |n⟩ are evaluated ¯ (t) along the trajectory. Since x ¯ changes with time, only locally, i.e. for a given position x the MD states |n(t)⟩ are time-dependent (rather than x-dependent, as we have assumed previously in the derivation of Eq. (1)). Substitution of the wavefunction in Eq. (2), but with |n⟩ being t-dependent rather than x-dependent, again into the time-dependent Schrödinger equation gives another effective Hamiltonian, which we term as the “length gauge” Hamiltonian,

ˆ lg = H

∑ pˆ2µ µ

2Mµ

+

∑

Vnn′ (x, t)|n′ (t)⟩⟨n(t)| ,

(3)

n,n′

with ˆ e (x)|n′ (t)⟩ . Vnn′ (x, t) = i⟨n(t)|∂t n′ (t)⟩ + ⟨n(t)|H

(4)

Local Approximation The potential energy in Eq. (4) can be put in a more transparent form if we assume that the ¯ (t). Then molecular wavefunction is sufficiently localized in the x-space around a position x ˆ e as we expand H ˆ e (x) = H ˆ e (¯ ˆ e (¯ ¯ ] · [x − x ¯ ] + ... . H x) + [∂ H x)/∂ x 5


(5)

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


When substituted into Eqs. (3, 4), the second and the higher order terms in Eq. (5) are assumed to be relatively small due to the expected localized nature of the nuclear states ¯ . One can show that the neglect of the offψn (x), centered around the average coordinates x diagonal forces (and the higher order terms) in a certain time-dependent basis leads to the celebrated Ehrenfest method, where the nuclei are propagated along an average trajectory according to an average (classical) force. We will see in a moment, however, that the linear terms in Eq. (5) account for important electron-nucleus correlations and therefore need to be retained for accurate description of non-adiabatic dynamics. The length gauge Hamiltonian in Eqs. (3, 4) can be set in a more transparent form ¯ (t)). Then we can readily if choose |n(t)⟩ to be a local adiabatic basis (i.e. at point x evaluate the matrix elements in Eq. (4). For n = n′ the second matrix element in Eq ¯ ), where fn is the classical electron-nucleus force for the n’s (4) gives En (¯ x) − fn (¯ x)(x − x PES, En (¯ x). For n ̸= n′ , by virtue of Hellmann-Feynman theorem, this matrix element is ¯ ), where Ann′ is the n, n′ element of the NACV matrix (Eq. 1) and ∆Enn′ (¯ x)Ann′ (¯ x)(x − x ∆Enn′ = En −En′ . Furthermore, in the absence of magnetic field, the electronic Hamiltonian ˆ e is real and therefore the states |n(t)⟩ can also be chosen real. Then, for diagonal n = n′ H terms, the first matrix element in Eq. (4) vanishes, while the off-diagonal Vnn′ can be written as ¯ + ∆Enn′ Ann′ · (x − x ¯ ), n ̸= n′ . Vnn′ = iAnn′ · v

(6)

¯ , x, 0)|n′ (0)⟩, Let us assume that at time t = 0 the wavefunction is given by |Ψ(0)⟩ = g(¯ x, p ˆ x(t)]+i¯ p(t)·[x−¯ x(t)] ¯ , x, t) = N ei[x−¯x(t)]·α(t)·[x−¯ g(¯ x, p .

(7)

Re Im where α ˆ is a complex matrix, αµν = αµν + iαµν and N = [23N det(ˆ αℑ )/π 3N ]1/4 is the normal-

ization coefficient. The Gaussian in Eq. (7) describes a nuclear wavepacket centered around ¯ and having momentum p ¯ . Application of the Hamiltonian in Eq. (3), classical position x with the off-diagonal term in Eq. (6), on |Ψ(0)⟩ would create two wavepackets “on” state 6



Page 10 of 30

¯ ). Repeated applications of |n(t)⟩, one a Gaussian and the other a Gaussian scaled by (x − x the interaction will create higher order polynomial-scaled Gaussians, see adiabatic expansion diagram in Fig. 1-a. Thus it becomes difficult to maintain a minimal basis, e.g. one Gaussian on each state. Many wavepacket methods, 64,68,79 e.g. ab-initio multiple spawning and Coupled Coherent States, essentially project these polynomial-scaled Gaussians back onto a complete set of unscaled Gaussians. However, in practice, the Gaussian basis is rarely truly 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

complete. Thus capturing the best possible result with a minimal basis is desirable. 0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

a)

0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.)

x (a.u.)

1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.0 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

x (a.u.)

1.0 0.8 0.6 '(-%0.4. 0 +'+ , 0.2 0.0 -0.2 -0.4 -0.60 0 x (a.u.) x (a.u.) -0.8 % -1.0

1.01.0 0.80.8 0.60.6 0.40.4 0.20.2 0.00.0 -0.2-0.2 -0.4-0.4 -0.6-0.6 -0.8-0.8 -1.0-1.0

+

&

0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

& '(-% . 0 +'+ ,

+ 0

!" ⋅ $0 % + " ⋅ 2 − 2% Δ50 0 !" ⋅ $ + " ⋅ 2 −0 2% Δ5 x (a.u.)

x (a.u.)

(a.u.) x x(a.u.)

x (a.u.) 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.41.0 -0.60.8 -0.80.6 0.4 -1.00.2

0

0

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 1.0 0.8 0.8

1.0 0.8 0.6 0.4 1.0 0.2 0.8 0.6 0.0 0.4 -0.2 0.2 -0.4 0.0 -0.6 -0.2 -0.8 -0.4 -1.0 -0.6

0.6 , 0.6 & '(-% . / +'+ 0.4 0.4 0.2 0.2 0.0 0.0 -0.2 -0.2

b)

0

x (a.u.)

!"x (a.u.) ⋅ $% & '()* 0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 1.0 -0.4 -0.4 0.8 0.8 -0.6 -0.6 0.6 0.6 -0.8 -0.8 0.4 0.4 -1.0 0.2 0.2 -1.0 0.0 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

+'+ ,

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

-0.8 -1.0

x (a.u.)

0

-0.6 -0.8 -1.0

0

1.0 0.8 1.0 0.6 0.8 0.4 0.6 0.2 0.4 0.0 0.2 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

x (a.u.)

00

x x(a.u.) (a.u.)

0

1.0 0.8

1.0 0.6 0.8 0.6 0.4 0.4 0.2 0.2 0.0 0.0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0

1.0 0.8 0.6 0.41.0 0.20.8 0.00.6 -0.20.4 -0.40.2 0.0 -0.6 -0.2 -0.8 -0.4 -1.0 -0.6

1.01.0 0.80.8

0.20.2 0.00.0

+ +

0 (a.u.) xx(a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

x (a.u.)

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 0 0

0

x (a.u.) x (a.u.)

x (a.u.)

& '(-% . / +'+ ,

0

0

x (a.u.) x (a.u.)

&

&

+'+ ,

|8 ⟩ 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.) 0

1.4 1.2 1.0 0.8 0.6 0.4 1.0 0.2 0.8 1.0 0.8 0.0 0.6 0.6 -0.2 0.4 -0.4 0.4 0.2 -0.6 0.0 0.2 -0.2 -0.8 -0.4 -1.0 0.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2

'(-% . / +'+0.0 ,

0

x (a.u.) -0.2 -0.4 -0.6 -0.8 -1.0

0

-0.6

-0.2 -0.8 x (a.u.) -1.0 1.0 -0.41.0 0 x 0.8 (a.u.) -0.60.8 1.0 0.6 0.6 1.0 0.8 -0.8 1.0 0.8 0.4 0.6 0.4 0.8 -1.00.2 0.6 +'+0.2 0.6 0.4 '(-%. 0 , 0.40 0.2 0.0 0.4 0.0 0.2 0.2 0.0 -0.2 0.0 -0.2 -0.2 -0.4 -0.4 -0.4 -0.6 -0.6-0.6 -0.8 -0.8-0.8 -1.0 -1.0

0.0 x (a.u.) -0.2 -0.2 -0.4-0.4 -0.6-0.6 -0.8-0.8 -1.0 0 -1.0 0 x (a.u.)0 x (a.u.)

!" ⋅ $% & '()* -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 1.0 -0.6 0.8 -0.8 0.6 -1.0

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 0 x (a.u.) -0.4 -0.6 -0.8 -1.0

1.0 0.8 0.6 0.4 0.2 0.0 0.4 -0.2 0.2 -0.4 0.0 -0.6 -0.2 -0.8 -0.4 -1.0 -0.6 -0.8 0

+ +

-1.0 x (a.u.)

x (a.u.) 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

|8′⟩

time 0

0 x (a.u.)

00

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.) 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

|8 ⟩ 0

x (a.u.) x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

& '(-% . / +'+ ,

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

x (a.u.)

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

0

0

x (a.u.)

x (a.u.)

1.0 0.8 1.0 0.6 0.8 0.4 0.6 0.2 0.4 0.0 0.2 -0.2 0.0 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1.0 -1.0 0

+'+ ,

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.)

x (a.u.)

1.0 0.8 0.6 0.41.0 0.20.8 0.00.6 -0.20.4 0.2 -0.4 0.0 -0.6 -0.2 -0.8 -0.4 -1.0 -0.6

& '(-%. 0 +'+ ,

!" ⋅ $% & '()*

00

0

0

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

+ + +

!" ⋅ $% + " ⋅ 2 − 2% Δ5

0 -0.8 0 -0.2 -0.2x (a.u.) -0.8 -1.0 x (a.u.) x (a.u.) -1.0 1.0 1.0 -0.4 -0.4 0 0 0 0.8 0.8 x-0.6 (a.u.) x (a.u.) x (a.u.) -0.6 1.0 0.6 1.0 0.6 0.81.0 0.8 -0.8 -0.8 1.0 0.4 0.6 0.4 0.60.8 0.8 -1.0 0.40.6 0.2 0.40.6 0.2 -1.0 0 0 0.20.4 0.0 0.20.4 0.0 0.00.2 0.00.2 (a.u.) x (a.u.) -0.2 x -0.2 -0.20.0 -0.20.0 -0.2 -0.2 -0.4 -0.4 -0.4 -0.4 -0.4 -0.4 -0.6 -0.6 -0.6 -0.6 -0.6 -0.6 -0.8 -0.8 -0.8 -0.8 -0.8 -0.8 -1.0 -1.0 -1.0 -1.0 0 -1.0 00 0 -1.0 0 0 0 (a.u.) x (a.u.) xx0(a.u.) x (a.u.) (a.u.) xx(a.u.)

x (a.u.) x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

x (a.u.)

x (a.u.)

+'+ , 0.60.6 & '(-% . /0.4 0.4

& '(-% . / +'+ ,

x (a.u.)

1.4 1.2 1.0 0.8 0.6 0.4 0.2 1.0 0.8 0.0 0.6 -0.2 0.4 -0.4 0.2 -0.6 0.0 -0.2 -0.8 -1.0 -0.4

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

0

x (a.u.)

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0

|8′⟩

0

0

0

0

0

0

x (a.u.)

x (a.u.)

x (a.u.)

x (a.u.)

x (a.u.)

x (a.u.)

Figure 1: (a) Wavepackets generated and propagated by application of Hamiltonian (3) with Eq. (6); (b) Wavepackets generated and propagated by “Momentum-Jump” Hamiltonian (9).

Momentum-Jump Approximation and Hamiltonian It is possible to avoid the creation of polynomial-scaled Gaussians by transforming the offdiagonal term in Eq. (6) to: ¯ ei∆pnn′ ·(x−¯x) , n ̸= n′ Vnn′ ≃ iAnn′ · v

(8)

¯ ) and v ¯ ≡ x ¯˙ . The equality of Eq. 6 and Eq. 8 where ∆pnn′ = −∆Enn′ Ann′ /(Ann′ · v ¯ ) ≪ 1, which is the case only if the molecular state |Ψ⟩ is implies that ∆pnn′ · (x − x ¯ , as we have already assumed in Eq. (5). The physical sufficiently localized around position x significance of the phases in the off-diagonal coupling coefficients in Eq. (8) can be readily 7


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


understood. Upon application of the off-diagonal interaction term, Eq. (8), the electronic ¯ changes to p ¯ + ∆pnn′ , state |n′ (t)⟩ again switches to state |n(t)⟩, the initial momentum p and the phase changes accordingly. The momenta of the old and the new wavepackets approximately satisfy classical energy conservation (i.e. for p ¯ ≫ ∆¯ pnn′ ). This is equivalent to the momentum-jump approximation typically applied in the Liouville formualtion of mixed quantum-classical dynamics and utilized in numerous ad-hoc numerical approaches, such as surface hopping. 24,80,81 Thus Eq. (3) can be approximately written as a “momentum-jump” Hamiltonian:

ˆ mj ≃ H

∑ pˆ2µ µ

2Mµ +

+

∑

∑

¯ )]|n⟩⟨n| [En (¯ x) − fn (¯ x) · (x − x

n

¯ ei∆pnn′ ·(x−¯x) |n⟩⟨n′ | , iAnn′ (¯ x) · v

(9)

n̸=n′

where we have used a shorthand notation |n⟩ ≡ |n(t)⟩. We emphasize that the “energy conservation” is a direct consequence of the choice of adiabatic basis set |n(t)⟩ in Eq. (9). If ˆ e (¯ |n(t)⟩ is not to be the eigenstates of H x(t)), the off-diagonal couplings in Eqs. (6) and (9) ¯ ). The role of energy conservation cannot have the exponential form with phases ∆pnn′ ·(x− x and momentum-jumps in NAMD simulations is still questioned, 29,42,81 but here we have shown that it plays an important role in optimizing transfer between wavepackets “on” different electronic states, without significantly altering their Gaussian form, see MomentumJump Hamiltonian diagram in Fig. 1-b. Thus with Eq. (9), it is possible to formulate a Gaussian wavepacket based NAMD scheme using a minimal basis.

Numerical Simulation Approach To simplify, we will assume that there are only two relevant time-dependent locally-adiabatic electronic states, |1(t)⟩ and |2(t)⟩. Generalization to a higher number of electronic states is straightforward. The nuclear wavefunction has two-components and the corresponding

8



Page 12 of 30

Schrödinger equation, using Eq. (9), reads

mj mj iψ˙1 (x, t) = H11 (¯ x)ψ1 (x, t) + H12 (¯ x)ψ2 (x, t)

(10)

mj mj iψ˙2 (x, t) = H22 (¯ x)ψ2 (x, t) + H21 (¯ x)ψ1 (x, t)

¯ (t) to be the mean coordinate, We take x ∫

¯ (t) = x

dx (ψ1∗ xψ1 + ψ2∗ xψ2 ) ,

(11)

ˆ µ−1 . dx (ψ1∗ pˆµ ψ1 + ψ2∗ pˆµ ψ2 )M

(12)

¯ (t) to be the mean velocity, and v ∫

v¯µ (t) =

We assume that initially the system is a local adiabatic state, |Ψ(0)⟩ = ψ1 (x, 0)|1(0)⟩, with nuclear state ψ1 being a Gaussian, i.e. ψ1 (x, 0) = g(x1 , p1 , x, 0) from Eq. (7). Note that, more generally, any initial ψ(x, 0) can be represented by a sum of Gaussians. In the absence of coupling, A12 = 0, ψ1 retains it’s Gaussian form during propagation by Eq. 10 with Eq. (9). The Gaussian coefficients satisfy equations of motions: 82

x˙ 1µ ≡ x˙ 0µ (t) = p˙1µ /Mµ , p˙ 1 = f1 (x1 ), α˙ µν = −2

∑

(13)

αµλ αλν /Mλ .

λ

In the presence of coupling, A12 ̸= 0 we expect that, for short times,

ψ1 (t) = c1 (t) g1 (x), ψ2 (t) = c2 (t) g2 (x) ,

9


(14)

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


( gn (x) ≡ g(xn , pn , x) ) is a good solution, provided that g1 ( x) ∝ ei∆p12 ·(x−¯x) g2 (x),

(15)

during the course of evolution. Condition (15) breaks down when the wavepackets in Eq. (14) spatially separate due to the difference in forces f1 and f2 in Eq. (9). However, if the wavepackets traverse the non-adiabatic regions rapidly, the condition (15) holds approximately while the time-dependent coupling A12 (¯ x) · v ¯ is non-zero. To find the coefficients c1 and c2 we project Eq. (10) for ψ1 onto state ψ2 and vice versa. Then we find that ∫

c˙n (t) =

dx − gn∗ (x)g˙ n (x)cn (t)−i

∑ n′

mj ¯ )gn′ (x)cn′ (t) . gn∗ (x)Hnn ′ (x, x

(16)

Since g1 and g2 are Gaussian functions, e.g. Eq. (7), the matrix elements can be calculated analytically; the explicit expressions are presented in the Supporting Information. Equations of motion for the parameter xn (t), pn (t) and αn,µν (t) defining the Gaussians can be found by relating these quantities to the expectation values of various operators. xn (t) =

∫

dx ψn∗ xψn /|cn |2 and pn (t) =

∫

dx ψn∗ p ˆ ψn /|cn |2 are the state dependent expectation

values of coordinate and momentum operators. Using Eqs. (10) and (14), after some algebra, one arrives at cn c∗n′ ∫ x˙ n = vn − 2Im{ dxgn∗ ′ (x)Hnmj′ n (¯ x)[x − xn ]gn (x)} , 2 |cn | n′ ̸=n ∑

cn c∗n′ ∫ x)[ˆ p − pn ]gn (x)} , p˙ n = fn (¯ x) − dxgn∗ ′ (x)Hnmj′ n (¯ 2Im{ 2 |cn | n′ ̸=n ∑

(17)

with n′ ̸= n. Furthermore, the quantities αnµν (t) can be expressed in terms of the expectation values of the products of coordinate and momentum operators. Specifically, elements of the

10



Page 14 of 30

inverse matrix α−1 is related to the square deviations of x as ∫

and, similarly,

∫

1 Im −1 dx gn (x)[x − xn ] ⊗ [x − xn ]gn (x) = [ˆ α ] , 4 n

(18)

dx gn (x)[x − xn ] ⊗ [ˆ p − pn ]gn (x) = α ˆnRe · [ˆ αnIm ]−1 .

(19)

Then, after some algebra we obtain ˆ −1 ·ˆ α ˆ˙ n = −2ˆ αn ·M αn −

∑ { 2c∗n cn′ ∫ n′ ̸=n

|cn |2

}

mj dx gn∗ (x)Hn,n x)×{4ˆ αnIm ·[x−xn ]⊗[x−xn ]·ˆ αnIm −ˆ αnIm }gn′ (x) . ′ (¯

(20) Eqs. (16, 17) and (20) complete the approximate evolution of the wavefunction according to the Hamiltonian in Eq. (10). While individual nuclear configurations, xn ’s, are propagated, electronic structure calculations are only required for the x ¯ nuclear configuration, similar to traditional Ehrenfest method. Thus the x ¯ and v ¯ pair constitute a single “trajectory”, with the set of Gaussian variables (xn , vn , etc) describing multiple wavepackets. This is a significant departure from traditional Gaussian wavepacket methods, which typically calculate electronic structure at the center of each Gaussian, 64,68,83–89 i.e. one trajectory per wavepacket. Thus, we call this new approach Ehrenfest-Plus (EP), as it incorporates the equation of motion for the state dependent Gaussian variables, in addition to the mean-field variables, while requiring electronic structure only for the mean-field position of the nuclei. Note that since we assume that c1 (0) = 1 and c2 (0) = 0, in order to avoid divergence in the second terms in the rhs of Eqs. 17 and 20, one needs to choose the initial parameters by

x2 (0) = x1 (0) , p2 (0) = p1 (0) + ∆p12 (0) , α ˆ 2 (0) = α ˆ 1 (0).

11


(21)

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


Eqs. (21) uniquely determine the initial values of the second wavepacket (with zero amplitude at t = 0) by guaranteeing that the second terms in Eqs. (17) and (20) are finite. The approximations presented are only valid for a finite time. As the wavepackets separate in phase space, (¯ x ̸= x1 ̸= x2 and p1 ̸= p2 − ∆p12 ), the two wavepacket picture will become incomplete. The equations of motion of the wavepackets must “separate”, allowing new wavepackets to be “spawned”, similar to the ab-initio multiple spawning 64 or decoherence induced surface hopping 41 approaches. Thus, at the spawning point we begin to propagate multiple and independent sets of Eq. (10). This effectively resets the variables, (¯ x = x1 = x2 and p1 = p2 − ∆p12 ), at the cost of propagating an additional set of wavepackets. For example, a single trajectory comprised of two wavepackets on two electronic PES becomes two trajectories with four wavepackets. Those two trajectories are independent after branching. While reducing the time between “spawns” can control the accuracy, this will exponentially increase the computational cost of simulations. The natural spawning criteria are:

|

∫

|cn (t ̸= 0)| ≤ cmin ,

(22)

|Amn (¯ x) · v ¯| ≤ Avmin ,

(23)

∗ dx gm (x)ei∆pmn ·(x−¯x) gn (x)| ≤ Omin .

(24)

The first criterion results from the divergence in Eqs. (17) and (20) when |cn |2 = 0. However, “spawning” at this point does not increase the number of wavepackets as there is no need to propagate a wavepacket with no amplitude. The second criterion comes from the divergence of ∆pmn when Amn (¯ x) · v ¯ = 0. The phase of the coupling term, Eq. (9), changes rapidly at this point, quickly invalidating approximation (14). The third criterion is a direct measure of approximation (14). While taking Omin close to one will negate the need for the other criterions, it can create an undesirably high “spawning” rate. 84 In practice, for the following simulation results, we only require the second criterion: The spawning occurs only when

12



Amn (¯ x) · v ¯ drops below a threshold (or changes sign in one timetep). The stability and reliability of the algorithm can also be improved by limiting non-adiabatic dynamics to regions where |Amn (¯ x) · ∆p| >> |Amn (¯ x) · p ¯ | is not true, which is consistent with the Massy criterion for nonadiabatic transitions. 24

Numerical Simulation Results E2

A12/100

E1

%

60

60

a) 20

20

2

Difference from exact (%)

1

A12

Refl. 2

E2

Refl. 1

E1

b)

30

k (a.u.) Trans. 1

% Trans. 1

30

10

Trans. 2

40

E1

% Trans. 1

30

A12/15

E2

90

90

20

30

c) 40

k (a.u.)

%

15

Trans. 1

30

45

k (a.u.)

%

90

90

V33

V11

60 V12

0 Momentum (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

Page 16 of 30

10

15

20

25

30

30 0

15

0

10

V13

V22

60 E2

-0.04

30

E1

-0.08 -0.12

d)

5

-4

0

4

-4

0 10

12

14

16

k (a.u.)

18

0

e)

4

20

f) 0

30

k (a.u.)

1000

2000

3000

time (a.u.)

Figure 2: a-c) Scattering probabilities for different incoming momentum (k) for Tully models, a) SAC, b) DAC, c) and ECR PESs. EP results shown by red markers, Ehrenfest by dashed blue line, and exact by solid/short-dashed/dash-dotted black line. Insets show adiabatic PES for state 1 (red dot) and 2 (Blue dashed) with scaled NACV (black solid). d) Detailed SAC results: Upper panel: difference between the exact scattering probability on PES 1, definitions same as in panel a. Lower panel: outgoing momenta, exact on PES [1/2] (black [solid/dashed] line ), and EP on surface [1/2] (red [squares/diamonds]), and Ehrenfest (blue dashed line) vs k. e) Scattering Results for two-dimensional PES, definitions same as in panel a with additional SCMC result shown by green long-dashed line. Inset shows PESs. f) Time dependent populations for diabatic states 1 (thick lines) and 3 (thin lines) for threelevel photodissociation model, exact shown as black solid, Ehrenfest as blue short-dashed, EP as red long-dashed.

13


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


We compare our EP results to the exact Schrödinger equation and Ehrenfest method for five model problems, which are standard tests of electron-nuclear correlation in non-adiabatic transitions, 24 see Fig. 2. The details of these scattering potentials and the initial conditions of the wavepacket are given in the Supporting Information. The first model is a single avoided crossing (SAC) that does not require any branching. Here, Ehrenfest (blue dashed line) and the Ehrenfest-Plus (red circles) both quantitatively agree with the exact solution for scattering probabilities, see Fig. 2-a, with a slight increase in accuracy for EP, see upper panel of Fig. 2-d. However the EP correctly calculates the outgoing momenta of the wavepackets on the two surfaces, while the Ehrenfest method results in an “average” momentum, see lower panel of Fig. 2-d. For the second model, a double avoided crossing (DAC), quantum interference between two pathways (crossing x = 0 on either surface 1 or 2) leads to “Stueckelberg” oscillations in the scattering probability, see Fig. 2-b. By neglecting difference in forces on PES 1 and 2, the Ehrenfest method (blue dashed line) shifts the phase of the oscillations at low k. This interference is correctly captured, or significantly improved, by the EP method with a single branching at the point where the NACV changes sign (red circles), or no branching (green square), respectively. This demonstrates that the equations of motion, derived here from first principles, naturally account for a “phase correction”. However, at least a single branch is required to recover the qualitative structure of the outgoing wavepackets, as depicted in Ref. 86. Full multiple spawning (FMS) can recover these oscillations as well but, according to Ref. 90, ten basis functions are required. The third model, Fig. 2-c, includes an extended coupling region and a reflection (ECR) on PES 2 for low momentum. At low momentum, after initial crossing of the finite NACV region, the wavepacket on PES 2 will reflect and re-enter the NACV region. This will cause probability to be transfered from state 2 back to state 1. Unlike the DAC model, that pathway (leading to reflection on PES 1) does not interfere with the wavepacket that transmitted on PES 1 after the first crossing. This is accurately represented in the exact

14



solution and EP method. However, by only propagating an “average” momentum, the Ehrenfest method misses the reflection entirely. The EP method results in a single branching at the reflection point, or when passing a threshold on A, depending on initial momentum. In addition to the Tully problems, we have simulated the scattering by a non-separable two-dimensional two-level potential well (Fig. 2-e), originally proposed by Shenvi, Subotnik, and Yang. 91 Like the DAC problem, interference effects are not accurately captured by the Ehrenfest method, at low momenta. The EP produces oscillations with the phase and amplitude in good agreement with the exact TDSE. Additionally we show the results from semiclassical Monte-Carlo (SCMC), which is a stochastic Monte-Carlo integration of the integral probabilities, where branching is allowed at every time point. 86,92 The SCMC requires 75000 trajectories for convergence. The EP results in only a single spawning, at the zero of the NACV, as in the one-dimensional DAC model. This suggests that, for some cases, the EP does not require increased spawning with higher dimensions. Finally, since our target application is relaxation dynamics from high excited states, we have simulated photodissociation dynamics in a three-level one-dimensional models proposed by Coronado, Xing, and Miller. 93 In Fig. 2-e, we show the time dependence of the diabatic state populations, though dynamics are run using the local adiabatic basis for EP and Ehrenfest, for the third model of Ref. 93 . The final scattering probabilities of the EP and TDSE agree to within ∼1%, while Ehrenfest results are significantly shifted after crossing the second region of significant NAC. Additionally the temporal behavior of the Ehrenfest populations shows large fluctuations, suggesting the adiabatic coherences are too large in the region of NAC. This is substantially reduced in the EP method. In the EP simulation, a single branching occurs when leaving the first region of NAC, where V12 ̸= 0. For the EP result we have not sampled initial conditions, thus our transitions are significantly sharper than the exact result. This is expected behavior for any mixed quantum-classical dynamics which starts from a single trajectory or wavepacket. The computational method developed in this paper is based on expansion of the molec-

15


Page 18 of 30

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


ular wavefunction into the basis of “local”, time-dependent, electronic wavefunctions, i.e. dependent on the nuclear positions at a particular time. Therefore, it can be used for numerical simulations of NAMD using “on-fly” calculations of the electronic structure. We envision its application as an efficient and accurate solver in multiple-spawning like computational schemes, 64 being applied to evaluate parameters of the nuclear wavepackets during their passing through a level crossing region. In our numerical studies we have demonstrated that, for realistic parameters, the new approach allows for accurate determination of the parameters of the outgoing wavepackets with the same efficiency as the Ehrenfest method. The Gaussian parameters, i.e. the mean coordinates, mean momenta and the wavepacket widths are determined by numerical solution of new nonadiabtic equations of motion, derived using the momentum jump Hamiltonian formalism developed in this paper. In contrast to the Ehrenfest method, the Ehrenfest plus method accurately describes the electron-nuclear correlation effects. In particular, it properly accounts for the separation of the wavepackets propagating along different PES, which, in its turn, leads to a more accurate description of electronic coherences and occupation numbers. Since the EP method and ab initio multiple spawning both involve the dynamic generation of semi-classical wavepackets, a few words about the differences of the two approaches is beneficial. FMS has been applied to the Tully models in Ref. 75. EP, with no branches, achieves the same level of accuracy as FMS in the SAC and DAC models. Additionally, in EP the electronic structure is only calculated at the mean-position of the coupled wavepackets. Depending on implementation AIMS requires electronic structure at the center of each wavepacket, and possibly the midpoint between them. 66 In EP, after branching occurs the previously coupled wavepackets propagate totally independently. Thus, EP does not fail even if the basis becomes linearly dependent, a significant headache in AIMS implementation. Finally, the non-classical equations of motion, Eqs. 17, help the EP method “self start”, there is no need to propagate backwards in time, as is the case in some AIMS approaches. 66 We note that the effects of quantum tunneling are not accounted for within present

16



method. While one can potentially attempt this by representing the initial states as superpositions and utilizing similar approach to Refs 62,94,95, this might not be feasible for systems with multiple nuclear degree of freedom and therefore is out of scope of this paper. Similarly, the effects of zero-point fluctuations can presumably be, at least partially, accounted for by branching at the turning points of the well.

Acknowledgement This work was supported by the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218NCA000001). We also acknowledge the LANL Institutional Computing (IC) Program provided computational resources.

Supporting Information Available Derivation details, explicit equations of motion, and model potential parameters.

References (1) Doltsinis, N. L.; Marx, D. First Principles Molecular Dynamics Involving Excited States and Nonadiabatic Transitions. J. Theor. Comput. Chem. 2002, 01, 319–349. (2) Curchod, B. F. E.; Mart´ınez, T. J. Ab Initio Nonadiabatic Quantum Molecular Dynamics. Chem. Rev. 2018, 118, 3305–3336. (3) Tavernelli, I. Nonadiabatic Molecular Dynamics Simulations: Synergies between Theory and Experiments. Acc. Chem. Res. 2015, 48, 792–800.

17


Page 20 of 30

Page 21 of 30 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


(4) Sapunar, M.; Ponzi, A.; Chaiwongwattana, S.; Maliˇs, M.; Prlj, A.; Decleva, P.; Doˇslić, N. Timescales of N–H Bond Dissociation in Pyrrole: A Nonadiabatic Dynamics Study. Phys. Chem. Chem. Phys. 2015, 17, 19012–19020. (5) Dawes, R.; Jiang, B.; Guo, H. UV Absorption Spectrum and Photodissociation Channels of the Simplest Criegee Intermediate (CH2OO). J. Am. Chem. Soc. 2015, 137, 50–53. (6) Jonas, D. M. Vibrational and Nonadiabatic Coherence in 2D Electronic Spectroscopy, the Jahn–Teller Effect, and Energy Transfer. Annu. Rev. Phys. Chem. 2018, 69, 327– 352. (7) Fazzi, D.; Barbatti, M.; Thiel, W. Unveiling the Role of Hot Charge-Transfer States in Molecular Aggregates via Nonadiabatic Dynamics. J. Am. Chem. Soc. 2016, 138, 4502–4511. (8) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic ExcitedState Molecular Dynamics: Modeling Photophysics in Organic Conjugated Materials. Acc. Chem. Res. 2014, 47, 1155–1164. (9) Oberhofer, H.; Reuter, K.; Blumberger, J. Charge Transport in Molecular Materials: An Assessment of Computational Methods. Chem. Rev. 2017, 117, 10319–10357. (10) Mitrić, R.; Werner, U.; Bonaˇcić-Kouteck´ y, V. Nonadiabatic Dynamics and Simulation of Time Resolved Photoelectron Spectra Within Time-Dependent Density Functional Theory: Ultrafast Photoswitching in Benzylideneaniline. J. Chem. Phys. 2008, 129, 164118. (11) Fazzi, D.; Barbatti, M.; Thiel, W. Modeling Ultrafast Exciton Deactivation in Oligothiophenes via Nonadiabatic Dynamics. Phys. Chem. Chem. Phys. 2015, 17, 7787–7799.

18



(12) Jankowska, J.; Prezhdo, O. V. Ferroelectric Alignment of Organic Cations Inhibits Nonradiative Electron–Hole Recombination in Hybrid Perovskites: Ab Initio Nonadiabatic Molecular Dynamics. J. Phys. Chem. Lett. 2017, 8, 812–818. (13) Liu, J.; Neukirch, A. J.; Prezhdo, O. V. Non-Radiative Electron–Hole Recombination in Silicon Clusters: Ab Initio Non-Adiabatic Molecular Dynamics. J. Phys. Chem. C 2014, 118, 20702–20709. (14) Correa, A. A.; Kohanoff, J.; Artacho, E.; Sánchez-Portal, D.; Caro, A. Nonadiabatic Forces in Ion-Solid Interactions: The Initial Stages of Radiation Damage. Phys. Rev. Lett. 2012, 108, 213201. (15) Zeb, M. A.; Kohanoff, J.; Sánchez-Portal, D.; Arnau, A.; Juaristi, J. I.; Artacho, E. Electronic Stopping Power in Gold: The Role of d Electrons and the H/He Anomaly. Phys. Rev. Lett. 2012, 108, 225504. (16) Yost, D. C.; Yao, Y.; Kanai, Y. Examining Real-Time Time-Dependent Density Functional Theory Nonequilibrium Simulations for The Calculation of Electronic Stopping Power. Phys. Rev. B 2017, 96, 115134. (17) Kendrick, B. K. Non-Adiabatic Quantum Reactive Scattering in Hyperspherical Coordinates. J. Chem. Phys. 2018, 148, 044116. (18) Xie, C.; Malbon, C.; Yarkony, D. R.; Guo, H. Nonadiabatic Photodissociation Dynamics Of The Hydroxymethyl Radical via The 22A(3S) Rydberg State: A Four-Dimensional Quantum Study. J. Chem. Phys. 2017, 146, 224306. (19) Xie, C.; Ma, J.; Zhu, X.; Yarkony, D. R.; Xie, D.; Guo, H. Nonadiabatic Tunneling in Photodissociation of Phenol. J. Am. Chem. Soc. 2016, 138, 7828–7831. (20) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106. 19


Page 22 of 30

Page 23 of 30 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


(21) Tomotaka, K.; Hiroshi, U.; Koichi, Y. A New Implementation of Ab Initio Ehrenfest Dynamics using Electronic Configuration Basis: Exact Formulation With Molecular Orbital Connection and Effective Propagation Scheme with Locally Quasi-Diabatic Representation. Int. J. Quantum Chem. 2016, 116, 1205–1213. (22) Makhov, D. V.; Symonds, C.; Fernandez-Alberti, S.; Shalashilin, D. V. Ab Initio Quantum Direct Dynamics Simulations of Ultrafast Photochemistry with Multiconfigurational Ehrenfest Approach. Chem. Phys. 2017, 493, 200 – 218. (23) Saita, K.; Shalashilin, D. V. On-The-Fly Ab Initio Molecular Dynamics with Multiconfigurational Ehrenfest Method. J. Chem. Phys. 2012, 137, 22A506. (24) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061–1071. (25) Subotnik, J. E. Fewest-Switches Surface Hopping and Decoherence in Multiple Dimensions. J. Phys. Chem. A 2011, 115, 12083–12096. (26) Subotnik, J. E.; Shenvi, N. A New Approach to Decoherence and Momentum Rescaling in the Surface Hopping Algorithm. J. Chem. Phys. 2011, 134, 024105. (27) Subotnik, J. E.; Shenvi, N. Decoherence and Surface Hopping: When Can Averaging Over Initial Conditions Help Capture the Effects of Wave Packet Separation? J. Chem. Phys. 2011, 134, 244114. (28) Subotnik, J. E. Augmented Ehrenfest Dynamics Yields a Rate for Surface Hopping. J. Chem. Phys. 2010, 132, 134112. (29) Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N. Understanding the Surface Hopping View of Electronic Transitions and Decoherence. Annu. Rev. Phys. Chem. 2016, 67, 387–417.

20



(30) Shenvi, N.; Yang, W. Achieving Partial Decoherence in Surface Hopping Through Phase Correction. J. Chem. Phys. 2012, 137, 22A528. (31) Shenvi, N.; Subotnik, J. E.; Yang, W. Simultaneous-Trajectory Surface Hopping: A Parameter-Free Algorithm for Implementing Decoherence in Nonadiabatic Dynamics. J. Chem. Phys. 2011, 134, 144102. (32) Bittner, E. R.; Rossky, P. J. Quantum Decoherence in Mixed Quantum-Classical Systems: Nonadiabatic Processes. J. Chem. Phys. 1995, 103, 8130–8143. (33) Zhu, C.; Jasper, A. W.; Truhlar, D. G. Non-Born Oppenheimer Trajectories with SelfConsistent Decay of Mixing. J. Chem. Phys. 2004, 120, 5543–5557. (34) Zhu, C.; Nangia, S.; Jasper, A. W.; Truhlar, D. G. Coherent Switching with Decay of Mixing:

An Improved Treatment of Electronic Coherence for Non-

¨ ıOppenheimer Trajectories. J. Chem. Phys. 2004, 121, 7658–7670. Born?A` (35) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic ExcitedState Molecular Dynamics: Treatment Of Electronic Decoherence. J. Chem. Phys. 2013, 138, 224111. (36) Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum–Classical Dynamics. Chem. Rev. 2018, 118, 7026–7068. (37) Akimov, A. V.; Prezhdo, O. V. Second-Quantized Surface Hopping. Phys. Rev. Lett. 2014, 113, 153003. (38) Wang, L.; Sifain, A. E.; Prezhdo, O. V. Fewest Switches Surface Hopping in Liouville Space. J. Phys. Chem. Lett. 2015, 6, 3827–3833. (39) Wang, L.; Trivedi, D.; Prezhdo, O. V. Global Flux Surface Hopping Approach for Mixed Quantum-Classical Dynamics. J. Chem. Theory Comput. 2014, 10, 3598–3605.

21


Page 24 of 30

Page 25 of 30 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


(40) Prezhdo, O. V.; Rossky, P. J. Mean-Field Molecular Dynamics with Surface Hopping. J. Chem. Phys. 1997, 107, 825–834. (41) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-Induced Surface Hopping. J. Chem. Phys. 2012, 137, 22A545. (42) Martens, C. C. Surface Hopping by Consensus. J. Phys. Chem. Lett. 2016, 7, 2610– 2615. (43) 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. (44) Menzeleev, A. R.; Bell, F.; Miller, T. F. Kinetically Constrained Ring-Polymer Molecular Dynamics for Non-Adiabatic Chemical Reactions. J. Chem. Phys. 2014, 140, 064103. (45) Boucher, W.; Traschen, J. Semiclassical Physics and Quantum Fluctuations. Phys. Rev. D 1988, 37, 3522–3532. (46) Doll, J.; Freeman, D.; Gillan, M. Stationary Phase Monte Carlo Methods: An Exact Formulation. Chem. Phys. Lett. 1988, 143, 277 – 283. (47) Herman, M. F. Nonadiabatic Semiclassical Scattering. I. Analysis of Generalized Surface Hopping Procedures. J. Chem. Phys. 1984, 81, 754–763. (48) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919–8929. (49) Kapral, R. Progress in the Theory of Mixed Quantum-Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129–157.

22



(50) Mac Kernan, D.; Ciccotti, G.; Kapral, R. Trotter-Based Simulation of QuantumClassical Dynamics. J. Phys. Chem. B 2008, 112, 424–432. (51) Makri, N. Time-Dependent Quantum Methods for Large Systems. Annu. Rev. of Phys. Chem. 1999, 50, 167–191. (52) Martens, C. C.; Fang, J.-Y. Semiclassical-Limit Molecular Dynamics on Multiple Electronic Surfaces. J. Chem. Phys. 1997, 106, 4918–4930. (53) Meyer, H.; Miller, W. H. A Classical Analog for Electronic Degrees of Freedom in Nonadiabatic Collision Processes. J. Chem. Phys. 1979, 70, 3214–3223. (54) Meyer, H.; Miller, W. H. Analysis and Extension of Some Recently Proposed Classical Models for Electronic Degrees of Freedom. J. Chem. Phys. 1980, 72, 2272–2281. (55) Miller, W. H. Classical-Limit Quantum Mechanics and the Theory of Molecular Collisions; John Wiley and Sons, 2007. (56) Pechukas, P. Time-Dependent Semiclassical Scattering Theory. I. Potential Scattering. Phys. Rev. 1969, 181, 166–174. (57) Pechukas, P. Time-Dependent Semiclassical Scattering Theory. II. Atomic Collisions. Phys. Rev. 1969, 181, 174–185. (58) Stock, G.; Thoss, M. Semiclassical Description of Nonadiabatic Quantum Dynamics. Phys. Rev. Lett. 1997, 78, 578–581. (59) 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. (60) Makri, N. Quantum-Classical Path Integral: A Rigorous Approach to Condensed Phase Dynamics. Int. J. Quantum Chem. 2015, 115, 1209–1214. 23


Page 26 of 30

Page 27 of 30 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


(61) Martens, C. C. Fully Coherent Quantum State Hopping. J. Chem. Phys. 2015, 143, 141101. (62) Chen, X.; Batista, V. S. Matching-Pursuit/Split-Operator-Fourier-Transform Simulations of Excited-State Nonadiabatic Quantum Dynamics in Pyrazine. J. Chem. Phys. 2006, 125, 124313. (63) Coalson, R. D. A Wavepacket Path Integral Method for Curve-Crossing Dynamics. J. Phys. Chem. 1996, 100, 7896–7902. (64) Ben-Nun, M.; Quenneville, J.; Mart´ınez, T. J. Ab Initio Multiple Spawning: Photochemistry from First Principles Quantum Molecular Dynamics. J. Phys. Chem. A 2000, 104, 5161–5175. (65) Grossmann, F. Semiclassical Real-Time Tunneling by Multiple Spawning of Classical Trajectories. Phys. Rev. Lett. 2000, 85, 903–907. (66) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of Ab Initio Multiple Spawning in the Molpro Quantum Chemistry Package. Chem. Phys. 2008, 347, 3 – 16. (67) Martinez, T. J. Insights for Light-Driven Molecular Devices from Ab Initio Multiple Spawning Excited-State Dynamics of Organic and Biological Chromophores. Acc. Chem. Res. 2006, 39, 119–126. (68) 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. (69) 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. 24



(70) Liu, L.; Liu, J.; Martinez, T. J. Dynamical Correlation Effects on Photoisomerization: Ab Initio Multiple Spawning Dynamics with MS-CASPT2 for a Model trans-Protonated Schiff Base. J. Phys. Chem. B 2016, 120, 1940–1949. (71) Mori, T.; Glover, W. J.; Schuurman, M. S.; Martinez, T. J. Role of Rydberg States in the Photochemical Dynamics of Ethylene. J. Phys. Chem. A 2012, 116, 2808–2818. (72) Mignolet, B.; Curchod, B. F. E. A Walk Through the Approximations of Ab Initio Multiple Spawning. J. Chem. Phys. 2018, 148, 134110. (73) Thompson, A. L.; Punwong, C.; Mart´ınez, T. J. Optimization of Width Parameters for Quantum Dynamics with Frozen Gaussian Basis Sets. Chem. Phys. 2010, 370, 70 – 77. (74) Yang, S.; Coe, J. D.; Kaduk, B.; Mart´ınez, T. J. An “Optimal” Spawning Algorithm for Adaptive Basis Set Expansion in Nonadiabatic Dynamics. J. Chem. Phys. 2009, 130, 134113. (75) Martinez, T. J.; Ben-Nun, M.; Levine, R. D. Multi-Electronic-State Molecular Dynamics: A Wave Function Approach with Applications. J. Phys. Chem. 1996, 100, 7884–7895. (76) Barbatti, M. Nonadiabatic Dynamics with Trajectory Surface Hopping Method. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 620–633. (77) Baer, M. Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections; John Wiley & Sons, 2006. (78) Gherib, R.; Ye, L.; Ryabinkin, I. G.; Izmaylov, A. F. On the Inclusion of the Diagonal Born-Oppenheimer Correction in Surface Hopping Methods. J. Chem. Phys. 2016, 144, 154103. (79) Child, M. S.; Shalashilin, D. V. Locally Coupled Coherent States and Herman–Kluk Dynamics. J. Chem. Phys. 2003, 118, 2061–2071. 25


Page 28 of 30

Page 29 of 30 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


(80) Kapral, R. Progress in the Theory of Mixed Quantum-Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129–157. (81) Ando, K.; Santer, M. Mixed Quantum-Classical Liouville Molecular Dynamics Without Momentum Jump. J. Chem. Phys. 2003, 118, 10399–10406. (82) Heller, E. J. Time?Dependent Approach to Semiclassical Dynamics. J. Chem. Phys. 1975, 62, 1544–1555. (83) Richings, G.; Polyak, I.; Spinlove, K.; Worth, G.; Burghardt, I.; Lasorne, B. Quantum Dynamics Simulations Using Gaussian Wavepackets: The VMCG Method. Int. Rev. Phys. Chem. 2015, 34, 269–308. (84) 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. (85) White, A. J.; Gorshkov, V. N.; Tretiak, S.; Mozyrsky, D. Non-Adiabatic Molecular Dynamics by Accelerated Semiclassical Monte Carlo. J. Chem. Phys. 2015, 143, 014115. (86) White, A. J.; Gorshkov, V. N.; Wang, R.; Tretiak, S.; Mozyrsky, D. Semiclassical Monte Carlo: A First Principles Approach to Non-Adiabatic Molecular Dynamics. J. Chem. Phys. 2014, 141, 184101. (87) Meek, G. A.; Levine, B. G. The Best of Both Reps—Diabatized Gaussians on Adiabatic Surfaces. J. Chem. Phys. 2016, 145, 184103. (88) Vacher, M.; Bearpark, M. J.; Robb, M. A. Direct Methods for Non-Adiabatic Dynamics: Connecting the Single-Set Variational Multi-Configuration Gaussian (VMCG) and Ehrenfest Perspectives. Theo. Chem. Acc. 2016, 135, 187. (89) Makhov, D. V.; Martinez, T. J.; Shalashilin, D. V. Toward Fully Quantum Modelling of

26



Ultrafast Photodissociation Imaging Experiments. Treating Tunnelling in the Ab Initio Multiple Cloning Approach. Faraday Discuss. 2016, 194, 81–94. (90) Ben-Nun, M.; Mart´ınez, T. J. Ab Initio Quantum Molecular Dynamics; John Wiley and Sons, 2002. (91) Shenvi, N.; Subotnik, J. E.; Yang, W. Phase-Corrected Surface Hopping: Correcting the Phase Evolution Of The Electronic Wavefunction. J. Chem. Phys. 2011, 135, 024101. (92) Gorshkov, V. N.; Tretiak, S.; Mozyrsky, D. Semiclassical Monte-Carlo Approach for Modelling Non-Adiabatic Dynamics in Extended Molecules. Nat. Comm. 2013, 4, 2144. (93) 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. (94) Kong, X.; Markmann, A.; Batista, V. S. Time-Sliced Thawed Gaussian Propagation Method for Simulations of Quantum Dynamics. J. Phys. Chem. A 2016, 120, 3260– 3269. (95) Wu, Y.; Batista, V. S. Matching-Pursuit for Simulations of Quantum Processes. J. Chem. Phys. 2003, 118, 6720–6724.

27


Page 30 of 30

An Improved Ehrenfest Approach to Model Correlated Electron

Recommend Documents