Geometric Phase Effects in Nonadiabatic Dynamics near Conical

Article pubs.acs.org/accounts

Geometric Phase Effects in Nonadiabatic Dynamics near Conical Intersections Ilya G. Ryabinkin,†,‡ Loïc Joubert-Doriol,†,‡ and Artur F. Izmaylov*,†,‡ †

Department of Physical and Environmental Sciences, University of Toronto Scarborough, Toronto, Ontario M1C 1A4, Canada Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, Ontario M5S 3H6, Canada

‡

CONSPECTUS: Dynamical consideration that goes beyond the common Born−Oppenheimer approximation (BOA) becomes necessary when energy differences between electronic potential energy surfaces become small or vanish. One of the typical scenarios of the BOA breakdown in molecules beyond diatomics is a conical intersection (CI) of electronic potential energy surfaces. CIs provide an efficient mechanism for radiationless electronic transitions: acting as “funnels” for the nuclear wave function, they enable rapid conversion of the excessive electronic energy into the nuclear motion. In addition, CIs introduce nontrivial geometric phases (GPs) for both electronic and nuclear wave functions. These phases manifest themselves in change of the wave function signs if one considers an evolution of the system around the CI. This sign change is independent of the shape of the encircling contour and thus has a topological character. How these extra phases affect nonadiabatic dynamics is the main question that is addressed in this Account. We start by considering the simplest model providing the CI topology: two-dimensional two-state linear vibronic coupling model. Selecting this model instead of a real molecule has the advantage that various dynamical regimes can be easily modeled in the model by varying parameters, whereas any fixed molecule provides the system specific behavior that may not be very illustrative. After demonstrating when GP effects are important and how they modify the dynamics for two sets of initial conditions (starting from the ground and excited electronic states), we give examples of molecular systems where the described GP effects are crucial for adequate description of nonadiabatic dynamics. Interestingly, although the GP has a topological character, the extent to which accounting for GPs affect nuclear dynamics profoundly depends on topography of potential energy surfaces. Understanding an extent of changes introduced by the GP in chemical dynamics poses a problem of capturing GP effects by approximate methods of simulating nonadiabatic dynamics that can go beyond simple models. We assess the performance of both fully quantum (wave packet dynamics) and quantum-classical (surface-hopping, Ehrenfest, and quantum-classical Liouville equation) approaches in various cases where GP effects are important. It has been identified that the key to success in approximate methods is a method organization that prevents the quantum nuclear kinetic energy operator to act directly on adiabatic electronic wave functions.

■

INTRODUCTION

conceptual advantage of separating energies of electronic states maximally for static nuclei. This separation allows one to think about nuclear dynamics as if it evolves primarily under forces from a single uncoupled potential energy surface (PES). Due to parametric nuclear dependence of the electronic problem, there are regions where PESs can either come close to each other or even intersect. Although the BO approximation fails in these regions, the AR only needs to incorporate all electronic states whose PESs are energetically close. The big advantage of the AR for the intuitive interpretation of dynamics on multiple PESs (i.e., nonadiabatic dynamics) is that couplings between electronic states are usually localized in the space of nuclear configurations, thus, dynamics can be thought intuitively as shown in Figure 1. The locality of couplings also helps in

Most of the first-principles investigations of dynamics in molecules and solids usually start with solving the electronic structure problem or introducing the Born−Oppenheimer (BO) separation of electronic and nuclear variables.1,2 Compared to other approaches, the BO picture has a clear advantage of general rigor. Here we distinguish the BO separation procedure from the BO approximation, the former introduces the adiabatic representation (AR) for electronic wave functions with parametric dependence on nuclear coordinates and is potentially exact, whereas the latter ignores residual electronic state couplings and constitutes an approximation. The AR forms the backbone of our understanding of molecular dynamics. Compared to another frequently used diabatic representation, it has a practical advantage of being delivered by standard electronic structure programs and a © 2017 American Chemical Society

Received: May 1, 2017 Published: June 30, 2017 1785

DOI: 10.1021/acs.accounts.7b00220 Acc. Chem. Res. 2017, 50, 1785−1793

Article

Accounts of Chemical Research

where Φk(r|R) are the electronic adiabatic wave functions parametrically dependent on R. They are defined as eigenstates of the electronic Hamiltonian Ĥ e Φk (r|R) = Ek (R)Φk (r|R)

(2)

Projection of the TDSE onto {Φk(r|R)} allows one to define equations for the nuclear components, χk(R, t)

∑ [(Tn̂ + Ek(R))δkj + τkĵ ]χj (R, t ) = i∂tχk (R, t ) (3)

j

where τ̂kj = −⟨Φk(R)|∇2Φj(R)⟩/2 − ⟨Φk(R)|∇Φj(R)⟩∇ are the nonadiabatic couplings (NACs) with ∇ representing gradient over all mass-weighted nuclear coordinates R. We will refer to the operator on the left-hand side of eq 3 as Ĥ adi. Existence and finiteness of the Φj(r|R) derivatives with respect to the nuclear coordinates in NACs requires Φj(r|R) to be at least continuous with R. On the other hand, Φk(r|R) can be multiplied by any arbitrary R dependent phase factor, eiφk(R), because this does not affect the electronic eigenvalue problem

Figure 1. Molecular dynamics on multiple adiabatic electronic surfaces.

quantitative modeling using mixed quantum-classical (MQC) approaches.3 However, for the fully quantum nuclear dynamics the locality of couplings frequently leads to divergencies in case of PES crossings and numerical difficulties in simulations. Moreover, for the conical intersections (CIs) that are particularly ubiquitous in molecules beyond diatomics,4−6 there is another complication arising from nontrivial geometric phases (GPs) appearing in the adiabatic electronic and nuclear wave functions.7−9 These GPs were always somewhat mysterious features of the AR for many practitioners of nonadiabatic dynamics. And questions like what do these phases do to dynamics, when are they important, or how does one account for them properly in practical schemes were not covered in the literature in adequate depth. Mainly for numerical reasons, a preferable representation for the fully quantum nuclear dynamics has been predominantly diabatic. In the diabatic representation we sacrifice the locality of couplings for simplicity of involved operators. The diabatic representation is free of divergencies and nontrivial GPs, but it requires a system dependent diabatization procedure that can introduce additional errors; also this representation is somewhat inferior for the popular fewest-switches surface hopping (FSSH) approach.3,10 In this Account, we will highlight what GPs add to the intuitive picture of nonadiabatic dynamics in the AR, how important these contributions can be, and when they occur. Also, it was found that a proper account for the GP can be used as an exact constraint in developing approximate methods in nonadiabatic dynamics. Thus, the second major theme will be introducing the GP in practical simulation techniques and developing better methods for nonadiabatic dynamics within the AR.

■

Ĥ e eiφk(R)Φk (r|R) = Ek (R) eiφk(R) Φk (r|R)

Electronic structure codes have algorithms to select this phase factor using information within a single nuclear configuration. To perform nuclear dynamics, this phase freedom can be used to make Φk(r|R) continuous and differentiable with respect to R by performing a phase matching procedure.12 Of course, even if the differentiability condition is imposed, the phase is still defined only up to any differentiable function of R, which is an inherent gauge freedom of the BH expansion. In all following discussions, to mirror results obtained for model electronic wave functions, their electronic structure counterparts need to undergo the phase matching procedure ensuring differentiability with R. When NACs are small, nuclear dynamics of χk(R, t) simply follows a single electronic state. It is easy to show that generally NACs are inversely proportional to electronic energy gaps, τ̂jk ∼ (Ej(R) − Ek(R))−1.13 Thus, regions of small or zero gap become especially interesting. Conical Intersections

In the absence of symmetric constraints, the most general form of PES intersections is conical. Let us consider a simple example of 2-fold degeneracy, when two adiabatic PESs have the same energy E0 at nuclear configuration R0. Assume that in the vicinity of R0 all the other PESs are well-separated and perform a Taylor expansion of a 2 × 2 submatrix of the electronic Hamiltonian, He, with respect to a displacement Q = R − R0: ⎛⟨Φ (R )|Ĥ (R)|Φ (R )⟩ ⟨Φ (R )|Ĥ (R)|Φ (R )⟩ ⎞ 1 0 1 0 e 2 0 ⎜ 1 0 e ⎟ ⎜ ⎟ ̂ ̂ ⎝⟨Φ2(R 0)|He(R)|Φ1(R 0)⟩ ⟨Φ2(R 0)|He(R)|Φ2(R 0)⟩⎠

THEORY

Adiabatic Representation

Fully quantum molecular dynamics of an isolated system in a pure state, Ψ, is given by the time-dependent Schrödinger equation (TDSE) Ĥ Ψ(r, R, t) = i∂tΨ(r, R, t), where r and R are electronic and nuclear degrees of freedom (DOF), respectively. The total Hamiltonian, Ĥ , consists of the kinetic energy of nuclei, T̂ n, and the electronic Hamiltonian, Ĥ e. Solving the TDSE can be initiated by expanding the total wave function in the Born−Huang (BH) expansion11 Ψ(r, R, t ) =

∑ χk (R, t )Φk(r|R) k

(4)

⎛ καQ α cαQ α ⎞ ⎟⎟ + O(Q 2) α κ c Q Q ̃ ⎝ α α α=1 α α⎠ N

= E012 +

∑ ⎜⎜

(5)

where 12 is the unit 2 × 2 matrix, N = dim(R), and components of vectors κ = (κ1, ..., κN), κ̃̃ = (κ̃1, ..., κ̃N), and c = (c1, ..., cN) are the first-order partial derivatives of the corresponding matrix elements at R0. When the displacement Q is small, one can neglect higherorder terms in eq 5 and introduce new variables κX = (κ − κ̃)· Q/2 and cY = c·Q, where κ = ∥κ − κ̃∥/2 and c = ∥c∥. In these

(1) 1786


Article


compensate the double-valuedness of Φk(r|R). However, modeling χk(R, t) with double-valued boundary conditions is generally difficult. Also, note that this consideration will become even more complicated for cases with multiple CIs. To avoid problems with double-valued boundary conditions and yet have a proper representation, Mead and Truhlar15 suggested to exploit the gauge freedom in the BH expansion by introducing a resolution of the identity (RI), 1 = exp(iθ) exp(−iθ)

new coordinates, the gap between two eigenvalues of the matrix (adiabatic PESs) is E+(R) − E−(R) = 2 κ 2X2 + c 2Y 2

(6)

5

and has a topography of a CI. (X,Y)-subspace is known as the branching space, whereas its (N − 2)-dimensional complement constitutes the CI seam. In the vicinity of the CI, the two PESs remain degenerate in the seam subspace X = Y = 0. When N ≫ 2, the seam conditions are easily satisfied, which makes CIs almost unavoidable in large molecules.14 However, unless nuclear wave functions have appreciable probability densities in the vicinity of a CI, CIs existence may not affect molecular dynamics. Thus, cases where a CI affects dynamics in spite of negligible probability density at the CI are especially interesting and will be highlighted below.

Ψ(r, R, t ) =

∑ χk (R, t ) eiθ e−iθ Φk(r|R) k

(10)

Note that the most general RI would involve mixing electronic and nuclear wave functions of different states, 1 = U†U, but that would ruin the AR. Due to double-valuedness of exp(±iθ) factors (Figure 3), one can define new single-valued nuclear, χ̃k

Geometric Phase

Adiabatic PESs and electronic wave functions are determined by diagonalization of the electronic Hamiltonian matrix, He, using the unitary transformation ⎛ cos θ −sin θ ⎞ U(θ ) = ⎜ ⎟ ⎝ sin θ cos θ ⎠

(7)

where θ (X , Y ) =

⎛ γY ⎞ 1 arctan⎜ ⎟ ⎝X⎠ 2

(8)

is a mixing angle between basis states |Φk=1,2(R0)⟩, and γ = c/κ is a coupling strength constant. Columns of U are the adiabatic wave functions |Φ1(R)⟩ = cos θ|Φ1(R 0)⟩ + sin θ|Φ2(R 0)⟩ |Φ2(R)⟩ = −sin θ|Φ1(R 0)⟩ + cos θ|Φ2(R 0)⟩

Figure 3. Phase factor exp(iθ).

(9)

= χk exp(iθ), and electronic, Φ̃k = exp(−iθ) Φk, wave functions. Subsequently, one can rewrite the nuclear TDSE (eq 3) as

Encircling a single CI point, placed at the origin for convenience, the geometric angle ϕ = arctan(Y/X) grows from 0 to 2π, whereas θ evolves only up to π (Figure 2);

∑ [eiθ Ĥ

adi −iθ

e

]jk χk̃ = i∂tχj̃

(11)

k

and introduce the modified nuclear Hamiltonian GP adi Ĥ = eiθ Ĥ e−iθ

(12)

Owing to χ̃k(R, t) single-valuedness, variational techniques can be applied to eq 11 with an ordinary single-valued boundary conditions, employing for example, Gaussians or harmonicoscillator eigenfunctions as a basis set. The single-valued electronic wave functions, Φ̃k(r|R), can define the GP according to Berry’s expression8 Figure 2. Mixing angle (θ) as a function of the geometric angle (ϕ) for γ = 0.3: when the CI is encircled (red), and when the CI is not encircled (blue).

θ(2π ) − θ(0) = i

∫C ⟨Φ̃k(R)|∇Φ̃k(R)⟩ dl = π

(13)

where the integral is taken over a contour C encircling the CI. The integral eq 13 with the double-valued functions Φk(r|R) gives zero, because Φk(r|R) already have a sign change in them. There is still a freedom to augment θ(R) with any continuous function of R that acquires multiples of 2π when the CI is encircled. In all our model studies,16−19 we used θ(R) from eq 8. This choice has an advantage that all terms in eq 11 have the same functional form as in eq 3, which not only is computationally convenient but also allows us to consider an interplay between regular nonadiabatic and GP-related terms.

therefore, both adiabatic wave functions flip their signs. This sign change can be thought of as accumulation of the GP, θ(R), during a continuous change of the Hamiltonian parameters.8,9 The sign change due to the GP makes the electronic wave functions double-valued with respect to the nuclear variables. Accounting for the double-valuedness of Φk(r|R) becomes indispensable in solving for the total Ψ(r, R, t). To keep the total molecular wave function Ψ(r, R, t) single-valued, one has to impose double-valued boundary conditions for χk(R, t) to 1787


Article


Figure 4. (a) Low-energy dynamics where SP1 and SP2 label saddle points between the two potential minima; (b) excited-state dynamics.

■

GEOMETRIC PHASE EFFECTS IN DYNAMICS To avoid a complex interplay between GP and non-GP effects (e.g., quality of electronic structure data and quantum dynamics approximations), which is inevitably present in studies of real molecules and in comparison of theory to experiment,20−22 our approach was to consider GP effects in the simplest possible models for CIs. Such models allow one to establish clearly what GP effects are and when they can affect dynamics significantly.

̂ = τ21 ̂† = τ12

The dimensionality of the branching space naturally suggests a minimal CI model to be the two-dimensional (2D) linear vibronic coupling (LVC) model in the diabatic representation

(14)

1 (ωX 2X2 + ωY 2Y 2 + Δ) + κX 2

(16)

V22 =

1 (ωX 2X2 + ωY 2Y 2 − Δ) − κX 2

(17)

V12 = cY

̂GP = τ22 ̂ GP = 2τ11 ̂ + τ12 ̂ τ11

(23)

̂GP = τ21 ̂ GP † = τ12 ̂ − 2τ11 ̂ τ12

(24)

Low-Energy Dynamics

(15)

V11 =

(22)

Comparing these equations with eqs 21 and 22 reveals that adding GP makes both diagonal and off-diagonal terms to become linear combinations of potential-like and L̂ z-dependent terms.

where 2 ⎞ ⎛ 2 ̂ = − 1 ⎜ ∂ + ∂ ⎟ ≡ − 1 ∇·∇ T2D 2 ⎝ ∂X2 2 ∂Y 2 ⎠

4(γ −1X2 + γY 2)

where L̂ z = −i(X∂/∂Y − Y∂/∂X) is the angular momentum operator, and the overhead arrows indicate the operator’s directionality. Without GP, diagonal and off-diagonal NACs were clearly split in two types: potential-like terms, also known as the diagonal Born−Oppenheimer correction (DBOC),24,25 and L̂ z-dependent terms. The corresponding adiabatic Hamiltonian with GP, Ĥ GP 2D , has the same form as Ĥ adi but different NACs 2D

Two-Dimensional Linear Vibronic Coupling Model

⎛V V ⎞ ̂ ⊗ 12 + ⎜ 11 12 ⎟ Ĥ 2D = T2D ⎝ V12 V22 ⎠

→ ← Lẑ − Lẑ

GP profoundly affects quantum dynamics even at low energies, when nuclei move close to minima of the lowest adiabatic PES and may not even reach the CI (Figure 4a). To illustrate GP effects, we considered dynamics of a localized wave packet corresponding to the ground state of the uncoupled diabatic Hamiltonian T̂ 2D + V11 using three Hamiltonians: (1) BO, Ĥ BO GP ̂ = T̂ 2D + E1(X, Y); (2) BO + GP, Ĥ GP BO = T2D + τ11 + E1(X, Y); 26 ̂ (3) exact, H2D in eq 14. Comparing dynamics of the initial adiabatic well population for Ĥ BO and Ĥ GP BO (Figure 5 left)

(18)

Here, Δ is the energy difference between minima of the V11 and V22 paraboloids. We start from the diabatic model because its transformation to the AR is exact, while the inverse transformation is generally not.23 Diabatic-to-adiabatic transformation (eq 7) brings Ĥ 2D to the form ⎛T̂ + τ ̂ τ ̂ ⎞ ⎛E 0 ⎞ 2D 11 12 adi ⎟+⎜ 1 ⎟ Ĥ 2D = ⎜⎜ ⎟ ̂ + τ22 ̂ ̂ ⎠ ⎝ 0 E2 ⎠ T2D ⎝ τ21

(19)

Figure 5. Probability (P) to find the system in the initial well as a function of time for two sets of parameters.

(20)

shows that tunnelling between the two wells taking place in the BO approach can be quenched by including the GP. The exact dynamics of H2D does not have qualitative differences with that of Ĥ GP BO, which confirms that the GP is the main feature separating the BO approximation from the exact result. The reason for the GP tunnelling freeze is that the initial localized wave packet is almost exactly equal to the ground vibronic ̂ eigenstate of the Ĥ GP BO and H2D Hamiltonians. The origin of a

where E1,2 =

1 1 (V11 + V22) ∓ (V11 − V22)2 + 4V12 2 2 2

are the adiabatic PESs with the minus (plus) sign for E1 (E2), see Figure 4. The NACs for the adiabatic states defined by eq 9 are ̂ = τ22 ̂ = τ11

X2 + Y 2 8(γ −1X2 + γY 2)2

(21) 1788


Article


Figure 6. Nuclear probability density at t = 1.6 au for Δ = 0, γ = 1 in exact (the nodal line is denoted by the dashed line) and BO simulations.

drastic difference between eigenstates of the Hamiltonians accounting for the GP and Ĥ BO can be traced to the symmetry of the problem.26 The eigenstate localization is robust with respect to variations in Δ and γ.26 This robustness can be rationalized from a dynamical point of view. If one puts a localized wave packet in one minimum, it will split into two parts that will try to reach the other minimum through symmetry-equivalent saddle points SP1 and SP2, (Figure 4a). Along these paths, fragments will acquire the same dynamical but opposite geometric phases. Even when localization is removed by lowering the barriers (Figure 5 right), GP induced destructive interference causes a nodal line to appear and to reduce the nuclear probability density in the acceptor well (Figure 6). Interestingly, the destructive interference observed with the node-less initial nuclear wave function can be turned into constructive interference if one considers a nuclear wave function with a nodal line along the x-direction.27 GP induces the constructive interference for the latter case because the two halves of the wave packet encircling the CI have different signs from the beginning. A symmetry breaking that makes energies of SP1 and SP2 different can remove the localization for a node-less nuclear wave function by energetic suppression of one of the paths around the CI, thus destroying interference. Can this GP effect survive in large molecules? This question was addressed by considering N-dimensional LVC models.16 Initially, Kelly and Kapral28 found that the nodal line signature of destructive interference disappeared in the N-dimensional case. However, a more systematic analysis established that the nodal line can simply move within a space of N nuclear DOF depending on how the branching plane coordinates are coupled to the remaining DOF, and its disappearance does not generally mean the absence of GP effects.16 A general class of molecules where GP effects can play a role in determining the structure and dynamics is Jahn−Teller (JT) distorted compounds. Bis(methylene) adamantyl (BMA) carbocation is the system where the minimum of the CI seam has the D2d symmetry. Thus, the system undergoes a JT distortion into two conformations of the C2v symmetry (Figure 7). Adiabatic barriers between two conformers are relatively small so that BO dynamics starting from a localized distribution in one well will oscillate between the two wells. In more accurate simulations accounting for GP effects dynamics do not produce significant oscillations and are qualitatively very similar to that in Figure 5 left.29,30 Interestingly, even for non-JT compounds similar GP effects have been found recently to be

Figure 7. Jahn−Teller distortions in the BMA carbocation.

responsible for reducing the rate of photoinduced dissociation of the phenol molecule in gas phase.31 Excited-State Dynamics

To model a typical photoinduced dynamics, one starts with an initial nuclear wave packet on the excited PES. Due to the funnel-like shape of the excited state PES near the CI (Figure 4b), some part of a nuclear wave packet will pass directly through the CI in the radiationless deactivation process.17 This dynamics is clearly different from that considered for the lowenergy case not only by the wave packet motion with respect to the CI but also by involvement of both electronic states; thus one can expect different GP effects in this case. In the vicinity of the CI (X ≈ −Δ/κ, Y ≈ 0), nonadiabatic matrix elements are strongly singular (eqs 21 and 22); therefore, the details of PESs become irrelevant. Analysis of NACs close to the CI revealed that GP-related terms in the Hamiltonian of eq 12 play a 2-fold role.17 First, they compensate for purely repulsive DBOC terms in eq 21. DBOC repulsion affects dynamics especially strongly in systems with small γ’s, where the DBOC rises as an extended wall on a way of the wave packet to the CI (Figure 8). In the adiabatic representation, DBOC extension can be correlated with deviation of ratio ∥h∥/∥g∥ from one, where h = (E1 − E2)⟨Φ1|∇Φ2⟩ and g = ∇(E1 − E2) are orthogonal vectors often evaluated in CI studies.32 Second, GP-related terms enhance the probability of nonadiabatic transitions by modifying the off-diagonal term (eq 22). Consideration of a “Mexican hat” (γ = 1, ωX = ωY) case reveals that τ̂12 (no GP) can be written in the polar (r,ϕ) coordinates centered at the CI as τ̂12 = L̂ z/(2r2). Expanding the nuclear wave packet, χ(X,Y,t), as a Fourier series of the L̂ z eigenfunctions, ∞

χ (X , Y , t ) =

∑ m =−∞

Cm(r , t ) e−imϕ

(25)

helps to analyze efficiency of nonadiabatic transitions without GP by applying τ̂12 to the χ(X,Y,t) expansion. The m = 0 term in eq 25 does not undergo transfer, because L̂ z e−imϕ|m=0 = 0, 1789


Article


Figure 8. DBOC for different γ’s.

and thus the action of τ̂12 on this term is 0. Once the GP is included, the transfer of the m = 0 component becomes C0(r,t) = C0(r,t)/(8r2). Therefore, large differpossible τ̂(GP) 12 ences between dynamics with and without GP should be expected when the nuclear wave packet has a significant portion 2 of the m = 0 component, w0(t) = ∫ ∞ 0 |C0(r,t)| r dr, on its arrival at the vicinity of the CI. To illustrate GP effects in excited-state dynamics, we simulated the adiabatic excited-state population for the butatriene and BMA carbocations using three Hamiltonians: (a) exact, Ĥ 2D eq 14; (b) “noGP”, Ĥ adi 2D with eqs 21 and 22; (c) “noGP noDBOC”, Ĥ adi with eq 22 and τ̂ii = 0. Parameters of 2D 2D-LVC models for these systems were obtained using ab initio data and detailed elsewhere.17 The initial nuclear wave function was assumed to be vertically excited from the minimum of the ground state and evolved from the Franck−Condon point toward the CI. For the butatriene carbocation removing the DBOC does not affect dynamics appreciably because this system has γ = 0.7, which makes the DBOC relatively compact and unnoticeable for the wave packet (Figure 9). Here, GP effects are related mainly to the transfer enhancement due to a large m = 0 weight, w0 = 88%. BMA presents a somewhat opposite case, γ = 0.1 and w0 = 42%, where the DBOC compensation is the main GP contribution, and without GP,

the population dynamics is slower and smaller in its amplitude compared to that in the exact simulation (Figure 9).

■

METHODS Provided examples unambiguously illustrated that ignoring the GP in the AR can lead to inadequate representation of nonadiabatic dynamics. The GP factors [exp(±iθ)] allow one to account for GP effects in model problems. Building a diabatic model where θ(R) can be obtained analytically has been extended to molecular systems recently.30,33 The next question becomes whether it is possible to account for the GP without evaluating the θ(R) function. We identify the root of GP-associated problems in the projection of the kinetic nuclear operator, T̂ n, onto adiabatic electronic wave functions. This projection produces an operator that inherits the double-valued character of the electronic wave functions. Therefore, to eliminate complications with GP, such projections should be avoided without abandoning the AR. We illustrate how this can be achieved in frameworks of a few nonadiabatic dynamics modeling techniques. Ehrenfest and Surface Hopping Methods

Ehrenfest34 and FSSH34,35 schemes treat nuclear dynamics classically, allowing only the electronic DOF to evolve according to the TDSE. Thus, even if electronic wave functions of different trajectories go around the CI from different sides acquiring opposite GPs, they will not be able to interfere because classical particles undergo independent dynamics. Therefore, there is no mechanism for getting destructive interference patterns in low-energy dynamics around a CI in these MQC methods. On the other hand, it was found that Ehrenfest and FSSH methods are capable of mimicking GP effects for the excitedstate dynamics (Figure 9).19,36 First, the DBOC does not appear in the working expressions of these MQC approaches, since T̂ n is treated classically from the beginning.37 Thus, there is no need for the GP-induced DBOC compensation. Moreover, adding the DBOC as an extra potential to PESs in the FSSH method can be detrimental.36 Second, both MQC methods have enhanced transition probabilities compare to those in the “noGP” quantum approach (Figure 9). However, this enhancement is not related to the GP but rather is a result of a continuous nature of the classical nuclear angular momentum, Lz. Off-diagonal matrix element (eq 22) of the adiabatic Hamiltonian (eq 19) has its counterpart in the MQC case but with the classical Lz function. Thus, the absence of the m = 0 component transfer in the quantum “noGP” case is paralleled in MQC only by a negligible fraction of trajectories

Figure 9. Adiabatic excited-state population (Pex) dynamics in the butatriene and BMA carbocations: quantum (solid) and quantumclassical (dashed) methods. 1790


Article


Figure 10. Nuclear probability densities at t = 1.6 au (Δ = 0, γ = 1) with the two QCLE approaches.

with Lz = 0, while the majority of trajectories within the range 0 < |Lz| < 1 do transfer on average as if they had |Lz|= 1/2.19 Therefore, MQC transition probabilities are enhanced, the amount being fortuitously equal to the corresponding GP contribution in quantum dynamics.

considering GP effects has distinguished clearly the correct path in introducing the key steps in the QCL formalism. Time-Dependent Variational Principle

One of the most popular approaches to quantum dynamics with on-the-fly evaluation of PESs is based on the timedependent variational principle (TDVP)42 and involves using frozen-width Gaussians to represent nuclear wave function.43−49 In the nonadiabatic case within the AR, one can use the following two forms for the total nonstationary wave function44,48

Quantum-Classical Liouville Approach

The QCL approach is a more advanced MQC method, which is capable of including more quantum nuclear effects such as interference and decoherence.38 The QCL approach starts with the full Liouville equation, i∂tρ = [Ĥ , ρ], where ρ is the total electron−nuclear density operator ρ(r,R,r′,R′,t) = Ψ*(r,R, t)Ψ(r′,R′,t). In the AR, the QCL approach involves two steps: (1) Wigner transformation (WT) of the nuclear DOF and (2) projection of the electronic DOF onto adiabatic electronic wave functions Φk(r|R). Because Φk(r|R) depend on nuclear coordinates, these two steps do not commute and give rise to two formalisms: Wigner-then-Adiabatic (WA)39 and Adiabaticthen-Wigner (AW).40 Results of these two formalisms for model nonadiabatic problems with avoided crossings were found to be very similar, and it was not clear which of them is generally preferable.41 Investigation of GP effects in CI problems revealed that the AW formalism encounters mathematical inconsistencies in its formulation leading to inability to account for GP effects in the low-energy dynamics.18 In contrast, the WA formalism is mathematically consistent and accounts for GP induced destructive interference (Figure 10). The reason for the WA success is 2-fold. First, when the WT is applied to the Liouville equation, it transforms the T̂ n operator into its classical counterpart, which does not give rise to NACs when projected onto the adiabatic basis. Second, the Wigner transformed density, ρW(r,r′,Rc,P,t), contains only a single nuclear coordinate, Rc. Therefore, when the projection onto adiabatic states is done ρW,ij(R c, P, t ) = ⟨Φi(R c)|ρW (R c, P, t )|Φj(R c)⟩

Ψ(r, R, t ) =

k ,s

(28)

or Ψ(r, R, t ) =

∑ Ck(s)(t )gk(R|p(ks), q(ks))Φs(r|q(ks)) k ,s

(29)

(s) where gk(R|p(s) k , qk ) are frozen-width Gaussians, which are parametrized by the set of time-dependent position {q(s) k } and (s) momenta {p(s) k } functions of time, and Ck are linear time(s) dependent coefficients. All time-dependent quantities (C(s) k , pk , (s) qk ) are evaluated based on equations arising from application of the TDVP, ⟨δΨ|Ĥ − i∂t|Ψ⟩ = 0.48 The key difference between the two representations is whether the electronic wave functions have a global nuclear coordinate dependence, Φs(r| R), or they are evaluated only at the Gaussian centers q(s) k , Φs(r |q(s) k ). The former choice when used in the TDVP generates NACs because T̂ n is projected onto Φs(r|R). Therefore, unless the double-valuedness of Φs(r|R) is compensated by putting in the GP factors, this form fails to produce the destructive interference pattern in low-energy dynamics.48 A possible solution for this representation has been suggested by introducing a local time-dependent diabatization.46 On the other hand, if one uses eq 29, application of the TDVP does ̂ not lead to NACs because Φs(r |q(s) k ) are not affected by the Tn operator even though they are obtained solving the electronic structure problem. We refer to this problem-free representation as a moving crude adiabatic (MCA) basis.48 In the MCA basis, each nuclear Gaussian carries its own electronic wave function. The latter gathers phase information via continuous evolution of the qk(s) parameters. The accumulated phase of Φs(r |q(s) k ) is a sum of dynamic and geometric parts. Numerical simulations of a symmetric 2DLVC model in the MCA basis confirmed the appearance of the nodal line in low-energy dynamics as in Figure 6.48

(26)

the projected density matrix elements are single-valued, because the double-valuedness of bra and ket components compensate each other. Both of these features of WA are not present in AW. Since the projection onto adiabatic states is done first in AW, it creates double-valued matrix functions of nuclear coordinates Xij(R, R′, t ) = ⟨Φi(R)|ρ(R, R′, t )|Φj(R′)⟩

∑ Ck(s)(t )gk(R|p(ks), q(ks))Φs(r|R)

(27)

In AW, the WT is applied on Xij, which is mathematically equivalent to a Fourier transformation of double-valued functions, a mathematically ill-defined procedure. Thus, 1791



■

CONCLUSIONS

■

ABBREVIATIONS

■

REFERENCES

AR, adiabatic representation; AW, Adiabatic-then-Wigner; BH, Born−Huang; BMA, bis(methylene) adamantyl; BO, Born− Oppenheimer; BOA, Born−Oppenheimer approximation; CI, conical intersection; DBOC, diagonal Born−Oppenheimer correction; DOF, degrees of freedom; FSSH, fewest-switches surface hopping; GP, geometric phase; JT, Jahn−Teller; LVC, linear vibronic coupling; MCA, moving crude adiabatic; MQC, mixed quantum-classical; NAC, nonadiabatic coupling; QCL, quantum-classical Liouville; PSE, potential energy surface; RI, resolution of the identity; TDSE, time-dependent Schrödinger equation; TDVP, time-dependent variational principle; WA, Wigner-then-Adiabatic; WT, Wigner transform

Although the GP is an entity arising only in the adiabatic representation, its consideration provides great qualitative insights in nonadiabatic dynamics near CIs and facilitates development of accurate simulation techniques for a quantitative description of nonadiabatic dynamics. Dependence of the GP phenomenon on the chosen representation cannot be disregarded; thus, the GP absence in the recently revived exact-factorization representation50,51 can be compared to the GP disappearance in the diabatic representation and is not very surprising. This does not reduce the importance of the GP in the widely used adiabatic representation, where it can be considered as a molecular analogue of topological phases in condensed matter physics of topological insulators.52 Indeed, the GP is a topological phase because its value does not depend on a path shape encircling a CI seam. However, the importance of GP effects not only depends on the CI topology but also is strongly affected by the topography of PESs. From a practical point of view, in this Account, we provided criteria when the GP is important for nonadiabatic simulations and what methods can be used to account for it properly. Fundamentally, understanding GP effects provides an opportunity for an unambiguous detection of CIs by observing associated GP signatures in the nuclear density, which are unavailable in other cases such as avoided crossings.

■

Article

(1) Akimov, A. V.; Prezhdo, O. V. Large-Scale Computations in Chemistry: A Bird’s Eye View of a Vibrant Field. Chem. Rev. 2015, 115, 5797−5890. (2) Subotnik, J. E.; Alguire, E. C.; Ou, Q.; Landry, B. R.; Fatehi, S. The Requisite Electronic Structure Theory To Describe Photoexcited Nonadiabatic Dynamics: Nonadiabatic Derivative Couplings and Diabatic Electronic Couplings. Acc. Chem. Res. 2015, 48, 1340−1350. (3) Tully, J. C. Perspective: Nonadiabatic dynamics theory. J. Chem. Phys. 2012, 137, 22A301. (4) Migani, A.; Olivucci, M. In Conical Intersection Electronic Structure, Dynamics and Spectroscopy; Domcke, W., Yarkony, D. R., Köppel, H., Eds.; World Scientific: Hackensack, NJ, 2004; p 271. (5) Domcke, W.; Yarkony, D.; Köppel, H. Conical Intersections: Electronic Structure, Dynamics & Spectroscopy; Advanced series in physical chemistry; World Scientific: Singapore, 2004. (6) Yarkony, D. R. Conical Intersections: Diabolical and Often Misunderstood. Acc. Chem. Res. 1998, 31, 511. (7) Longuet-Higgins, H. C.; Opik, U.; Pryce, M. H. L.; Sack, R. A. Studies of the Jahn−Teller Effect. II. The Dynamical Problem. Proc. R. Soc. London, Ser. A 1958, 244, 1. (8) Berry, M. V. Quantal Phase Factors Accompanying Adiabatic Changes. Proc. R. Soc. London, Ser. A 1984, 392, 45. (9) Mead, C. A. The geometric phase in molecular systems. Rev. Mod. Phys. 1992, 64, 51−85. (10) Hack, M. D.; Jasper, A. W.; Volobuev, Y. L.; Schwenke, D. W.; Truhlar, D. G. Do Semiclassical Trajectory Theories Provide an Accurate Picture of Radiationless Decay for Systems with Accessible Surface Crossings? J. Phys. Chem. A 2000, 104, 217−232. (11) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; International series of monographs on physics; Clarendon Press: Oxford, 1954. (12) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. Fast Numerical Evaluation of Time-Derivative Non-Adiabatic Couplings for Mixed Quantum-Classical Methods. J. Phys. Chem. Lett. 2015, 6, 4200−4203. (13) Cederbaum, L. In Conical Intersection Electronic Structure, Dynamics and Spectroscopy; Domcke, W., Yarkony, D. R., Köppel, H., Eds.; World Scientific: Hackensack, NJ, 2004; p 3. (14) Truhlar, D. G.; Mead, A. C. Relative likelihood of encountering conical intersections and avoided intersections on the potential energy surfaces of polyatomic molecules. Phys. Rev. A: At., Mol., Opt. Phys. 2003, 68, 032501. (15) Mead, C. A.; Truhlar, D. G. On the determination of Born− Oppenheimer nuclear motion wave functions including complications due to conical intersections and identical nuclei. J. Chem. Phys. 1979, 70, 2284. (16) Joubert-Doriol, L.; Ryabinkin, I. G.; Izmaylov, A. F. Geometric phase effects in low-energy dynamics near conical intersections: A study of the multidimensional linear vibronic coupling model. J. Chem. Phys. 2013, 139, 234103. (17) Ryabinkin, I. G.; Joubert-Doriol, L.; Izmaylov, A. F. When do we need to account for the geometric phase in excited state dynamics? J. Chem. Phys. 2014, 140, 214116.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Artur F. Izmaylov: 0000-0001-8035-6020 Funding

This work was supported by NSERC and Sloan Fellowship. Notes

The authors declare no competing financial interest. Biographies Ilya Ryabinkin earned his Ph.D. in physical chemistry from Moscow State University in 2006. As a postdoctoral fellow, he worked on fundamentals of density functional theory with Viktor Staroverov and solved the long-standing problem of constructing exchange-correlation potentials from wavefunctions. In 2012, he joined Izmaylov’s group and spearheaded the project on the role of the geometric phase in nonadiabatic dynamics. Loı ̈c Joubert-Doriol obtained his Ph.D. in 2012 from the University of Montpellier under supervision of Fabien Gatti developing models for nonadiabatic quantum dynamics of molecules interacting with UV/vis light. As a Marie Curie postdoctoral fellow, he worked with Artur Izmaylov and Massimo Olivucci on nonadiabatic quantum effects in biological systems. Since 2016, he develops quantum variational approaches for nonadiabatic dynamics in Izmaylov’s group. Artur Izmaylov is an Assistant Professor at the University of Toronto. He received a Ph.D. with Gustavo Scuseria at Rice University in 2008 and worked as a postdoctoral fellow with John Tully (Yale University) and Michael Frisch (Gaussian Inc.). The main efforts of his group are directed toward developing simulation techniques for quantum molecular dynamics involving multiple electronic states to investigate energy and charge transfer in molecules and nanostructures. 1792


Article


(42) Broeckhove, J.; Lathouwers, L.; Kesteloot, E.; Van Leuven, P. On the equivalence of time-dependent variational principles. Chem. Phys. Lett. 1988, 149, 547−550. (43) Ben-Nun, M.; Martinez, T. J. Ab Initio Quantum Molecular Dynamics. Adv. Chem. Phys. 2002, 121, 439−512. (44) 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. (45) Worth, G. A.; Robb, M. A.; Burghardt, I. A Novel Algorithm for Non-Adiabatic Direct Dynamics Using Variational Gaussian Wavepackets. Faraday Discuss. 2004, 127, 307−323. (46) Meek, G. A.; Levine, B. G. The Best of Both Reps−Diabatized Gaussians on Adiabatic Surfaces. J. Chem. Phys. 2016, 145, 184103. (47) Izmaylov, A. F. Perturbative wave-packet spawning procedure for non-adiabatic dynamics in diabatic representation. J. Chem. Phys. 2013, 138, 104115. (48) Joubert-Doriol, L.; Sivasubramanium, J.; Ryabinkin, I. G.; Izmaylov, A. F. Topologically Correct Quantum Nonadiabatic Formalism for On-the-Fly Dynamics. J. Phys. Chem. Lett. 2017, 8, 452−456. (49) Izmaylov, A. F.; Joubert-Doriol, L. Quantum Nonadiabatic Cloning of Entangled Coherent States. J. Phys. Chem. Lett. 2017, 8, 1793−1797. (50) Min, S. K.; Abedi, A.; Kim, K. S.; Gross, E. K. U. Is the Molecular Berry Phase an Artifact of the Born-Oppenheimer Approximation? Phys. Rev. Lett. 2014, 113, 263004. (51) Curchod, B. F. E.; Agostini, F. On the Dynamics through a Conical Intersection. J. Phys. Chem. Lett. 2017, 8, 831−837. (52) Moore, J. E. The birth of topological insulators. Nature 2010, 464, 194−198.

(18) Ryabinkin, I. G.; Hsieh, C.-Y.; Kapral, R.; Izmaylov, A. F. Analysis of geometric phase effects in the quantum-classical Liouville formalism. J. Chem. Phys. 2014, 140, 084104. (19) Gherib, R.; Ryabinkin, I. G.; Izmaylov, A. F. Why do Mixed Quantum-Classical Methods Describe Short-Time Dynamics through Conical Intersections So Well? Analysis of Geometric Phase Effects. J. Chem. Theory Comput. 2015, 11, 1375. (20) Kendrick, B. K. In Conical Intersections. Electronic Structure, Dynamics and Spectroscopy; Domcke, W., Yarkony, D. R., Köppel, H., Eds.; Advanced Series in Physical Chemistry; World Scientific: Singapore, 2003; Vol. 15, Chapter 12, pp 521−553. (21) Juanes-Marcos, J. C.; Althorpe, S. C.; Wrede, E. Theoretical study of geometric phase effects in the hydrogen-exchange reaction. Science 2005, 309, 1227. (22) Hazra, J.; Balakrishnan, N.; Kendrick, B. K. The geometric phase controls ultracold chemistry. Nat. Commun. 2015, 6, 7918. (23) Mead, C. A.; Truhlar, D. G. Conditions for the definition of a strictly diabatic electronic basis for molecular systems. J. Chem. Phys. 1982, 77, 6090−6098. (24) Handy, N. C.; Lee, A. M. The adiabatic approximation. Chem. Phys. Lett. 1996, 252, 425. (25) Meek, G. A.; Levine, B. G. Wave Function Continuity and the Diagonal Born-Oppenheimer Correction at Conical Intersections. J. Chem. Phys. 2016, 144, 184109. (26) Ryabinkin, I. G.; Izmaylov, A. F. Geometric Phase Effects in Dynamics Near Conical Intersections: Symmetry Breaking and Spatial Localization. Phys. Rev. Lett. 2013, 111, 220406. (27) Xie, C.; Kendrick, B. K.; Yarkony, D. R.; Guo, H. Constructive and Destructive Interference in Nonadiabatic Tunneling via Conical Intersections. J. Chem. Theory Comput. 2017, 13, 1902−1910. (28) Kelly, A.; Kapral, R. Quantum-classical description of environmental effects on electronic dynamics at co nical intersections. J. Chem. Phys. 2010, 133, 084502. (29) Endicott, J. S.; Joubert-Doriol, L.; Izmaylov, A. F. A perturbative formalism for electronic transitions through conical intersections in a fully quadratic vibronic model. J. Chem. Phys. 2014, 141, 034104. (30) Joubert-Doriol, L.; Izmaylov, A. Molecular “topological insulators”: A case study of electron transfer in the bis(methylene) adamantyl carbocation. Chem. Commun. 2017, 53, 7365−7368. (31) 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. (32) Yarkony, D. R. Marching along ridges. An extrapolatable approach to locating conical intersections. Faraday Discuss. 2004, 127, 325−336. (33) Malbon, C. L.; Zhu, X.; Guo, H.; Yarkony, D. R. On the incorporation of the geometric phase in general single potential energy surface dynamics: A removable approximation to ab initio data. J. Chem. Phys. 2016, 145, 234111. (34) Tully, J. C. Mixed quantum-classical dynamics. Faraday Discuss. 1998, 110, 407−419. (35) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011−2015. J. Phys. Chem. Lett. 2016, 7, 2100− 2112. (36) 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. (37) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061. (38) Kapral, R. Progress in the theory of mixed quantum-classical dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129−157. (39) Kapral, R.; Ciccotti, G. Mixed quantum-classical dynamics. J. Chem. Phys. 1999, 110, 8919. (40) Horenko, I.; Salzmann, C.; Schmidt, B.; Schütte, C. Quantumclassical Liouville approach to molecular dynamics: Surface hopping Gaussian phase-space packets. J. Chem. Phys. 2002, 117, 11075. (41) Ando, K.; Santer, M. Mixed quantum-classical Liouville molecular dynamics without momentum jump. J. Chem. Phys. 2003, 118, 10399. 1793


Geometric Phase Effects in Nonadiabatic Dynamics near Conical

Recommend Documents