The damped Ehrenfest (D-Eh) method_ Application to non-adiabatic reaction paths

An implementation of the Ehrenfest method with damped velocity is discussed. The method is then applied to study the non-adiabatic reaction paths for two simple chemical systems: the isomerization of the allene radical cation in its excited state and the channel 3 photochemical transformation of benzene to benzvalene. For both systems the initial conditions for the Ehrenfest trajectory were either an adiabatic eigenstate with the geometry close to a conical intersection, or a superposition of eigenstates at the geometry close to a conical intersection. In allene we were able to show that the adiabatic reaction, which passes through a conical intersection, stimulates electron dynamics. In benzene we were able to show the importance of the phase at the conical intersection for the generation of the benzevalene intermediate.


Introduction
The concept of a reaction path is fundamental in mechanistic chemistry. Intuitively it corresponds to a sequence of chemical structures that connects reactants and products. As we discuss below the concept has a precise definition if one has a transition state. However, in photochemistry, the reaction path has 2 segments: one on the excited state, the other on the ground state. The transition between the 2 segments takes place at a conical intersection (CoIn) [1][2][3][4][5]. [6]. Here the potential surface is discontinuous and the wavefunction may become a superposition of eigenstates. The objective of this paper is to show how reaction paths can be computed even when the electronic wavefunction is not an eigenstate.
One method for following the reaction path is the intrinsic reaction coordinate (IRC). The IRC approach has been used extensively in quantum chemistry to follow the reaction path from a transition structure to local minima representing the reactants and products. This allows complicated reaction mechanisms involving multiple nuclear rearrangements to be understood easily. Starting at a transition state (TS), the initial step of an IRC path is in the direction of the normal coordinate corresponding to the imaginary frequency at the TS, denoted v. The IRC proceeds [7] via integration of the differential equation, where q are the mass-weighted Cartesian coordinates and s the coordinate along the IRC. For subsequent steps, v becomes the massweighted gradient vector. A full review of IRC and reaction dynamics is available [8] so we do not go into detail here. Alternatively the PES can be explored with classical trajectories. Dynamic reaction path (DRP) methods use trajectories to approximately follow the IRC. The IRC represents the vibrationless-rotationless path, therefore removing all kinetic energy in a DRP trajectory begun at the TS and heading in the direction of v will result in DRP following the IRC path exactly. Damped velocity-Verlet (DVV) [9] also approximates the IRC path using classical trajectories but instead of removing all kinetic energy, the magnitude of the velocity is fixed at a small non-zero value, approximating the IRC. The DVV scheme has been shown to be stable and particularly efficient for large systems [9].
When a system is in an electronically excited state, for example, after vertical excitation, the system may follow different paths on the PES that may take it to excited state products, or to ground state minima via passage through a CoIn. When trying to connect an excited state transition structure with reactants and (ground state) products, one cannot rely on IRC because an IRC is only strictly defined if the initial step is taken along the normal co-ordinate corresponding to the imaginary frequency at a TS, this breaks down at the CoIn. At a CoIn the initial reaction path has been defined using a constrained optimization procedure [10]. Schapiro et al. [11] follow the reaction path through a CoIn by using a surface hopping algorithm to hop between adiabatic states, the system only ever populating one adiabatic state at a time.
However, passing through a CoIn can give a superposition of states. An alternative is to use the Ehrenfest method since motion can occur on an adiabatic surface or on a superposition of surfaces. In the present work we use a simple analogue/variant of the DVV algorithm of Schlegel [9] in which we damp the magnitude of the velocity vector to give an overall magnitude of between 0.01 and 0.04 Bohr/fs. We demonstrate the utility of the damped Ehrenfest method by studying the isomerization of the allene radical cation in its excited state and the channel 3 photochemical transformation of benzene to benzvalene.

Theoretical method: the damped Ehrenfest method (D-Eh)
The Ehrenfest method is now an established technique for nonadiabatic dynamics [12-16 17-19 20,21]. In this work we use a complete active-space self-consistent field (CASSCF) formalism [22] with a configuration interaction (CI) Hamiltonian t H ( ) e n represented in a basis of configuration state functions. The orbitals are not optimized at each step so we refer to the method as CAS-CI. Otherwise all of the usual concepts of CASSCF apply (e.g. choice of active space). The CAS-CI wavefunction is propagated as a solution of the time dependent Schrodinger equation rather than being found as a stationary eigenstate. We have developed the main theoretical ideas elsewhere [22][23][24] as well as many applications some of which are reviewed in a research volume [25]. The Ehrenfest CAS-CI wavefunction propagation corresponds to electron dynamics. The nuclei are propagated using standard on-the fly dynamics where the gradient is computed from the Ehrenfest wavefunction. In this case the wavefunction is not a stationary state so the full expression for the gradient and second derivative must be used [26]. In this subsection we give only the essential details.
The most important point is that the Ehrenfest method yields a single (i.e. Ehrenfest) state irrespective of how many adiabatic states are involved. As a consequence there is only a single time-dependent potential surface. Because of this fact there are no discontinuities at a CoIn. Thus we can adapt the method of Schlegel [9] where the velocity is scaled to produce an approximate reaction path rather than a trajectory.
We now proceed to give the main ideas associated with the Ehrenfest propagation and the computation of the derivatives. At any time the Ehrenfest state corresponds to a superposition of eigenstates of the electronic Hamiltonian at that particular geometry. Thus neither the wavefunction or the orbitals are optimized. Thus we are building a reaction path based upon a non-stationary state wavefunction. The only difficulty with this is that one must use the general expressions for the gradient [26]. Many criticisms of the Ehrenfest method stem merely from the fact that the gradient is computed in an approximate manner rather than using these general expressions (described below). In other work we have shown that [23] the Ehrenfest method can be derived as an approximation to full quantum dynamics methods. But we emphasise that we are not concerned with dynamics here except that we use damped dynamics to constrain the kinetic energy to be approximately zero so that we generate reaction paths. In the method of Schlegel [9], the electronic motion was always constrained to be a solution that was an eigenstate. In this work, we allow electron dynamics via the Ehrenfest propagation as well.
We should also point out that are two possible approaches to damped dynamics. In the first, one could follow the gradient of the Ehrenfest state (as we have done in this work). If the Ehrenfest state has collapsed to an adiabatic state then this is the same as the method proposed by Schlegel [9]. In a second, alternative, algorithm one could use the Ehrenfest wavefunction to determine a surface hop as proposed by Schapiro et al. [11] in the context of a conventional IRC. In this case the gradient would always correspond to an adiabatic state. In this work we choose the first algorithm because we are interested in the effect of following the gradient of a superposition of states, corresponding to laser excitation or created after passage through a conical intersection.
In our current implementation the time-dependent electronic wavefunction is expanded in the basis of (the configuration state functions) CSFs. By gathering the expansion coefficients c t ( ) s n at time t n in the vector C t ( ) n : (3) Using matrix notation we have Thus, the CI vector is being propagated in discretized steps, with the orbital coefficients remaining constant.
The wavefunction given in Eq. (4) does not correspond to a stationary state and is not an optimized wavefunction. As a consequence one needs the general expression for the gradients and second derivatives as we now discuss.

ii) Derivatives
The nuclei are propagated under the gradient derived from the Ehrenfest wavefunction r t ( , ).
= dP dt r t H r R r t ( , )| ( ; )| ( , ) As discussed previously, the Ehrenfest wavefunction is not a stationary state and is not optimized. Thus the full equations for the gradient and second derivatives must be used [26]. In these equations we need to distinguish 3 types of variables: 1) X the orthogonal rotation of the orbitals among themselves as one displaces the geometry with the derivative X RI with respect to the nuclear co-ordinateR I 2) C the orthogonal rotation of the CI eigenvectors among themselves as one displaces the geometry with the derivative C RI with respect to the nuclear co-ordinate and 3) Y the re-orthogonalization of the orbitals as one displaces the geometry with the derivative S RI .
The gradient then involves E C , the gradient of the energy with respect to the rotation of the CI eigenvectors, E X the gradient of the energy with respect to the rotation of the orbitals and E Y the gradient of the energy with respect to the re-orthogonalization of the orbitals. The leading term in the gradient E R I is the gradient of the energy due to the change in the molecular Hamiltonian with nuclear geometry (The Hellmann-Feynman term). The gradient [26] is given compactly as.
For an optimized wavefunction that corresponds to a stationary state, E X and E C vanish because the orbitals and CI coefficients have been optimized. For an Ehrenfest state we must determine [27] and C X R R I I by solution of the coupled perturbed MCSCF equations The second derivatives have the general form given in Eq. (7), where the notation is the same as in Eq. (6), e.g. E CRJ is the mixed second derivative with respect to CI vector rotation and nuclear displacement.
Again these equations are simplified dramatically if the orbitals and CI coefficients are optimized i.e. E X and E C are zero.
Quantities such as X RI have to be obtained from the coupled perturbed equations which have to be solved for each nuclear displacement. These equations are similar to the CASSCF coupled perturbed equations [27] except that the CI vector rotation coefficients correspond to rotations between the Ehrenfest vector and its orthogonal complement.
iii) Implementation Our implementation of this method, in a development version of Gaussian [28], uses our CASSCF formulation of the Ehrenfest method described in Vacher et al [23,22]. A review of this method has also been given [24]. iv) Scaling the velocity At each step of the Ehrenfest propagation, the velocity vector is scaled such that the overall magnitude is constant. The value of this constant varies between calculations but is chosen to be between 0.01 and 0.04 Bohr/fs, a similar range to that used in DVV. This factor was chosen so that a large step would be possible without compromising energy conservation. For the computations reported in the next section on allene we used 0.01 and in 0.04 in Benzene.

Results and discussion
In this section, we report results of our damped Ehrenfest method just described. We must emphasize that we are not performing quantum dynamics per se. Rather we carry out non-adiabatic reaction path studies that are different from the usual IRC method [7] and are more general. On the one had we can begin at an arbitrary geometry. The gradient is either the gradient of an adiabatic state or the gradient of a superposition of such states depending on the initial conditions. We use the general expression for the gradient, which is correct whether or not the wavefunction is optimized or not. The propagation of the reaction path depends on three factors (i) the nature of the Ehrenfest wavefunction at that time step (ii) the gradient of the Ehrenfest wavefunction at that time step and (iii) the propagation of the nuclei using the electronic gradient and the usual classical trajectory equations. At each step, the kinetic energy is constrained to be almost zero by damping. Thus we have a damped trajectory on a time-dependent potential surface (corresponding to the Ehrenfest state which is, in general, a superposition of eigenstates). In more approximate Ehrenfest implementations, the gradient might have been chosen to be an average over the component states of the Ehrenfest superposition. We do not make such an approximation here. Further, we use a CAS-CI formulation of the Ehrenfest method. Thus we use all the CAS electronic states rather than including a small number of eigenstates.
In our results, which we will now discuss, we will see several types of non-adiabatic behavior that cannot be represented in a conventional IRC. A reaction path through a conical intersection is one example. However, in addition, in the examples we will now discuss, we will see situations where the electron and nuclear motion become asynchronous. This interplay between electronic coherence and the control of the passage through conical intersections has been discussed previously with different methods and under several points of view [29][30][31][32].   The allene radical cation is of interest because there are two almost degenerate cationic states with positive charges at alternate ends of the molecule. (We have also studied this system with the Quantum Ehrenfest method [12]). There have been several quantum dynamics studies on this system which is a simple "model" for intramolecular charge transfer [33][34][35].
The geometry is shown in Fig. 1 together with the torsional angle, which is the main geometric motion on this potential surface. The ground state minimum (the D 0 symmetry broken solution) has a torsional angle of approximately 45°. Upon excitation from D 0 to D 1 , one reaches an adiabatic potential surface which has a minimum at a CoIn with an approximately 90°twisted geometry. We have carried out two simulations on the system. Firstly we discuss an Ehrenfest trajectory starting on the excited D 1 adiabatic state at a 45°twisted geometry. Secondly, we have carried out a simulation where we begin, again at a 45°twisted geometry, but this time with the electron hole localized on one of the carbons. Thus, in this latter case, we start with a superposition of the two adiabatic states. For our simulations we used a CAS space with three electrons associated with the two allene CeC double bonds as shown in Fig. 2. Thus there are 2 doublet states where the single electron is localized in either the orbital shown in Fig. 2 or the equivalent orbital localized on the opposite carbon.
Let us consider the adiabatic simulation beginning on D 1 first. The main results are summarized in Fig. 3 and Fig. 4. In Fig. 3 we can see that the torsional angle increases passing through 90°to the conformer with a torsional angle of approximately 135°. (note: the small "wobbles" at longer times are a result of the fact that the kinetic energy is not exactly zero and do not indicate an instability in the Ehrenfest wavefunction). In other words the system appears to do a non-adiabatic isomerization passing through a CoIn at approximately 90°and ending up on the ground state. However, on closer inspection, we have a surprising effect shown in Fig. 4. In the upper half of Fig. 4 we show the absolute value of D 0 -D 1 energy gap, which indicates that we have passed through the CoIn at  5. Trajectory initiated on a superposition of D 0 and D 1 states of allene cation (non-stationary non-adiabatic state). (The Ehrenfest state has collapsed to D 0 by 50 fs, so the system has returned to the ground state before reaching the conical intersection). Upper panel (Fig. 5a): spin density on atoms B and C (see Fig. 1). Lower panel (Fig. 5b): Dihedral angle ABCD approximately 50 fs. Now, in the lower panel of Fig. 4, we show the spin density on the two terminal methylenes. As expected, for an isomerization of a radical cation, taking place initially on a adiabatic surface, we see that the odd electron is initially delocalized equally on both of the terminal methylene centers. However as we pass through the CoIn, we see that motion through the region of strong interaction at the CoIn, stimulates electron dynamics. Thus we have a situation where the charge distribution between the two centers is instantaneously asymmetric (i.e. electron dynamics where the electron and nuclear motion are asynchronous). Obviously, such a situation could not be described with a classical trajectory or classical reaction path involving a surface hop. However, as we have shown in other work [36] this dynamics may be washed out when the width of the initial wavepacket is taken into account. In Fig. 5 we show the trajectory where we start with a superposition of D 0 and D 1 (where the odd electron is localized on one of the 2 terminal methylenes). Such a state might be created in attosecond laser experiment [25]. In the top panel of Fig. 5 (showing the spin density on the terminal methylenes), we see rapid electron dynamics, which eventually dies out longer-term. The corresponding nuclear geometry is shown in the bottom panel of Fig. 5. In this example, the Ehrenfest state has collapsed to D 0 by 50 fs, so the system has returned to the ground state before reaching the conical intersection. We can see rotation towards the CoIn at 90°, and then the system relaxes on the ground state D 0 back to the 45°twisted geometry.
These two simulations on allene nicely illustrate the range of applicability of the study of reaction paths using the damped Ehrenfest method. In the adiabatic reaction path we can observe stimulated electron dynamics as one passes through the CoIn region. In this instance, electron and nuclear motion are asynchronous, accompanied by smooth torsional motion. In the simulation, where we start from a superposition of D 0 and D 1 , the reaction path is slightly different, since the transition to the ground state takes place before isomerization is complete.
i) Benzene S 1 S 0 adiabatic and non-adiabatic reaction paths Now let us consider a second example, the prefulvene CoIn in the Channel 3 benzene-benzvalene photochemistry (see for example [37]), which is much more complex. We are interested in the non-adiabatic reaction paths from the "channel 3" region, which we shall call prebenzvalene (PBV in Fig. 6). (A review of the experimental and theoretical results has recently been published [38]). A schematic view of the potential surface for S 1 and S 0 is shown in Fig. 6 (which we have adapted from the early work of Palmer [37]).
The low-lying excited states of benzene are combinations of VB resonance structures as indicated in Fig. 6. The S 1 /S 0 conical intersection has recently been reviewed using VB concepts [1]. Thus it is convenient to describe the electronic structure along possible reaction paths using such VB ideas. Further, it is convenient to localize the active orbitals of the CAS-CI space. In this case the configuration state functions (CSF) of the CI expansion correspond to VB configurations built from localized atomic orbitals. Thus, referring to Fig. 6, S1 in the Franck-Condon region is the antisymmetric combination of the two standard Kekulé structures while the benzvalene region is a combination of Kekulé and Dewar structures as we shall elaborate subsequently.
We will consider two types of simulation: adiabatic simulations initiated near the PBV-TS region (Fig. 6) and non-adiabatic simulations starting at the CoIn P-BV-CoIn. For the adiabatic simulations, we will begin the simulations in the forward direction towards P-BV-CoIn and in the backward direction towards Benzene-S1. In this simulation, via the CoIn P-BV-CoIn, we will have a change of potential surface (at the CoIn) to the region of BV-1, BV-2 and BV-M on S 0 . For the non-adiabatic simulations beginning at the P-BV-CoIn we will examine the effect of a change of phase in the initial superposition of states (i.e. we change the mixing of the two degenerate states that form the initial wavefunction  Table 1 Benzene S 1 S 0 adiabatic and non-adiabatic reaction paths. Initial geometry (Fig. 6) Initial wavefunction (Fig. 8) Short-time a reaction path (Fig. 6) Long time b reaction path (Fig. 6) PBV-TS → P-BV-CoIn C 2 eC 6 = 1. for the Ehrenfest dynamics; for a full discussion of this approach see Meisner et al. [32]). The results are summarized in Table 1 and we now explore some of the implications. All of our computations were run with a 6 orbital, 6 electron CAS space in a 6-31g * basis. The active orbitals were localized onto C atom sites as discussed below. In Table 1 we distinguish (column 1) the two types of simulation we have carried out: "adiabatic" where the initial conditions corresponded to an electronic eigenfunction and non-adiabatic where the computation was initiated from a state that was a superposition. In the case where we begin at a conical intersection, this division is not meaningful because any linear combination is also an eigenstate. In this case we are investigating the effect of changing the phase. In column 2 we give the chemical nature of the initial wavefunction in terms of the VB configurations of Fig. 9. Finally in columns 3 and 4 we summarise the nature of the electronic structure at short-time and long-time. These definitions are intended merely to indicate the qualitative nature of the reaction path (after a few fs.) at the beginning of the reaction path the initial characteristics of the reaction path are established and later when the nature of the products has been identified.

Benzene S 1 S 0 adiabatic reaction paths from the prefulvene-like S 1 transition state
We begin with the example that is closest to a traditional trajectory. We start at the PBV-TS where the geometry has been displaced slightly along the normal co-ordinate corresponding to the imaginary frequency in the direction of P-BV-CoIn (this geometry is indicated PBV-TS → P-BV-CoIn in Table 1). Motion along this co-ordinate in the positive direction PBV-TS → P-BV-CoIn corresponds to compression of the C 2 eC 6 and C 1 eC 3 co-ordinates. The main results are summarized in Fig. 7. In the upper panel we have plotted the energy difference S 1 and S 0 . This passes through zero quite quickly after the initiation of the trajectory.
In the lower panel we show the energy. Notice that the kinetic energy is approximately zero (because we have a damped trajectory). Also notice that the potential energy goes down monotonically. In fact this trajectory ultimately ends up in the region of benzene itself and does not pass near the P-BV-M. In contrast, the trajectory run in the reverse direction from the PBV-TS (PBV-TS → Benzene S 1 Table 1) goes smoothly to the region of benzene-S 1 . In order to understand this result, we must look at the effect of phase of the CoIn and we now proceed to discuss these simulations.

Benzene non-adiabatic reaction paths from the S 1 /S 0 CoIn
In this section we discuss the Ehrenfest trajectories starting from the P-BV-CoIn changing the phase of the mixing (see Fig. 8 and Fig. 9) of the two degenerate eigenstates. We have discussed this effect earlier [32] and noted how this effect may be used to control reactivity. The results of these simulations are summarized in the last four rows of Table 1 and Figs. 8-11. In order to interpret the results we have carried out the simulations using a localized basis of carbon p 2 orbitals, localized on each of the carbon centers. (These orbitals are merely the localized orbitals of the 6 orbital CAS space). This localized active space gives us a valence bond [39] like representation of the wavefunction and enables a physical interpretation of the angle of mixing between the two eigenvectors. However, the interpretation of the results is complicated by the fact that the many electron structures of traditional Valence Bond theory shown in row a of Fig. 9 are not orthogonal and do not correspond to the basis functions usually used in electronic structure computations. The orthogonal basis functions used in electronic structure computations are adapted to the standard canonical order based upon the symmetric group. These structures, shown as non-canonical Rumer structures, are collected in row b of Fig. 9 and are obviously linear combinations of the structures shown in row a. The valence bond structure of benzvalene is shown in row c (and we have discussed the complete CoIn seam using this methodology elsewhere [1,40,41]). In this localized AO basis, at the transition state, and the CoIn, the wavefunction is dominated by configurations U and X shown in the row d of Fig. 9. The prebenzvalene VB structures are shown in Fig. 9c. It must be emphasized that these structures are intended for qualitative interpretation only. Nevertheless is should be clear that U-X contains the pre-benzvalene component A + C plus a Dewar benzene component D. Thus U-X is the phase that is consistent with one of the P-BV structure. Now let us turn to Fig. 8 where we illustrate the effect, in a schematic way, of the mixing of the two eigenstates (which are U + X and U-X) of the CoIn. The black-white pattern indicates the relative weights of U and X. (If one proceeds through the apex of the cone (from the upper sheet to the lower one) then the configuration does not change. On the other hand if one remains on the same sheet then the phase changes from a + b to a-b [42,1,4]). The S 1 transition vector at P-BV-TS is the sum of U + X and corresponds to = 0 in Fig. 8. In contrast, for  Fig. 9c and 9d) and Dewar benzene structures ( Fig. 9c and 9d). Thus U-X = 90 leads towards the prebenzvalene region of the potential surface (Fig. 6) In Fig. 10, we show the results for = 90 , here we have plotted the C 2 -C 6 bond length (which is the main change associated with motion towards BV). One can see that it passes through a minimum around A 1.90 (the CoIn occurs at A 1.94 ) and the corresponding molecular structure is shown in Fig. 11.
Thus in order to reach the BV region of the surface one needs to enter the CoIn region with the correct phase = 90 . We've not discussed the simulations for the other phase angles (summarized in Table 1) because they all ultimately end up in ground state benzene. The schematic version of the potential surface gives one some indication as to why this occurs. The BV1, BV2 and BVM region of the surface is a flat floppy biradicaloid region that is very localized in the sense that there are only small geometric changes involved. Thus looking at Fig. 6 one can see that the region labeled by BV-1, BV-2 and BV-M is an almost continuous set of structures with an unpaired electron on C 1 loosely  Bond length (Angstrom) Time (fs) Fig. 10. C 2 eC 6 bond length for a trajectory from P-BV-CoIn with an initial state mixing of 90 coupled to C 5 . Only the 90°reaction path passes close enough to this region, and only then very briefly before falling off this ridge and returning to benzene itself. Of course this is completely consistent with the vanishingly small quantum yield for this reaction observed experimentally [43]

Conclusions
In this paper we have shown how the damped Ehrenfest method can be used to define and study reaction paths that are non-adiabatic. We illustrate two different situations. In one of these, the initial wavefunction conditions correspond to an adiabatic state and one observes non-adiabatic behaviour in the region where the adiabatic surfaces cross at a CoIn. The other situation is characterized by choosing the initial conditions to be a superposition of two eigenstates with a variety of phase relationships. We have chosen two simple chemical systems, allene and the prefulvene region of benzene to illustrate the main ideas. The most interesting points are the stimulation of electron dynamics while passing through the CoIn in allene, and the phase requirement to get the interesting chemical reaction path in the benzene channel 3 problem. Fig. 11. Geometry after 14 fs of a trajectory from P-BV-CoIn with an initial state mixing of 90 . At this geometry the C 2 eC 6 bond length is at its shortest (as seen in Fig. 10).