A Tractable Numerical Model for Exploring ... - ACS Publications

Article pubs.acs.org/jchemeduc

A Tractable Numerical Model for Exploring Nonadiabatic Quantum Dynamics Evan Camrud and Daniel B. Turner* Department of Chemistry, New York University, 100 Washington Square East, New York, New York 10003, United States S Supporting Information *

ABSTRACT: Numerous computational and spectroscopic studies have demonstrated the decisive role played by nonadiabatic coupling in photochemical reactions. Nonadiabatic coupling drives photochemistry when potential energy surfaces are nearly degenerate at avoided crossings or truly degenerate at unavoided crossings. The dynamics induced by nonadiabatic coupling are challenging to comprehend; here we describe a versatile one-dimensional model that is numerically tractable and illustrates fundamental aspects of nonadiabatic quantum dynamics. This model reinforces and builds on concepts taught in graduate-level quantum mechanics and therefore should be accessible to advanced undergraduate and graduate students. We use the model to demonstrate how the local topography of an unavoided crossing can affect a photochemical quantum yield. KEYWORDS: Upper-Division Undergraduate, Graduate Education/Research, Physical Chemistry, Computer-Based Learning, Molecular Properties/Structure, Photochemistry, Quantum Chemistry

■

esoteric1,2 but are now theorized to be ubiquitous in chemical reactions involving excited states. Indeed, conical intersections have become a major focus in photochemistry and related fields since the 1990s,3−12 and they are now thought to play a key role in processes as diverse as charge transfer,13−15 human vision,16−18 DNA photoprotection,19 phototropism,20 and photosynthesis.21 In regions of large nonadiabatic coupling, the nuclear and electronic motions become coupled. This results in a collapse of the Born−Oppenheimer approximation and thereby makes computations and simulations extremely challenging. Much like a conventional transition state, avoided and unavoided crossings are transient molecular configurations created during a chemical reaction.11 In these regions of strong nonadiabatic coupling, the adiabatic states that arise from the Born− Oppenheimer approximation change character quickly, but in contrast, diabatic states by definition do not depend on the nuclear coordinate. Because of their importance in physical chemistry, the topics of nonadiabatic coupling, the adiabatic approximation, and conical intersections are now described in textbooks22−24 and taught in advanced undergraduate and graduate courses on quantum mechanics. Impressively, synthetic chemistry25 and photochemistry26 textbooks introduce the conical intersection as a mechanistic concept. Nevertheless, these topics confuse many students. Part of the challenge in understanding nonadiabatic dynamics undoubtedly

INTRODUCTION Scientists envision molecules and chemical reactions through potential energy surfaces, which are representations of how atomic motions modulate the energies of molecular orbitals. The Born−Oppenheimer approximation leads to a set of solutions known as adiabatic potential energy surfaces, which are usually adequate for chemistry restricted to ground states because the reaction proceeds along one electronic eigenstate that changes energy and character as the atoms rearrange. However, the Born−Oppenheimer approximation cannot be used to describe photochemical transformations because it ignores nonadiabatic coupling. Photochemical reactions inherently involve multiple adiabatic potential energy surfaces, and nonadiabatic coupling allows the system to switch between these surfaces and thereby mediate photochemistry. Figure 1 illustrates that the strength of the electronic * ) strongly component of the nonadiabatic coupling (VNADC depends on the energetic difference between the two potential energy surfaces (ΔE): V*NADC ∝ 1/ΔE. Because the adiabatic potentials change energy as a function of the nuclear coordinate, X, the electronic component of the nonadiabatic * (X). The critical photocoupling is position-dependent, VNADC chemical transformation occurs at the molecular configuration that brings two adiabatic potentials close in energy. These configurations are known as avoided crossings. In certain systems the two adiabatic potentials are truly degenerate at specific molecular configurations, and the nonadiabatic coupling becomes infinite. Such unavoided crossings are known as conical intersections when depicted along two coordinates because the topology is that of a cone. These features were once considered © XXXX American Chemical Society and Division of Chemical Education, Inc.

Received: November 8, 2016 Revised: February 21, 2017

A

DOI: 10.1021/acs.jchemed.6b00861 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Article

Figure 1. Photoexcitation (green) promotes the system from the ground adiabatic state (black) to the excited adiabatic state (gray). Nonadiabatic * is finite and smooth. (right) For an unavoided coupling returns the system to the ground adiabatic state. (left) For an avoided crossing (AC), VNADC crossing (UC), V*NADC is infinite at the point where the two surfaces are degenerate. In two dimensions, these points of degeneracy are known as conical intersections.

is that scientists use advanced ab initio techniques to simulate and study specific molecules. A tractable numerical model that eschews specificity in favor of generality and accessibility would ease the illustration and comprehension of nonadiabatic coupling. The model presented herein is most appropriate for students who have knowledge of the harmonic oscillator, twolevel system, bra−ket notation, time dependence, and basis transformations. Prior use of software environments such as Python or MATLAB is not required, although it would enable a student to grasp more quickly the scientific principles. In the model described here, by varying just three parameters, one can tune the topography of the unavoided crossing from peaked to sloped27−29 or induce an avoided crossing. In what follows, we begin by describing the model and the numerical method, and then we discuss the results for varying topographies. We conclude by describing avenues for expanding this model.

■

Figure 2. (top) Diabatic potentials |Gd⟩ and |Ed⟩ (solid blue and red lines) are coupled via Vd(X) to yield adiabatic potentials |Ga⟩ and |Ea⟩ (dashed cyan and yellow lines). At the position of the avoided crossing, XAC, the nonadiabatic coupling, VNADC(X), is nonzero. (bottom) The diabatic states do not change character as a function of the nuclear coordinate X. The adiabatic states change character significantly near the avoided crossing.

METHOD The model depicted in Figure 2 begins with a Hilbert space formed by two harmonic oscillators representing diabatic electronic states of a molecule. Diabatic states have distinct chemical identities, but their energetic ordering may change as a function of nuclear displacement, as described by van Voorhis in an exceptionally clear manner using the example of NaCl.8 In the top panel of Figure 2, the solid red and solid blue lines are the energies of the diabatic states, which form a complete orthonormal basis for the electronic system. At large values of X, the ground diabatic state, |Gd⟩, has a higher energy than the excited diabatic state, |Ed⟩. In contrast, the adiabatic states, |Ea⟩ and |Ga⟩, do not have distinct chemical identities, but their energetic ordering is always preserved. The adiabatic states also form a complete orthonormal basis for the electronic system. Chemists are intimately familiar with ground adiabatic potential energy surfaces, as these are a key component of the ubiquitous “reaction coordinate diagrams”. Mathematically, the change in chemical identity means that each adiabatic state becomes a linear combination of the diabatic states in the region of the avoided crossing. The amount of each diabatic state used to create each adiabatic state changes as a function of the nuclear coordinate: for example, |Ga⟩ = c1(X)|Gd⟩ + c2(X)|Ed⟩ and |Ea⟩ = d1(X)|Gd⟩ + d2(X)|Ed⟩. This is illustrated in the bottom panel of Figure 2, where the dashed and dotted traces depict d2(X) and c2(X), respectively.

A key aspect of this model is that the coupling between the diabatic potentialsnormally a constant in electron-transfer models such as Marcus theoryis a function of the nuclear coordinate X. We couple the two diabatic states to yield analytic expressions for the adiabatic potentials. To solve for the vibrational eigenstates of the adiabatic potentials, we use a numerical approximation technique taught in graduate quantum courses known as the discrete variable representation (DVR).30,31 Using the adiabatic potentials and their vibrational eigenfunctions, we compute the nonadiabatic coupling and then propagate a wavepacket that would result from optical excitation in the Franck−Condon region. The workflow is as follows: 1. postulate diabatic states and their coupling; 2. derive the adiabatic states; 3. solve numerically for the adiabatic vibrational eigenstates using DVR; 4. compute the nonadiabatic coupling matrix elements; 5. construct a total Hamiltonian; B



Article

⎛ ω 2X ⎞ 0 ̂ e (X ̂ ) = 1 ⎜ Vadia (2X̂ − X 0) + Ee⎟ 2⎝ 2 ⎠

6. propagate a Gaussian wavepacket under the total Hamiltonian. The subsections below detail each of these steps. Deriving the Adiabatic Potentials and Computing Their Vibrational Eigenstates

+

First we postulate diabatic potentials and then derive analytic expressions for the adiabatic potentials. We then build the full Hamiltonian by adding kinetic energy. The key elements in the Hamiltonian are the two diabatic potential energy surfaces, V̂ gdiab and V̂ ediab, each described as a harmonic oscillator, and their coupling, V̂ d. These operators can be applied to position basis states X̂ . The excited state is displaced both vertically and horizontally from the ground state: ̂ g (X ̂ ) Vdiab

1 2 = ω 2 X̂ 2

e

̂ (X ̂ ) = Vdiab

1 2 ̂ ω (X − X 0)2 + Ee 2

⎡ c1(X ) ⎤ ⎥ |Ga⟩ = ⎢ ⎢⎣ c 2(X )⎥⎦

and ⎡ d1(X ) ⎤ ⎥ |Ea⟩ = ⎢ ⎢⎣ d 2(X )⎥⎦

(1) (2)

As the notation suggests, each adiabatic state is a positiondependent linear combination of the diabatic states. The Supporting Information contains the lengthy expressions for c1(X), c2(X), d1(X), and d2(X). The bottom panel of Figure 2 illustrates how the adiabatic states depend on position in one example of an avoided crossing. One goal of this model is to represent both avoided and unavoided crossings. Most coupling values will generate an avoided crossing. With conventional Condon coupling (μ = 0), an avoided crossing will always be created for any nonzero value of β. In fact, this situation leads to the well-known noncrossing rule for diatomic molecules: When there is just one internal coordinate, two adiabatic states of the same symmetry cannot be degenerate. Under Condon coupling, the only possible unavoided crossing is for the trivial coupling V̂ d = 0; there is no distinction between the diabatic and adiabatic potentials. The non-Condon coupling used here, where V̂ d is linear in X̂ , can produce nontrivial unavoided crossings. To represent an unavoided crossing, we must choose

(3)

To solve the system of coupled equations, we construct a 2 × 2 matrix describing purely the potential energy: ⎡ V ̂ g (X ̂ ) V ̂ (X ̂ ) ⎤ d ⎥ ⎢ diab ⎢ V̂ (X̂ ) V̂ e (X̂ )⎥ ⎦ ⎣ d diab

where the diabatic potentials are the diagonal matrix elements and the off-diagonal matrix elements are equal to the coupling. This is equivalent to constructing a two-level Hilbert space for the diabatic potentials, where the diabatic states are

⎡1 ⎤ |Gd⟩ = ⎢ ⎥ ⎣0⎦

⎛X E ⎞ β = −μ⎜ 0 + 2 e ⎟ ω X0 ⎠ ⎝ 2

and ⎡0⎤ |Ed⟩ = ⎢ ⎥ ⎣1 ⎦

⎡ c ( X ̂ ) d ( X ̂ ) ⎤ ⎡ V ̂ g ( X ̂ ) V ̂ ( X ̂ ) ⎤ ⎡ c (X ̂ ) c (X ̂ ) ⎤ d 1 2 ⎥⎢ 1 ⎢ 1 ⎥⎢ diab ⎥ ⎢⎣ c (X̂ ) d (X̂ )⎥⎦⎢ V̂ (X̂ ) V̂ e (X̂ )⎥⎢⎣ d (X̂ ) d (X̂ )⎥⎦ ⎣ d ⎦ 1 2 2 2 diab

(4)

to yield general analytic expressions for the adiabatic potential energy surfaces: ⎛ ω 2X ⎞ 0 ̂ g (X ̂ ) = 1 ⎜ Vadia (2X̂ − X 0) + Ee⎟ 2⎝ 2 ⎠ −

⎛ ω 2X ⎞2 1 0 (μX̂ + β) + ⎜ (2X̂ − X 0) − Ee⎟ 4⎝ 2 ⎠

(6)

which makes the two surfaces degenerate at a single point. Any other value of β will produce an avoided crossing, wherein the two potentials are not degenerate. Although we have derived analytic expressions for the adiabatic potentials, we cannot solve for the vibrational eigenfunctions analytically except in certain cases such as zero coupling or perfect symmetry along the position coordinate. The discrete variable representation (DVR) method allows us to determine the vibrational eigenfunctions of an arbitrary potential. The method, detailed in the Supporting Information, works by diagonalizing the position matrix of a known set of basis functions, in this case the analytic solutions to the quantum harmonic oscillator. This diagonalized position matrix creates what are known as DVR “grid points”, which one uses to evaluate functions such as the potential energy surfaces because operations involving diagonal matrices are straightforward. Once evaluated, the function can be transformed back to the initial basis, which in this case is the quantum harmonic oscillator. This process allows for the creation of a Hamiltonian of an arbitrary potential energy function within the basis of the harmonic oscillator. By diagonalizing the Hamiltonian, one

As diabatic states, the elements of these vectors do not depend on the nuclear coordinate X̂ . Then we diagonalize the system:

⎡ V ̂ g (X ̂ ) 0 ⎤ adia ⎥ =⎢ e ⎢ 0 ⎥ ̂ ̂ V ( X ) ⎣ ⎦ adia

(5b)

The diagonalization process yields the eigenvalues presented in eqs 5 as well as the adiabatic eigenfunctions

where ω is the normal vibrational frequency, Ee is the energetic shift of the excited state, and X0 is its displacement along the position coordinate X̂ . We define the linear non-Condon electronic coupling as Vd̂ (X̂ ) = μX̂ + β

2 ⎞2 1 ⎛ ω X0 2 ̂ ̂ (μX + β) + ⎜ (2X − X 0) − Ee⎟ 4⎝ 2 ⎠

2

(5a) C



Article

eigenfunctions of the two adiabatic potentials for −0.5 < X < 1.5.

obtains the eigenvalues and eigenvectors needed to numerically approximate the solutions to the arbitrary potential energy. Each eigenvector is a linear combination of the original basis functions. To compute the vibrational eigenfunctions, we use the DVR of the harmonic oscillator (HO). We perform DVR separately for the ground and excited states and hold the transformation matrices for later use. Following the DVR procedure, we truncate the harmonic oscillator position matrix, X̂ HO, at a large value k. We denote the transformation matrix and eigenvalues formed by diagonalization as ̂ † XHO ̂ ̂ UDVR ̂ UDVR = XDVR

(7)

After evaluating the adiabatic potentials, V̂ adia(X̂ DVR), we transform back to the harmonic oscillator position basis: ̂ Vadia ̂ g (XDVR ̂ † = Vadia ̂ g (XHO ̂ )UDVR ̂ ) UDVR

(8a)

̂ Vadia ̂ e (XDVR ̂ † = Vadia ̂ e (XHO ̂ )UDVR ̂ ) UDVR

(8b)

Subsequent addition of these matrices to the harmonic oscillator kinetic energy matrix produces the adiabatic Hamiltonians g ̂ g (XHO ̂ ) + 1 PHO ̂ 2 Ĥ adia = Vadia 2

Figure 3. Selected vibrational eigenfunctions Ψgm(X) and Ψen(X) of the ground and excited adiabatic potentials, respectively, for a model of a conical interaction. Shifted energy scales aid in viewing. (bottom) Despite the pronounced cusp in the ground adiabatic potential, the vibrational eigenfunctions are smooth. (top) Nearly all of eigenenergies of the excited adiabatic potential are higher than the potential energy of the cusp.

(9a)

e ̂ e (XHO ̂ ) + 1 PHO ̂ 2 Ĥ adia = Vadia 2

(9b) 1

2

̂ , is one-half of the where the kinetic-energy operator, 2 PHO square of the momentum operator in the harmonic oscillator basis (the mass is set to 1 in this unit system). Diagonalizing each Hamiltonian produces eigenvalues corresponding to the vibrational energy levels and the eigenvectors that are used to build the vibrational eigenfunctions by taking linear combinations of the eigenfunctions of the quantum harmonic oscillator. Each eigenvector is a linear combination that is an approximate solution to the Schrödinger equation with its corresponding eigenvalue. The solutions for the ground and excited adiabatic states satisfy g Ĥ adiaϕmg = Emg ϕmg

Having derived the adiabatic potentials and solved numerically for their vibrational eigenfunctions, we can then form a Hilbert space by construction of ⎛ k ⎞ Ĥ adia = |Ga⟩⟨Ga| ⊗ ⎜⎜ ∑ Emg |m⟩⟨m|⎟⎟ + |Ea⟩⟨Ea| ⎝m=0 ⎠ ⎛ k ⎞ ⊗ ⎜⎜ ∑ Ene|n⟩⟨n|⎟⎟ ⎝n=0 ⎠

The size of this adiabatic Hamiltonian matrix is 2k × 2k, and it spans a Hilbert space sufficient to include nonadiabatic coupling. The adiabatic Hamiltonian matrix is diagonal.

(10a)

k

Ψmg (X ) =

∑ ϕmg(i)ψi(X ) i=1

Computing the Nonadiabatic Coupling

(10b)

The physical basis behind the Born−Oppenheimer approximation is that the nuclei are much more massive than the electrons. Thus, the electrons are more mobile than the nuclei, and at a given moment, the electrons’ positions depend on the positions of the nuclei parametrically. In quantum chemistry, this is implemented by solving a Schrödinger equation that involves only the electronic degrees of freedom, where the Hamiltonianand its eigenstates and eigenvaluesdepend on the nuclear configuration. This procedure neglects the nuclear momentum, and in fact the Born−Oppenheimer approximation is valid whenever the electronic states change character slowly along the nuclear coordinate. At any given moment, the electrons are in an eigenstate of the electronic Hamiltonian; the nuclei do not move fast enough to violate this adiabatic assumption. Thus, chemists generally visualize dynamics in the adiabatic basis, because chemical reactions usually proceed along the ground adiabatic state. The Born−Oppenheimer

and e Ĥ adiaϕne = Eneϕne

(11a)

k

Ψ en(X )

=

∑ ϕne(i)ψi(X ) i=1

(12)

(11b)

ϕjn(i)

respectively, where is the ith element of the eigenvector for the jth adiabatic electronic state and ψi is the ith analytic solution to the quantum harmonic oscillator. The mth vibrational eigenfunction of the ground adiabatic state is Ψgm(X), and the nth vibrational eigenfunction of the excited adiabatic state is Ψen(X). In these expressions, the variable X is distinct from the operator X̂ in that it can be a new coordinate vector used only to visualizenot computethe eigenfunctions and wavepackets. Figure 3 shows selected vibrational D



Article

cne(t0) = ⟨Ψ en(X )|g (X )⟩

approximation is invalid when the electronic wave functions change character quickly as a function of nuclear motion. The distinction between “quickly” and “slowly” depends on the slopes of the adiabatic potentials as well as their energetic separation. To determine when the Born−Oppenheimer approximation breaks down, we must compute the nonadiabatic coupling, V̂ NADC. Here we use the conventional expression24 NADC

Vm̂ , n

=

∫ dX Ψen(X )∇Ψmg (X )·⟨Ea(X )|∇|Ga(X )⟩

which yields the coefficients of the excited-state vibrational eigenfunctions at time t0. In this Hilbert space, the ground-state vibrational eigenfunctions must also be included even though their initial coefficients are zero (i.e., cgn(t0) = 0). To propagate the wavepacket, we must compute the coefficients at each time step, cji(t), using the full Hamiltonian that includes the nonadiabatic coupling, Ĥ tot adia, as written in eq 15. Diagonalization produces a new set of eigenvalues, Ê *, and their corresponding eigenvectors that function as a basis transformation matrix, Û :

(13)

where |Ga⟩ and |Ea⟩ are the electronic eigenfunctions first introduced in eq 4 and detailed in the Supporting Information. We have neglected the higher-order term that depends on ⟨Ea(X)|∇2|Ga(X)⟩. The electronic component of the nonadiabatic coupling is precisely what was depicted in Figure 1, V*NADC(X) = ⟨Ea(X)|∇|Ga(X)⟩. We compute and use the lengthy analytic expression for VNADC * (X) (see the Supporting Information). Since we use finite matrices, the integral in eq 13 becomes a sum over the y subdivisions of the nuclear coordinate vector, which, as noted above, is distinct from the number of basis vectors used in the DVR procedure. In discrete form, we can write the nonadiabatic coupling as NADC

Vm̂ , n

† tot Û Ĥ adiaÛ = E*̂

* (X i ) ∑ ΔXi Ψen(Xi)∇Ψmg (Xi)·V NADC i=1

tot Ĥ adia|Ul⟩ = El̂ *|Ul⟩

= Ĥ adia +

NADC Vm̂ , n

⎡ c1g(t0) ⎤ ⎡ c1*(t0) ⎤ ⎥ ⎢ ⎥ ⎢ ⎢ ⋮ ⎥ ⎢ ⋮ ⎥ ⎥ ⎢ g ⎥ ⎢ * ⎥ ⎢ ( ) c ( t ) c t ⎢ ⎥ k 0 k 0 † ⎥ Û ⎢ e ⎥=⎢ ⎢ c1 (t0) ⎥ ⎢ ck*+ 1(t0) ⎥ ⎥ ⎢ ⎥ ⎢ ⎢ ⋮ ⎥ ⎢ ⋮ ⎥ ⎥ ⎢ c e(t ) ⎥ ⎢ * ⎣ k 0 ⎦ ⎢⎣ c 2k(t0) ⎥⎦

(14)

(15)

⎡ * −iE1*(t 0 +Δt ) ⎤ ⎤ ⎡ ⎢ c1 (t0)e ⎥ ⎢ c1*(t0 + Δt ) ⎥ ⎢ ⎥ ⎢ ⋮ ⎥ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ * −iEk*(t 0 +Δt ) ⎥ * ⎥ ⎢ + Δ ( ) c ( t )e c t t ⎢ k 0 ⎥ ⎢ k 0 ⎥ = ⎢ −iEk*+ 1(t 0 +Δt ) ⎥ ⎢ ck*+ 1(t0 + Δt ) ⎥ * c ( t )e ⎢ k+1 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ * *k(t0 + Δt ) ⎥⎦ −iE 2*k (t 0 +Δt ) ⎥ c ⎢ 2 ⎣ ⎣ c 2k(t0)e ⎦

Creating and Propagating a Wavepacket

Because the purpose of this exercise is to represent quantum dynamics due to nonadiabatic coupling in, for example, spectroscopic measurements, we choose to evaluate the result of lifting a wavepacket in the Franck−Condon region to the excited adiabatic state. It is therefore applicable to produce a wavepacket in the excited state directly over the minimum of the ground state, here taken as X = 0. We begin with the general equation for a time-dependent wavepacket created from a linear combination of the ground and excited eigenfunctions:

i=1

⎡ c *(t + Δt ) ⎤ ⎡ c g(t + Δt ) ⎤ ⎢ 1 0 ⎥ ⎢ 1 0 ⎥ ⎢ ⎥ ⎢ ⋮ ⎥ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ ck*(t0 + Δt ) ⎥ ⎢ ckg(t0 + Δt ) ⎥ ⎥=⎢ Û ⎢ ⎥ ⎢ ck*+ 1(t0 + Δt ) ⎥ ⎢ c1e(t0 + Δt ) ⎥ ⎢ ⎥ ⎢ ⎥ ⋮ ⎢ ⎥ ⎢ ⋮ ⎥ ⎢ ⎥ ⎢ e ⎥ + Δ c ( t t ) * ⎦ ⎢⎣ c 2k(t0 + Δt ) ⎥⎦ ⎣ k 0

(16)

We use an initial wavepacket that has the form of a Gaussian function because such a wavepacket will vary smoothly with time in a harmonic potential. This means that a Gaussian wavepacket will remain in the shape of a Gaussian at all times until the wavepacket encounters anharmonic regions of the potential.24 To create a Gaussian wavepacket induced by excitation in the Franck−Condon region, we first define a Gaussian function centered at X = 0 with width σ, g (X ) = e − X

2

/σ 2

(22)

Then we transform the coefficients back to the adiabatic basis,

k

∑ [cig(t )Ψig(X ) + cie(t )Ψie(X )]

(21)

and the phase factorsbased on the corresponding eigenvalues32 and the time step Δtare included with the coefficients:

We will discuss the diagonalization of this Hamiltonian and its use to propagate the wavepacket in the next subsection.

Ψ(X , t ) =

(20)

where |Ul⟩ is the lth column of Û . The coefficients from the initial wavepacket are transformed into this new basis,

Hence, it is helpful at this point to use a finely spaced X vector. The quantity V̂ NADC gives the coupling value between states m,n |Ea; n⟩ and |Ga; m⟩. We can produce a final Hamiltonian that will couple the necessary elements of the ground and excited adiabatic states: tot Ĥ adia

(19)

or in vector form,

y

=

(18)

(23)

and we use eq 16 to produce the wavepacket at time step t0 + Δt. It is this mixing of adiabatic coefficients by the nonadiabatic coupling that induces nontrivial quantum dynamics.

(17)

Computational Methods

and then take the inner product of this Gaussian with each excited-state eigenfunction Ψen,

We performed simulations in MATLAB 2013b because this software is available through many universities for student use. E



Article

We first evaluate a set of parameters (Ee = 500/18, X0 = 20/ 3, β = −μ(X0/2 + Ee/(ω2X0)) + 36.549) that should induce a large avoided crossing. Indeed, the adiabatic states presented in Figure 4 confirm that the avoided crossing involves an energetic

The three parameters that allow for exploration of possible topologies and topographies are Ee, μ, and X0. For most of the simulations, the parameters were Ee = 500/18 (0.076 hartree), X0 = 20/3 (0.6449 Å), ω = 1 (0.0027 hartree), σ = 1 (0.097 Å), and μ = 1/4 (0.0071 hartree/Å). We used 1500 basis states and a minimum of 75 Hermite polynomials to construct the eigenfunctions to make the results smoothly varying. For dynamics simulations, we took 1 fs steps for a total duration of 200 fs. These computations required about 15 min on a personal laptop computer (Intel Core i7 processor), demonstrating the feasibility of the method for students. The quantum harmonic oscillator used for DVR had the potential energy curve equal to the ground-state diabatic potential (eq 1). We implemented three noteworthy protocols to reduce computing time and increase the accuracy of the numerical results: 1. The first regards differentiation during computation of the nonadiabatic coupling. In general, vector differentiation returns a new vector that is reduced in length by one element. Such a vector then cannot be added or multiplied with an undifferentiated vector because the grid points will be distinct. To ameliorate this problem, a vector containing one extra element must be used for differentiation. In this model, we exploit the fact that by using the eigenvectors ϕgm and ϕen, we can compute the vibrational eigenfunctions on any arbitrary grid points. Hence, if the undifferentiated vector Ψen(X) is evaluated on grid points from 0 to 1 in steps of 0.1, we would evaluate the Ψgm(X) vector on grid points −0.05 to 1.05 in steps of 0.1. Differentiation of this Ψgm(X) returns a vector evaluated on the grid points corresponding to Ψen(X). 2. For an unavoided crossing, the adiabatic potentials must be displaced to intersect perfectly at one of the grid points. Because there are a finite number of points, we found it most effective to use a trial and error method to determine the displacement of X0 = 20/3 (0.6449 Å) noted above. 3. An additional concern is that the nonadiabatic coupling reaches a singularity at the point of the unavoided crossing. This renders numerical solutions useless in considering the true coupling at the intersection. Here we manually replace the not-a-number (NaN) element with a value of 20 (0.055 hartree) at the unavoided crossing. This may seem small, but it is many orders of magnitude larger than the nonadiabatic coupling at all other points (∼10−17 Hartree). The most important reason for using this number, however, is that at low coupling, the wavepacket will oscillate normally, as one would imagine without coupling. At very high couplings such as 1000, the wavepacket will reflect at the intersection, just as would occur if there were an extremely large (or infinite) potential barrier. We posit that when the wavepacket tunnels through the coupling potential, it proceeds onto the ground adiabatic state.

Figure 4. Adiabatic potentials and resultant quantum dynamics of a large avoided crossing. The wavepacket is restricted to the upper adiabatic potential (excited state) for the duration of the simulation. No amplitude develops on the lower adiabatic potential (ground state).

separation of about 0.1 hartree, and the initial wavepacket prepared on the excited adiabatic potential remains there throughout the entire duration of the simulation (200 fs). The periodicity of this wavepacket is approximately 55 fs. The wavepacket oscillates over a distance of approximately 0.8 Å, similar in magnitude to nuclear vibrations. The wavepacket decreases in amplitude as it begins to propagate. Even as it peaks again halfway through its period, the amplitude is not as large as it is at the beginning of each period. This is due to the shape of the potential energy surface at this point. The nonCondon coupling forces the surface to pitch higher on one end. When the wavepacket is on the higher end, it must spread its position to compensate for the higher potential. The decreased amplitude at the bottom of the potential is a manifestation of the uncertainty principle, as the wavepacket has more momentum. The increased momentum forces the position uncertainty to increase, which expands the wavepacket and lowers its amplitude. The reverse of this occurs at each end of the wavepacket oscillations, wherein the correspondence principle gives an idea of the wavepacket “slowing down” as if it were a pendulum at its peak. This provides extreme certainty of the position and thus creates a large uncertainty in the momentum. As the potential energy surfaces are brought closer together, the nonadiabatic coupling between the adiabatic states is no longer negligible. Because of this, when a wavepacket nears the avoided crossing, it imparts some of its amplitude to the lower adiabatic state. As shown in Figure 5, for an avoided crossing with a separation of about 0.01 hartree, one can see the slight transfer of the wavepacket onto the ground state. The amplitude of the ground-state wavepacket is so small that it does not appear on this color scale until the fourth pass through the crossing region. The periodicity of the wavepacket has decreased to approximately 45 fs. This is due to the narrowed potential surface relative to the large avoided crossing. The shape of the potential also contributes to the increased

■

RESULTS AND DISCUSSION This model is intended to provide qualitative insights into the quantum dynamics induced by nonadiabatic coupling. As a reference, we begin the discussion with the dynamics that result from avoided crossings. We evaluate two such cases before turning our attention to three types of unavoided crossings. F



Article

Figure 5. Adiabatic potentials and resultant quantum dynamics of a close avoided crossing. The wavepacket is mostly retained on the upper adiabatic potential, although a small amplitude develops on the lower adiabatic potential after about 100 fs, which corresponds to the fourth pass through the crossing region.

oscillation distance, now approximately 0.95 Å. One can also see that over time, the wavepacket will eventually dephase because of its loss of amplitude to the ground state. This happens on a relatively long time scale compared with other events, however. This avoided crossing simulation allows us to comment on the “paradigm of photochemistry of vision” described by Shank and colleagues,33 as coherent wavepacket transfer from one adiabatic potential to another is a signature of a conical intersection. Figure 5 suggests that close avoided crossings can lead to nonzero wavepacket amplitude transfer between adiabatic states, supporting what others have discussed.34,35 It is important to mention, however, that a one-dimensional model does not fully describe the character of adiabatic potentials or conical intersections in molecules. It is unclear how far this model and its results can be interpreted. The avoided crossing depicted in Figure 5 presents an opportunity to enhance the conceptual understanding of the nonadiabatic coupling matrix elements. On the basis of eq 14, each element of the nonadiabatic coupling matrix is the product of three functions: the mth vibrational eigenfunction of the ground adiabatic state (Ψgm(X)), the nth vibrational eigenfunction of the excited adiabatic state (Ψen(X)), and the electronic component (VNADC * (X)), for which we have an analytic expression. For completeness, one of the vibrational eigenfunctions will require differentiation; this derivative will be zero wherever the amplitude is also zero. Conceptually, then, a nonadiabatic coupling matrix element is nonzero only when all three of these functions are simultaneously nonzero over a significant portion of the nuclear coordinate axis. The electronic component is independent of the chosen vibrational eigenfunctions, and it therefore provides a spatial filter. This means that the only portion of a vibrational eigenfunction that is relevant to the nonadiabatic coupling is the portion that is nonzero at the same position value(s) where VNADC * (X) is also nonzero. In Figure 6 we present all of the nonadiabatic coupling matrix elements for the specific conditions of the close avoided crossing in Figure 5. The indices for the vibrational eigenfunctions of the ground (excited) adiabatic potential are listed vertically (horizontally). Generally, the amplitude of the nonadiabatic coupling increases for higher-energy vibrational

Figure 6. (top) Two vibrational eigenfunctions of each adiabatic potential, as marked. The electronic component of the nonadiabatic * (X), sets a spatial filter (red-shaded region) over which coupling, VNADC the amplitude of both vibrational eigenfunctions must be nonzero. (bottom) Selected matrix elements V̂ NADC m,n . The 20 lowest-energy vibrational eigenfuctions of the ground adiabatic state (m < 20) because they have very little amplitude contribute negligibly to V̂ NADC m,n in the red-shaded region.

eigenfunctions in both electronic states. The 20 lowest-energy ground vibrational eigenfunctions contribute negiligibly to the nonadiabatic coupling because these vibrational eigenfunctions have essentially zero amplitude where V*NADC(X) is nonzero. Substantially different dynamics results from unavoided crossings. Figure 7 shows that significant wavepacket amplitude transfers to the ground adiabatic potential even on the first pass

Figure 7. Adiabatic potentials and resultant quantum dynamics of an unavoided crossing. The wavepacket splits upon the first encounter with the degeneracy, and significant amplitude develops on the lower adiabatic potential in just 25 fs. G



Article

through the unavoided crossing. One may notice the effect of the nonadiabatic coupling mentioned earlier: the excited-state wavepacket reflects. This can be seen in the sharp cusp at the unavoided crossing in the position versus time plot of the excited state. Because of this reflection, the wavepacket has a period of just 30 fs. The portion of the wavepacket transferred to the ground adiabatic potential shares a similar period, though it becomes difficult to observe past the first 100 fs, as too much amplitude has been imparted into it. The distance of the total process is certainly larger than 2.0 Å, which is quite interesting. This is likely due to imparting a high-energy wavepacket into a low-energy potential. Figure 8 shows two potential curves that are displaced in space and energy such that the lower adiabatic potential has

Figure 9. Adiabatic potentials and resultant quantum dynamics of a sloped13 unavoided crossing. This topography leads to negligible transfer from the excited state to the ground state.

lower state, supporting the notion that the wavepacket must actually reach the unavoided crossing to transfer. Because the upper adiabat is not displaced as far as other curves, the distance over which the wavepacket propagates is much smaller, a mere 0.45 Å in this case. There are differences and similarities between the 1D unavoided crossings discussed in this work and true conical intersections, which can only be described in two dimensions. A 1D unavoided crossing and a conical intersection both have an electronic degeneracy between two adiabatic potentials, and they also both have an infinite nonadiabatic coupling at one position value. A 1D unavoided crossing arises from nonCondon type coupling, here V̂ d(X̂ ), that is known as a “derivative” coupling in the nomenclature of conical intersections. However, a 1D unavoided crossing does not have a “coupling” coordinate, which is a second spatial coordinate that serves to lift the degeneracy. The model does not include this second coordinate because it would significantly increase the complexity of all aspects of these calculations and make the results more conceptually challenging. As a consequence, the Berry phasethe phase shift that occurs as the wavepacket encounters a conical intersectionis not induced by the 1D unavoided crossing. The 1D unavoided crossing model does adequately represent special cases of conical intersections wherein the tuning coordinate either does not lift the degeneracy or does so only slightly. Examples include some models of electron and energy transfer.

Figure 8. Adiabatic potentials and resultant quantum dynamics of an unavoided crossing meant to mimic the conical intersection proposed for rhodopsin.33 This topography causes the fastest wavepacket transfer from the excited state to the ground state.

two minima. The values X0 = 6.56 (0.6346 Å) and Ee = 2.01 (0.055 hartree) yielded this topography. Potential surfaces such as these are similar to that proposed for the biomolecule rhodopsin, the photopigment used to detect light in retinal rod cells.33 The wavepacket transfer still occurs through the unavoided crossing, but the momentum of the wavepacket and shape of the surface allow a large amount of the wavepacket transferred to the rightmost minimum to remain there instead of returning to the ground-state global minimum on the lefthand side. The shape of this surface also shows other interesting effects, most notably the amplitude loss of the excited-state wavepacket. The oscillation distance visibly decreases with each consecutive pass through the intersection. Because the unavoided crossing is at the minimum of the upper adiabatic potential, the wavepacket exists near the unavoided crossing for more extended periods of time and thus transfers more quickly. This is shown through the dephasing of the wavepacket, which takes just 55 fs to visibly begin to scramble the amplitude. Comparing this time with the time required for the previous unavoided crossing to dephase (110 fs) supports the conclusion that the topology of the potential energy surface affects the rates and yields of wavepacket transfer. Figure 9 shows two potential curves that are displaced enough to create an unavoided crossing at a much higher energy than the wavepacket in the Franck−Condon region. As expected, the wavepacket oscillates freely in the excited state with a period of 55 fs. This wavepacket never transfers to the

■

CONCLUSIONS Beginning in general chemistry, students develop a conceptual understanding of reaction coordinate diagrams and their utility in describing reactions such as catalysis. Students of quantum mechanics are able to quantify and describe the ground-state potential energy surface of diatomic molecules such as H2 using the Born−Oppenheimer approximation. This Journal has addressed important topics aimed to increase both the conceptual and quantitative understanding of quantum mechanics and quantum chemistry. Specifically, McCoy36 addressed an unexpected feature of modern quantum-chemical methods, and Kaliakin et al.37 showed how computations and 3D printing can be combined to enhance the visualization of transition states and reaction pathways. Castet and co-workers H



Article

ORCID

described a computational experiment for the photochromism of naphthopyrans, which is one of the few examples of a chemistry laboratory exercise involving excited-state potential energy surfaces.38 We have developed a minimum working example to help students visualize and study nonadiabatic quantum dynamics, thereby providing another example of chemical reactivity involving excited states. The model is general, does not require quantum-chemical computations, and is within the grasp of advanced undergraduate or graduate students of quantum mechanics. Even in its one-dimensional form, the model shows significant wavepacket transfer through the avoided and unavoided crossings. This occurs through nonadiabatic coupling and not through population relaxation, which could be introduced to allow the wavepacket to remain on the ground state after it has transferred through the nonadiabatic region. The method not only increases student understanding of nonadiabatic quantum dynamics but also provides a thorough reinforcement of many basic quantum-mechanical concepts. The important ideas of the harmonic oscillator, two-level system, time dependence, basis transformations, the Born− Oppenheimer approximation, and wavepackets are all strengthened by means of this exercise. Instructors could incorporate this model into the quantum-mechanics curriculum in a multitude of ways. An instructor could use the results as a visual aid during a discussion of adiabatic and diabatic bases, DVR, or quantum dynamics. Going further, the model could be usedin whole or in partas a problem set or semester-long project. The model may find use by undergraduate research groups interested in studying these topics, especially nonadiabatic dynamics. There is also room for expansion. Addition of more electronic states and vibrational modes would improve the model, although the present model is already sufficient to enhance student understanding of the material. These additions would provide further insight into what happens when the nonadiabatic coupling value is too large as well as whether an avoided crossing still merits wavepacket transfer. Altering the model to include dephasing and conventional population relaxation would allow for other interesting phenomena to be observed. Solutions to the Schrödinger equation also allow for the computation of a large range of observables. These could include but are not limited to average position, average momentum, and spectroscopic signals such as the transient absorption spectrum. Finally, a change in the value of σ should adjust the wavepacket momentum and induce novel dynamics, as we have shown that the nonadiabatic coupling behaves as a potential barrier.

■

Daniel B. Turner: 0000-0002-3148-5317 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The National Science Foundation supported this work through CAREER Grant CHE-1552235.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available on the ACS Publications website at DOI: 10.1021/acs.jchemed.6b00861. Detailed protocol, technical details, intermediate results, and derivations of the diabatic Hamiltonian and the electronic matrix elements for the nonadiabatic coupling (PDF) Sample MATLAB script (.m inside ZIP)

■

REFERENCES

(1) Teller, E. The crossing of potential energy surfaces. J. Phys. Chem. 1937, 41, 109−116. (2) Herzberg, G.; Longuet-Higgins, H. C. Intersection of Potential Energy Surfaces in Polyatomic Molecules. Discuss. Faraday Soc. 1963, 35, 77−82. (3) Stock, G. Classical description of nonadiabatic photoisomerization processes and their real-time detection via femtosecond spectroscopy. J. Chem. Phys. 1995, 103, 10015−10029. (4) Bernardi, F.; Olivucci, M.; Robb, M. A. Potential energy surface crossings in organic photochemistry. Chem. Soc. Rev. 1996, 25, 321− 328. (5) Domcke, W.; Stock, G. Theory of ultrafast nonadiabatic excitedstate processes and their spectroscopic detection in real time. Adv. Chem. Phys. 1997, 100, 1−169. (6) Yarkony, D. R. Conical intersections: The new conventional wisdom. J. Phys. Chem. A 2001, 105, 6277−6293. (7) Worth, G. A.; Cederbaum, L. S. Beyond Born-Oppenheimer: Molecular dynamics through a conical intersection. Annu. Rev. Phys. Chem. 2004, 55, 127−158. (8) van Voorhis, T.; Kowalczyk, T.; Kaduk, B.; Wang, L.-P.; Cheng, C.-L.; Wu, Q. The diabatic picture of electron transfer, reaction barriers, and molecular dynamics. Annu. Rev. Phys. Chem. 2010, 61, 149−170. (9) Tully, J. Perspective: Nonadiabatic dynamics theory. J. Chem. Phys. 2012, 137, 22A301. (10) Domcke, W.; Yarkony, D. Role of conical intersections in molecular spectroscopy and photoinduced chemical dynamics. Annu. Rev. Phys. Chem. 2012, 63, 325−352. (11) Schapiro, I.; Melaccio, F.; Laricheva, E. N.; Olivucci, M. Using the computer to understand the chemistry of conical intersections. Photochem. Photobiol. Sci. 2011, 10, 867−886. (12) Bernardi, F.; De, S.; Olivucci, M.; Robb, M. A. Mechanism of Ground-State-Forbidden Photochemical Pericyclic Reactions: Evidence for Real Conical Intersections. J. Am. Chem. Soc. 1990, 112, 1737−1744. (13) Coe, J. D.; Martínez, T. J. Competitive decay at two- and threestate conical intersections in excited-state intramolecular proton transfer. J. Am. Chem. Soc. 2005, 127, 4560−4561. (14) Hammes-Schiffer, S.; Stuchebrukhov, A. A. Theory of coupled electron and proton transfer reactions. Chem. Rev. 2010, 110, 6939− 6960. (15) Brazard, J.; Bizimana, L. A.; Gellen, T.; Carbery, W. P.; Turner, D. B. Experimental detection of branching at a conical intersection in a highly fluorescent molecule. J. Phys. Chem. Lett. 2016, 7, 14−19. (16) Levine, B. G.; Martínez, T. J. Isomerization through conical intersections. Annu. Rev. Phys. Chem. 2007, 58, 613−634. (17) Frutos, L. M.; Andruniów, T.; Santoro, F.; Ferré, N.; Olivucci, M. Tracking the excited-state time evolution of the visual pigment with multiconfigurational quantum chemistry. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 7764−7769. (18) Polli, D.; Altoe, P.; Weingart, O.; Spillane, K. M.; Manzoni, C.; Brida, D.; Tomasello, G.; Orlandi, G.; Kukura, P.; Mathies, R. A.; Garavelli, M.; Cerullo, G. Conical intersection dynamics of the primary photoisomerization event in vision. Nature 2010, 467, 440−443. (19) Middleton, C. T.; de La Harpe, K.; Su, C.; Law, Y. K.; CrespoHernández, C. E.; Kohler, B. DNA excited-state dynamics: From single bases to the double helix. Annu. Rev. Phys. Chem. 2009, 60, 217−239.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. I



Article

(20) Kim, P. W.; Rockwell, N. C.; Freer, L. H.; Chang, C.-W.; Martin, S. S.; Lagarias, J. C.; Larsen, D. S. Unraveling the primary isomerization dynamics in cyanobacterial phytochrome Cph1 with multipulse manipulations. J. Phys. Chem. Lett. 2013, 4, 2605−2609. (21) Tiwari, V.; Peters, W. K.; Jonas, D. M. Electronic resonance with anticorrelated pigment vibrations drives photosynthetic energy transfer outside the adiabatic framework. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 1203−1208. (22) Sakurai, J. J. In Modern Quantum Mechanics, revised ed.; Tuan, S. F., Ed.; Addison-Wesley: Reading, MA, 1994; pp 464−480. (23) Griffiths, D. J. Introduction to Quantum Mechanics; Prentice Hall: Upper Saddle River, NJ, 1995; pp 323−351. (24) Tannor, D. J. Introduction to Quantum Mechanics: A TimeDependent Perspective; University Science Books: Sausalito, CA, 2007; pp 23 and 335−348. (25) Carey, F. A.; Sundberg, R. J. Advanced Organic Chemistry, Part A: Structure and Mechanisms, 5th ed.; Springer: New York, 2008; pp 1079−1080. (26) Turro, N. J. Modern Molecular Photochemistry; University Science Books: Sausalito, CA, 1991; pp 52−75. (27) Robb, M. A.; Olivucci, M. Photochemical processes: Potential energy surface topology and rationalization using VB arguments. J. Photochem. Photobiol., A 2001, 144, 237−243. (28) Malhado, J. P.; Hynes, J. T. Photoisomerization for a model protonated Schiff base in solution: Sloped/peaked conical intersection perspective. J. Chem. Phys. 2012, 137, 22A543. (29) Virshup, A. M.; Chen, J.; Martínez, T. J. Nonlinear dimensionality reduction for nonadiabatic dynamics: The influence of conical intersection topography on population transfer rates. J. Chem. Phys. 2012, 137, 22A519. (30) Light, J. C.; Hamilton, I. P.; Lill, J. V. Generalized discrete variable approximation in quantum mechanics. J. Chem. Phys. 1985, 82, 1400. (31) Littlejohn, R. G.; Cargo, M.; Carrington, T., Jr.; Mitchell, K. A.; Poirier, B. A general framework for discrete variable representation basis sets. J. Chem. Phys. 2002, 116, 8691. (32) Because these states are now, finally, eigenfunctions of the total Hamiltonian, propagation under the time-dependent Schrödinger equation involves straightforward exponentiation. (33) Wang, Q.; Schoenlein, R. W.; Peteanu, L. A.; Mathies, R. A.; Shank, C. V. Vibrationally coherent photochemistry in the femtosecond primary event of vision. Science 1994, 266, 422−424. (34) Balzer, B.; Hahn, S.; Stock, G. Mechanism of a photochemical funnel: a dissipative wave-packet dynamics study. Chem. Phys. Lett. 2003, 379, 351−358. (35) Gelman, D.; Katz, G.; Kosloff, R.; Ratner, M. A. Dissipative dynamics of a system passing through a conical intersection: Ultrafast pump-probe observables. J. Chem. Phys. 2005, 123, 134112. (36) McCoy, A. B. Investigating the Adiabatic Approximation in Quantum Mechanics through the Analysis of Two Coupled Harmonic Oscillators. J. Chem. Educ. 2001, 78, 401−404. (37) Kaliakin, D. S.; Zaari, R. R.; Varganov, S. A. 3D Printed Potential and Free Energy Surfaces for Teaching Fundamental Concepts in Physical Chemistry. J. Chem. Educ. 2015, 92, 2106−2112. (38) Castet, F.; Méreau, R.; Liotard, D. Photochemistry of Naphthopyrans and Derivatives: A Computational Experiment for Upper-Level Undergraduates or Graduate Students. J. Chem. Educ. 2014, 91, 924−928.

J


A Tractable Numerical Model for Exploring ... - ACS Publications

Recommend Documents