Ab Initio Multiple Spawning Method for Intersystem ... - ACS Publications

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Ab Initio Multiple Spawning Method for Intersystem Crossing Dynamics: Spin-Forbidden Transitions between B and A States of GeH 3

1

1

1

2

Dmitry A Fedorov, Spencer R. Pruitt, Kristopher Keipert, Mark S. Gordon, and Sergey A Varganov J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b01406 • Publication Date (Web): 11 Apr 2016 Downloaded from http://pubs.acs.org on April 16, 2016

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

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

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

Ab Initio Multiple Spawning Method for Intersystem Crossing Dynamics: Spin-Forbidden Transitions between 3B1 and 1A1 States of GeH2

Dmitry A. Fedorova, Spencer R. Pruittb, Kristopher Keipertc, Mark S. Gordonc and Sergey A. Varganova*

a

Department of Chemistry, University of Nevada, Reno, 1664 N. Virginia St., Reno, NV

89557-0216, United States b

Argonne Leadership Computing Facility, Argonne National Laboratory, 9700 S. Cass

Avenue, Lemont, IL 60439, USA c

Ames Laboratory, Department of Chemistry, Iowa State University, Ames, IA 50010,

United States

*E-mail: [email protected] Phone: +1 (775) 784-1406

1

ACS Paragon Plus Environment


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

Abstract Dynamics at intersystem crossings are fundamental to many processes in chemistry, physics and biology. The ab initio multiple spawning (AIMS) method was originally developed to describe internal conversion dynamics at conical intersections where derivative coupling is responsible for nonadiabatic transitions between electronic states with the same spin multiplicity. Here, the applicability of the AIMS method is extended to intersystem crossing dynamics in which transitions between electronic states with different spin multiplicities are mediated by relativistic spin-orbit coupling. In the direct AIMS dynamics, the nuclear wave function is expanded in the basis of frozen multidimensional Gaussians propagating on the coupled electronic potential energy surfaces calculated on-the-fly. The AIMS method for intersystem crossing is used to describe the nonadiabatic transitions between the 3B1 and 1A1 states of GeH2. The potential energies and gradients were obtained at the CASSCF(6,6)/6-31G(d) level of theory. The spin-orbit coupling matrix elements were calculated with the configuration interaction method using the two-electron Breit-Pauli Hamiltonian. The excited 3B1 state lifetime and intersystem crossing rate constants were estimated by fitting the AIMS state population with the first-order kinetics equation for a reversible unimolecular reaction. The obtained rate constants are compared with the values predicted by the statistical nonadiabatic transition state theory with transition probabilities calculated using the Landau-Zener and weak coupling formulas.

2


Page 2 of 35

Page 3 of 35

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


I. Introduction Intersystem crossings (ISC), radiationless transitions between electronic states with different spin multiplicities, play an important role in many chemical, physical and biological processes. These transitions happen mostly due to the spin-orbit coupling (SOC) between electronic states, a relativistic effect arising from the interaction between the orbital angular momentum and spin of the electron. Since the SOC increases with nuclear charge, historically, ISC processes were believed to be important only in systems that contain heavy atoms. In organic molecules, the ISC rates were assumed to be much lower than the rates of internal conversion (IC) and radiative processes. However, the rate of ISC depends not only on the magnitude of the SOC, but also on the shapes of the potential energy surfaces (PESs) and the accessibility of the intersection between PESs. With advances in ultrafast sub-picosecond spectroscopy it has been shown that ISC can compete with IC and radiative processes even in organic molecules with relatively small SOCs.1-4 Therefore, fast ISC, along with spontaneous emission and internal conversion, plays an important role in the photochemistry and photophysics of many different molecular systems.5-12 The ISC rates can be calculated using time-independent (statistical) or timedependent (molecular dynamics) approaches. In the time-independent approach, it is assumed that the ISC occurs at a minimum energy crossing point (MECP), the lowest energy point on the crossing seam of two spin-diabatic PESs. The simplest way to calculate the probability of transition between two states at the MECP is by using the Landau-Zener (LZ) formula.13-15 Since the LZ formula is derived assuming the classical movement of nuclei, it predicts a nonzero probability of transitions between spin-diabatic

3



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

PESs only if the internal energy of a system is higher than the MECP energy. Hence, the LZ formula does not account for quantum effects such as tunneling and interference between different nuclear trajectories. In contrast, the weak coupling (WC) formula16,17 takes into account these quantum effects and predicts nonzero transition probabilities even for internal energies below the MECP energy. Knowing the transition probabilities calculated with the LZ or WC formulas and the density of rovibrational states for the reactants and MECP, the ISC rate constant can be calculated using the nonadiabatic transition state theory (NA-TST).18,19 The time-dependent approach to the calculation of the ISC rates is based on the sampling of different reaction pathways through molecular dynamics simulations. Quantum and semiclassical simulations, in which the time-dependent Schrödinger equation for nuclear motion is solved on pre-computed PESs, are restricted to relatively small systems and short-time nuclear dynamics. Examples of ISC dynamics on precomputed PESs include the studies of ISC in sulfur dioxide20 and the collision of atomic oxygen with the deuterium molecule.21 Alternatively, small parts of the electronic PESs along the nuclear trajectories can be calculated on the fly, as in direct molecular dynamics. In Ehrenfest mean-field molecular dynamics,22,23 nuclear trajectories evolve on an average PES corresponding to multiple electronic states. The Ehrenfest dynamics can be justified for regions of strong coupling, where the electronic wave function can be accurately represented as a mixture of state specific wave functions. However, due to a lack of decoherence, Ehrenfest dynamics predicts unphysical behavior in the asymptotic regions (e.g. dissociation, bond breaking), where the wave function splits into multiple components localized in different regions of space. Another quantum-classical approach

4


Page 4 of 35

Page 5 of 35

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


to direct dynamics is the trajectory surface hopping (TSH) method24 with a fewestswitches criterion.25 The original TSH method was formulated for IC dynamics at a conical intersection and later was adapted for ISC dynamics.4,12,26-31 In the TSH method, the nuclear wave packet is represented by a set of independent classical trajectories. The changes of electronic states (hops between PESs) occur stochastically, with the probability of hopping determined based on the coupling between the electronic states and state populations. Since for most of the time nuclear trajectories propagate independently, quantum coherence effects are not accounted for. An example of quantum coherence effects is a frustrated hop, in which the energy of a trajectory cannot be adjusted to conserve the classical energy. Frustrated hops are classically forbidden, however, they are quantum mechanically allowed. Ignoring the frustrated hops leads to decoherence. To resolve this issue a number of improvements to the fewest-switches TSH method was introduced, such as the time-uncertainty32 and the gradV33 methods. However, the major drawback of both TSH and Ehrenfest methods is that they do not provide a scheme to systematically improve the description of quantum coherence effects. Therefore, in general, the convergence of these methods to the exact solution of the Schrödinger equation is not guaranteed.34,35 In the ab initio multiple spawning (AIMS) method,36 the total wave function is written as a product of the electronic and nuclear wave functions. The nuclear wave function is expanded in a basis of frozen multidimensional Gaussians. The Gaussian basis functions are used due to their well-known connections to classical mechanics.37 Each Gaussian basis function propagates on a single electronic state. The positions, momenta and phases of the Gaussians are propagated using classical equations of motion. The

5



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

expansion coefficients (complex amplitudes) are propagated by solving the timedependent Schrödinger equation, thereby accounting for quantum coherence effects in a natural way. The time-dependent Gaussian basis set is expanded when needed by adding new basis functions on neighboring electronic states. This basis set expansion, called spawning, happens in the proximity of the intersections between electronic states, if the coupling between the states is large. The advantage of the AIMS method over TSH is that AIMS converges to the exact solution of the Schrödinger equation with increasing numbers of basis functions.38,39 The AIMS method has been used to study IC dynamics at conical intersections and has been interfaced with several electronic structure codes,40-43 including MOLPRO44,45 and GAMESS.46,47 The current work extends the AIMS method to ISC dynamics, in which electronic states with different spin multiplicities are coupled through spin-orbit interactions. In addition, a direct interface between the AIMS and GAMESS programs is developed. The ISC between the first excited 3B1 state and the ground 1A1 state of GeH2 is used to demonstrate how the ISC rates and electronic state lifetimes can be calculated. The ISC rate constants and excited state lifetime calculated from the AIMS dynamics are compared with the corresponding values obtained using the statistical NA-TST where transition probabilities between electronic states were computed using the LZ and WC formulas. While this paper was under review a similar implementation of the AIMS method for ISC dynamics based on the AIMS-MOLPRO interface was reported.48

6


Page 6 of 35

Page 7 of 35

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


II. Theoretical Methods A. Nonadiabatic Transition State Theory According to the NA-TST, for a unimolecular ISC reaction with internal energy E, the rate constant can be expressed as18,19,49,50

( )

k E =

E

1

ρ ( )∫

hρ R E

MECP

(E − ε ⊥ )Ptrans (ε ⊥ )d ε ⊥ .

(1)

0

Here h is the Plank’s constant and ε ⊥ is the energy associated with the reaction coordinate orthogonal to the crossing seam between two spin-diabatic PESs. The density of rovibrational states at the MECP, ρ MECP (E − ε ⊥ ) , is convoluted with the probability of transition between two PESs, Ptrans (ε ⊥ ) , and divided by the density of rovibrational states

( )

of the reactants, ρ R E . We use the double passage version of the LZ formula,17,51 where the probability of transition between two spin-diabatic states is expressed as

( ) (

)(

)

Ptrans ε ⊥ = 1− P LZ 1+ P LZ ,

(2)

(3)

In (3), H SO and ∆G = G1 - G2 are the spin-orbit coupling and the difference of the gradients of two PESs at the MECP, respectively; µ ⊥ is the reduced mass of the degree of freedom orthogonal to the crossing seam. The MECP energy, EMECP, is defined with respect to the energy of the reactants. The LZ formula does not account for quantum tunneling and interference between different reaction pathways. However, the following WC double passage formula does take these quantum effects into account.16,17,52 7



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

Page 8 of 35

(4)

12

Here G = ( G1G 2 ) is the geometric mean of the PES gradients at the MECP. The Airy function, which can be defined as

Ai(x) =

1

∞

∫ cos (t π

3

3+ xt ) dt,

(5)

0

is responsible for the nonzero transition probability due to tunneling at energies below EMECP, as well as for the oscillations from the interference of the forward and backward reaction pathways at energies above EMECP. The density of rovibrational states at the MECP, ρ MECP (E − ε ⊥ ) , can be obtained from the harmonic oscillator and classical rotor models. Since the MECP is a stationary point on the 3N-7 dimensional crossing seam (N is the number of nuclei), but not on any PES, the effective Hessian defined in terms of the PES gradients and hessians (H1, H2) is introduced,19

H′eff =

G1 H 2 ± G 2 H1 ∆G

.

(6)

In (6), the plus and minus signs are used for peaked ( G1 ⋅ G 2 < 0 ) and sloped ( G1 ⋅ G 2 > 0 ) intersections, respectively. The reaction coordinate orthogonal to the seam is projected out of H′eff leaving the reduced-dimension Hessian19,53

 ∆G   ∆G  H eff =  1 H ′  eff  1 ,  ∆G   ∆G 

(7)

where 1 is the unit matrix. The effective Hessian (7) is diagonalized to obtain the 3 N-7 harmonic frequencies of the crossing seam. These frequencies are used to calculate the

8


Page 9 of 35

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


density of vibrational states at the MECP. The rovibrational densities of states at the MECP and at the reactants were obtained by convoluting the densities of the harmonic vibrational states calculated though direct counting54 with the densities of rotational states from the asymmetric top model.55

B. Ab Initio Multiple Spawning for ISC Dynamics The full description of the AIMS method originally developed to describe IC dynamics can be found elsewhere.36,38,56,57 Here, only the essential equations are presented, focusing on the details relevant to ISC dynamics. In the AIMS method, the total wave function is expanded in terms of the products of electronic and nuclear wave functions:

ψ = ∑ χ I ( R;t )φ I ( r; R) .

(8)

I

In (8), χ I and φ I are, respectively, the nuclear and electronic wave functions associated with the electronic state I. The nuclear wave function is a superposition of frozen Gaussian basis functions with time-dependent coefficients: N I (t )

χ ( R;t ) = ∑ CkI ( t ) χ kI ( R; R kI ( t ), P kI ( t ) , γ kI ( t ), α ρI ) , I

(9)

k

k

where k labels the nuclear basis functions and NI(t) is the number of nuclear basis functions on the electronic state I. The frozen Gaussian basis functions χ kI are parameterized with the time-independent width α kI , and the time-dependent average

()

()

()

position R kI t , momentum P kI t and phase γ kI t . The multidimensional Gaussians

χ kI are products of one-dimensional Gaussian functions

9



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

χ

I k

χ

( R; R ( t ), P (t ),γ (t ),α ) = e I k

I

ρk

I k

 2α ρk I  =   π 

I

I k

1

4

iγ k I ( t )t

k

Page 10 of 35

3N

χ ρ ( R; R ρ ( t ) , P ρ ( t ) , α ρ ) , ∏ ρ I

=1

(

I

k

)

I

k

I

k

(

k

)

2 exp  −α ρk Rρk − Rρk ( t ) + iP ρI k ( t ) Rρk − Rρk ( t )  ,  

(10)

(11)

where ρk labels the 3N Cartesian coordinates of the nuclei for the basis function

()

()

(trajectory) k. The positions, R kI t , and momenta P kI t are propagated according to the classical Hamilton equations; the phase γ kI is propagated using the semi-classical Lagrangian,

∂Rρk I ∂t = P ρI k M ρ ,

(12)

∂P ρI k ( t ) ∂t = −  ∂EI ( R ) ∂R  R(t ) , ρk

(13)

3N

(

)

2 ∂γ k ( t ) ∂t = ∑  P ρI k ( t ) − 2α ρk I  2M ρ − EI ( Rk I ( t )) .   ρ =1 I

(14)

Here Mρ is the mass of the nuclei associated with the ρ-th coordinate and EI is the potential energy of the electronic state I. The time-independent Gaussian width, α ρk , depends only on the type of the nucleus and does not depend on either the basis function (trajectory) index k or the electronic state index I.58 The set of coupled differential equations, obtained by substituting (6) into the timedependent Schrödinger equation, is solved to find the complex amplitudes CkI , (15) Eq 15 is written in matrix notation for clarity. Since the nuclear basis functions χ kI are time-dependent and nonorthogonal, propagation of the amplitudes requires the overlap

10


Page 11 of 35

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


( )

matrix elements S II

kl

and the time derivative of the overlap matrix elements

. In

this notation, the electronic state indices I and J label different matrices, while the trajectory indices k and l label elements of the matrices.

( SII ) kl = χ kI χ lI

,

(16) (17)

In (16) and (17), the integration is performed only over the nuclear coordinates, since the overlap of the electronic wave functions, φ I φ J , is zero for I≠J and unity for I=J. In general, the full Hamiltonian matrix H IJ is the sum of the electronic Hamiltonian and the nuclear kinetic energy:

( H IJ ) kl = χ kI φ I

Hˆ e + TˆR φ J χ lJ .

(18)

The integration in (18) is performed over both electronic and nuclear degrees of freedom. The diagonal matrix element, H II , is the sum of the kinetic and potential energies of the electronic state I. In the spin-diabatic representation of two electronic states I and J with different spin multiplicities, the off-diagonal matrix elements of the kinetic energy,

χ kI φ I TˆR φ J χ lJ , are zero due to orthogonality of the spin parts of the electronic wave functions. Thus the off-diagonal elements of the Hamiltonian (18) can be written as:

( H IJ ) kl = χ kI φI

Hˆ e φ J χ lJ = χ kI φ I Hˆ SO φ J χ lJ ,

(19)

where Hˆ SO is the spin-orbit operator. Here we use the Breit-Pauli spin-orbit operator with one- and two-electron terms59 that can be expressed in atomic units as

11



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

Page 12 of 35

Zα ˆ 1  e2 ˆ ˆ H SO = 2  ∑ ∑ 3 l iα ⋅ sˆ i − ∑ ∑ 3 l ij ⋅ sˆ i + 2sˆ j 2c  i α r iα i j r ij

(



) .

(20)



Roman and Greek subscripts refer to electrons and nuclei, respectively. Space and spin angular momentum operators are defined in terms of the electronic momentum, p, and

(

)

the distances between the particles, ˆliα = ( ri − Rα ) × pˆ i , ˆl ij = ri - r j × pˆ i . The real positive spin-orbit coupling (SOC) matrix element between two electronic states with spins S and S´ averaged over magnetic quantum numbers M S and M S′ is calculated as

φ I Hˆ SO φ J

2

=

S

S′

∑ ∑

M S =− S M S′ =− S ′

2

SM S Hˆ SO S ′M S′ .

(21)

Integration over the nuclear degrees of freedom in eq 19 can be done using the first-order saddle-point approximation,60

χ kI φ I Hˆ SO φ J χ lJ ≈ φ I Hˆ SO ( RC ) φ J

χ kI χ lJ ,

(22)

where RC is the location of the centroid of the product of the nuclear Gaussian basis functions χ kI and χ lJ . In the AIMS method, the nuclear motion on several PESs is described with an adaptive basis set for improved computational efficiency. At each time step, the SOC matrix elements are calculated for each nuclear basis function (trajectory). If the SOC between electronic states I and J is large, then a population transfer from the basis function χ kI propagating on the electronic state I is likely. To facilitate this population transfer a new basis function, χ lJ , is added (spawned) on the state J, if the dimensionless effective coupling parameter

12


Page 13 of 35

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


Λ eff =

φ I Hˆ SO ( RC ) φ J

χ kI χ lJ

EI − EJ

(23)

is larger than the predefined threshold.61 In (23) EI and EJ are the potential energies of the electronic states I and J. Another parameter that controls spawning is the overlap between the initial (parent) and the spawned (child) basis functions. If the overlap is smaller than a certain threshold, a new basis function is not spawned even if the effective coupling is large. It is important to note that new basis functions are created with no initial population. The population transfer between different states occurs during the propagation of the basis function coefficients, CkI , according to eq 15. An important distinction between IC dynamics at a conical intersection and ISC dynamics is the dependence of the coupling magnitude on the molecular geometry. In IC dynamics, the vector derivative coupling between two adiabatic states with the same spin peaks at the conical intersection and decreases rapidly with distance from the conical intersection (however, if the diabatic basis is used, the derivative coupling is often smooth). In contrast, the scalar SOC between spin-diabatic states usually changes slowly with molecular geometry. Coupling (22) decreases significantly only if the overlap between nuclear basis functions, χ kI χ lJ , is small. As a result, this coupling can be non-negligible even at relatively large distances from the crossing seam. Hence, to fully capture the ISC population transfer, new nuclear basis functions have to be spawned farther away from the crossing seam than in IC dynamics. Convergence of the population transfer rate with respect to the effective coupling parameter Λ eff is shown in Supporting Information (Figure S1). Population transfer between singlet and triplet electronic states can be described by 13



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

Page 14 of 35

the first-order kinetics equation for a reversible unimolecular reaction: (24) where F3 and F1 are the populations of the triplet and singlet states; k3 and k1 are the ISC rate constants. The solution of the forward rate equation for (24) can be written as:62

 t − t0  F3 t = F3 t0 − k1τ exp  − + k1τ ,  τ 

() ( ( )

)

(25)

where t0 is the initial time, F3(t0)=1 is the initial population of the triplet states, and

τ = 1 ( k1 + k3 ) is the triplet state lifetime. The individual state population is a sum of the populations from each trajectory propagated on this state. NI

FI = ∑ FkI ,

(26)

k =1

In (26), FkI is the population of the individual basis function k on the electronic state I, NI

( )

F = ∑ CkI* S II kiCiI . I k

(27)

i=1

By fitting the triplet state population with eq 25 the rate constants k1, k3 and the initial time t0 can be obtained.

III. Computational Details The AIMS code is responsible for the propagation of the nuclear basis functions and their complex amplitudes. At each time step, the AIMS code provides a new molecular geometry and requests the electronic structure properties from the GAMESS suite of programs. These properties include: a) state-specific energies and gradients for

14


Page 15 of 35

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


all electronic states considered in the AIMS dynamics, b) vectors of molecular orbital and configuration interaction (CI) coefficients, and c) spin-orbit coupling between electronic states. The direct interface used in this study was developed independently from a previous interface.42 The potential energy, energy gradient, Hessian and geometry optimization calculations were performed using the state specific complete active space self-consistent field (CASSCF) method63 with the 6-31G(d) basis set. The full valence active space contained 6 electrons and 6 orbitals. Harmonic frequencies at the MECP were computed by diagonalizing the effective Hessian, after projecting out the reaction coordinate orthogonal to the crossing seam. Spin-orbit coupling was computed with the CAS-CI method using the state-averaged CASSCF orbitals and the partial 2-electron method (HSO2P).64 The AIMS molecular dynamics simulations on the lowest energy singlet 1A1 and triplet 3B1 electronic states of GeH2 were done without any symmetry restrictions. Integration of the classical equations of motion was carried out using the velocity Verlet algorithm with a variable time step. To avoid any numerical issues with the integration of the equations of motion, small time steps of 0.1 fs in the region of strong coupling (

Λ eff > 0.1) and 0.2 fs everywhere else were used. For spawning to occur, in addition to Λ eff > 0.1, the overlap between the parent and the child basis functions must be larger than 0.5. The widths of the Gaussian basis functions were 30 bohr-2 and 4.7 bohr-2 for Ge and H atoms, respectively.58 Each AIMS simulation was initialized with a single nuclear basis function propagated on the excited 3B1 state. The initial conditions were generated with a Wigner distribution42,65 using the ground state normal modes and equilibrium

15



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

geometry. The zero-point energy (ZPE) of the ground state was equally distributed between the kinetic and potential energy. A direct interface between the GAMESS quantum chemistry program and the AIMS dynamics program was developed. Rather than rely on a third program to facilitate data transfer through files and scripts, the AIMS program is compiled as a library, linked into a single executable with GAMESS, and called directly. The existing file I/O library within GAMESS is used for shared file access. All calculations are started from a standard GAMESS CAS-CI input file, as well as the standard AIMS input file. Initial vibrational frequencies and coupling between states are calculated with GAMESS. The AIMS program is called directly, and reads the energy and gradient vectors, converged configuration interaction (CI) vectors, MO vectors, transition dipole moments, and SOC between electronic states. All dynamics and spawning are performed with the AIMS program in serial, while all computationally expensive electronic structure quantities are calculated using the full parallelization of the GAMESS program.

16


Page 16 of 35

Page 17 of 35

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


Table 1. Bond length (r), angle (θ), energy relative to the singlet minimum (∆E), harmonic vibrational frequencies (ω) and spin-orbit coupling (SOC) for the 1A1 minimum, 3

B1 minimum and MECP geometries. 1

r, Å θ, deg ∆E, cm

-1

ω, cm-1

SOC, cm-1

3

A1

B1

MECP

1.621

1.564

1.541

92.4

118.7

132.8

0

7484

8201

924

806

1839

1991

1844

2052

374

367

1795 2151

352

IV. Results and Discussion The geometries, relative energies, vibrational frequencies and SOC for the lowest energy singlet (1A1) and triplet (3B1) electronic states and the MECP of GeH2 are reported in Table 1. The contour plots of the two PESs and SOC, φ I Hˆ SO φ J , are shown in Figure 1, with the black line indicating the crossing seam. There is a conical intersection between the 3B1 state and another triplet state in the region where the bond lengths are longer than 1.8 Å and the angles are smaller than 80°. The chosen initial conditions do not allow the AIMS trajectories to reach this conical intersection, which is approximately 24000 cm-1 higher in energy than the 1A1 minimum.

17



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

Figure 1. (a, b) Contour plots of the GeH2 PES cross sections (r = rGe-H1 = rGe-H2) for the 1

A1 ground state and 3B1 excited state, respectively. The black line represents the crossing

seam. Energies are shown relative to the 1A1 minimum. (c) Contour plot of the SOC between the 1A1 and 3B1 states. The white rectangle corresponds to the high-energy region where another triplet state is lower in energy than the 3B1 state. Molecular dynamics trajectories never reach this region.

A total of 48 AIMS simulations were run for 150 fs, each starting with a single nuclear basis function (trajectory) on the excited 3B1 state. On average, 14 trajectories were spawned from the excited 3B1 to the ground 1A1 state and 7 trajectories spawned back to

18


Page 18 of 35

Page 19 of 35

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 3B1 from the 1A1 state during each simulation. Figure 2 shows the distributions of bond lengths and angles at all 1024 spawning events. The bond lengths and angles were grouped into 0.05 Å and 2° bins, respectively. The spawning events occur at bond lengths between 1.3 Å and 1.9 Å. The maximum of the distribution (1.55 Å -1.60 Å) is close to the MECP bond length of 1.54 Å. The distribution of the angles has a maximum at (135°137°), which is also close to the MECP value of 133°. Since the two bond lengths, r1 and r2, are grouped in a single distribution the total number of spawns is two times larger for the bond length distribution than for the angle distribution. In the region where spawning events occur, the SOC matrix element varies only between 300 cm-1 and 380 cm-1 (Figure 1c). Thus the overlap between nuclear basis functions, χ kI χ lJ , and the energy gap between electronic states, EI - EJ, are the major factors responsible for a large effective coupling (eq 23).

Figure 2. Distributions of the bond lengths r1 and r2 (a) and angles (b) for all 1024 spawning events. The MECP bond length and angle are 1.54 Å and 133°, respectively.

19



The population decay from the 3B1 excited state to the 1A1 ground state is plotted in Figure 3. The population of the 3B1 state at each given time was calculated by averaging the populations from all 48 simulations; the values were grouped into 5 fs bins; and a mean value was obtained for each bin. After 150 fs approximately 45% of the population was transferred from the excited 3B1 to the ground 1A1 state. A first-order decay fit with eq 25 predicts the rate constants k3=5.3×1012 s-1 (4.7×1012 s-1, 5.9×1012 s-1) for the 3

B1→1A1 transitions and k1=3.1×1012 s-1 (1.5×1012 s-1, 4.7×1012 s-1) for the 1A1→3B1

transitions. The lower and upper bounds with 95% confidence are reported in parentheses; the R2 value for the fit is 0.988. The corresponding lifetime of the 3B1 state, τ, is 119 fs (94 fs, 164 fs). 1 0.9

3

1

Population of B state

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

0.8 0.7 0.6 0

50

Time, fs

100

150

Figure 3. Population of the 3B1 state as a function of time calculated by averaging the populations from all 48 simulations (black dots). The first-order decay fit using eq 23 is shown with a red line. The rate constants obtained from the fit are k3=5.3×1012 s-1 for the 3

B1→1A1 transitions, and k1=3.1×1012 s-1 for the reverse 1A1→3B1 transitions. The

corresponding 3B1 state lifetime is 119 fs.

The ISC transition probabilities calculated with the LZ and WC formulas, as well as

20


Page 21 of 35

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 NA-TST rate constants are shown in Figure 4. If the energy orthogonal to the crossing seam, ε ⊥ , is lower than the MECP energy (EMECP = 717 cm-1 relative to the 3B1 minimum, and EMECP =8201 cm-1 relative to the 1A1 minimum) then the LZ transition probability is zero (Figure 4a). Since the density of vibrational states at the MECP for the internal energy below the first vibrational level, ZPEMECP=(ω1+ω2)/2=1973 cm-1, is zero, the convolution of the density of states with the transition probability, as in eq 1, produces a threshold below which the rate constant is zero. If the transition probability is calculated using the LZ formula, then the rate constant is nonzero only for the internal energies E > EMECP + ZPEMECP. This corresponds to E > 2688 cm-1 for the 3B1→1A1 transition, and to E > 10173cm-1 for the 1A1→3B1 transition (Figure 4b). In contrast, the rates calculated with the WC probability, which accounts for tunneling, can be nonzero even for internal energies below EMECP + ZPEMECP, as can be clearly seen for the k1 rate constant. Since GeH2 has only three vibrational modes, its density of rovibrational states has prominent peaks at the energies corresponding to the vibrational states at the MECP. These peaks are responsible for the small oscillations in the rate constants close to the MECP energy. The average internal energy of all 48 AIMS trajectories is 5733 cm-1, with respect to the 3B1 minimum, and 13217 cm-1, with respect to the 1A1 minimum. These energies give NA-TST rate constants of k3=3.7×1012 s-1, k1=6.3×1011 s-1 and the 3B1 state lifetime of 231 fs, for the WC transition probability (Figure 5b). For the LZ transition probability, the corresponding values are k3=3.8×1012 s-1, k1=4.9×1011 s-1 and 233 fs. The rate constants and the lifetime obtained from the AIMS dynamics and calculated with the NATST are reported in Table 2. The lifetime of the 3B1 state predicted by the AIMS

21



molecular dynamics is about two times shorter than the lifetime calculated with the NATST. This discrepancy can be explained by dynamic effects. In the NA-TST, for a transition between the two states to occur the reaction must proceed through the MECP, whereas in the AIMS molecular dynamics, transitions can occur anywhere close to the intersection seam (see eq 23 for spawning criterion). From Figure 2, it is clear that a significant number of the spawning events happen far from the MECP, which results in more population transfer and ultimately in a shorter lifetime of the 3B1 state.

Transition probability

1

a)

0.8 0.6 0.4 0.2 0

0

1013

k3

1012

k1

1011 1010

5,000 10,000 15,000 20,000 25,000 30,000 Internal energy, cm -1 b)

Rate constant, s-1

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

0

5,000 10,000 15,000 20,000 25,000 30,000 Internal energy, cm -1

Figure 4. (a) The 3B1→1A1 transition probabilities calculated using the LZ formula with double passage (solid blue line) and the WC formula (solid green line). The vertical dashed red line represents the MECP energy with respect to the 3B1 minimum. (b) The 3

B1→1A1 rate constants k3 (red); and the 1A1→3B1 rate constants k1 (blue) as functions of

the internal energy obtained with LZ transition probability (solid lines) and WC transition probability (dashed lines). The internal energy corresponding to k3 is relative to the 3B1

22


Page 22 of 35

Page 23 of 35

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


minimum; the internal energy corresponding to k1 is relative to the 1A1 minimum.

Table 2. The rate constants k3, k1 and the 3B1 state lifetime τ obtained with the AIMS molecular dynamics, as well as with the NA-TST using the LZ and WC transition probabilities. The lower and upper bounds with 95% confidence are reported in parentheses. NA-TST

AIMS

LZ

WC

k3, s-1

5.3×1012 (4.7×1012, 5.9×1012)

3.8×1012

3.7×1012

k1, s-1

3.1×1012 (1.5×1012, 4.7×1012)

4.9×1011

6.3×1011

τ, fs

119 (94, 164)

233

231

V. Conclusions An extension of the ab initio multiple spawning molecular dynamics method to study the ISC between electronic states with different spin multiplicities coupled by spinorbit interaction was presented. The new method was applied to population decay from the excited 3B1 to the ground 1A1 state of GeH2. The rate constants for the 3B1→1A1 and reverse 1A1→3B1 transitions, as well as the lifetime of the 3B1 state were calculated by fitting the excited state population with the first-order kinetics equation. The ISC rate constants obtained from the ab initio multiple spawning dynamics were compared with the values calculated using the statistical nonadiabatic transition state theory. In the statistical theory, the probability of transitions between the 3B1 and 1A1 states was

23



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

calculated with the Landau-Zener and weak coupling formulas. The ab initio multiple spawning dynamics predicts a 3B1→1A1 rate constant of k3=5.3×1012 s-1 and a reverse 1A1→3B1 rate constant of k1=3.1×1012 s-1. These rate constants correspond to a 3B1 state lifetime of 119 fs. The k3 values calculated with nonadiabatic transition state theory using the Landau-Zener and weak coupling probability formulas are about a factor of 1.4 smaller than the value predicted by the ab initio multiple spawning dynamics; the k1 values are about 5 times smaller. This leads to a longer, compared with multiple spawning dynamics, 3B1 state lifetime of 231 fs for the weak coupling transition probability, and 233 fs for the Landau-Zener probability. The longer 3B1 state lifetime predicted by the NA-TST can be explained by the local nature of the statistical theory, which assumes that transitions between electronic states always occur at the MECP. In contrast, in the ab initio multiple spawning dynamics, these transitions can occur at any point along the intersection seam between the potential energy surfaces of two states.

Acknowledgements DAF and SAV thank Profs. Benjamin Levine and Todd Martinez for fruitful discussions, and the University of Nevada, Reno for financial support. MSG and KK acknowledge support from an Air Force Office of Scientific Research CoDesign grant, under AFOSR Award No. FA9550-12-1-0476. SRP was supported by the Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357.

Supporting Information Population of 3B1 state as a function of time calculated with different values of the spawning threshold Λ eff (Figure S1).

24


Page 24 of 35

Page 25 of 35

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


References (1)

Cavaleri, J. J.; Prater, K.; Bowman, R. M. An Investigation of the Solvent Dependence on the Ultrafast Intersystem Crossing Kinetics of Xanthone. Chem. Phys. Lett. 1996, 13, 495–502.

(2)

Satzger, H.; Schmidt, B.; Root, C.; Zinth, W. Ultrafast Quenching of the Xanthone Triplet by Energy Transfer: New Insight into the Intersystem Crossing Kinetics. J. Phys. Chem. A. 2004, 108, 10072–10079.

(3)

Minns, R. S.; Parker, D. S. N.; Penfold, T. J.; Worth, G. A.; Fielding, H. H. Competing Ultrafast Intersystem Crossing and Internal Conversion in the "Channel 3" Region of Benzene. Phys. Chem. Chem. Phys. 2010, 12, 15607– 15615.

(4)

Balucani, N.; Leonori, F.; Casavecchia, P.; Fu, B.; Bowman, J. M. Crossed Molecular Beams and Quasiclassical Trajectory Surface Hopping Studies of the Multichannel Nonadiabatic O(3P) + Ethylene Reaction at High Collision Energy. J. Phys. Chem. A. 2015, 119, 12498–12511.

(5)

Griesbeck, A. G.; Abe, M.; Bondock, S. Selectivity Control in Electron Spin Inversion Processes: Regio- and Stereochemistry of Paternò-Büchi Photocycloadditions as a Powerful Tool for Mapping Intersystem Crossing Processes. Acc. Chem. Res. 2004, 37, 919–928.

(6)

McCusker, J. K.; Vlček, A., Jr. Ultrafast Excited-State Processes in Inorganic Systems. Acc. Chem. Res. 2015, 48, 1207–1208.

(7)

Auböck, G.; Chergui, M. Sub-50-fs Photoinduced spin crossover in [Fe(bpy)3]2+ Nat. Chem. 2015, 7, 629–633.

(8)

Maeda, S.; Taketsugu, T.; Ohno, K.; Morokuma, K. From Roaming Atoms to Hopping Surfaces: Mapping Out Global Reaction Routes in Photochemistry. J. Am. Chem. Soc. 2015, 137, 3433–3445.

(9)

Marian, C. M. Spin-orbit Coupling and Intersystem Crossing in Molecules. WIREs Comput. Mol. Sci. 2011, 2, 187–203.

(10)

Zhang, W.; Alonso-Mori, R.; Bergmann, U.; Bressler, C. Tracking Excited-State Charge and Spin Dynamics in Iron Coordination Complexes. Nature. 2014, 509, 345–348.

(11)

Yson, R. L.; Gilgor, J. L.; Guberman, B. A.; Varganov, S. A. Protein Induced Singlet-Triplet Quasidegeneracy in the Active Site of [NiFe]-hydrogenase. Chem. Phys. Lett. 2013, 577, 138–141.

(12)

Zaari, R. R.; Varganov, S. A. Nonadiabatic Transition State Theory and

25



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

Trajectory Surface Hopping Dynamics: Intersystem Crossing Between 3B1 and 1 A1 States of SiH2. J. Phys. Chem. A. 2015, 119, 1332–1338. (13)

Zener, C. Nonadiabatic Crossing of Energy Levels. Proc. R. Soc. A. 1932, 137, 696–702.

(14)

Wittig, C. The Landau-Zener Formula. J. Phys. Chem. B. 2005, 109, 8428–8430.

(15)

Landau, L. D. A Theory of Energy Transfer. II. In Collected Papers of L. D. Landau; ter Haar, D., Ed.; Elsevier, 1965; pp. 63–66.

(16)

Nikitin, E. E. Nonadiabatic Transition Near the Turning Point in Atomic Collisions. Opt. Spectrosc. 1961, 11, 246–248.

(17)

Delos, J. B. On the Reactions of N2 with O. J. Chem. Phys. 1973, 59, 2365–2369.

(18)

Lorquet, J. C.; Leyh-Nihant, B. Nonadiabatic Unimolecular Reactions. 1. A Statistical Formulation for the Rate Constants. J. Phys. Chem. 1988, 92, 4778– 4783.

(19)

Harvey, J. N.; Aschi, M. Spin-Forbidden Dehydrogenation of Methoxy Cation: a Statistical View. Phys. Chem. Chem. Phys. 1999, 1, 5555–5563.

(20)

Lévêque, C.; Taïeb, R.; Köppel, H. Theoretical Prediction of the Importance of the 3B2 State in the Dynamics of Sulfur Dioxide. J. Chem. Phys. 2014, 140, 091101.

(21)

Zhao, J. A Quantum Time-Dependent Wave-Packet Study of Intersystem Crossing Effects in the O(3P0, 1, 2) + D2(ν=0, j=0) Reaction. J. Chem. Phys. 2013, 138, 134309.

(22)

Billing, G. D. On the Use of Ehrenfest's Theorem in Molecular Scattering. Chem. Phys. Lett. 1983, 100, 535–539.

(23)

Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106.

(24)

Tully, J. C. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H+ with D2. J. Chem. Phys. 1971, 55, 562.

(25)

Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061.

(26)

Mai, S.; Marquetand, P.; González, L. A General Method to Describe Intersystem Crossing Dynamics in Trajectory Surface Hopping. Int. J. Quant. Chem. 2015, 115, 1215–1231.

26


Page 26 of 35

Page 27 of 35

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


(27)

Mai, S.; Marquetand, P.; González, L. Non-adiabatic and Intersystem Crossing Dynamics in SO2. II. The Role of Triplet States in the Bound State Dynamics Studied by Surface-Hopping Simulations. J. Chem. Phys. 2014, 140, 204302.

(28)

Richter, M.; Marquetand, P.; González-Vázquez, J.; Sola, I.; González, L. SHARC: Ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. J. Chem. Theory Comput. 2011, 7, 1253–1258.

(29)

Hu, W.; Lendvay, G.; Maiti, B.; Schatz, G. C. Trajectory Surface Hopping Study of the O(3P) + Ethylene Reaction Dynamics. J. Phys. Chem. A. 2008, 112, 2093– 2103.

(30)

Jasper, A. W.; Dawes, R. Non-Born–Oppenheimer Molecular Dynamics of the Spin-Forbidden Reaction O(3P) + CO(X1Σ+) → CO2( X1Σg+ ). J. Chem. Phys. 2013, 139, 154313–154314.

(31)

Maiti, B.; Schatz, G. C.; Lendvay, G. Importance of Intersystem Crossing in the S(3P,1D) +H2 → SH + H. J. Phys. Chem. A. 2004, 108, 8772–8781.

(32)

Jasper, A. W.; Stechmann, S. N.; Truhlar, D. G. Fewest-switches with Time Uncertainty: A Modified Trajectory Surface-Hopping Algorithm with Better Accuracy for Classically Forbidden Electronic Transitions. J. Chem. Phys. 2002, 116, 5424.

(33)

Jasper, A. W.; Truhlar, D. G. Improved Treatment of Momentum at Classically Forbidden Electronic Transitions in Trajectory Surface Hopping Calculations. Chem. Phys. Lett. 2003, 369, 60–67.

(34)

Schmidt, J. R.; Parandekar, P. V.; Tully, J. C. Mixed Quantum-Classical Equilibrium: Surface Hopping. J. Chem. Phys. 2008, 129, 044104–044107.

(35)

Parandekar, P. V.; Tully, J. C. Detailed Balance in Ehrenfest Mixed QuantumClassical Dynamics. J. Chem. Theory Comput. 2006, 2, 229–235.

(36)

Ben-Nun, M.; Martinez, T. J. Ab Initio Quantum Molecular Dynamics. Adv. Chem. Phys. 2002, 121, 439–512.

(37)

Heller, E. J. Frozen Gaussians: A Very Simple Semiclassical Approximation. J. Chem. Phys. 1981, 75, 2923–2931.

(38)

Ben-Nun, M.; Martı́nez, T. J. Nonadiabatic Molecular Dynamics: Validation of the Multiple Spawning Method for a Multidimensional Problem. J. Chem. Phys. 1998, 108, 7244.

(39)

Ben-Nun, M.; Martinez, T. J. A Continuous Spawning Method for Nonadiabatic Dynamics and Validation for the Zero-Temperature Spin-Boson Problem. Isr. J. Chem. 2007, 47, 75–88. 27



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)

Ben-Nun, M.; Martı́nez, T. J. Ab Initio Molecular Dynamics Study of cis–trans Photoisomerization in Ethylene. Chem. Phys. Lett. 1998, 298, 57–65.

(41)

Levine, B. G.; Martínez, T. J. Ab Initio Multiple Spawning Dynamics of Excited Butadiene: Role of Charge Transfer. J. Phys. Chem. A. 2009, 113, 12815–12824.

(42)

Gaenko, A.; DeFusco, A.; Varganov, S. A.; Martínez, T. J.; Gordon, M. S. Interfacing the Ab Initio Multiple Spawning Method with Electronic Structure Methods in GAMESS: Photodecay of trans-Azomethane. J. Phys. Chem. A. 2014, 118, 10902–10908.

(43)

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.

(44)

Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G. et al. MOLPRO, version 2012.1, a package of ab initio programs, 2012; see http://www.molpro.net.

(45)

Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a General-Purpose Quantum Chemistry Program. WIREs Comput. Mol. Sci. 2011, 2, 242–253.

(46)

Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: GAMESS a Decade Later. In Theory and Applications of Computational Chemistry; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier, Amsterdam, 2005; pp 1167-1189.

(47)

Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347–1363.

(48)

Curchod, B. F. E.; Rauer, C.; Marquetand, P.; González, L.; Martínez, T. J. Communication: GAIMS — Generalized Ab Initio Multiple Spawning for Both Internal Conversion and Intersystem Crossing Processes. J. Chem. Phys. 2016, 144, 101102–101107.

(49)

Jasper, A. W. Multidimensional Effects in Nonadiabatic Statistical Theories of Spin- Forbidden Kinetics: A Case Study of 3O + CO → CO2 . J. Phys. Chem. A. 2015, 119, 7339–7351.

(50)

Lykhin, A. O.; Kaliakin, D. S.; dePolo, G. E.; Kuzubov, A. A.; Varganov, S. A. Nonadiabatic Transition State Theory: Application to Intersystem Crossings in the Active Sites of Metal-Sulfur Proteins. Int. J. Quant. Chem. 2016, 116, 750761.

(51)

Harvey, J. N. Understanding the Kinetics of Spin-Forbidden Chemical Reactions. Phys. Chem. Chem. Phys. 2007, 9, 331.

28


Page 28 of 35

Page 29 of 35

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


(52)

Nikitin, E. E.; Umanskii, S. Y. Theory of Slow Atomic Collisions; SpringerVerlag Berlin Heidelberg, New York, 1984.

(53)

Miller, W. H.; Handy, N. C.; Adams, J. E. Reaction path Hamiltonian for Polyatomic Molecules. J. Chem. Phys. 1980, 72, 99.

(54)

Baer, T.; Hase, W. L. Unimolecular Reaction Dynamics: Theory and Experiments; Oxford University Press: New York, 1995.

(55)

Green, N. J. B. Unimolecular Kinetics, Part 1. The Reaction Step; Green N. J. B., Ed.; Elsevier Science: Amsterdam, 2003.

(56)

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.

(57)

Martinez, T. J.; Ben-Nun, M.; Levine, R. D. Multi-Electronic-State Molecular Dynamics: A Wave Function Approach with Applications. J. Phys. Chem. A. 1996, 100, 7884–7895.

(58)

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.

(59)

Bethe, H. A.; Salpeter, E. E. Quantum Mechanics of the One and Two Electron Atoms; Plenum, New York, 1977.

(60)

Mathews, J.; Walker, R. L. Mathematical Methods of Physics; W. A. Benjamin, New York, 1980.

(61)

Martinez, T. J.; Ben-Nun, M.; Levine, R. D. Molecular Collision Dynamics on Several Electronic States. J. Phys. Chem. A. 1997, 101, 6389–6402.

(62)

Granucci, G.; Persico, M. Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping. J. Chem. Phys. 2007, 126, 134114.

(63)

Roos, B. O. Advances in Chemical Physics; Wiley Interscience, New York, 1987; Vol. 69, pp 339–445.

(64)

Fedorov, D. G.; Gordon, M. S. A Study of the Relative Importance of One- and Two-Electron Contributions to Spin–Orbit Coupling. J. Chem. Phys. 2000, 112, 5611-5623.

(65)

Hillery, M.; O'connell, R. F.; Scully, M. O.; Wigner, E. P. Distribution Functions in Physics: Fundamentals. Phys. Rep. 1984, 106, 121–167.

29



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

30


Page 30 of 35

Page 31 of 35


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


Number of spawns

Number of spawns

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 32 of 35

400 a) 300 200 100 0 1.3

1.4

1.5

1.6 1.7 r, angstrom

1.8

1.9

200 b) 150 100 50 0 125

130

135 140 θ, degrees


145

150

Page 33 of 35

Population of B 1 state

1

3

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


0.9 0.8 0.7 0.6 0

50

Time, fs


100

150


Transition probability

1

a)

0.8 0.6 0.4 0.2 0

0

5,000 10,000 15,000 20,000 25,000 30,000 -1 Internal energy, cm

1013 b)

Rate constant, s-1

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 34 of 35

10

k

12

k1

1011 1010

3

0

5,000 10,000 15,000 20,000 25,000 30,000 -1 Internal energy, cm


Page 35 of 35

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


MECP$

S=1$ S=0$


Ab Initio Multiple Spawning Method for Intersystem ... - ACS Publications

Recommend Documents