A Straightforward Method of Analysis for Direct Quantum Dynamics

J. Phys. Chem. A 2010, 114, 8713–8729

8713

A Straightforward Method of Analysis for Direct Quantum Dynamics: Application to the Photochemistry of a Model Cyanine† Charlotte S. M. Allan,‡ Benjamin Lasorne,*,‡,§ Graham A. Worth,| and Michael A. Robb*,‡ Department of Chemistry, Imperial College London, London SW7 2AZ, United Kingdom, School of Chemistry, UniVersity of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom, and Institut Charles Gerhardt Montpellier, UMR 5253, CNRS-UM2-ENSCM-UM1, CTMM, UniVersite´ Montpellier 2, CC 1501, Place Euge`ne Bataillon, 34095 Montpellier, France ReceiVed: February 22, 2010; ReVised Manuscript ReceiVed: May 4, 2010

We present a new way of analyzing direct quantum dynamics simulations based on a Mulliken-type population analysis. This provides a straightforward interpretation of the wavepacket in much the same way as semiclassical trajectories are usually analyzed. The result can be seen as a coupled set of quantum trajectories. We apply this to the study of the photochemistry of a 12-atom model cyanine to explore possibilities for intelligent optimal control. The work presented here builds on previous semiclassical dynamics simulations [Hunt, P. A.; Robb, M. A. J. Am. Chem. Soc. 2005 127, 5720]. Those calculations suggested that, by controlling the distribution of momentum components in the initial wavepacket, it should be possible to drive the system to a specific region of the conical intersection seam and ultimately control the product distribution. This was confirmed experimentally by optimal control methods [Dietzek, B.; Bru¨ggemann, B.; Pascher, T.; Yartsev, A. J. Am. Chem. Soc. 2007 129, 13014]. This paper aims to demonstrate this in a quantum dynamics context and give further insight into the conditions required for control. Our results show that directly addressing the trans-cis torsional modes is not efficient. Instead, one needs to decrease the momentum in the skeletal deformation coordinates to prompt radiationless decay near the minimum conical intersection at large twist angles. 1. Introduction Radiationless decay is an intrinsically quantum effect that involves the breakdown of the Born-Oppenheimer approximation. Quantum dynamics simulations have thus proved to be invaluable for understanding the complexity of photochemical processes of organic molecules. However, such accurate calculations become increasingly expensive as the size of the system grows, and semiclassical trajectory methods have been implemented in quantum chemistry programs as cheaper, but approximate, alternatives. In addition, in quantum dynamics, the time evolution of the photoexcited molecule is represented by a global wavepacket, a complicated function in the space of the nuclear coordinates. In the method used here (direct dynamics variational multiconfiguration Gaussian wavepackets (DD-vMCG)1-3) the global wavepacket is a collection of coupled Gaussian functions that follow “quantum trajectories”. We present a new way of analyzing the simulations based on a Mulliken-type population analysis. This provides a straightforward interpretation of the wavepacket. The method allows direct analysis of the quantum trajectories in much the same way as semiclassical trajectories are usually analyzed. Focus on the trans-cis photoisomerization of cyanines stems not only from interest in the small cyanines used as fluorescent dyes for biomedical imaging but from the desire to understand the process in larger protonated Schiff base (PSB) systems such †

Part of the “Klaus Ruedenberg Festschrift”. * Corresponding authors. E-mail: B.L., [email protected]; M.A.R., [email protected]. ‡ Imperial College London. § Universite´ Montpellier 2. | University of Birmingham.

CHART 1: Dihedral N1-C2-C4-N5 Measures the Extent of Trans-Cis Isomerization and C2-C3-C4-H8 and C3-C4-N5-H9 Measure Out-of-Plane Deformation

as the rhodopsin chromophore, responsible for vision in humans4 and photosynthesis in halobacteria.5 Thus, a great deal of experimental and theoretical work has focused on cyanines and protonated Schiff bases. These have employed both small model systems and larger systems that allow for the influence of the protein environment in the reaction. Although it has been shown that the length of the carbon chain can have a significant effect on the shape of the potential energy surface (PES)6 and therefore the progress of the reaction, small cyanine models continue to be used to explore the possibilities for intelligent optimal control. In this, the aim of theoretical investigations is not the calculation of a laser pulse that will achieve the target, but rather the identification of the mechanistic key features that can be utilized to control the photoreactivity in laser-driven experiments. In this work we have chosen the 12-atom model system trimethine cyanine (1), Chart 1. This molecule is small enough that quantum dynamics simulations are not prohibitively expensive and is the same model system used for previous semiclassical dynamics simulations so that meaningful comparisons can be drawn. The total dihedral angle N1-C2-C4-N5 indicates the degree of trans-cis isomerization, a twist angle

10.1021/jp101574b  2010 American Chemical Society Published on Web 05/25/2010

8714

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Figure 1. Cartoon of the two radiationless decay mechanisms for the cyanine model. The black arrow starts at the Franck-Condon point, stays close to the S1 reaction path (an asynchronous, concerted mixture of two sequential branches (the two perpendicular black lines): first, an in-plane relaxation, second, a trans-cis torsion), crosses at a twisted conical intersection, and continues on S0 to form the cis product. The blue arrow starts at the Franck-Condon point, follows the initial relaxation path (with more available energy) up to a quasi-planar conical intersection, and continues on S0 to regenerate the trans reactant. The white line represents the seam of intersection between the two potential energy surfaces and is approximately parallel to the second branch of the minimum energy path on S1(torsion).

of 0° in the trans starting structure and 180° in the cis isomerization product. The quantum dynamics study presented here builds on previous semiclassical dynamics simulations.7 Those calculations suggested that, by controlling the distribution of momentum components in the initial wavepacket, it should be possible to drive the system to a specific region of the conical intersection seam and ultimately control the product distribution. This has been shown to be possible experimentally by optimal control methods.8 The aim of this paper is to support the validity of this idea using quantum dynamics simulations. Previous work and this work confirm how we can control reactivity. The key feature is the topological relationship between the hyperline of conical intersection between the S1 and S0 PES (the “seam”, of dimension 3N - 8 in the space of the 3N - 6 nuclear coordinates) and the S1 minimum energy path (MEP) originated from the Franck-Condon (FC) point. Dynamics simulations and systematic exploration of the seam show that the most accessible part of the seam in terms of energy is approximately parallel to the MEP along the torsional coordinate on S1. The terminus of the MEP, the S1 minimum, lies near the minimum energy point on the seam. Therefore, decay can occur at a range of torsional angles if one can vary the excitation of vibrational motions orthogonal to the trans-cis reaction path (see Figure 1). Accordingly, the product distribution is controlled by the partition of the momentum between the torsional coordinate and the coordinates of the branching space:9-12 namely, in-plane skeletal deformation. It is worth noting that the reaction coordinate on S1 is a complicated mixture of concerted but asynchronous deformations.6 Nevertheless, the key features can be understood with two general coordinates as shown in Figure 1. In the following, we will refer to the beginning of the MEP as “the initial relaxation” (skeletal deformation coordinate, also involved in the branching space) and to the end as “the reaction path” (trans-cis torsional coordinate). We present one of the first calculations performed in the full space of the 3N - 6 nuclear coordinates using the DD-vMCG

Allan et al.

Figure 2. Three successive snapshots of a 3-GBF wavepacket example. The black arrowed curves represent the quantum trajectories followed by the centers of the three GBF in the space of the nuclear coordinates.

method.1-3 Direct dynamics approaches calculate the PES “onthe-fly” using quantum chemistry calculations, as it is required along a trajectory representing the motion of the nuclei in time. The DD-vMCG method is a direct quantum dynamics approach. The PES is still calculated on-the-fly but along trajectories followed by the centers of coupled Gaussian basis functions (GBF). This method is distinguished from the closely related ab initio multiple spawning (AIMS) method13 by the fact that the positions of the GBF are determined variationally from the timedependent Schro¨dinger equation applied to the motion of the nuclei. Our approach ensures that the basis functions cover the most important parts of the space. In addition, the DD-vMCG equations of motion for the GBF are quantum rather than classical as they are in AIMS, which ensures a faster convergence and thus fewer basis functions needed. The result can be seen as a set of quantum trajectories (see Figure 2) “steered” by the occupation of the various basis functions. Complicated quantum mechanical effects are easily understood by consideration of the populations of the individual basis functions that make up the global wavepacket, using our new Mulliken-type population analysis, which we will discuss subsequently. We present two representative calculations to illustrate radiationless decay that are distinguished by the region of the seam where internal conversion from S1 to S0 takes place (see Figure 1). The first occurs late in the reaction coordinate, at twisted geometries, near the minimum of the reaction path; and the second early in the reaction coordinate, namely at “planar” geometries (we will refer to this as a planar geometry although, as will be discussed, the geometries of structures on the seam are only planar with respect to the N-C-C-C-N backbone). We show how it is possible to control the reactivity of the model trimethine cyanine, forcing decay at planar geometries (expected to lead to regeneration of the molecule in the trans geometry) or at twisted geometries (expected to lead to the cis product). The important control feature for directing reactivity appears to be the motion of the wavepacket biased along the skeletal deformation coordinate. For the largest magnitude of the momentum the wavepacket moves quickly to the seam in the region of quasi-planar geometries. As the excess energy is decreased systematically, the angle at which the nonadiabatic event occurs changes continuously toward the fully twisted minimum of the conical intersection seam. So, perhaps counterintuitively, one needs to decrease the momentum in the skeletal deformation coordinates to prompt decay near the minimum conical intersection at large twist angles. This has

Method of Analysis for Direct Quantum Dynamics

Figure 3. Cartoon representation of the minimum energy path showing initial skeletal stretching, followed by conrotatory torsion and finally disrotatory torsion.

also been observed for fulvene photochemistry, and it has been shown theoretically that this damping is possible for that system.14 2. Background Previous theoretical studies of model cyanines have focused largely on mapping the MEP from the FC point to the minimum on S1.6 For the model system discussed here, the MEP shows a stepwise breaking of symmetry leading from the C2V FC point toward C2-like geometries and finally the C1 minimum on the excited state (see Figure 3). There are no stationary points on the MEP. Symmetry-restricted minima (C2V and C2, strictly saddle points within the complete coordinate space) can be located on the S1 surface although the MEP does not pass through them directly: it traces a path toward them, but the symmetry is broken before they are reached. Initial relaxation consists of in-plane symmetric stretching (toward the C2V saddle point). This is followed by conrotatory torsional motion (with respect to the two central CC bonds) leading toward C2 geometries. Finally, before the C2-symmetryrestricted minimum is reached, disrotatory torsion takes over, the C2 symmetry is broken, and the path reaches the S1 C1 minimum. A C1 minimum conical intersection lies nearby. The end of the reaction path is dominated by skeletal deformation (single-double bond inversion) together with a barrierless rotation around the central bond (trans-cis isomerization).6 The minimum conical intersection is a structure that is twisted by approximately 100° about one of the central CC bonds. Subsequently, the energetic accessibility of an extended conical intersection seam in the model cyanine trimethine was demonstrated using semiclassical dynamics trajectories.7 The accessible seam was shown to cover the complete range of dihedral angles N1-C2-C4-N5. This in turn suggested that the geometry of the conical intersection reached should determine the S0 product observed: molecules that reach the seam at twisted geometries (N1-C2-C4-N5 > 90°) will decay to the cis product and those that reach the seam at higher energies and at less twisted geometries (N1-C2-C4-N5 < 90°) will regenerate the trans reactant (see Figure 1). The dynamics showed that nonadiabatic events should occur at all twist angles, suggesting that the cis or trans product might be targeted by controlling the distribution of momentum components in the initial wavepacket. The branching-space vectors9-12 are in-plane skeletal deformations at each point on the seam. Thus, it was suggested that excitation of these modes should direct the wavepacket toward the seam (to flat geometries where N1-C2-C4-N5 < 90°) and excitation of torsional modes should direct it down the MEP to the twisted minimum conical intersection.

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8715 Investigation of the reaction starting at the excited cis product15 provided additional information about the S1 surface of a similar set of molecules. The protonated Schiff base tZt 1-iminium-2,4-pentadiene (H2CdCHsCHdCHsCHdNH2+) was shown to exhibit a local minimum on S1 in which the terminal NH2 group is 90° twisted. However, no conical intersection was found near this S1 minimum (although the other features of the S1 surface and conical intersection seam are similar). Experimental work on cyanines and protonated Schiff bases has employed, among others, the molecule 1,1′-diethyl-2,2′cyanine iodide, which has been shown to behave in accordance with the theoretical results described above. Control of product formation has been demonstrated using optimal control methods.8,16 It was found that specific stimulation of the wavepacket prior to excitation to S1 caused a reduction in the excited-state population. This was explained by assuming that S1-S0 decay through a conical intersection occurs early in the reaction path (i.e., at less twisted geometries). The mode excited by this laser pulse was assigned (through comparison with previous theoretical studies and density functional theory calculations) as a symmetric stretch. The dump pulse required to efficiently dump the remaining excited molecules later in the reaction was composed of low-frequency torsional modes. These modes were identified as modes that should prevent early decay at “flat” geometries and induce twisting and therefore decay near the minimum conical intersection. It has become increasingly clear that a one-dimensional model is not sufficient for describing the isomerization process in these systems and that control is achieved by specific excitation of modes orthogonal to the torsional coordinate rather than directly addressing the torsional modes themselves.8,16-18 In support of this, solvent studies have shown that torsion (but not orthogonal motion) is damped by viscous solvents. In viscous solvents the isomerization yield is reduced, implying that the wavepacket remains in the vicinity of the higher energy region of the conical intersection seam for longer and there is increased decay through flat geometries, with regeneration of the trans reactant.19 This paper aims to demonstrate this in a quantum dynamics context and give further insight into the conditions required for control. 3. Direct Dynamics Methods and Their Analysis as Coupled Trajectories Previous work20-22 has described the results of quantum dynamics calculations using the DD-vMCG method1-3 to simulate radiationless decay. This is a quantum dynamics algorithm based on the propagation of coupled Gaussian functions. Each function follows a quantum trajectory along which the PES are calculated “on-the-fly” by a quantum chemistry program. We present here a new method of analyzing the coupled trajectories based on the Mulliken population analysis.23 First, we recall a few basic principles involved in quantum dynamics beyond the Born-Oppenheimer approximation and we give a brief discussion of the DD-vMCG method. Excited-State Quantum Dynamics. Radiationless decay happens at molecular geometries where the Born-Oppenheimer approximation is not valid. In quantum dynamics, this implies that the molecular wavepacket must be expanded in a basis set of coupled electronic states,

|Ψ(Q,t)〉 )

∑ ψ(s)(Q,t)|s;Q〉 s

(1)

8716


Allan et al.

where s is the label of the electronic state, and Q is the collection of nuclear coordinates that define the molecular geometry (e.g., a set of 3N - 6 rectilinear displacements along normal modes, or normal coordinates, for an N-atom nonlinear molecule). The notation |s;Q〉 means that the electronic state labeled s depends parametrically on Q. Each individual component in this expansion is a nuclear wavepacket. The time evolution of the collection of coupled nuclear wavepackets is obtained by solving the time-dependent Schro¨dinger equation for the molecular Hamiltonian represented in the electronic basis set,

∂ ip ψ(s)(Q,t) ) ∂t

∑

(ss') (s') ˆQ H ψ (Q,t)

(2)

s'

(ss′) ˆQ acts in the Q-space and connects components where H corresponding to the electronic states s and s′. Most photochemical applications involve two coupled singlet states, S1 and S0, and start with an initial condition where the photon absorption creates a Franck-Condon wavepacket on S1. This corresponds to a sudden projection of the vibrational ground state of S0, χ0(0)(Q), to S1, such that ψ(0)(Q,0) ) 0 and ψ(1)(Q,0) ) χ0(0)(Q). This hypothesis is largely used in quantum dynamics applications. In particular, this allows for the direct calculation of the electronic absorption spectrum at all wavelengths (see, e.g., ref 24). This also is the reason why semiclassical simulations use a random sampling of the zero-point-energy vibrations in S0 to generate the distribution of initial positions and momenta for the swarm of trajectories that will start on S1 (see, e.g., ref 7). This creates a Wigner distribution in the classical phase space that is the best statistical representation of the quantum nuclear wave function χ0(0)(Q). Under the harmonic approximation, χ0(0)(Q) is a multidimensional Gaussian function of the normal coordinates centered at the S0 equilibrium geometry. This Gaussian wave function is equivalent to the corresponding Wigner distribution used in semiclassical harmonic sampling. The molecular wavepacket is normalized, but the electronic components, ψ(s)(Q,t), are not. The corresponding integrated densities of probability are the electronic populations, which read

P(s)(t) )

∫Q |ψ(s)(Q,t)|2 dQ

(3)

where integration is performed over the Q-space. Consequently, the normalization condition reads

∫Q 〈Ψ(Q,t)|Ψ(Q,t)〉 dQ ) ∑ P(s)(t) ) 1

(4)

s

where 〈|〉 denotes integration over the space of the electronic states at each Q. The Franck-Condon initial conditions detailed above imply that the electronic populations at time t ) 0 are P(0)(t) ) 0 and P(1)(t) ) 1 on S0 and S1, respectively. The time evolution of the individual electronic populations is due to (01) (10)† ˆQ ˆQ nonzero off-diagonal elements, H ) H . These nonBorn-Oppenheimer terms called nonadiabatic couplings account for the nonadiabatic event: the decay of population from S1 to S0 observed through internal conversion. They take various expressions depending on the definition of the electronic states (see, e.g., ref 2 and references therein for more details; see also Appendices A and B).

Direct Quantum Dynamics. In the DD-vMCG method, each nuclear wavepacket (i.e., each component on a given electronic state) is expanded in a basis set of time-dependent parametrized Gaussian functions, the GBF, gj(Q,t). In this work, we used the single-set formalism, where the same GBF are used for all electronic states. The corresponding ansatz reads

ψ(s)(Q,t) )

∑ Aj(s)(t) gj(Q,t)

(5)

j

where j is the label of the GBF. Note that the multiset formalism was used in previous work (in this, GBF are indexed by s as well). The single-set ansatz for the molecular wavepacket can be recast as

|Ψ(Q,t)〉 )

∑ gj(Q,t) ∑ Aj(s)(t)|s;Q〉 j

(6)

s

which means that each GBF has a different coefficient on each electronic states. It can be noted that the one-GBF limit is related to the Ehrenfest approach for semiclassical dynamics, where the trajectories Q(t) are driven by a mixture of several electronic states, ΣsA(s)(t)|s;Q(t)〉 (see, e.g., refs 25 and 26). The multiset DD-vMCG ansatz is the same as in the ab initio multiple spawning (AIMS) method.13 Both are called Gaussian wavepacket (GWP) methods. However, they correspond to different equations of motion. GWP basis functions, i.e., GBF, are localized in space so that a local harmonic approximation (LHA) is a reasonable representation of the PES around their center.27 Both DD-vMCG and AIMS have thus been implemented as direct quantum dynamics methods, where the potential energy and its first and second derivatives are calculated “on-the-fly” by a quantum chemistry program along the trajectories followed by the centers of the GBF in much the same way as in classical or semiclassical dynamics algorithms. Note that frozen-width GBF were used in this work, which requires more functions for convergence but makes the calculations more stable. Indeed, “thawed” GBF can spread and make the LHA approximation more likely to break down. The DD-vMCG approach provides variational equations of motion not only for the expansion coefficients, Aj, but also for the GBF, gj [eq 5]. In brief, the center of each GBF follows a “quantum trajectory” determined from the time-dependent Schro¨dinger equation, whereas it follows a classical trajectory in the AIMS approach. In both cases, the trajectories are coupled through their coefficients, as opposed to semiclassical methods such as trajectory surface hopping (TSH) where trajectories are independent (see, e.g., ref 28). This coupling between the trajectories makes GWP methods capable of directly accounting for quantum effects such as tunneling and radiationless decay. Further, quantum rather than classical trajectories allow for a faster convergence of the solution because they are variational and accumulate where the exact quantum density is large so that fewer are needed. For a recent review of direct semiclassical and quantum dynamics methods, see ref 2. The DD-vMCG method is related to the multiconfiguration timedependent Hartree (MCTDH) quantum dynamics method29 and has been implemented in a development version of the Heidelberg MCTDH package.30 In this work, it was interfaced with a development version of the Gaussian quantum chemistry program.31 Mulliken-Type Population Analysis of the Quantum Dynamics. Semiclassical dynamics simulations approximate the quantum wavepacket as a swarm of classical trajectories. Each


J. Phys. Chem. A, Vol. 114, No. 33, 2010 8717

trajectory can be seen as a moving picture of the molecular geometry with a statistical weight. Quantum dynamics approaches such as MCTDH generate delocalized wavepackets, the interpretation of which in classical terms is not straightforward. In contrast, the Gaussian functions used in the DD-vMCG approach are localized around their centers and can be seen as following quantum trajectories. Each one will contribute to some extent to the global wavepacket. The problem here is that the basis set is not orthogonal, which makes the definition of individual contributions ambiguous. We propose here a solution in the form of a new analysis derived from the Mulliken analysis of orbital populations.23 For a basis set made of Gaussian functions, the pseudo Gaussian populations (PGP),

electronic state. The global wavepacket is expanded in a basis set of four GBF, 4

ψ(x,t) )

∑ Aj(t) gj(x,t)

(13)

j)1

where x is the coordinate. In Figure 4, the black curve with four maxima represents a snapshot at a given time t of the global density of probability, 4

|ψ(x,t)| ) 2

Aj'(t) g*(x,t) gj'(x,t) ∑ A*(t) j j

(14)

j, j')1

PGPj(s)(t) ) |Aj(s)(t)| 2

(7)

cannot be used directly as a measure of the GBF weights. Indeed, the normalization condition,

∑ Sjj'(t) ∑ Aj(s)*(t) Aj'(s)(t) ) 1 j,j'

(8)

s

In Figure 4, the four gray Gaussian curves are the weighted densities of probability of each Gaussian functions, |Aj(t) g(x,t)|2. The global density of probability, |ψ(x,t)|2, is not the sum of these, since nonzero overlaps create interference terms (as when overlapping atomic orbitals create bonding or antibonding molecular orbitals). As a matter of simplicity, we assume here that each Gaussian function is normalized and real with a width equal to 1/2,

involves nonzero overlaps,

Sjj'(t) )

gj'(Q,t) dQ ∫Q g*(Q,t) j

gj(x,t) )

(9)

and thus contributions due to cross-terms, Aj(s)*(t) Aj′(s)(t) with j * j′. This is a problem similar to the estimation of how much each atomic orbital contributes to the global electronic wave function. An imperfect but operational solution is the Mulliken analysis of orbital populations.23 We extend this formalism to the DD-vMCG method by introducing gross Gaussian populations (GGP),

GGPj(s)(t) ) R

∑ Sj'j(t) Aj'(s)*(t) Aj(s)(t)

(10)

j'

which satisfy

∑ GGPj(s)(t)

(11)

∑ ∑ GGPj(s)(t) ) 1

(12)

P(s)(t) )

e-1/2[x-xj(t)] 4

√π

2

(15)

and all coefficients are real. The width is the standard deviation of the corresponding density of probability, |gj(x,t)|2. In this example, we chose the centers {xj(t)} ) {3, 6, 9, 12} and the coefficients {Aj(t)} ∝ {1, 3, 5, 2} (corrected by a constant normalization factor). The overlap between two neighboring Gaussian basis functions is approximately 11% and is virtually zero between functions further apart. The PGP and the GGP are {PGPj(t)} ) {2%, 20%, 56%, 9%} and {GGPj(t)} ) {3%, 24%, 62%, 11%}, respectively. The former do not add up to one but 87% (as expected due to nonzero overlaps), while the latter add up to one by construction. Note that the example illustrated here is a case where all interferences are constructive (coefficients and overlaps are all real and positive). The sum of the PGP is less than one, and interfering contributions account for the rest (13% of the global density of probability). In practice, the coefficients and the overlaps are complex numbers,

j

and

s

j

In eq 10, R stands for the real part. The imaginary parts can be ignored as they sum up to zero. In our calculations, the individual imaginary parts were always small. This method calculates the contribution made by each GBF to the total wavepacket. Each GBF can be examined separately, providing a detailed breakdown of the total wavepacket. The separate trajectories describe the space explored in the simulation and the populations of the individual basis functions indicate when and where electronic population transfers occur. These concepts will be illustrated with an example (see Figure 4). Let us consider a one-dimensional case with a single

Figure 4. Snapshot of the density of probability of a global wavepacket (black curve) and of its weighted Gaussian components (gray curves). The vertical gray bars indicate the gross Gaussian populations (the corresponding numbers are given in %). The horizontal axis corresponds to the spatial coordinate x.

8718


and the sum of the PGP can be much greater than one when there are destructive interferences. In Figure 4, the four vertical gray bars represent the four GGP (their respective heights are indicated) positioned at the centers of the Gaussian basis functions. They can be considered as a reduction of the quantum distributionsthe global density of probabilitysto a statistical sample of four representative positions with their respective weights. Following these four positions and weights along time can be seen as following four “quantum trajectories” that best represent the global wavepacket. To plot the four GGP on the same graph as the densities of probabilities and have quantities that can all be compared with respect to the vertical axis, the GGP are scaled by a factor of 1/3, and an error bar added of width 3 ((1.5 around each central position) corresponding to the spacing between each neighbor centers. Doing so, the surface area of each rectangle is thus equal to the GGP calculated at its center, which can be seen as an estimate of the integrated probability of finding x within the boundaries of the corresponding domain. The GGP are thus effective indicators of the actual weight of the Gaussian basis functions within the global density of probability. They account for the individual contributions (as the PGP do) but also for the interferences due to overlapping basis functions. As for Mulliken’s gross orbital populations, the formula used to decide how much of the interfering contribution is to be attributed to a particular orbital is somewhat arbitrary and can be criticized. Nevertheless, this represents an operational procedure that helps the interpretation of the whole quantum distribution in terms of a few weighted quantum trajectories in much the same way as delocalized molecular orbitals can be seen as superpositions of localized atomic orbitals. In treating nuclear motion quantum mechanically, one allows for the fact that the nonadiabatic event can be nonlocal and only partially complete. In contrast, one can see nonlocality in a semiclassical trajectory calculation only in the distribution of the results, which are independent of each other. Semiclassical nonlocality is related to the statistical spreading of the points where surface hopping occurs. The population transfer occurs at different times for different trajectories. Of course, this also means different geometries, but we will refer to such results as nonlocal in time. Quantum nonlocality has to do with the fact that quantum trajectories communicate with each other at all times. Results will be called nonlocal in space if the population transfer involves a correlated decrease of GGP(1) j (t) and increase of GGPj′(0)(t) with j * j′, i.e., between two GBF gj(Q,t) and gj′(Q,t) on different trajectories (see Figure 5). The results presented in this paper (albeit with a limited number of basis functions) suggest that the nonadiabatic event is rather local in space but not local in time. This implies that semiclassical surface hopping trajectories, if one performs enough of them, are a reliable way of treating the problem. This thus justifies the validity of the conclusions drawn in previous work.7 However, it is clear that quantum trajectories are an efficient way of getting mechanistic information in photochemical processes. Direct quantum dynamics simulations need far fewer functions to reach convergence, as illustrated in ref 2. In addition, semiclassical trajectories usually give good estimates of the yield and time scale of a single-event decay but may fail to describe the subsequent time evolution of the system (groundstate dynamics or further recrossings between the two states). 4. Computational Details Direct dynamics calculations employed the DD-vMCG method1-3 as implemented in a development version of the

Allan et al.

Figure 5. Comparison of local transfer (left panel) vs nonlocal transfer (right panel). The upper panel shows a 2-GBF case before any transfer and the lower panel the same case after the transfer from S1 to S0. The horizontal axis corresponds to a generic nuclear coordinate x.

Heidelberg MCTDH package.30 The LHA representation of the PES requires the calculation of the energy and first and second derivatives of the S1 and S0 adiabatic surfaces at the center of each GBF in the wavepacket expansion. These are evaluated on-the-fly using quantum chemistry calculations made by a development version of the Gaussian program31 interfaced with the quantum dynamics code. As the GBF propagates along the quantum trajectory, each new complete active space selfconsistent field (CASSCF) calculation uses the active space of the previous calculation as its guess. The 3-21G basis set was used with a CASSCF wave function. The active space was the set of five π-molecular orbitals and six electrons, obtained at the ground-state C2V minimum geometry. The same stationary points and optimized structures of the seam were found as when the 6-31G(d) basis set was used (see Supporting Information). The surfaces topology was found not to differ significantly between the two levels of calculation. The maximum error on the relative energy is no more than 8 kcal mol-1. Note that we are not interested in any feature that depends quantitatively on the energetics. Our reason for using the 3-21G basis set was merely to save computer time in direct quantum dynamics simulations. The pair of coupled electronic states was represented in a diabatic picture to make the off-diagonal nonadiabatic coupling, (01) ˆQ [eq 2], enter the Hamiltonian matrix as a local potentialH like function of Q. A simple form of the regularized diabaticstates method of Ko¨ppel et al.32,33 (see also Appendix B for more details) was used to transform the adiabatic information from the quantum chemistry calculations into the LHA expansion of both diabatic PES and the coupling function. In brief, this approach uses the gradient difference and derivative coupling vectors calculated at a particular conical intersection point, called the reference point. These branchingspace vectors9-12 are the two directions that lift degeneracy around the apex of the double cone at first order. They are used here to linearly extrapolate the diabatic Hamiltonian matrix from the reference point to any current point, so that the two eigenvalues coincide with the adiabatic energies calculated with the quantum chemistry program at the current point. Various reference points were used, according to where the wavepacket is expected to cross the seam. In addition, the two branchingspace vectors were rotated within the plane they span to ensure that the off-diagonal coupling is zero at the FC geometry and

Method of Analysis for Direct Quantum Dynamics thus make the adiabatic and diabatic states coincide at the starting point of all calculations. The difference between the adiabatic and diabatic representations in our simulations may need some discussion here. In brief, a naive picture would assume diabatic states and adiabatic states coincide before the crossing and swap but still coincide after the crossing. If so, assigning diabatic populations [eq 3] to adiabatic populations is simply a matter of knowing the energetic order of the two diabatic states: the initial diabatic state is first S1 and becomes S0 once its energy is lower than that of the other diabatic state, and populations change accordingly. This is, however, only correct if the state crossing is not accompanied by a change in the coordinates that provide the nonadiabatic coupling. The nonadiabatic coupling responsible for the transfer of population drives the system to a region where the coupling between the diabatic states increases, so that the diabatic states are different from the adiabatic states. The 2 × 2 electronic Hamiltonian matrix in the diabatic representation used here [eq 2] has nonzero off-diagonal terms (couplings), which implies that the diagonal elements (diabatic PES) are different from the eigenvalues (adiabatic PES). This deviation along the direction of the nonadiabatic coupling is why semiclassical TSH trajectories gain extra kinetic energy along this direction when they hop down at finite energy gap to account for the loss of potential energy. Despite this, in the present study the ordered diabaticstate populations are taken as approximate adiabatic-state populations. The nonzero nonadiabatic coupling, however, leads to Rabi-like population oscillations (see Appendix A). In all simulations, the initial wavepacket was formed as a product of 30 one-dimensional Gaussian functions representing a good approximation of the vibrational ground state in S0, χ0(0)(Q) (see section 3). The frequencies were obtained from CASSCF calculations performed at the S0 equilibrium geometry. Normal coordinates Q were scaled to frequency-mass-weighted dimensionless displacements from the FC point, so that χ0(0)(Q) corresponds to a width equal to 1/2 in each of the 30 coordinates. All initial wavepackets were started at the FC point on the upper diabatic state, which coincides there with S1. To control the reactivity, we modulated the stimulation of selected modes in the initial conditions. To this end, the initial wavepacket was given a nonzero group velocity in the form of an extra average momentum, P0, by multiplying the real-valued multidimensional Gaussian function by the spatial phase factor exp (iP0 · Q/p), where P0 is a vector of components with chosen magnitudes along the normal coordinates. Classically, this corresponds to shifting the center of the momentum distribution of the sample of initial conditions from zero momentum to a finite P0. The resulting initial wavepacket at time t ) 0 [eq 6] is such that ψ(0)(Q,0) ) 0 and ψ(1)(Q,0) ) g1(Q,0) ) χ0(0)(Q) exp(iP0 · Q/p). All A coefficients are zero except A(1) 1 (0) ) 1, which corresponds to the first GBF, g1(Q,0). The remaining nonpopulated GBF are arbitrarily shifted around the first one at time t ) 0. Then, from the first step of the simulation on, they start gaining population and moving where they are most needed, according to the variational solution of the equations of motion. We have run simulations with between four and 16 basis functions. Only 4-GBF calculations will be shown here. Some 1-GBF and 16-GBF results are shown in Supporting Information. Simulations with many basis functions are difficult to converge with respect to CASSCF calculations due to GBF with small populations reaching very distorted geometries. The limiting case of a 1-GBF calculation is close to the Ehrenfest

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8719 model where the motion of the trajectory is determined by a quantum average of the PES of the two states, which induces very large oscillations in the populations. This unphysical behavior is known to happen because the overconstrained wavepacket cannot split into distinct components on S0 and S1.26 Another reason for using extra initial momentum is related to symmetry breaking and the method used to generate diabatic states from adiabatic states. The starting point corresponds to a C2V geometry. The three kinds of nontotally symmetric representations (A2, B1, B2) are present in the list of normal modes in this system. For each three, there are thus two equivalent sides (clockwise-anticlockwise, up-down, left-right), which makes eight equivalent ways of breaking C2V into C1. In other words, each C1 structure corresponds to eight equivalent structures with the same energy but distinguishable coordinates. Such structures can be enantiomers or, more generally, enantiotops that are distinct only because identical nuclei have different labels. In terms of normal coordinates, this means that the space they span can be split into eight equivalent domains all having in common the same C2V points (e.g., the FC point) but distinct C1 points that are mirror images of each other. In particular, all conical intersection points used as references for the adiabatic-to-diabatic transformation have C1 symmetry. To make this locally valid transformation fit in the region actually explored by the wavepacket, we thus added residual momentum along some selected nontotally symmetric modes to break symmetry in a controlled way. Thus, we could drive the system to the specific space domain among the eight equivalent ones that contains the C1 conical intersection chosen as a reference. In addition, this bias helps to reduce the number of basis functions. A wavepacket sent with no momentum along nontotally symmetric modes cannot break symmetry. It would thus spread symmetrically along such modes along its course. For example, in Figure 2, the vertical axis could be seen as a nontotally symmetric coordinate along which the wavepacket spreads symmetrically. The function on the negative side is the mirror image of the one on the positive side. Here, this implies that eight GBF (with equal coefficients and trajectories that are mirror images of each other) would be needed instead of a single representative one when C1 regions are explored in the absence of any momentum bias. 5. Results and Discussion In this section, we will first use knowledge of the S1 MEP and S1/S0 conical intersection seam to define potential targets for controlling the branching ratio between regeneration of the trans reactant and formation of the cis product. Then, we will present the construction of the initial momenta necessary to reach the two limiting targets (twisted vs quasiplanar conical intersections), and we will show the geometries visited along the time evolution of the corresponding wavepackets. Finally, the quantum trajectories will be analyzed in semiclassical terms, and the local character of the population transfer will be discussed. This will be followed by a systematic investigation whereby we show that a range of intermediate targets can be reached by reducing the magnitude of the momentum used for targeting the quasi-planar conical intersection. Definition of the Targets. Knowledge of the shape of the S1 PES and the conical intersection seam provides information on the static mechanism and a starting point for defining a target to use for controlling trans-cis photoisomerization. The MEP (unconstrained relaxation path leading from the FC point to the minimum energy point on S1) of the model trimethine cyanine with CASSCF and a 6-31G(d) basis set has

8720


Allan et al.

Figure 6. Direction of the geometrical displacement induced by the initial momentum vector used to direct the wavepacket to the conical intersection seam at large twist angles (upper panel) vs low twist angles, N1-C2-C4-N5 ∼ 0° (lower panel).

been fully described elsewhere.6 In this work we have used the 3-21G basis set. Additionally, our dynamics calculations employ state-averaged orbitals at all times during the simulation. We therefore located the symmetry restricted critical points (saddle points) identified previously with state-averaged orbitals and at the 3-21G basis set level to confirm that the PES involved in our simulations is essentially the same as that already described (see Supporting Information). As has already been discussed (see Figure 3), the initial stage of the MEP is described by initial stretching motions coupled with conrotatory torsional motion (with respect to the two central carbon-carbon bonds) forming structures with pseudo-C2 symmetry. In contrast, disrotatory motion would induce C2Vto-Cs symmetry breaking. The S1 surface in the FC region is convex (ridge-like) for the conrotatory torsional coordinate.6 This implies that the initial relaxation (totally symmetric stretching) on the ridge will rapidly couple to the unstable torsion down one side of the ridge, thus inducing symmetry breaking. Following this, the conrotatory torsion couples to the disrotatory one to become a one-sided torsion, and the path continues to a minimum on S1. The minimum conical intersection lies nearby. Previous work7 suggested that there is an energetically accessible segment of the seam that is approximately parallel to the S1 MEP at the end of it when plotted in the space of the torsional coordinate and skeletal deformations involved in the initial relaxation and in the branching space (see Figure 1). This feature was proposed to be the key element for controlling cyanine trans-cis photoisomerizations. The theoretical prediction was confirmed in recent optimal control experiments. Dietzek et al. have shown how stimulation of a particular vibrational mode can reduce the photoisomerization yield of 1,1′-diethyl-2,2′-cyanine iodide.8 The results were consistent with a picture in which the conical intersection seam lies parallel to the MEP where excitation of modes orthogonal to the reaction coordinate (torsion) can induce decay to the trans product. Intuitively, the key to control would be to hinder/excite twisting to access a planar/twisted crossing and thus form corresponding trans/cis ground-state products. Targeting a twisted crossing is not necessary. As has already been stated, after the initial in-plane relaxation has taken place, twisting dominates the reaction coordinate so that in the absence of

control the system will naturally decay down the MEP to a twisted conical intersection. Thus, the only way to avoid twisting is to excite orthogonal motions such that the dynamical pathway can deviate from the MEP. This was confirmed by a large number of unsuccessful attempts at targeting a flat crossing. To illustrate this, we present one simulation where the system decays down the MEP and reaches the seam at twisted geometries. In the second example, an initial momentum is used to target a region of the seam where the molecule is essentially flat (N1-C2-C4-N5 ∼ 0°). This was achieved by using an initial momentum vector parallel to the geometrical displacement from the FC point to a point on the seam that had a total dihedral angle of 0°. In an additional attempt to confirm that twisting will occur unless specifically prevented, we systematically reduced the magnitude of the momentum vector that had led to the “flat” region of the seam. This resulted in a larger average twist angle of the crossing geometries. Initial Momentum Construction. The two limiting cases presented below correspond to different sets of initial conditions. These correspond to an extra momentum, P0, given to the wavepacket at time t ) 0, such that ψ(0)(Q,0) ) 0 and ψ(1)(Q,0) ) χ0(0)(Q) exp(iP0 · Q/p) (see section 4). The first case corresponds to an initial momentum chosen to reach the seam at geometries with large twist angle, N1-C2-C4-N5 (see Figure 6, upper panel). For the second case, we gave an initial momentum designed to direct the wavepacket to a region of the conical intersection seam where the total dihedral angle is close to 0° (see Figure 6, lower panel). This approach aims to simulate experiment whereby the laser pulse is chosen to target a twisted vs quasi-planar point on the seam, with the aim of forming the cis product vs regenerating the trans reactant (see Figure 1). (Note that the arrows in Figure 6 are not the momentum vector itself but correspond to the geometrical displacement, parallel to the initial velocity, induced by the initial momentum.) In spite of the nonzero initial momentum, the first case is representative of simulations that are given no initial momentum: after initial in-plane relaxation, the wavepacket quickly twists and reaches the seam at geometries with large twist angles. Due to the possibility of breaking symmetry in different directions,

Method of Analysis for Direct Quantum Dynamics we had to give a small initial momentum to ensure that all the basis functions evolved to the same conical intersection enantiomer (see section 4). This was necessary to describe the nonadiabatic event properly. The initial momentum contained components of three torsional normal modes (Q1, Q3, Q7) and a stretching mode (Q12). The total excess energy was 2.5 eV distributed equally between the four modes in terms of momentum components conjugated to frequency-mass-weighted coordinates. For the second case (quasi-planar conical intersection target), Figure 1 suggests that motions involved in the initial relaxation direction (mostly in-plane skeletal deformation) might be an appropriate choice to excite. However, this strategy was unsuccessful (due to the fact that enantiomer control is lost if no out-of-plane motion is given in the initial momentum vector as well as requiring very long simulation times). Instead, construction of an appropriate initial momentum vector was achieved by targeting a specific geometry on the conical intersection seam. A quasi-planar structure with a dihedral angle of 0° was chosen. It was the result of an energy gap minimization on a geometry taken from a dynamics simulation that was given no initial momentum. The geometry has one of the NH2 groups twisted and the adjacent methine hydrogen and the central methine hydrogen projecting out of the plane of the N-C-C-N backbone. The geometry of this point was converted to frequency-mass-weighted normal coordinates and used as the direction of the initial momentum vector. The vector has components along every normal mode including torsional modes. These torsional modes are necessary not only to ensure the diabatization is correct but also to force the relevant hydrogen atoms out of the N1-C2-C4-N5 plane. To force the wavepacket close to the seam, a large excess energy (6 eV) was required. This suggests that even though the energy of the geometry chosen as the target lies only 1.8 kcal mol-1 above the FC point, there is a significant barrier separating the two along the path followed by the wavepacket (such a large excess energy may be seen as acting against the entropic difficulty of finding the optimal path to the target: the wavepacket is sent along a direct path that cuts through energetic barriers). This was confirmed by examining the energetic path of the S1 components of the basis functions: all trace pathways exhibiting a barrier of more than 3 eV (69 kcal mol-1) early in the simulation, before the crossing is reached. (This suggests that a better way of reaching this point on the seam might be to inject momentum later in the simulation after the initial relaxation. This strategy has indeed been demonstrated experimentally for the molecule 1,1′-diethyl-2,2′-cyanine iodide8 and suggested from DD-vMCG simulations applied to benzene photochemistry.20) The point targeted on the seam is structurally related to a local minimum on the S1 surface that was located for longer Schiff bases.15 A stationary point (transition state) corresponding to this minimum can also be found for the trimethine cyanine examined here and is 17.8 kcal mol-1 below the FC point (see Supporting Information) adding weight to the suggestion that delaying stimulation until after this point is reached might provide an alternative route to this crossing point. The key structural difference between the local S1 stationary point and the targeted seam point is the position of the central methine hydrogen. In the S1 transition state it is in plane while on the seam it extends out of the plane of the carbon backbone. The importance of the relative positions of the two hydrogen atoms attached to the isomerizing bond has attracted much interest in


Figure 7. Twisted conical intersection target: variation of the total dihedral angle (N1-C2-C4-N5) over time for each basis function (blue, 1; red, 2; orange, 3; green, 4).

the study of longer Schiff base models (see ref 34 and references therein). It has been suggested that hydrogen-out-of-plane (HOOP) motion can induce decay to S0 in the absence of large backbone twist angles (i.e., causing decay back to the reactant). This is exactly what is observed here: the (energetically accessible) target geometry exhibits a large HOOP deformation and a small twist angle (0°). Crucially, this deformation is not induced by the gradient if the wavepacket is simply placed at the FC point. Geometries Visited. The four basis functions that make up the global wavepacket can be viewed as four coupled trajectories, each of which contribute varying amounts to the global wavepacket and reach the seam at a range of geometries and over a range of times. In the examples shown here the time does not correspond to any measurable time because the wavepackets in the simulation have been given a fairly large initial momentum. However, it is the relative times that different basis functions reach the seam that we wish to highlight. In the first case (target: twisted conical intersection), the geometries at which S1-S0 transfer occurs span a region around a twisted conical intersection geometry. This is illustrated in Figure 7, which shows how the total dihedral angle of the four quantum trajectories (centers of each GBF) varies over the duration of the simulation. Each basis function reaches the seam at a similar (but not identical) geometry, as illustrated in Figure 8. After reaching the seam, the spread increases. Figure 7 suggests that GBF1, -3, and -4 continue twisting after reaching the conical intersection seam, giving the cis product. In contrast, GBF2 returns to the trans reactant after reaching the seam. Any further interpretation requires determining the weight of each GBF in the wavepacket expansion. To this end, we will use the GGP of each basis function (i.e., each quantum trajectory) for each electronic state over time [eq 10]. For now, we will quote the corresponding numbers, but a more detailed analysis will be given in the next subsection. GBF2 accounts for approximately 9% of the total wavepacket and on reaching the seam almost all of this is transferred to the ground state (7% of the total wavepacket). This implies that the progress of GBF2 is a good description of the trajectory of molecules that decay to the ground state at this time and geometry. Again, for the second case (target: quasi-planar conical intersection), the geometries of the structures at the time that transfer occurs span a range of geometries but this time around a seam point that has a N1-C2-C4-N5 dihedral angle of 0°. Three basis functions exhibit this behavior and progress to geometries with large twist angles (>120°), whereas GBF2 crosses later, at an already twisted geometry (see Figure 9). As expected, after the crossings, when the initial stretching that leads the wavepacket to the “flat”-geometry region of the seam has ended, the wave function begins to follow the PES

8722


Allan et al.

Figure 8. Twisted conical intersection target: selected geometries of the four quantum trajectories over time.

downhill. As has been discussed previously, the end of the MEP involves torsion of the total dihedral bond and in the absence of any countering motion the wavepacket will twist. By the end of the simulation, three of the four functions have geometries much like that of the twisted S1 minimum while only the fourth returns to the flat trans conformer. It must be stressed here that each quantum trajectory “feels” both PES at the same time in proportion to their weights on each electronic state. The population transfer is not 100%, so that quantum trajectories after the crossing still feel the influence of S1 as well as S0 (similar to a quantum superposition of Ehrenfest semiclassical trajectories). Again, for further interpretation, we need the weight of each GBF in the wavepacket expansion given in terms of GGP [eq 10]. The time evolution of the GGP will be analyzed in more detail in the next subsection, but we quote some significant values here. The S1 components of the twisted functions account for more than half of the contributions of those functions to the global wavepacket so the trajectories of those functions are mostly driven by the S1 PES. It is also worth noting that GBF4, which returns to the trans conformer, although only accounting for 7% of the global wavepacket, transfers more than half of its population to S0. Thus it is this function that is most representative of the dynamics on S0. We present here calculations with four basis functions simply for ease of visualization. Similar simulations were performed with up to 16 basis functions, and the results do not differ significantly. When the number of basis functions is increased to 8, the results are similar although the range of geometries sampled by the individual basis functions is of course larger. A

total of 27% of the global wavepacket is accounted for by a single basis function, most of which (23% of the global wavepacket) is transferred to S0 on reaching the seam. As in the 4-GBF example presented above, this function retains the trans configuration, suggesting that this function is most representative of the dynamics on S0. The total S1-S0 transfer is 54%, confirming that convergence had not been reached with four basis functions. Some results produced with 16 basis functions are shown in Supporting Information. As the number of basis functions is increased, CASSCF convergence becomes increasingly difficult. This is due to a greater chance of exploring geometries where the active space, which was chosen to be stable around the reaction path, may become unsuitable. The initial momentum vectors had to be slightly modified to produce simulations of adequate length. This was done in an ad hoc way, by visual inspection of the geometries and reduction of components that appeared to be causing the convergence problems. Note that the initial conditions and subsequent results of calculations with different numbers of GBF are similar enough qualitatively to lead to consistent conclusions. Analysis of the Quantum Trajectories in Semiclassical Terms. As illustrated above, we need to analyze how much each quantum trajectory contributes to the total wavepacket to understand the quantum simulations in semiclassical terms. To this end, we computed the GGP of each basis function (i.e., each quantum trajectory) for each electronic state over time [eq 10] and analyzed their time evolution. We start with the case corresponding to the twisted conical intersection target. On the approach to the seam (where the



Figure 9. Quasi-planar conical intersection target: selected geometries of the four quantum trajectories over time. The time/geometry at which each basis function reaches the seam and transfers population to S0 is highlighted.

Figure 11. Quasi-planar conical intersection target: gross Gaussian populations of the S0 components of each of the four basis functions and total S0 population over time (blue, 1; red, 2; orange, 3; green, 4; purple, total).

Figure 10. Twisted conical intersection target: (upper panel) contributions (gross Gaussian populations) of each basis function to the global wavepacket over time; (lower panel) gross Gaussian populations of the S0 components of each of the four basis functions and total S0 population over time (blue, 1; red, 2; orange, 3; green, 4; purple, total).

functions overlap significantly) the contributions of the four basis functions fluctuate (see Figure 10; upper panel). After the functions reach the seam and then spread further, the contributions become stable. By the end of the simulation, basis function 4 accounts for 42% of the global wavepacket, basis functions 1 and 3 for 24% each, and basis function 2 for 10%.

As the four basis functions evolve over the surface, each reaches the seam at a different time and at a different geometry. The total population transfer to S0 is not immediate but occurs over a range of time and gradually reaches a maximum value of 68%: 14% via GBF1, 7% via GBF2, 17% via GBF3, and 30% via GBF4 (shown in Figure 10; lower panel). The same analysis can be done for the second case corresponding to the quasi-planar conical intersection target. The total S1 to S0 transfer is approximately 28% (see Figure 11). The first basis function approaches the seam at 16 fs and there is a population transfer of ca. 12% from the component on S1 to the component on S0. Similarly, GBF2 is responsible for a 5% transfer starting at 23 fs, GBF3 for a 5% transfer at 29 fs, and GBF4 for a 6% transfer at 14 fs. The Rabi-like oscillations seen in the electronic populations are due to the diabatic populations as approximate adiabatic populations. The actual population of the adiabatic S0 state

8724


Allan et al.

Figure 12. Twisted conical intersection target: S1/S0 adiabatic energy differences (∆V ) V1 - V0) for the trajectory followed by the center of each basis function (top-left, 1; top-right, 2; bottom-left, 3; bottom-right, 4) over time (red); total contributions (gross Gaussian populations) of each individual basis function to the global wavepacket (green), and corresponding S0 and S1 populations (blue and orange, respectively) over time.

should be a smoother function of time but with a similar magnitude. This is discussed in more detail in section 4 and Appendix A. Such oscillations occur because of a nonzero offdiagonal coupling between the two diabatic electronic states. The oscillation period is proportional to the inverse of the energy difference between the corresponding adiabatic energies (∆V ) V1 - V0). For example, an energy difference ∆V ) 3.0 eV corresponds to an oscillation period equal to 1.4 fs. For a given difference between the two diagonal elements of the diabatic Hamiltonian matrix, the amplitude of the oscillations increases monotonically with an increase of the off-diagonal element. In simulations with fewer basis functions, we noticed that quantum trajectories explored regions with larger values of the off(01) ˆQ , hence producing larger amplitude diagonal coupling H oscillations. Local Character of the Population Transfer. The time evolution of the GGP can be analyzed also to determine whether the transfer of population is local or nonlocal. As already discussed in section 3, quantum dynamics is nonlocal in space, but a small nonlocality means that a semiclassical approach is a valid approximation. The extent of radiationless decay from S1 to S0 is correlated to the difference in adiabatic energy, ∆V, between the two components of each quantum trajectory. Population transfer from components of the basis functions on S1 to components of basis functions on S0 occurs approximately at a minimum in ∆V where ∆V is close to zero, indicating at what time each basis function reaches the seam. Figure 12 shows the evolution of the ∆V and the associated population transfer to S0 for the first case (twisted conical intersection target). As the difference in adiabatic energy ∆V between the two components of each basis function approaches zero, the basis function nears the conical intersection seam where the states are degenerate. In this case the first basis function approaches the seam at 17 fs, the second at 18 fs, the third at 23 fs, and the fourth also at 17 fs. The values of ∆V at which population transfer to S0 occurs are 2.4, 1.4, 3.5, and 7.3 kcal mol-1 for GBF1-4, respectively. That ∆V does not reach exactly zero indicates that the centers of the basis functions come very close to the seam but do not actually cross it. Nevertheless, the quantum description of the system allows population transfer to occur.

A more subtle feature of Figure 12 is the relationship between the S1, S0, and total populations of the individual basis functions. Where the populations of the S1 and S0 components of an individual basis function are mirror images of each other, the total contribution of that function to the global wavepacket is necessarily constant. Population transfer is “local” in space (population is only transferred between S1 and S0 between components of the same basis function; see Figure 5). When the populations of the S0 and S1 components of a basis function are not mirror images, there may be population transfer to components of one or more of the other basis functions (even if the total contribution of the basis function in question remains constant). A method of analysis such as that presented in this work is required to assess the quantum effect described by transfer of population from a component of a basis function on S1 to a component of another basis function on S0 (nonlocal transfer in space). Analysis of Figure 12 indicates limited nonlocal transfer. In this example the mirror imaging is always accompanied by a constant total population. An increase (decrease) in the total population of a basis function and the populations of both the components indicates nonlocal transfer (in space) to (from) this basis function. This is observed between 18 and 20 fs for GBF1 and GBF4. During this time both the S0 and S1 components of GBF1 (GBF4) and the total population are increasing (decreasing). The S1 component on GBF4 is decreasing at a faster rate than the GBF4 component on S0, suggesting transfer from GBF4 S1 to GBF1 S0. This quantum effect cannot be described using classical trajectories. However, the effect is small in this case, which explains why semiclassical surface hopping trajectories have been used successfully to describe this system in the past. Similar conclusions hold for the second case corresponding to the quasi-planar conical intersection target (see Figure 13). Intermediate Targets: Systematic Analysis. Finally, to demonstrate that control can be modulated over a range of twisted geometries on the seam (intermediate targets), we repeated the simulation but systematically reduced the excess energy given to the initial wavepacket. In all other respects the simulations were identical to the second example presented above.



Figure 13. Quasi-planar conical intersection target: S1/S0 adiabatic energy differences (∆V ) V1 - V0) for the trajectory followed by the center of each basis function (top-left, 1; top-right, 2; bottom-left, 3; bottom-right, 4) over time (red); total contributions (gross Gaussian populations) of each individual basis function to the global wavepacket (green), and corresponding S0 and S1 populations (lavender and orange, respectively) over time.

Figure 14. Upper panel: variation of the total dihedral angle (N1-C2-C4-N5) of crossing geometries of individual basis functions with the initial excess energy. Crossing geometries correspond to the geometries at which population is transferred to S0 and ∆V exhibits a minimum. Bar widths correspond to the amount of population transferred to S0 by a particular basis function (as a percentage of the global wavepacket). Lower panel: variation of the weighted mean total dihedral angle (N1-C2-C4-N5) of crossing geometries averaged over individual basis functions with the initial excess energy.

Figure 14 indicates that an initial excess energy of 5 eV is in fact the most efficient value for inducing decay at low twist angles (with the initial momentum vector used here). In that simulation, 44% of the wavepacket is transferred to S0 (compared with 28%, 43%, 36%, 34%, and 27% when the excess energy is 6, 4.5, 3, 2.5, and 2 eV, respectively). When the initial excess energy is lower than 5 eV, the twist angle of the basis function responsible for the greatest population transfer (within

each wavepacket) to S0 increases (from 3.7° at 5 eV to 4.0° (4.5 eV), 31.6° (3.5 eV), 55.4° (2.5 eV), and 88.8° (2 eV)) as the initial excess energy is reduced. The weighted mean twist angle at the crossing geometry, shown in Figure 14 (lower panel), is not smoothly increasing with decreasing initial excess energy but does indicate that lower initial excess energies result in crossings at more twisted geometries. This comparison may be limited by the moderate number (4) of basis functions employed. Additionally, the reference structure used as the origin for extrapolation of the diabatic states was the same for all simulations in the series (to provide consistency) but will inevitably become less accurate in cases where the crossing geometries are highly twisted. As expected, as the initial energy is reduced, the time taken for basis functions to reach the seam increases: the mean crossing times increase from 32 fs with 5 eV initial excess energy to 63 fs with 2 eV (weighted mean crossing times are very similar). The range of twist angles of the basis functions in each simulation varies between 64.5° and 95.0°, confirming that the basis functions spread over a large area of the space even when a large initial momentum is used to direct the reaction outcome. This reflects the quantum spreading of the wavepacket along an unstable direction of the PES orthogonally to its main direction of travel. 6. Conclusions Direct quantum dynamics is an efficient approach for understanding the photochemistry of organic molecules. In the DD-vMCG method used here, the global wavepacket that represents the time evolution of the photoexcited molecule is made up of a collection of coupled Gaussian functions. This is a complicated function that is much delocalized in the space of the nuclear coordinates. We have presented a new way of analyzing the simulations based on a Mulliken-type population analysis, which provides a straightforward interpretation of the wavepacket in terms of weighted, coupled quantum trajectories. The method allows direct analysis of these in much the same way as semiclassical trajectories are usually analyzed. We have

8726


introduced a new quantity, the GGP, for gross Gaussian population, inspired from Mulliken’s gross orbital population. For each basis function, the GGP proved to be an effective indicator of its actual weight within the global density of probability. Contributions on each electronic state could also be distinguished. We applied this method to the 12-atom model system trimethine cyanine. Studying the trans-cis photoisomerization of cyanines raises interest not only for their potential use as fluorescent dyes for biomedical imaging but also because they are models of larger chromophores involved in light-induced processes in biological systems. This has prompted much theoretical and experimental work to explore possibilities for intelligent optimal control of cyanines photodynamics. In this work, we have shown that trans-cis isomerization will occur unless specifically prevented. The actual key to control is thus the stimulation of the in-plane skeletal deformations, which forces radiationless decay at flat geometries earlier on the reaction coordinate (expected to lead to regeneration of the trans reactant). In the absence of this, decay occurs at the twisted conical intersection (expected to lead to formation of the cis product). Thus, one needs to decrease the momentum in the skeletal deformation coordinates to prompt decay near the minimum conical intersection at large twist angles. In this work, we did not simulate a control experiment, but we identified (i) a tunable target somewhere between reactant and product (a range of molecular geometries where the reaction path can meet the seam of conical intersection) and (ii) the underlying mechanism able to tune between two extreme behaviors (reactant regeneration vs product formation). We thus suggested, or confirmed, that experiments should play with the amount of excitation given to the in-plane skeletal deformations to change the cis:trans branching ratio. Appendix A: Rabi-Like Oscillations in a Two-Level System This appendix gives some insight into the physical origin of the Rabi-like oscillations that can be seen in the diabatic electronic populations. At the conical intersection the adiabatic states have near infinite coupling through the nuclear kinetic energy operator. However, the derivative coupling is quite localized around the crossing and decreases fast. In contrast, this term is zero (small) everywhere by construction in the (quasi)diabatic representation, but the potential-like coupling term (off-diagonal term of the electronic Hamiltonian matrix) can take finite values. As a crude approximation, we can neglect the influence of the nuclear coordinates far from the crossing and consider that the adiabatic populations correspond to almost uncoupled electronic states, while the diabatic states have an almost constant coupling. Transformation between Eigenstates and Noneigenstates. ˆ . {|0〉,|1〉} is We consider a two-level system of Hamiltonian H a generic basis set assumed to be real-valued for simplicity. The corresponding Hamiltonian matrix reads

[

H00 H01 H) H10 H11

]

where H10 ) H*01 ) H01. The eigenstates noted {|0a〉,|1a〉} are chosen real-valued and correspond to

Allan et al.

Ha )

[

V0 0 0 V1

]

The unitary transformation between eigenstates and noneigenstates,

{

|0a〉 ) cos θ|0〉 + sin θ|1〉 |1a〉 ) -sin θ|0〉 + cos θ|1〉

can be represented as a rotation matrix

[

cos θ -sin θ sin θ cos θ

U)

]

such that the two Hamiltonian matrices are related through

Ha ) U†HU This leads to

[

]

H00 + H11 V0 0 ) 1+ 0 V1 2 H11 - H00 H11 cos 2θ + H01 sin 2θ 2 H11 H11 - H00 sin 2θ + H01 cos 2θ 2

[

- H00 sin 2θ + H01 cos 2θ 2 - H00 cos 2θ - H01 sin 2θ 2

]

To make this matrix diagonal, the rotation angle must satisfy

tan 2θ ) -

2H01 H11 - H00

unless H11 ) H00, in which case the angle must be set to π/4 (modulo π/2) if H01 * 0 or is undetermined if the off-diagonal term vanishes too. Further, we get

sin 2θ ) -ε cos 2θ ) ε

2H01

√(H11 - H00)2 + 4H012 H11 - H00

√(H11 - H00)2 + 4H012

with undetermined absolute signs (ε ) (1). The two corresponding eigenvalues read

{

H00 + H11 -ε 2 H00 + H11 V1 ) +ε 2

V0 )

( (

H11 - H00 2 H11 - H00 2

) )

2

2

+ H012 + H012

A phase convention is set according to V0 < V1, so that ε ) 1. Time Evolution of the Populations. Let us consider a wavepacket expanded in the eigenbasis set,

|ψ(t)〉 ) aa0(t)|0a〉 + aa1(t)|1a〉

Method of Analysis for Direct Quantum Dynamics The time-dependent Schro¨dinger equation reads

[ ] [ ]

a aa0(t) d a0(t) ) Ha a a dt a1(t) a1(t)

ip

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8727 ˜ 2 ) (pΩ01)21. This leads to a simple form for the evolution H operator, ∞

˜ t/p -iH

e

)

1 t ∑ (-1)n (2n)! (p)

(pΩ01)2n1 -

n)0

∞

and gets trivial solutions,

i

As expected for eigenstates, the populations are stationary,

Pa0(t) ) |〈0a |ψ(t)〉| 2 ) |aa0(t)| 2 ) Pa0(0) ) |〈1 |ψ(t)〉| ) 2

|aa1(t)| 2

2n+1

) cos(Ω01t)1 - i sin(Ω01t)

aa1(t) ) aa1(0)e-iV1t/p

a

∑ (-1)n (2n +1 1)! ( pt )

˜ (pΩ01)2nH

n)0

aa0(t) ) aa0(0)e-iV0t/p

Pa1(t)

2n

)

Pa1(0)

˜ H pΩ01

In turn, the populations read P0(t) ) cos2(Ω01t)(1 - R) - sin(2Ω01t) sin(2θ)T β + sin2(Ω01t)[cos2(2θ)(1 - R) + sin2(2θ)R + sin(4θ)Rβ] P1(t) ) cos2(Ω01t)R + sin(2Ω01t) sin(2θ)T β +

Now, let us consider the general case where the Hamiltonian matrix is not diagonal. The wavepacket reads

|ψ(t)〉 ) a0(t)|0〉 + a1(t)|1〉 Intermediate coefficients corresponding to shifted energies are introduced,

sin2(Ω01t)[cos2(2θ)R + sin2(2θ)(1 - R) - sin(4θ)Rβ]

where R ) P1(0) is the initial population of state number 1, and β ) a0*(0)a1(0) is the initial coherence between the two states. Note that P0(0) ) 1 - R. R and T stand for the real part and imaginary part, respectively. In addition, we used (see above)

a0(t) ) b0(t)e-i[(H00+H11)/2p]t

H01 pΩ01 H11 - H00 cos 2θ ) 2pΩ01

sin 2θ ) -

-i[(H00+H11)/2p]t

a1(t) ) b1(t)e

as well as a traceless Hamiltonian matrix,

[

-

˜ ) H

H11 - H00 2

H01 H11 - H00 2

H01

]

It is trivial to show that the time-dependent Schro¨dinger equation can be recast as

ip

[ ] [ ]

b (t) d b0(t) ˜ 0 )H b b1(t) (t) dt 1

which also reads

[ ]

P0(t) )

[ ]

b0(t) b0(0) ˜ ) e-iHt/p b1(t) b1(0)

)

˜, We introduce -(pΩ01) as the determinant of H 2

|

-

H11 - H00 2 H01

H01 H11 - H00 2

|

(

The populations of the noneigenstates oscillate at twice the Rabilike angular frequency, i.e., ω01 ) 2Ω01 ) (V1 - V0)/p, the Bohr angular frequency between the two levels V1 and V0. For example, V1 - V0 ) 3.0 eV corresponds to an oscillation period equal to 1.4 fs. All other things being equal, the amplitude of the oscillations is determined by the strength of the coupling (the limiting case of uncoupled states corresponds to eigenstates with stationary populations). To make this obvious, let us assume R ) 1 and β ) 0. The transfer of population from state number 1 to state number 0 reads

H11 - H00 )2

)

2

( ) H01 pΩ01

(

H11

2

sin2(Ω01t)

H012 sin2(Ω01t) - H00 2 + H012 2

)

Appendix B: Diabatization Scheme

-

H012 ) -(pΩ01)2

Ω01 is the Rabi-like angular frequency of the system. Note that pΩ01 ) (V1 - V0)/2. The Cayley-Hamilton theorem shows that

This appendix provides some details on the diabatization scheme used in our direct quantum dynamics simulations, which corresponds to a simple form of the regularized diabatic states method of Ko¨ppel et al.32,33 This method makes use of the fact that the singularities of the adiabatic representation are only due to the linear derivative coupling terms. In view of this, the regularized transformation from the adiabatic basis set to the quasi-diabatic one is obtained

8728


Allan et al.

by only removing the first-order terms in the Taylor expansion of the energy difference, and the rest are neglected because they are small. The corresponding two-level molecular Hamiltonian operator in the regularized quasi-diabatic representation reads

where the value of the mixing angle φ will be set below. The current point corresponding to Q is represented by the positionvector b R (in the same frame as b R0 and in mass-weighted Cartesian coordinates). Doing so, the elements of W(1)(Q) are defined as linear extrapolations of the quasi-diabatic energies from the reference conical intersection point,

∆(Q) (1) (ss') ˆQ W (Q) [H ] ) TˆQ1 + Σ(Q)1 + (1) ∆ (Q) where TˆQ is the kinetic energy operator acting in the space of the nuclear coordinates Q, Σ(Q) ) [V0(Q) + V1(Q)]/2, ∆(Q) ) [V1(Q) - V0(Q)]/2, V1(Q) and V0(Q) are the adiabatic potential energy functions calculated by the quantum chemistry program,

W(1)(Q) )

[

-k(Q)X1 λ(Q)X2 λ(Q)X2 k(Q)X1

]

The significance of k(Q)X1 and λ(Q)X2 will be given later. The first-order mixing angle (see Appendix A) is defined from

λ(Q)X2 k(Q)X1

such that

b0〉 ) cos φ|00;R b0〉 + sin φ|10;R b0 〉 |0;R 0 0 0 b 〉 ) -sin φ|0 ;R b 〉 + cos φ|10;R b0〉 |1;R

This is precisely why only the branching plane spanned by the two linearly independent branching-space vectors matters physically. In our method, we make use of this to ensure that the quasi-diabatic states coincide with the adiabatic states at the FC point. For this, the angle φ is chosen so as to make the off-diagonal term of W(1)(Q) vanish at Q ) QFC,

bFC - b λ(QFC)XFC x 2 · (R R0) ) 0 2 ) b

[

]

0 -∆(1)(Q) ) U(1)†(Q) W(1)(Q) U(1)(Q) (1) 0 ∆ (Q)

which implies

tan 2φ ) -

and

[

The two degenerate electronic states at the conical intersection are defined up to an arbitrary mixing angle that plays the role b0〉 and |10;R b0〉 of a gauge for the diabatization. If we label, |00;R the two adiabatic states that were actually used to produce b x10 and b x20 in the quantum chemistry calculation, then b x1 and b x2 as defined above correspond to another equivalent pair of adiabatic b0〉, such that b0〉 and |1;R states, |0;R

{

∆(1)(Q) ) √[k(Q)X1]2 + [λ(Q)X2]2

tan 2θ(1)(Q) ) -

b-b k(Q)X1 ) b x 1 · (R R 0) b-b λ(Q)X2 ) b x 2 · (R R0)

]

[

]

V0(Q) 0 ∆(Q) (1) W (Q) U(1)(Q) ) U(1)†(Q) Σ(Q)1 + (1) V1(Q) 0 ∆ (Q)

Thus, the regularized scheme is equivalent to using the exact adiabatic potential energy functions V1(Q) and V0(Q) of the twolevel molecular Hamiltonian operator in the adiabatic representation. However, the effective nonadiabatic coupling derived from the action of TˆQ on θ(1)(Q) is not exact, and this defines the quality of the approximation. In W(1)(Q), X1 and X2 are rectilinear displacements along the tuning and coupling modes, respectively. The latter two are chosen as linearly independent mixtures of the two branchingspace vectors9-12 calculated at a given conical intersection chosen as a reference point (represented by the position-vector b R0 in mass-weighted Cartesian coordinates). The original gradient-half-difference and interstate-coupling vectors given x10 and b x20, respectively (in in the same frame as b R0 are noted b mass-weighted Cartesian gradient components, where the square-root of the masses multiply the denominators rather than the numerators as in mass-weighted Cartesian coordinates). The rotated tuning and coupling modes are expressed as

{

b x 1 ) cos 2φx b01 - sin 2φx b02 b x 2 ) sin

2φx b01

+ cos

2φx b02

bFC - b b x 02 · (R R0) bFC - b b x 01 · (R R 0)

In the original method proposed by Ko¨ppel et al., the reference point b R0 used for the transformation can change and span the seam of conical intersection. In contrast, we use a simplified version with a unique reference that is set in the input file generated by the user of the code before any simulation. This may imply some “try and guess” because a specific choice may affect the results if the wavepacket does not pass anywhere near the reference point when it crosses the seam. In addition, giving the geometry of the reference point in terms of its coordinates means that nuclei are labeled. This means that equivalent geometries corresponding to a nuclear permutation (enantiomers or enantiotopes) will not be treated in the same way with respect to the diabatization scheme (as already discussed in the main text). Acknowledgment. Special thanks are due to David Mendive Tapia for collaboration in conceptual and technical aspects of this work. We are grateful to Dr. Michael J. Bearpark, Marta L. Arau´jo, Dr. David Asturiol Bofill, and Dr. Lluı`s Blancafort for helpful discussions. Support from the EPSRC-U.K. (Grant No. EP/F028296/1) is acknowledged. Supporting Information Available: Cartesian coordinates and energies of all optimized geometries; frequency-massweighted components of the initial momentum vectors and

Method of Analysis for Direct Quantum Dynamics corresponding Cartesian coordinates giving the directions of the induced geometrical displacements shown in Figure 6; plots of the time evolution of ∆V, individual S0 GGP, and total S0 population of a 1-GBF case crossing at a twisted geometry and of a 16-GBF case crossing at a flat geometry. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Lasorne, B.; Robb, M. A.; Worth, G. A. Phys. Chem. Chem. Phys. 2007, 9, 3210. (2) Worth, G. A.; Robb, M. A.; Lasorne, B. Mol. Phys. 2008, 106, 2077. (3) Lasorne, B.; Worth, G. A. Direct Dynamics with Quantum Nuclei. In Multidimensional Quantum Dynamics: MCTDH Theory and Applications; Meyer, H.-D., Gatti, F., Worth, G. A., Eds.; Wiley-VCH: Weinheim, 2009; p 113. (4) Schoenlein, R. W.; Peteanu, L. A.; Mathies, R. A.; Shank, C. V. Science 1991, 254, 412. (5) Gai, F.; Hasson, K. C.; McDonald, J. C.; Anfinrud, P. A. Science 1998, 279, 1886. (6) Sanchez-Galvez, A.; Hunt, P.; Robb, M. A.; Olivucci, M.; Vreven, T.; Schlegel, H. B. J. Am. Chem. Soc. 2000, 122, 2911. (7) Hunt, P. A.; Robb, M. A. J. Am. Chem. Soc. 2005, 127, 5720. (8) Dietzek, B.; Bru¨ggemann, B.; Pascher, T.; Yartsev, A. J. Am. Chem. Soc. 2007, 129, 13014. (9) Davidson, E. R. J. Am. Chem. Soc. 1977, 99, 397. (10) Desouter-Lecomte, M.; Dehareng, D.; Leyhnihant, B.; Praet, M. T.; Lorquet, A. J.; Lorquet, J. C. J. Phys. Chem. 1985, 89, 214. (11) Atchity, G. J.; Xantheas, S. S.; Ruedenberg, K. J. Chem. Phys. 1991, 95, 1862. (12) Yarkony, D. R. ReV. Mod. Phys. 1996, 68, 985. (13) Ben-Nun, M.; Martı´nez, T. J. AdV. Chem. Phys. 2002, 121, 439. (14) Belz, S.; Grohmann, T.; Leibscher, M. J. Chem. Phys. 2009, 131, 034305. (15) Garavelli, M.; Bernardi, F.; Celani, P.; Robb, M. A.; Olivucci, M. J. Photochem. Photobiol. A: Chem. 1998, 114, 109. (16) Dietzek, B.; Bru¨ggemann, B.; Pascher, T.; Yartsev, A. Phys. ReV. Lett. 2006, 97, 258301. (17) Dietzek, B.; Bru¨ggemann, B.; Persson, P.; Yartsev, A. Chem. Phys. Lett. 2008, 455, 13. (18) Improta, R.; Santoro, F. J. Chem. Theory Comput. 2005, 1, 215. (19) Dietzek, B.; Tarnovsky, A. N.; Yartsev, A. Chem 2009, 357, 54. Dietzek, B.; Christensson, N.; Pascher, T.; Pullerits, T.; Yartsev, A. J. Phys. Chem. B 2007, 111, 5396.

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8729 (20) Lasorne, B.; Bearpark, M. J.; Robb, M. A.; Worth, G. A. J. Phys. Chem. A 2008, 112, 13017. (21) Arau´jo, M.; Lasorne, B.; Magalha˜es, A. L.; Worth, G. A.; Bearpark, M. J.; Robb, M. A. J. Chem. Phys. 2009, 131, 144301. (22) Asturiol, D.; Lasorne, B.; Worth, G. A.; Robb, M. A.; Blancafort, L. Phys. Chem. Chem. Phys. 2010, 12, 4949. (23) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (24) Raab, A.; Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. J. Chem. Phys. 1999, 110, 936. (25) Klein, S.; Bearpark, M. J.; Smith, B. R.; Robb, M. A.; Olivucci, M.; Bernardi, F. Chem. Phys. Lett. 1998, 292, 259. (26) Amarouche, M.; Gadea, F. X.; Durup, J. Chem. Phys. 1989, 130, 145. (27) Heller, E. J. J. Chem. Phys. 1975, 62, 1544. (28) Hack, M. D.; Wensmann, A. M.; Truhlar, D. G.; Ben-Nun, M.; Martı´nez, T. J. J. Chem. Phys. 2001, 115, 1172. (29) Multidimensional Quantum Dynamics: MCTDH Theory and Applications; Meyer, H.-D., Gatti, F., Worth, G. A., Eds.; Wiley-VCH: Weinheim, 2009. (30) Worth, G. A.; Beck, M. H.; Ja¨ckle, A.; Burghardt, I.; Lasorne, B.; Meyer, H. D. The MCTDH Package, DeVelopment, Version 9.0; University of Heidelberg: Heidelberg, Germany, 2008. (31) Frisch, M. J. Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N. Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Parandekar, P. V.; Mayhall, N. J.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian DeVelopment Version, Revision H.01; Gaussian, Inc.: Wallingford, CT, 2009. (32) Ko¨ppel, H.; Gronki, J.; Mahapatra, S. J. Chem. Phys. 2001, 115, 2377. (33) Ko¨ppel, H.; Schubert, B. Mol. Phys. 2006, 104, 1069. (34) Sumita, M.; Ryazantsev, M. N.; Saito, K. Phys. Chem. Chem. Phys. 2009, 11, 6406.

JP101574B

A Straightforward Method of Analysis for Direct Quantum Dynamics

Recommend Documents