Dynamical Friction Effects on the Photoisomerization of a Model

ARTICLE pubs.acs.org/JPCA

Dynamical Friction Effects on the Photoisomerization of a Model Protonated Schiff Base in Solution Jo~ao Pedro Malhado,† Riccardo Spezia,‡ and James T. Hynes*,†,§ cole Normale Superieure, 24 rue Lhomond, 75005 Paris, France Departement de Chimie, CNRS UMR 8640 PASTEUR, E Laboratoire Analyse et Modelisation pour la Biologie et l’Environement, CNRS UMR 8587, Universite d’Evry Val d’Essonne, Bd. F. Mitterrand, 91025 Evry, France § Department of Chemistry and Biochemistry, University of Colorado, Boulder, Colorado 80309-0215, United States † ‡

bS Supporting Information ABSTRACT: Photoisomerization involving a conical intersection (CI) for a model protonated Schiff base (PSB) in modeled water and acetonitrile solvents is examined with the inclusion of energy- and momentum-transfer effects described via a generalized Langevin equation (GLE) frictional approach and surface-hopping dynamics. Short-time GLE frictional effects on the model’s three coordinates, the intramolecular bond length alternation and torsional PSB coordinates and a solvent coordinate, eliminate several unphysical features associated with a no-friction inertial description and have the general feature of accelerating nonadiabatic transitions to the ground electronic state. The inertial prediction of equal probability formation of groundstate trans and cis isomer products subsequent to the Franck-Condon excitation of the ground cis isomer is replaced by the GLE prediction of a preferential higher proportion of ground-state trans isomer, that is, a successful cis to trans photoreaction. This preference is solvent-dependent and is enhanced in water solvent with its higher friction intensity and short time scales. For the fast water solvent motion, the nonadiabatic transitions to the S0 ground state are centered around the CI seam (which is due to the solvent coordinate’s role as a tuning coordinate), facilitating direct transitions to the ground-state trans isomer. In contrast, for the slower acetonitrile solvent motion, the decay occurs, on average, away from the CI seam in regions with a finite free-energy gap between the excited and ground states, resulting in reduced trans isomer production. Some directions for the extension of the model description are also discussed.

1. INTRODUCTION The excited-state isomerization of protonated Schiff bases (PSBs), which involves a conical intersection (CI) facilitating the ultrafast photoreaction, has been the subject of a number of studies.1-13 The excited-state reaction in question is an isomerization, and PSBs are of special interest as they are often regarded as models for retinal, whose photoisomerization in rhodopsin is key in the vision of vertebrates.14-16 In the present work, we extend a recent series of theoretical studies of a simplified model, whose PSB coordinate and charge distribution characteristics were initially inspired1 by vacuum electronic structure calculations of Olivucci, Robb, and co-workers,2,17,18 on to this photoreaction in a polar solvent.1,3,19 In particular, we here include the important short-time dynamical friction effects of momentum and energy transfer for a model PSB photoisomerization in acetonitrile and water solvents. Surface-hopping trajectories are employed to examine the solvent effects on the excited-state decay rate and product isomer distribution, as well as details of the approach of the system to the CI region and the location and character of nonadiabatic transitions in this region. The PSB isomerization is strongly coupled to the surrounding polar solvent due to the marked charge transfer involved in the photoprocess. A series of gas-phase electronic structure r 2010 American Chemical Society

studies2,17,18,20-22 on short PSBs has indicated that there is a characteristic charge-transfer pattern in the photoisomerization; a net positive charge is transferred from one side of the PSB to the other side (across the isomerizing double bond) in the initial Franck-Condon (FC) transition to the excited state, and this charge transfer is reversed upon passage through the CI region to produce the ground-state product. Both the existence of a CI in this type of isomerizing systems and the charge character of the electronic states involved can be rationalized2 in terms of a two-electron, two-orbital model,12,23 where the electronic structure is described by valence bond (VB)-type states localized on the p orbitals involved in the isomerizing double bond. (A more detailed commentary on how the charge distribution is actually treated in the model is given in section 2.) This type of basis set with a charge character independent of the solute geometry (namely, it is conserved in the CI vicinity) proved to be very useful in the prior PSB model Special Issue: Graham R. Fleming Festschrift Received: July 1, 2010 Revised: September 14, 2010 Published: October 08, 2010 3720

dx.doi.org/10.1021/jp106096m | J. Phys. Chem. A 2011, 115, 3720–3735

The Journal of Physical Chemistry A description of the evolving electrostatic interaction of the solute with a polar/polarizable solvent,1,3 and its use in the electronic description for the PSB model is continued here. The model for the reacting PSB chromophore is a minimal one and also involves two molecular coordinates, a bond length alternation (BLA) coordinate (r), reflecting coupled double bond compression and single bond stretching, facilitating and accompanying the isomerization, and an angular coordinate (θ) characterizing the twisting, torsional isomerization motion of single double bond. This simplification is supported by the short PSB gas-phase electronic structure studies2,17,18,20,24 mentioned above, which have indicated that the first excited singlet state minimum energy path, leading to a CI and inducing photoisomerization, involves primarily, though not necessarily exclusively, these two molecular coordinates. The single torsional coordinate selected in the model requires further precision. Isomerization may occur at one of several double bonds of the conjugated carbon chain in a one bond flip mechanism4,5 (although the possibility of a hula-twist mechanism has been suggested25,26) in a process which depends on the charge distribution in the vicinity of the molecule.6,27 In solution, the predominant product results from the isomerization around the central double bond of a conjugated chain with an odd number of double bonds,6,28 and this will be the only case considered here. In the prior model studies,1,3 the above-described VB and twocoordinate representation was supplemented with a third solvent coordinate (z) characterizing the interaction of the PSB charge distribution with the nonequilibrium orientational polarization of the solvent, described at a dielectric continuum level. A nondissipative Hamitonian description for the PSB-solvent system was constructed using an r-, θ-, and z-dependent free energy, which functions as an effective potential energy. This description was employed to characterize, for example, the minimum energy path to a z-dependent seam of CIs. However, the dynamics per se of the model exhibited19 a number of nonphysical aspects, briefly described within, due to the absence of dissipative momentum and energy transfer, so that no meaningful study of the excited-state decay, reaction rate, or product distribution could be affected. This difficulty is here removed by the inclusion of dissipative effects via a generalized Langevin approach involving time-dependent frictions on all of the coordinates, modeled to be appropriate on the ultrafast time scale of the photoreaction. The outline of the remainder of this paper is as follows. In section 2, we collect aspects of the model characterization not involving frictional effects. In particular, we give only a brief summary of the essential features of the VB and nondissipative Hamiltonian description, whose details are already available in refs 3 and 19. Section 3 presents a detailed discussion of the dissipative short-time frictional description and its implementation in surface-hopping trajectory calculations. Section 4 is devoted to the model trajectory results for the photoisomerization and their interpretation. Concluding remarks are offered in section 5.

2. PSB MODEL DESCRIPTION 2.1. Electronic Valence Bond States. Following refs 1, 3, and 19, the PSB electronic structure is defined by two VB states {|B2æ,|ABæ}, schematically shown in Figure 1. In the representation {|B2æ,|ABæ}, the isolated model system electronic Hamiltonian defines the potential energy surfaces in

ARTICLE

Figure 1. (left) The localized states with the formal charges on each site are displayed for an example of a small PSB (the tZt-penta-2,4dieniminium cation). (right) The total charge on each side of the isomerizing double bond, which is, in general, extended over the left- or right-hand sides of the π-conjugated system for the |B2æ and |ABæ states, respectively. See text.

the nuclear coordinates, which has the nondiagonal matrix form ! VB2 ðr, θÞ γðθÞ ð1Þ He ¼ V ¼ γðθÞ VAB ðr, θÞ The analytical form and parametrization of the terms in eq 1 are detailed in the Appendix. In its original form, the two-electron, two-orbital model considers only the two p orbitals forming the isomerizing double bond.12,23 Although the correct qualitative picture of the important factors governing the electronic structure of isomerizing PSBs is obtained,2,17 this model gives rise to restricted charge characteristics for the states |ABæ and |B2æ represented on the left-hand side (LHS) of Figure 1. In our interpretation, they are considered single resonance structures of the two π systems extending on both sides of the molecule. In this way, a more accurate perception of the charge character of state |ABæ is with the positive charge extended over the nitrogencontaining moiety of the molecule, while for state |B2æ, the positive charge resides on the opposite side of the molecule with respect to the isomerizing double bond,17,23 as depicted on the right-hand site (RHS) of Figure 1. (A discussion of the model’s neglect of a third VB state can be found in ref 3.) Figure 2 shows the vacuum potential energy surfaces, together with their charge character, for the eigenstates S1 and S0 formed by diagonalization of eq 1. Finally, from a nuclear dynamics perspective, VB basis sets, such as {|B2æ,|ABæ} used here, are approximately diabatic.29,30 2.2. Solvent Electrostatic Effects: Solvent Coordinate. As outlined in the Introduction, there is a strong electrostatic solvent-PSB solute interaction. The formalism for this interaction when the solvent is polar and polarizable is now briefly sketched; further details can be found in ref 1. The solvent is approximated as a nonequilibrium dielectric polarizable continuum, characterized by a polarization field PB, with two components. The solvent electronic polarization PBel is assumed to always be equilibrated to the PSB solute charge distribution, while the nuclear orientational polarization PBor will generally be out of equilibrium with the rapidly varying PSB charge distribution. It is assumed (a) that the PSB electronic configuration in solution is still described by a linear combination of the VB states |B2æ and |ABæ introduced in section 2.1, which possess a fixed charge character independent of PBor, and (b) that PBor can be expressed as 3721

dx.doi.org/10.1021/jp106096m |J. Phys. Chem. A 2011, 115, 3720–3735

The Journal of Physical Chemistry A

ARTICLE

Figure 2. Vacuum PSB model ground (S0) and excited (S1) potential energy surfaces in the r and θ coordinates. Colors indicate the charge character according to the square coefficients of the projections of the ground and excited states on the VB basis states.134

a linear combination of polarization fields equilibrated to |ABæ and |B2æ, which constitute the extreme cases of the model’s electronic structure charge distributions eq

eq

PBor ¼ ð1 - zÞ PBor;AB þ z PBor, B2

ð2Þ

Here, z is a nondimensional solvent coordinate gauging the solvent polarization. When z = 0, PBor is in equilibrium with the |ABæ charge distribution, and for z = 1, it is equilibrated to that of state |B2æ. It can be shown1,31-33 that the solution free-energy analogue of the gas-phase potential energy matrix (eq 1) is 1 0 1 2 2 ðr, θÞ þ k ð1 zÞ γðθÞ V s B C B 2 C Gðr, θ, zÞ ¼ B @ 1 2A γðθÞ VAB ðr, θÞ þ ks z 2 ð3Þ which we now briefly discuss. First, the solvent “force constant” ks, which measures the ease of solvent coordinate departures away from the appropriate equilibrium, is Z eq eq ε ð PBor; AB ðr, θÞ - PBor, B2 ðr, θÞÞ2 d xB ks ¼ 4π ε¥ ðε - ε¥ Þ ð4Þ where ε¥ and ε are the solvent’s high-frequency and static dielectric constants, respectively, and the integration variable B x denotes the coordinates of three-dimensional space.34 Second, the diagonal terms VB2 and VAB in the free-energy matrix are the vacuum VB state potentials; in fact, each of these should be augmented by the equilibrium solvation free energy of the respective state. In the model, these are taken to be equal, as explained in ref 35, and result only in an unimportant shift in the free-energy zero and are thus suppressed. The ground (GS0) and excited (GS1) adiabatic free-energy surfaces correspond to the eigenvalues of the matrix (eq 3) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G 2 þ GAB 1 GS0, 1 ¼ B ð5Þ ðGB2 - GAB Þ2 þ 4γ2 2 2 which provide effective potential energy surfaces in the PSB and solvent coordinates. Note that comparison of eqs 3 and 5 shows

that the solvent coordinate z is linearly related to the solvent’s contribution to the VB state’s free-energy gap GB2 - GAB, which in turn influences the free-energy gap between the excited- and ground-state surfaces. A detailed analysis of the PSB model free-energy surfaces can be found in ref 1 and 3. A principal feature resulting from including the “tuning” solvent coordinate is that there is a seam of CIs, an aspect to which we return in section 4. This CI seam is located at θ = π/2; in order to reach it from the FC region, solvent coordinate motion is required, in addition to the solute’s internal coordinate displacements. Finally, as in refs 1, 3, and 19, the ks value used here, ks = 2.5 eV, is that obtained for acetonitrile using a two-sphere cavity model to treat an electron transfer in an aromatic radical anion.36 This value is employed to describe the nonequilibrium interaction of PSB molecules and the two solvents examined, water and acetonitrile. This is justified by noting that in eq 4, the difference in equilibrium polarizations depends mostly on the solute charge distribution and is thus similar for both solvents, as is the prefactor, since water and acetonitrile have a comparable ε¥ value which is significantly smaller than the respective static dielectric constant (about a factor of 20 for acetonitrile and a factor of 44 for water37). In this way, given that the solutesolvent electrostatic interaction is identical for the two solvents, their difference examined in the following sections is exclusively due to their different dynamical behavior, especially to the time scale of their molecular reorientation. 2.3. Inertial Dynamics. In order to consider dynamics of the PSB model system on each of the free-energy surfaces GS0,1, defined by eq 5, the kinetic energy K is required. Following ref 19, this is taken to be a sum of three independent terms involving the several conjugate momenta pu of each of the three system coordinates and the associated masses mu. In the inertial description of the dynamics of the system, the motion on the groundand excited-state surfaces is governed by classical mechanics determined by the Hamiltonian H ¼

pr 2 pθ 2 pz 2 þ þ þ GS0, 1 ðr, θ, zÞ 2mr 2mθ 2mz

ð6Þ

The Hamiltonians for the two surfaces are conservative, that is, there are no dissipative or frictional effects. We now discuss the masses in the kinetic energy K in eq 6. Focusing first on the internal coordinates r and θ, in general, the associated masses depend on the molecular geometry. As in previous studies,3,19 such kinetic coupling between internal degrees of freedom38-40 is not considered, and constant masses are taken. For the coordinate r, as in refs 3 and 19, the value mr = 6.0 g 3 mol-1 was used, corresponding to the reduced mass of a CC bond.41 This value yields, for the isolated solute potential surfaces defined in section 2.1, the harmonic frequencies 1100 and 1230 cm-1 for the equilibrium positions on the ground state and on the excited state for a planar configuration, respectively. Both values are lower than the predictions from electronic structure calculations in small PSBs18 (and experimental values for retinal in the protein environment42-44), which indicate frequencies higher than 1500 cm-1 for both instances (for the tZt-penta-2,4-dieniminium cation, a slightly smaller excited-state frequency value of 1447 cm-1 has been reported18). These frequency value differences should not introduce qualitative differences in the results since, as will be seen in section 4, the 3722



ARTICLE

dynamics in r is already significantly faster than the motion of the other model coordinates. For the mass (more precisely, moment of inertia) associated with the angular torsional coordinate, the value used is mθ = 3.8 g 3 Å2 3 mol-1. This yields a harmonic frequency in the ground-state equilibrium position of 543 cm-1, which compares well with the experimental value of 568 cm-1 for retinal within the rhodopsin protein environment.45 The excited-state FC frequency of 384 cm-1 is higher than that calculated for tZtpenta-2,4-dieniminium,18 which is just above 200 cm-1. Finally, a useful approximation for the solvent coordinate mass mz is required. (Note that the solvent coordinate z is adimensional (cf. eq 2), so that from eq 6, mz has the unusual units of mass 3 length2.) While considerations can be made connecting mz with inertial properties of solvent molecules, the moment of inertia along the dipole moment, in particular,31 here, its value will be taken from observed solvent dynamics in time-dependent fluorescence experiments. The spectral Stokes shift following excitation of a dissolved dye molecule can be related to the normalized time correlation function (tcf) of the solvent coordinate Æz(0) 3 z(t)æ/Æz(0)2æ.46,47 The short-time expansion around t = 0 coincides with that of a Gaussian, so that Æzð0Þ 3 zðtÞæ ωz 2 2 2 2 t e-ωz t =2 1 2 2 Æzð0Þ æ

ð7Þ

and for short times, the Gaussian behavior indeed dominates the solvent dynamics.47 If the free energy for the solvent coordinate motion is harmonic, the coefficient ωz in eq 7 is the frequency of the well in which the solvent coordinate moves. This frequency ωz and the solvent force constant ks (eq 4) provide a route to the solvent mass mz ks mz ¼ 2 ð8Þ ωz Two different mz values are required, using approximate solvation dynamics of acetonitrile and water. In the first case, a Gaussian fit to the first 120 fs of the simulated Stokes shift of the dye molecule coumarin 153 in acetonitrile48 yields ωz = 1.03 10-2 fs-1. For water, the short-time Gaussian fit to the experimental Stokes shift decay of coumarin 343 in water yields ωz = 3.85 10-2 fs-1.49 These values, together with the value of ks = 2.5 eV used in the definition of the model free-energy surfaces eq 3, define via eq 8 the solvent masses 227 and 16.3 g 3 Å2 3 mol-1, for acetonitrile and water respectively. These two constant solvent masses are assumed to apply over both ground- and excited-state surfaces.50

3. INCLUSION OF DISSIPATIVE EFFECTS FOR THE PSB AND THE SOLVENT In addition to the inclusion of electronic transitions associated with the PSB photoisomerization whose treatment will be described in section 3.3, two distinct treatments of the dynamics will be contrasted in the trajectory calculations to be presented in section 4. The first uses the inertial dynamics associated with the conservative Hamiltonian eq 6, the approach followed in ref 19. The second more realistic dynamical description, and our main focus, approximately accounts for energy and momentum flow between the PSB coordinates and the solvent coordinate and their respective “surroundings”. There are three components for our treatment involving generalized dissipative effects. The first two involve the dissipative influence of the solvent on the PSB

internal coordinates, that is, energy and momentum transfer between those coordinates and the solvent not accounted for by the coupling of these coordinates to the collective solvent coordinate z in eq 6. The third involves the dissipative influence of the remaining solvent degrees of freedom not included in the solvent coordinate z. In section 3.1, a generalized Langevin treatment accounting for dissipation is described. The model specification of the friction kernels entering this treatment is given in section 3.2, and the description of the surface hopping approach for the electronic transitions for the dissipative trajectories is given in section 3.3. 3.1. Generalized Langevin Equation Description. As noted above, excess energy deposited in the PSB via photoexcitation will flow out of the three degrees of freedom (r,θ,z) treated explicitly to the remaining ones forming an environment. In order to model such generalized dissipative or frictional effects arising from (a) the solvent on the PSB model internal coordinates r and θ and (b) the remaining solvent degrees of freedom (dof) on the solvent coordinate z, the generalized Langevin equations (GLEs) will be employed Z t DGS0, 1 d2 u du - ζu ðt - τÞ dτ þ Fr u ∈ fr, θ, zg mu 2 ¼ dt dτ Du 0 ð9Þ The first term here corresponds to a systematic force, that is, the negative gradient of the free-energy surfaces given by eq 5. So far, these are the equations of motion in the inertial description. The second term induces a dissipative contribution, in which ζu is a memory function or friction kernel ζu ðtÞ ¼

1 ÆFr ð0Þ 3 Fr ðtÞæ kB T

ð10Þ

where kB is Boltzmann’s constant and T is the absolute temperature.51 Finally, Fr is a time-dependent stochastic force representing fluctuations of the force imposed by the environment dof on the coordinate (not accounted for by the first two terms in eq 9) and has a zero average value. When the time scales of these forces are much smaller than the time scale of the velocity correlations of the coordinates of interest, the tcf eq 10 for Fr can be approximated by a delta function, and eq 9 reduces to the well-known Langevin equation (LE). The latter has been used in the past to interpret the dynamics of photoisomerization taking place through a CI.52,53 However, the PSB problem dynamics, involving isomerization and decay through a CI, occur on an ultrafast time scale on the order of some tens of fs, compared to which the correlation time of the friction kernels can certainly not be considered small, for any of the three coordinates. The LE description is particularly unrealistic for vibrational coordinates such as r and θ. In particular, vibrational energy transfer rates are overestimated by many orders of magnitude.54,55 A LE description for the solvent coordinate is a less drastic approximation than that for the PSB internal coordinates but, nonetheless, significantly distorts the dynamics of the solvent.46,54,56 Consequently, the GLE eq 9 must be used for each of the model coordinates. In the GLE (eq 9), the free-energy gradient force term couples all three coordinates since the free energy involves all of these coordinates; this coupling is already present in an inertial description. In general, there will also be a frictional coupling between these coordinates; instead of a single generalized friction 3723



ARTICLE

term in eq 9, a sum of integrals involving friction functions ζuv and the velocity in coordinate v ∈ {r,θ,z} will occur.57,58 This type of coupling will not be taken into account here. However, if the average initial velocities are zero, it is possible58 to recast such sums of frictional terms into a form equivalent to the GLE (eq 9). The GLE eq 9 is a phenomenological equation of motion. GLEs have been extensively used to study the dynamics of chemical systems59 and have been particularly successful for chemical reaction dynamics in solution.60-64 For example, reaction rate theory using a GLE agrees with MD simulation results within the error bars of the latter62 for the reaction transmission coefficients accounting for deviations from Transition-State Theory; by contrast, use of the LE can be in error by several orders of magnitude. Unfortunately, such a level of agreement is not to be expected for the PSB study since a GLE treatment will be used over a significant spatial extent of the coordinates, where the assumption of a single time-dependent friction (for each coordinate) is suspect.65,66 Nonetheless, this description will provide a first estimate of what will prove to be quite important dissipative effects for the PSB model problem. 3.2. Friction Kernels Description and Parametrization. We now turn to the form and parametrization of the GLE friction kernels eq 9 for the PSB model problem. While a simple Gaussian approximation for the time dependence is often useful when only short-time behavior is important,58 it cannot be employed for the PSB model intramolecular dynamics. This is because the BLA stretching and torsional coordinates of the model are bound, and vibrational energy relaxation is involved in their dynamics; MD simulations indicate that a Gaussian friction kernel can lead to an error of more than 10 orders of magnitude for vibrational energy relaxation rates.54 Thus, an alternative functional form is required for the model PSB molecular internal degrees of freedom. For this purpose, we employ the friction kernel functional form initially proposed by Douglass (as a velocity tcf)67,68 ζðtÞ ¼ ζ0

cosðRδtÞ coshðRtÞ

ð11Þ

It is more complex than a Gaussian and has three parameters (ζ0 setting the initial value; R controlling the decay time scale; and δ defining the shape/oscillatory behavior). The key point is that this friction kernel form was later used to study vibrational relaxation of diatomics with physically reasonable results.55,69 We now turn to its parametrization. The function eq 11 oscillates with the period τ = (2π/Rδ) and reaches zero at t = (π/2Rδ). Oscillations are small for small values of δ and become more pronounced for δ values on the order of unity or higher. A small value of δ = 0.4, for which the oscillatory character of eq 11 is strongly muted (Figure 3), gives a good qualitative account of the tcf of the forces of the environment on a vibrational coordinate55 and is adopted throughout this study. This yields a friction kernel which decays on a longer time scale than a Gaussian, as is appropriate, with the same curvature at time zero. With the parameter δ now selected, it remains to assign the time scale parameter R and the initial friction amplitude parameter ζ0 for the time-dependent friction eq 11 for the PSB model’s two internal degrees of freedom. To determine the first of these parameters, eq 11 is expanded up to second order in time and then matched to the corresponding expansion of ζ(t)/ζ0 for a MD-simulation-inferred result for the approximate translational

Figure 3. Normalized friction kernel functions. For the PSB internal coordinates (r, θ), the functional form eq 11 is used. For the two solvent cases studied, two different Gaussian functions are considered.

Table 1. Parameters Concerning Dynamics of the PSB Internal Coordinatesa u r θ a

mu 6.0 g 3 mol-1

3.8 g 3 Å2 3 mol-1

R/10-2 fs

(ζ0)/(mu)/10-3 fs-2

δ

3.7

1.8

0.4

3.7

3.2

0.4

The friction functional form is eq 11.

time-dependent friction for an ion in water, here, the chloride ion58 (for details on this approximate procedure, see also ref 70). This yields R = 3.7 10-2 fs-1, which is used for both r and θ coordinates. We pause to note that, because the frictional modeling is intended to be appropriate for the ultrashort time scales relevant for the CI problem, no long time (ps or longer) “tail” is included in the modeling; such tails contribute heavily to large friction constants, giving, for example, small rotational and translational diffusion constants in solution (see for example, refs 71 and 72), but have no relevance on the short time scales of interest here.73 We now deal with the parametrization for the initial friction kernel amplitude ζ0, which is proportional to the mean square value of external forces acting on the PSB model internal coordinates (see eq 10). For the BLA coordinate r, the same amplitude as that obtained via the above-mentioned calculation for chloride translation in water was used. For the θ coordinate, the ion translation value is expected to be a less good approximation since the “force” here is a torque and more sensitive to the molecular geometry. For the torsion coordinate, a value of ζ0 four times larger than that for r was adopted, reflecting the fact that, even in the case of small PSBs, opposition of the solvent to the end group motion will be more pronounced than that for the BLA motion, an aspect that will only be enhanced for larger PSBs. The friction parameter values for the internal coordinates are compiled in Table 1. Here, we add some comments concerning the simplifications involved in our choice of generalized friction description. In general, the time-dependent friction, in particular, its amplitude ζ0, should depend on the solute’s charge distribution.74 We have neglected this feature, consistent with the assumption in section 2.1 that the positive charge character, although not its location, in the present PSB model is the same for each VB state component of the ground and the excited states. Removal of that assumption would introduce discontinuities in the frictional time dependence due to the existence of hops between electronic states (to be discussed below), which implies an instantaneous change in the PSB charge distribution. The incorporation of such effects would at the very least require 3724



ARTICLE

the friction amplitude to be state-, r-, and θ-dependent. It is not known at present how to do this in any reasonably physically correct way, and this is not attempted here; the same friction kernels for the PSB internal coordinates (and the solvent coordinate; see below) are used for the ground and excited states. Turning to the friction kernel for the solvent coordinate z, a simple Gaussian form is adopted rffiffiffi 2 ζS -t2 =2σz 2 -t 2 =2σ z 2 ¼ e ð12Þ ζz ðtÞ ¼ ζ0 e π σz where ζS is the static friction constant and is the time integral from 0 to ¥ of the friction kernel Z ¥ ζz ðtÞ dt ð13Þ ζS ¼ 0

Parametrization of this kernel thus involves defining two parameters, the initial friction amplitude ζ0 (or, via eq 13, the friction constant ζS) and the time scale σz. To this end, the kernel is connected to solvation dynamics as determined in timedependent fluorescence Stokes shift experiments75 and MD simulations.56 Here, the basic perspective of refs 46, 47, and 56 is followed, in which such solvation dynamics is governed by the tcf of the solvent coordinate z. For this special purpose, we take the free-energy potential acting on the solvent coordinate in eq 9 to be harmonic, so that the free-energy derivative there is given by -

DG ¼ -mz ωz 2 z ¼ -ks z Dz

ð14Þ

involving the solvent mass, the square frequency, and the solvent force constant related by eq 8. The minimum of the free-energy curve has been taken at z = 0; this involves no loss of generality, only simplifying the notation to avoid introducing subtracted averages. Standard manipulations47,76 of eq 9 to form the solvent coordinate tcf can then be used to show that the full time integral of the normalized solvent coordinate tcf is Z τz ¼ 0

¥

Æzð0Þ 3 zðtÞæ ζS ζ dt ¼ ¼ S 2 2 mz ω z ks Æzð0Þ æ

ð15Þ

which defines the effective time τz. Since this relation is independent of the precise form of the friction kernel, it applies to the Gaussian eq 12. Thus, since the solvent force constant has been defined in section 2.2, when τz and σz are determined, the Gaussian friction is completely specified. As stated above and in section 2.3, within the applicability of linear response theory, it is possible to extract information about the solvent coordinate tcf from Stokes shift spectral data. For acetonitrile (with coumarin 153 dye), this tcf decay yields τz = 260 fs.75 This value and the solvent force constant ks used in the definition of the free-energy surfaces (section 2.2) together define through eq 15 the initial amplitude ζ0 of the solvent coordinate friction kernel for a solvent mass corresponding to acetonitrile. Information about the time scale of the memory function acting on the solvent coordinate, σz in (eq 12), was obtained from MD simulations of model dipolar solutes in model chloromethane,56 an aprotic polar solvent with a mass similar to that of acetonitrile. This yields the value σz = 45 fs used in the

Table 2. Parameters for Solvent Coordinate z Dynamics and Friction for the Two Solvent Mass mz Values Used mz/g 3 Å2 3 mol-1 CH3CN H2O

227 16.3

(ζ0)/(mz)/10-3 fs-2

τz/fs

45

0.49

260

16

4.1

55

σz/fs

simulations with the acetonitrile mz value. This then completes the friction parametrization for this solvent. For the model water solvent parametrization of the friction, an approximate route is taken. The above-mentioned MD study56 shows that both parameters τz and σz vary by changing the mass of the solvent molecules. Since mz is a measure of the inertia of the solvent in changing its polarization, in particular, by rotation, it should be directly related to the molecular mass. The fractional variation of τz and σz with fractional variation of molecular mass can thus be used to scale the values of these quantities, such that a fractional variation of 93% of mz going from acetonitrile to water will result in a variation of 79% on τz and 66% on σz, to yield τz = 55.0 fs and σz = 16 fs. The former value is lower than that of most current solvents75 since in a real system, the integral of the friction includes the effect of the long time tail mentioned above; it is nevertheless expected to give a good account of the shorttime solvent dynamics. The friction kernel time scales used for both solvents can be compared in Figure 3. Parameters employed for the friction on the solvent coordinate z are shown in Table 2. The initial solvent dynamics predicted by the GLE for acetonitrile agrees well with that calculated from MD simulations,48 while in the case of water, a slightly faster decay than that measured experimentally49 is obtained. 3.3. Frictional Force and Surface Hopping. The model PSB photoisomerization process involves a CI between the S1 and S0 electronic states. In order to account for nonadiabatic transitions between these states, the dynamical treatment adopted here, as in ref 19, follows a surface-hopping perspective,77-80 by which the system is described by an ensemble of independent trajectories which propagate classically on the adiabatic free-energy surfaces that result in the diagonalization of matrix eq 3, with electronic transitions (hops between free-energy surfaces) occurring along the trajectory according to a stochastic process. Our simulations employed Tully’s fewest switches scheme77,78,81 for the calculation of electronic transitions, also called “molecular dynamics with quantum transitions”. Comparison of the results of simulations using the fewest switches algorithm with full quantum dynamics simulations may range from good to fair,82-84 but this surface-hopping approach is found to give at least good qualitative agreement in general, and the fewest switches algorithm is usually the reference adopted for comparison with other classical nuclear nonadiabatic approaches. Surface-hopping simulations were carried out in which the dynamics on each surface was determined by the GLE eq 9 instead of the Hamiltonian equations derived from eq 6. Most details of these simulations are postponed until section 4, but here, we address special issues arising from the use of a GLE description. In particular, attention is required for the generation of the stochastic forces Fr for an arbitrary form of the friction kernel ζ(t) since their tcf must fulfill the condition eq 10. For this purpose, an algorithm treating the stochastic forces as an autoregressive time series85,86 was used (see Supporting Information for implementation details). For the nonadiabatic transitions per se, the introduction of friction presents no difficulties since in the fewest switches surface-hopping algorithm,77,78,81 although the 3725



ARTICLE

evolution of the quantum electronic system depends on the nuclear motion, there is no assumption on how their dynamics is propagated, either inertial or dissipative; thus, the hopping probability at each time step is calculated in the same way in both cases. Surface-hopping simulations have been reported for other problems where the classical degrees of freedom are evolved in time through either a GLE-related formalism87 (for interelectronic state transitions) or LE88 (for an internal conversion problem involving a CI). From the point of view of the GLE dynamics, the occurrence of hops entails a discontinuity of the velocities due to energy conservation criteria,77,78,81 which, in turn, will affect the dynamics after the hop has occurred, due to the convolution integral in eq 9 whose derivative will no longer be continuous.89,90 These effects will be small when hops occur in the vicinity of the CI seam where the S1-S0 free-energy gap is small and will have no effect on the dynamics taking place from the FC region and the intersection before any hop takes place.

4. PSB MODEL ISOMERIZATION RESULTS In this section, we present the surface-hopping trajectory results for the model PSB isomerization process for the two model solvents considered, water and acetonitrile. The focus is on the dynamics including the GLE frictional description, that is, with the equations of motion in eq 9 on the S1 and the S0 surfaces between nonadiabatic surface hops, but for perspective, we also include the inertial dynamics results, where the excited- and ground-state surface dynamics between hops is governed by the Hamiltonian eq 6. Before the presentation of these trajectory results, we provide a few further details on the calculations. Whenever electronic state populations are reported, these correspond to the fraction of trajectories that propagate on the respective S1 or S0 surface at a given moment in time. Simulations were performed using 10000 independent trajectories starting on the excited state with Gaussian distributions in both position and momentum for each of the model’s coordinates, the internal PSB BLA (r) and torsional (θ) coordinates, and the solvent coordinate (z). These distributions model a FC transition and are thus centered at the ground-state equilibrium positions and zero momenta, with widths consistent with the equipartition theorem for a thermalized system at 300 K with a harmonic approximation to the freeenergy surface around the equilibrium position. Trajectories propagating on the ground-state surface with a kinetic energy less than the free-energy gap between surfaces, a frustrated hop situation,78 were kept in the original surface with no change of velocity direction, as suggested in ref 84. The alternative prescription of inverting the velocity direction91 has been shown to give negligibly different results.19 Numerical integration of eq 9 was done by employing a modified velocity Verlet algorithm described in ref 92 (without the approximations proposed there to speed up the calculations). The computation procedure of stochastic forces associated with the GLE friction kernels has been discussed in section 3.3 and is detailed in the Supporting Information. 4.1. Population Density Evolution. Figure 4 presents the time evolution of the PSB model system state population density for the two model solvents considered, when the dynamics is inertial. Figure 5 represents the same two solvent cases in a GLE frictional description. This subsection gives an overview of the dynamics and contrasts the results obtained, neglecting and including dissipative effects, clearly indicating the importance of the latter.

Figure 4. Time evolution of the population density of the PSB model with water (A) and acetonitrile (B) solvents in the excited and ground states. No friction is applied. Bright color corresponds to higher population density.

There are two striking features observed in the inertial dynamics in Figure 4. First, the ground-state population distributions in the torsional coordinate (and, to a lesser extent, also in the excited state) evolve for continuously higher absolute values of θ. This reveals that, once arriving at the ground state, molecules possess enough rotational kinetic energy to surmount the double bond torsional barrier. The PSB molecule is thus performing free bond rotations in solution at a rate of 10 revolutions/ps, which is clearly unphysical. Associated with this behavior, and the fact that the ground-state torsional barrier region corresponds to a crossing or close crossing of the S0 and S1 surfaces, repetitive double bond rotations induce upward transitions from the ground to the excited state. Such repopulation of the excited state can be observed on the LHS of Figure 4 by the appearance of population density in the excited state at |θ| = 3π/2 and 5π/2 that does not result from the original distribution evolution in the excited state. A second striking feature of Figure 4 is that in both the excited and ground states, solvent coordinate oscillations are observed (this is more apparent in the H2O case in Figure 4A; for CH3CN, only half of an oscillation in z is observed on the simulation time scale). From eq 2, this corresponds to oscillations of the average polarization field surrounding the PSB molecule. From the point of view of time-dependent fluorescence Stokes shift experiments, 3726



Figure 5. Time evolution of the population density of the PSB model with water (A) and acetonitrile (B) solvents in the excited and ground states. Motion on each coordinate is governed by a generalized Langevin equation with friction kernels defined in section 3. Bright color corresponds to higher population density.

this would correspond to a time oscillation of the frequency of the observed fluorescence peak. Such a behavior has apparently never been observed; thus, such solvent coordinate oscillations are also unphysical. Both of these inertial simulation characteristics are related to the neglect of dissipative effects. In the case of the free bond rotations, energy conservation of the Hamiltonian eq 6 dictates that the excess energy due to photoexcitation be conserved within the model’s three degrees of freedom. However, as anticipated in section 3, in a realistic description, such energy would be transferred to the remaining degrees of freedom of the PSB molecule and to the solvent molecules in particular. In a similar fashion, while motion of the solvent molecules is dominated by inertial effects at short times, for longer times, dissipation arising from interactions between solvent molecules will dominate. Indeed, as shown in Figure 5, with the frictional effects of inclusion through the GLE formalism, both unphysical characteristics of the model inertial dynamics are suppressed. The frictional influence on the kinetic energies of the three model coordinates are discussed in the Supporting Information. We next turn to a more detailed analysis of the more significant results including friction in Figure 5. Starting with the BLA

ARTICLE

coordinate r dynamics, the population distribution remains compact on both the excited and ground states, with a motional time scale shorter than those of the other system coordinates and with a high-frequency oscillatory pattern, which is largely independent of the solvent identity. No solute-solvent vibration to vibration (VV) energy transfer is possible with the PSB internal coordinate frictions employed; such VV transfer would not however likely be effective on the time scales considered here (see, for example, ref 93). We next focus on the population distributions along the solvent coordinate z. As discussed in section 2.1, there is a charge transfer within the PSB molecule upon photoexcitation. In the FC region, the solvent is out of equilibrium, and the GLE-based Figure 5 shows a monotonic, rather than oscillatory, solvent coordinate motion toward the excited-state equilibrium solvation position, with the population distributions remaining monomodal. A transition to the ground state corresponds to a reverse charge transfer in the PSB solute, inducing an appropriate solvent adjustment toward the S0 charge distribution. For both solvents, the main characteristics of the z motion are similar and occur on clearly distinct time scales, faster for water (Figure 5A) than for acetonitrile (Figure 5B). While the inclusion of friction has an important effect on the solvent coordinate dynamics in suppressing unphysical oscillations, as discussed above, the difference in solvent motional time scale is mostly due to the solvent’s inertial properties (section 2.3) since it is already present in the inertial simulations (Figure 4). As will be seen in the following, the different solvent motion time scales, and, in particular, how they compare with the torsion time scale, have an important impact on the photoreaction outcome and mechanism. We now turn to torsional coordinate θ in the GLE results in Figure 5. The free-energy surfaces are periodic in θ,1,3 with the positions θ = 0 and π corresponding to stable configurations on the ground state (cis and trans isomers, respectively) and, except for the FC region, to maximum free-energy points on the excitedstate surfaces. S0 torsional coordinate evolution through a broadening of the initial Gaussian distribution is centered at θ = 0 (Figure 5). This is a relatively slow process, essentially due to two factors; the first factor is that in the FC region, the potential along θ is bound3,18 and only after an initial motion in r does it become unstable; and second, independent of the value of r, the gradient in θ at the population distribution’s center (θ = 0) is always zero, and thus, forces acting on the initial distribution are small. It is most important to note that in Figure 5, the appearance of population density on the ground state coincides with the excited-state population reaching the twisted geometry |θ| = π/2, which occurs about 50 fs after photoexcitation. In refs 1, 3, and 19, it was shown that the solvent coordinate motion is required in order to reach the CI seam from the FC region. While in the case of water (Figure 5A) the motions in θ and z occur on a similar time scale, for acetonitrile (Figure 5B), motion in θ is faster, such that it reaches |θ| = π/2 before the solvent coordinate brings the system on to the CI seam. Thus, transitions to the ground state occur for this case before the CI seam is reached. This point proves to be important and is discussed more thoroughly in the next subsection. A last point of interest in Figure 5 is that changes in the solvent induce changes in the population densities along the torsional coordinate. On one hand, for water solvent (Figure 5A), no excited-state population accumulation is observed at |θ| = π/2, suggesting that transitions to the ground state occur immediately upon reaching the twisted geometry. On the other hand, for the 3727



ARTICLE

Figure 6. Excited-state minimum free-energy paths for water (blue) and acetonitrile (red) solvents. Bright colors correspond to the paths in the three-dimensional space (r,θ,z), while darker colors correspond to projections of these paths on the planes r = 1.65 Å, θ = π/2, and z = 1.2. The black line indicates the CI seam on plane θ = π/2.

acetonitrile case (Figure 5B), a fraction of the trajectories remains in the excited state at |θ| = π/2 for a fraction of time, as is apparent from an oscillatory pattern in the S1 population density centered at this value of θ. This fact has significant consequences, as discussed in the following subsections. 4.2. Approach to CI Seam and Nonadiabatic Transition Locations. As was just emphasized, after the FC S0 to S1 transition, solvent motion on the excited-state surface, in addition to motion in the PSB internal coordinates, is required to reach the CI seam. Here, we examine aspects of the pathway of the CI seam and then turn to the location of the actual dynamic nonadiabatic transitions to the S0 ground state. A simple description of this pathway is a minimum (free) energy path (MEP). A MEP is of course not a real dynamical path. For the present system, since the masses associated with the different coordinates are independent of those coordinates, the MEP corresponds to the steepest descent path, starting from the FC region, in mass-weighted coordinates.31 The MEP can be viewed as a dynamical path with no inertia, with infinitely small velocity or, in the language of section 3, a strongly overdamped trajectory. Such paths, which are very familiar for gas-phase reactions,94,95 have been employed in several detailed electronic structure calculation studies of isolated molecule CIs2,17,18,20,21 and for ground electronic state CI model problems.36,58 Here, we will use MEPs to provide a comparison reference for the full GLE trajectory results for the dynamical effect of different solvents with different solvent masses. A previous study has detailed PSB model MEPs for different, more extreme solvent mass values.19 While we will focus on the GLE trajectory results, it is useful to begin with the MEP in all three of the model’s coordinates, shown in Figure 6. The first important things to note are that the paths differ, in how the CI seam is reached and where it is reached, for the two solvents. Both features reflect the different solvent masses and thus different time scales. After an initial, similar motion in the BLA coordinate for the two solvents, the paths diverge. In the fast solvent H2O, solvent coordinate evolution proceeds first to a location above the seam, and only afterward does the path, largely in the torsion coordinate but also including some further BLA motion, reach the CI seam. In the

Figure 7. Average excited-state trajectories on the (r,z) plane for both solvents (solid lines). Excited-state minimum free-energy paths (dashed lines). The seam is indicated by the solid black line, with X indicating the minimum free-energy point along the seam. The black dot corresponds to the FC position (note that it is located at θ = 0).

slower solvent CH3CN, the solvent coordinate motion is on a time scale slower than that of the torsion and BLA coordinates; the system first reaches the |θ| = π/2 plane and then proceeds to the CI seam by motion in the z coordinate. When projected on the (r,z) plane, the path to the CI seam (after an initial r motion portion very similar for the two solvents), the solvent coordinate motion precedes a final r motion for the fast H2O solvent but is the final motion for CH3CN. The dominance of the more rapid motions in the passage from the FC region to the CI seam over the slower motion is opposite to that occurring for a barriercrossing reaction.58 The MEPS and the average GLE trajectory results for the BLA and solvent coordinate for the water and acetonitrile solvents are shown in Figure 7. The excited-state distributions for these two coordinates shown in Figure 5 are fairly well localized, allowing a fair comparison with the MEPs. In the water case (Figure 7), the MEP reflects this rapid solvent’s feature that the solvent coordinate motion occurs before the final approach to the CI seam via the torsional coordinate and the accompanying BLA motion; in contrast, for the slower acetonitrile solvent (Figure 7), the final approach is in the solvent coordinate. (Note that the MEP is restricted to the CI seam once the seam is reached.) The excited-state GLE average trajectory results bear some similarity to these paths but exhibit 3728



Figure 8. Location of nonadiabatic transitions from the excited to ground state. Contours correspond to energy gaps between the excited and ground states for θ = π/2. The CI seam is indicated by the solid bold line, with X indicating the minimum free-energy point along the seam. The black dot corresponds to the FC position (located at θ = 0).

important differences. For water solvent, the most marked difference is the lightly damped character of the BLA coordinate motion, which continues (with less motion in the solvent coordinate) in the seam region. For acetonitrile solvent, the most marked differences from the MEP are first, the BLA coordinate oscillation accompanying the slow solvent motion and, second, the feature that the CI seam is not actually reached in the 290 fs time frame displayed. For the PSB isomerization per se, Figure 8 indicates, together with the free-energy gap between S1 and S0 surfaces, the positions in the r and z coordinates where the electronic transition from the excited to the ground state occurs for the first time for each trajectory. The most striking feature, compared to standard concepts, is the rather wide dispersion of the locations for the nonadiabatic transitions. Not only is there a significant transition probability away from the CI seam, there is also significant probability of transitions occurring away from the minimum free-energy point on the CI seam. The slower acetonitrile solvent case is particularly marked in these respects; Figure 8 shows that the probability of nonadiabatic transitions to the ground state is distorted away from the CI seam toward solvent coordinate values less than those appropriate for the seam, a clear reflection of this solvent’s slower time scale dynamics than that for water.

ARTICLE

Figure 9. Distribution of where transitions from the excited to the ground state occur for the first time for each trajectory. (uhop - uCI), where u is one of the coordinates r, θ, or z, corresponds to the distance to the CI seam on that coordinate, keeping other coordinates fixed at their uhop value.

These aspects are confirmed by the distributions displayed in Figure 9, where the transition position distributions are plotted with respect to the distance to the values (rCI,θCI,zCI) which correspond to the values of the respective coordinates on the CI seam, keeping the other coordinates at the values where the transition occurs. The width of the distributions indicates that transitions occur in the vicinity of the plane (r,|θ| = π/2,z) but not necessarily in the immediate vicinity of the seam. For the slower acetonitrile solvent, the transition position distributions in r and z are shifted with respect to the seam values, while such a shift does not occur for θ. This indicates that transitions occur in the vicinity of the plane (r,|θ| = π/2,z) but not necessarily in the immediate vicinity of the seam. Indeed, a reasonable description of the process in this solvent would be that transitions would take place each time the system reaches the twisted geometry vicinity while the solvent is moving comparatively slowly toward the seam value. 4.3. Excited-State Population Evolution. We next turn to the time decay of the excited S1 state evolution for the isomerizing PSB in the two solvents. Figure 10 for water and acetonitrile 3729



ARTICLE

Table 3. Parameters Resulting from a Fit of Equation 16 to the Excited-State Decay Profiles for H2O and CH3CN Solvents and with and without Friction along with the GroundState Population Ratio of the Different Ground-State PSB Product Isomeric Forms at Asymptotic Times t0/fs

k/fs-1

(PS0,trans/PS0,cis)(t = 290 fs)

H2O

no friction

39

0.004

1.0

CH3CN

with friction no friction

54 22

0.0151 0.004

3.0 1.0

with friction

71

0.0150

1.7

Figure 10. Excited-state population time evolution for the PSB model for water and acetonitrile solvents. Dashed lines correspond to simulations without friction.

solvents shows the S1 evolution in the GLE description for the trajectories, as well as the no-friction, purely inertial dynamics trajectory results. The latter suffer from the unphysical inertial dynamics characteristics described in detail in section 4.1, but they are nonetheless useful for perspective. In each case, there is an induction period, to which we return below, followed by an S1 population decay. The first striking point about the decay is that for each solvent, the inclusion of the generalized friction has accelerated the decay, a feature already apparent, albeit less clearly, in Figure 5 compared to Figure 4. Ordinarily, one would associate friction with a slowing in time; here, however, the acceleration is primarily a consequence of the suppression of the repeated rotations about the PSB model isomerizing bond found in the pure inertial simulations, which were causing repopulation of the excited state, mentioned in section 4.1, and which for the no-friction case in Figure 10 are responsible for the leveling (and occasional increase) of the excited-state population.96 In order to further discuss the results in Figure 10, excited-state population decay profiles were fit by ( 1 for t < t0 PS1 ðtÞ ¼ ð16Þ e-kðt - t0 Þ for t > t0 where t0 is taken to represent an induction time which can be seen as the time taken by the system to get from the FC region to the CI seam vicinity and k is a rate constant describing the speed of the decay. Equation 16 does not provide an exact description of the excited-state decays but, with some qualification to which we will return, gives a good qualitative account of the trends registered in Table 3. This applies even for the case of no friction, where the decay is nonexponential and the fit (eq 16) is very approximate. (Experimentally, dynamics is often monitored by loss of population from the FC region, which would not have a delay time; see the discussion below.) We consider first the induction time t0. As was just discussed, solvent coordinate motion is required to the reach the CI seam region from the FC location on the S1 surface, whatever the solvent’s mass. However, the time scales for this motion differ for water and acetonitrile solvents. For the more massive acetonitrile, its slowness compared to the PSB torsional motion leads to a lengthening of the induction time compared to that of the lighter water solvent. In contrast, and remarkably, the S1 population decay rate constant k in Table 3, although about three times larger than the no-friction rate constant, is the same for both solvents.97 This apparent equality of decay rate constants is somewhat misleading

Figure 11. Time evolution of the fraction of trajectories in the FC region for water and acetonitrile solvents. Dashed lines correspond to an exponential fit to the envelope of the curves.

since, as noted above, the finite friction exponential time fits used to determine them from Figure 10 is approximate and, in particular, is rather crude for CH3CN solvent, where a bimodal character of the decay (past the induction period) is clearly evident. A more detailed discussion of the decays in Figure 10 is useful here. For H2O solvent, for which the solvent and torsional time scales are comparable, the decay commences when the tail of the |θ| distribution reaches π/2; the exponential decay is governed by the angular distribution increasing in this neighborhood (see Figure 5A). For CH3CN solvent, the solvent coordinate is slow compared to the torsion, as discussed above, and the delay time is no longer governed exclusively by the |θ| distribution tail reaching π/2, and although the first transitions occur off of the CI seam, with a finite excited-ground-state gap, the decay will be slower on average since a fraction of trajectories does not make the transition. Subsequently, as the solvent coordinate approaches more closely to the seam, the gap closes, and the decay is even a bit more rapid than that in H2O since the torsional population at |π/2| is larger (see Figure 5). The differences in Figure 10, while visible, are not pronounced, and they are simply suppressed in the exponential decay fits. We noted above that the excited-state dynamics is typically monitored experimentally via loss of population in the FC region or from excited-state locales further along the path toward the CI region,98 rather than the global S1 population as in our PS1(t) in Figure 10. Figure 11 displays the calculated time-evolving fraction of trajectories located in the FC region. This region is defined by the parallelepiped in (r,θ,z) comprising a 3σ width for each side of the ground-state equilibrium location for each coordinate, with σ as the standard deviation of the ground-state Gaussian coordinate distribution.99 Aside from the high-frequency oscillations in the BLA coordinate, which bring the system out of and into the FC region, there is, on average, an overall (very) approximately exponential decay, with decay rates 3730



Figure 12. Ground-state population evolution of the two PSB isomeric forms for water and acetonitrile solvents.

of 0.020 and 0.027 fs-1 for CH3CN and H2O, respectively, a difference somewhat more pronounced than that for decay portions of the excited-state populations in Figure 10. The slower decay for CH3CN solvent is due to the torsional motion; the solvent coordinate motion is sufficiently slow as not to contribute significantly on the decay time scale. In contrast, the faster decay for H2O solvent is due to the fast solvent coordinate motion out of the FC region; in fact, if the z coordinate is not included in the FC region definition, there is no distinction between the two solvents’ behavior. 4.4. Ground-State Product Isomer Ratio. The excited-state population decay rate constant just discussed is of course moot on whether the PSB isomerization reaction has been “successful” or not, that is, whether the ground-state trans isomer has been produced via the photoreaction rather than the cis isomer. Table 3 shows that there is an important solvent effect here; the desired product ratio is higher, by almost a factor of 2, in water than that in acetonitrile (without friction, the trans/cis ratio is unity, a point to be returned to below). Figure 12 displays the time evolution of the ground-state PSB isomer distributions and can be used to gain insight into the origin of the greater trans product ratio in water solvent. It will be noticed that for each solvent, there is a delay in the rise of the cis curve compared to that of the trans curve. This suggests that ground-state cis isomer formation is related to trajectories in which an initial “direct” nonadiabatic transition from the excitedstate to ground-state trans has not occurred. Together with Figures 8 and 9 indicating the distributions of locations of the transitions to the ground state, this feature implies the following picture, schematically illustrated in the two panels of Figure 13 and now discussed. In the faster water solvent, transitions are concentrated about |θ| = π/2 with the solvent coordinate (and the BLA coordinate) concentrated on the CI seam. Thus, trajectories will have the dominant character not far from that shown in the top panel of Figure 13, where successful transitions to the trans branch of the

ARTICLE

Figure 13. Schematic for the general mechanism leading to the formation of the ground-state cis and trans isomers. In the top panel, on the time scale of the motion on the torsional coordinate θ, the excited and ground states are degenerate at θ = 0, leading to the formation of mainly the trans isomer. This is the case for the fast water solvent. In the bottom panel, motion in θ occurs with a finite energy gap between the excited and ground states, leading to formation of a significant amount of cis isomer. This result is observed for the slower acetonitrile solvent.

ground state occur on the first pass. In reality, the free-energy curves will often be slightly split at |θ| = π/2, but successful transitions to ground-state trans will remain in the majority. Nonetheless, there is finite probability that the solvent and BLA coordinates are off of the CI seam, such that there is a gap between the excited and ground states at and near |θ| = π/2, as depicted in the lower panel of Figure 13. Trajectories which do not make a transition to ground-state trans on the first pass will perform an oscillation in the excited-state well and, for example, decay on the second pass through the twisted geometry.100 For the slower acetonitrile solvent, at and in the neighborhood of |θ| = π/2, the situation differs from the water solvent case in that the solvent coordinate (and the BLA coordinate) are most likely off of the CI seam, and the more probable situation is like that depicted in the lower panel of Figure 13, with an S1-S0 energy gap and ground-state cis generation via excited-state torsional oscillation for those trajectories not making a transition to ground-state trans on the first pass. Again, trajectories remaining in S1 will tend to be stabilized in the S1 well by frictional dissipation (see Figure 5B) with comparable proportions of trans and cis ground-state isomer generation. The favoring of the desired ground-state trans product to the cis product in faster water (ratio = 3; see Table 3) compared to slower acetonitrile solvent (ratio = 1.7) is a marked solvent effect. Also shown in Table 3, the trans/cis product ratio for either solvent in an inertial description without friction is found to be very close to 1. This inertial description lack of discrimination is another consequence of the free double bond rotations described in section 4.1, with the population distribution spreading over several ground-state free-energy wells, yielding equal populations for both isomeric forms, a feature removed by the focusing effect of the friction. 3731



5. CONCLUDING REMARKS In this article, we have examined the energy- and momentumtransfer effects described via a generalized Langevin equation (GLE) frictional approach for a model protonated Schiff base (PSB) photoisomerization involving a conical intersection (CI) in modeled water and acetonitrile solvents. The inclusion of short-time GLE frictional effects on the three coordinates of the model, the intramolecular bond length alternation, torsional PSB coordinates, and a solvent coordinate, eliminated a number of unphysical features associated with a no-friction inertial description. The most striking GLE friction influence is the elimination of the inertial description of equal probability formation of ground-state trans and cis isomer products subsequent to the Franck-Condon excitation of the ground cis isomer. Somewhat similar effects are obtained in refs 101 and 102 in a reduced dimension photoisomerization model, although there are very considerable differences between the model and methods employed there and in the present work.103 Further, the importance of dissipative energy transfer to the surroundings has very recently been suggested for the photoisomerization of all-trans retinal.104 Among the most significant solvent-dependent effects found are the following. For the fast water solvent motion, the nonadiabatic transitions to the S0 ground state are centered around the CI seam (the seam itself having its origin in the solvent coordinate’s role as a tuning coordinate). In contrast, for the slower acetonitrile solvent motion, the decay occurs, on average, away from the CI seam in regions of space with a finite free-energy gap between the excited and ground states. These differences in the position of the decay were shown to translate to a solvent dependence of the distribution of photoproducts, with a preferential higher proportion of cis isomer in the latter case. The characteristics of the generalized friction, especially that acting on the torsional coordinate, also influence the distribution of photoproducts; higher friction intensity and short time scales favor the formation of ground-state trans isomer, that is, a successful cis to trans photoreaction.105 The numerical magnitude of differing solvent effects on the dynamics of the photoisomerization process found within via comparison of model water and acetonitrile solvents cannot be said to be pronounced. For example, the excited-state probability induction times differ (H2O compared to CH3CN) by a factor of 0.76, and the final decay rates are essentially the same, while the ground-state trans/cis isomer ratios differ by a factor of 1.8; the Franck-Condon decay “rates” differ by a factor of 1.4. In part, this could be due to the fact that both H2O and CH3CN are fairly polar solvents and, in a dielectric continuum description, do not differ much in their solvating power. Further, there are no “steric” short-range repulsive effects in the model free-energy description. While they are implicitly present in the modeled friction for the torsional motion, the present model does not consider the solvent dependence of that friction. In addition, there are no barriers along the pathway in the present model description, while the presence of an electronic barrier has been suggested for some longer PSBs in some polar solvents106 (see also ref 21). Such a barrier might account for the longer time scales and smaller trans product yields (compared to the present model) observed experimentally for some longer PSBs in polar solvents;28,104,106,107 however, as indicated below, the vacuum potential that we have employed may not be appropriate for such longer PSBs.

ARTICLE

However, in at least one aspect, the present model would tend to overemphasize the solute’s coupling to the solvent, namely, the modeled PSB-solvent interactions involve solvent interactions with fairly localized charges, which is certainly a simplification compared to real PSBs, especially longer ones. Recent experiments98,108 on rhodopsin photoisomerization show little solvent dependence for dynamical time scales whose character is most closely related to the Franck-Condon times calculated here. On the other hand, there are significant experimental differences between retinal photoisomerization in solution and in protein environments,28,106-110 aspects of which have been connected to the influence of charged groups in the local protein environment.6,27,111,112 However, any comparison of retinal and the present model, significant aspects of which are derived from the vacuum calculations of Olivucci, Robb, and co-workers,2,17,18,20,24 may not be appropriate; very recent vacuum calculations113 on various model retinal chromophores suggest pathways and potential surface characteristics, especially for larger PSBs, which differ rather significantly from those in that earlier work. Future model studies will need to account for this. There are a number of routes in which the work presented here could and should be extended. Here, some possibilities of future directions are briefly discussed. For convenience, these are addressed separately, although, clearly, their combination is desirable. One future direction is to improve the dynamical description by employing a quantum dynamical computational scheme.6,114-117 A possible choice of algorithm would be the multiconfiguration timedependent Hartree (MCTDH) wavepacket propagation method,117 which has been successfully employed in a number of isolated molecule conical intersection problems,118,119 as well as a version of the present model in an inertial description.120 A second important improvement would be to the model itself (beyond any amelioration of the vacuum potential model, vide supra) by addition of a generalized frictional contribution due to other internal molecular degrees of freedom not treated explicitly and intramolecular vibrational energy redistribution. This inclusion would be expected to reduce the importance of the friction related to the solvent. The IVR impact of solute internal degrees of freedom will assist and, in the limit, dominate the present model’s sole reliance on the solvent to affect the energy transfer necessary to “focus” trajectories into the CI region and to suppress repopulation of the excited state121 (see section 4.1). An illustration of the former effect is provided in ref 122, where IVR to the internal degrees of freedom of isolated butadiene is sufficient to suppress oscillations in a twisting coordinate. Such IVR, however, need not suppress, though it may reduce the probability of events such as those in the bottom panel of Figure 13, responsible for a solvent effect on the product yield and arising from slow solvent motion. Such a reduction would require, for example, the IVR to slow the torsional motion to such a great extent so as to render it slower than the slow solvent motion in order to avoid the presence of an excited-groundstate energy gap. In any event, these issues require detailed investigation, since IVR, particularly at lower energies and on very short time scales, need not be rapid or extensive.123-128 A future and significant model improvement would be to exploit the fact that molecular-level interactions can be recast in terms of free energies, leading to the definition of generalized, microscopic solvent coordinates to replace the present dielectric continuum one.52,53,129-132 A first step to exploit this possibility would be to retain the model description of the model PSB solute 3732



ARTICLE

Table 4. Parameters Used in the Definition of the PSB Model Matrix Elements of Equation 3a DAB

RAB

DB2

5.0 eV cθAB

4.0 eV cBθ2

1.64 Å

RB2

-1

kFC

2.05 Å

-1

rFC

req AB

rBeq2

Δr

λ

1.33 Å

1.48 Å

0.12 Å

0.53 Å

σFC r

σFC θ

ΔE

cγ

-1.09 eV 0.27 eV 2.18 eV 1.39 Å 0.075 Å 0.62 0.27 eV -1.77 eV a

Same as those reported in ref 19 (misprints in the latter reference are corrected in the present table).

but to replace the continuum solvent with a microscopic description. A final and more challenging direction for extension of the present work would be the description of the photoreaction in a protein environment. A primary goal would be the development of a solvent coordinate type of description adapted to the local environment in highly specific protein environments, whose spatially inhomogeneous character strongly contrasts with the more homogeneous solvent environment. A very first step would be to retain the model PSB description employed here but to recast the interaction with the environment and the environment itself via a formulation involving a dynamic, analytic generalization of a static Poisson-Boltzmann approach133 often employed in computational studies.

’ APPENDIX Analytic Forms of Terms in Equation 1. This appendix specifies the excited- and ground-state free-energy surfaces’ dependence on the PSB internal coordinates by specifying the functional forms of the terms in the matrix of eq 1. These terms were taken to have the same form as those described in ref 19, which have the same functional form as the ones in ref 3 (with a different parametrization), which, in turn, are a modification of the ones described in ref 1 to achieve a better description of the FrankCondon region. The diagonal elements of eq 1, the VB state potentials, are defined as morse shift θ ðrÞ þ VAB ðr, θÞ þ VAB ðθÞ VAB ðr, θÞ ¼ VAB θ ðrÞ þ VBshift VB2 ðr, θÞ ¼ VBmorse 2 2 ðr, θÞ þ VB2 ðθÞ þ

VBFC 2 ðr, θÞ þ ΔE

ð17Þ ð18Þ

whose individual components are eq Þ AB B2

-RAB, B2 ðr - r

morse VAB , B2 ðrÞ ¼ DAB, B2 ð1 - e

,

2 shift VAB , B2 ðr, θÞ ¼ -2DAB, B2 ðRAB, B2 Þ eq

-jr - r

Δrðr - rAB, B2 Þe

eq j=λ AB B2

,

ðsinðθÞÞ2

2 θ θ VAB , B2 ðθÞ ¼ cAB, B2 ðcosðθÞÞ FC VBFC ðsinðθÞÞ2 e-ðr - r 2 ðr, θÞ ¼ k

FC 2

Þ2

ð19Þ

ð20Þ

ð21Þ

2 2 FC 2 Þ =ðσFC r Þ -sinðθÞ =ðσ θ Þ

e

ð22Þ The off-diagonal VB state coupling terms of eq 1 have the simple form γðθÞ ¼ cγ cosðθÞ

ð23Þ

The parameter values used are collected in Table 4 and were chosen such as to approximately reproduce the main characteristics of the potential energy surfaces of the isolated tZt-penta2,4-dieniminium cation reported by vacuum electronic structure 17,18 The resulting free-energy surfaces from eq 3 are used studies. as a generic simplified PSB model in solution.

’ ASSOCIATED CONTENT

bS

Supporting Information. Description of the algorithm used to generate correlated stochastic forces in the GLE and of the kinetic energy evolution on the three model coordinates. This material is available free of charge via the Internet at http:// pubs.acs.org.

’ ACKNOWLEDGMENT The financial support of Fundac-~ao para a Ciência e a Tecnologia (FCT) in the form of a Ph.D. Grant to J.P.M. is acknowledged; the present work is based in part on the Ph.D. thesis of J.P. M.134 This work was supported in part (J.T.H.) by NSF Grant CHE-0750477. J.T.H. takes this opportunity to express his appreciation of Graham Fleming's scientific contributions. ’ REFERENCES (1) Burghardt, I.; Cederbaum, L. S.; Hynes, J. T. Faraday Discuss. 2004, 127, 395–411. (2) Gonzalez-Luque, R.; Garavelli, M.; Bernardi, F.; Merchan, M.; Robb, M. A.; Olivucci, M. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 9379– 9384. (3) Burghardt, I.; Hynes, J. T. J. Phys. Chem. A 2006, 110, 11411– 11423. (4) Szymczak, J. J.; Barbatti, M.; Lischka, H. J. Chem. Theory Comput. 2008, 4, 1189–1199. (5) Ruiz, D. S.; Cembran, A.; Garavelli, M.; Olivucci, M.; Fuss, W. Photochem. Photobiol. 2002, 76, 622–633. (6) Ben-Nun, M.; Molnar, F.; Schulten, K.; Martínez, T. J. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 1769–1773. (7) Weiss, R. M.; Warshel, A. J. Am. Chem. Soc. 1979, 101, 6131– 6133. (8) Bernardi, F.; Olivucci, M.; Robb, M. A. Chem. Soc. Rev. 1996, 25, 321–328. (9) Molnar, F.; Ben-Nun, M.; Martínez, T. J.; Schulten, K. J. Mol. Struct. 2000, 506, 169–178. (10) Levine, B. G.; Martínez, T. J. Annu. Rev. Phys. Chem. 2007, 58, 613–634. (11) Weingart, O.; Schapiro, I.; Buss, V. J. Phys. Chem. B 2006, 111, 3782–3788. (12) Bonacic-Koutecky, V.; Koutecky, J.; Michl, J. Angew. Chem., Int. Ed. Engl. 1987, 26, 170–189. (13) Yamazaki, S.; Kato, S. J. Chem. Phys. 2005, 123, 114510. (14) van der Horst, M. A.; Hellingwerf, K. J. Acc. Chem. Res. 2004, 37, 13–20. (15) Dugave, C.; Demange, L. Chem. Rev. 2003, 103, 2475–2532. (16) Kandori, H.; Schichida, Y.; Yoshizawa, T. Biochemistry (Moscow, Russ. Fed.) 2001, 66, 1197–1209. (17) Garavelli, M.; Celani, P.; Bernardi, F.; Robb, M. A.; Olivucci, M. J. Am. Chem. Soc. 1997, 119, 6891–6901. (18) Garavelli, M.; Bernardi, F.; Olivucci, M.; Vreven, T.; Klein, S.; Celani, P.; Robb, M. A. Faraday Discuss. 1998, 110, 51–70. (19) Spezia, R.; Burghardt, I.; Hynes, J. T. Mol. Phys. 2006, 104, 903– 914. (20) Cembran, A.; Gonzalez-Luque, R.; Serrano-Andres, L.; Merchan, M.; Garavelli, M. Theor. Chem. Acc. 2007, 118, 173–183. 3733


The Journal of Physical Chemistry A (21) Cembran, A.; Bernardi, F.; Olivucci, M.; Garavelli, M. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6255–6260. (22) Bonacic-Koutecky , V.; Sch€offel, K.; Michl, J. Theor. Chim. Acta 1987, 72, 459–474. (23) Michl, J.; Bonacic-Koutecky , V. Electronic Aspects of Organic Photochemistry; Wiley: New York, 1990. (24) Cembran, A.; Bernardi, F.; Olivucci, M.; Garavelli, M. J. Am. Chem. Soc. 2003, 125, 12509–12519. (25) Sumita, M.; Saito, K. Chem. Phys. Lett. 2006, 424, 374–378. (26) Liu, R. S. H. Acc. Chem. Res. 2001, 34, 555–562. (27) Cembran, A.; Bernardi, F.; Olivucci, M.; Garavelli, M. J. Am. Chem. Soc. 2004, 126, 16018–16037. (28) Freedman, K. A.; Becker, R. S. J. Am. Chem. Soc. 1986, 108, 1245–1251. (29) Garrett, B. C.; Truhlar, D. G. The Coupling of Electronically Adiabatic States in Atomic and Molecular Collisions. In Theoretical Chemistry: Advances and Perspectives; Henderson, D., Ed.; Academic Press: New York, 1981; Vol. 6A, pp215-289. (30) Sevin, A.; Giessner-Prettre, C.; Hiberty, P. C.; Noizet, E. J. Phys. Chem. 1991, 95, 8580–8584. (31) Lee, S.; Hynes, J. T. J. Chem. Phys. 1988, 88, 6853–6862. (32) Kim, H. J.; Hynes, J. T. J. Chem. Phys. 1992, 96, 5088–5110. (33) Bianco, R.; Hynes, J. T. J. Chem. Phys. 1995, 102, 7864–7884. (34) In general, this force constant is a function of the PSB internal coordinates, an aspect neglected here. See ref 31. (35) The diagonal vacuum VB state potentials in the matrix in eq 3 in general should be respectively supplemented by equilibrium solvation free energies ΔGB2eq ,AB(r,θ). The essential difference between the charge distribution of states |ABæ and |B2æ is, as indicated on the RHS of Figure 1, that the former has the overall positive charge on the RHS of the molecule and the latter on the LHS. Ignoring the details of the PSB molecule’s electronic structure and internal asymmetry, the states |ABæ and |B2æ can be seen to have an equivalent charge distribution, which is symmetric with respect to the isomerizing double bond. Thus, following refs 1, 3, and 19, this equilibrium stabilization free energy is taken to be identical for both states and independent of the internal coordinates. Then, this term merely shifts the free-energy surfaces of the ground and excited states and does not alter the free-energy gap between them, thus not introducing any deformation of the surfaces in the (r,θ) plane. For convenience, its value is simply set to zero ΔGB2eq(r,θ) = ΔGeq AB(r,θ) = 0. (36) Laage, D.; Burghardt, I.; Sommerfeld, T.; Hynes, J. T. J. Phys. Chem. A 2003, 107, 11271–11291. (37) Reichardt, C. Solvents and Solvent Effects in Organic Chemistry, 2nd ed.; VCH: New York, 1988. (38) Ayala, P. Y.; Schlegel, H. B. J. Chem. Phys. 1998, 108, 2314– 2325. (39) Decius, J. C. J. Chem. Phys. 1948, 16, 1025–1034. (40) Wilson, E. B., Jr.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover: Mineola, NY, 1980. (41) The appropriate value39,40 for a BLA mode on a linear chain of five sp2-hybridized atoms, with four carbons and a terminal nitrogen, is 4.2 g 3 mol-1. For a similar chain of nine atoms, this value would be 4.1 g 3 mol-1. These values will however vary when substituents are added to the molecular backbone. (42) Song, L.; El-Sayed, M. A. J. Am. Chem. Soc. 1998, 120, 8889–8890. (43) Smith, S. O.; Braiman, M. S.; Myers, A. B.; Pardoen, J. A.; Courtin, J. M. L.; Winkel, C.; Lugtenburg, J.; Mathies, R. A. J. Am. Chem. Soc. 1987, 109, 3108–3125. (44) Herbst, J.; Heyne, K.; Diller, R. Science 2002, 297, 822–825. (45) Lin, S. W.; Groesbeek, M.; van der Hoef, I.; Verdegem, O.; Lugtenburg, J.; Mathies, R. A. J. Phys. Chem. B 1998, 102, 2787–2806. (46) van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 1985, 89, 4181– 4188. (47) Carter, E. A.; Hynes, J. T. J. Chem. Phys. 1991, 94, 5961–5979. (48) Ingrosso, F.; Ladanyi, B. M.; Mennucci, B.; Elola, M. D.; Tomasi, J. J. Phys. Chem. B 2005, 109, 3553–3564. (49) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Nature 1994, 369, 471.

ARTICLE

(50) A previous study19 conducted similar simulations using an artificially light solvent mass (with a value of about half of that of water used in the present study). (51) Even though the system under study is out of equilibrium, eq 10 is taken to remain valid, and the temperature is taken to be T = 300 K throughout. (52) Warshel, A.; Chu, Z. T. J. Phys. Chem. B 2001, 105, 9857–9871. (53) Warshel, A.; Chu, Z. T.; Hwang, J.-K. Chem. Phys. 1991, 158, 303–314. (54) Bruehl, M.; Hynes, J. T. Chem. Phys. 1993, 175, 205–221. (55) Egorov, S. A.; Skinner, J. L. J. Chem. Phys. 1996, 105, 7047– 7058. (56) Smith, B. B.; Staib, A.; Hynes, J. T. Chem. Phys. 1993, 176, 521– 537. (57) Grote, R. F.; Hynes, J. T. J. Phys. Chem. 1981, 75, 2191–2198. (58) Burghardt, I.; Laage, D.; Hynes, J. T. J. Phys. Chem. A 2003, 107, 11292–11306. (59) Snook, I. The Langevin and Generalized Langevin Approach to the Dynamics of Atomic, Polymeric and Colloidal Systems; Elsevier: New York, 2007. (60) Tucker, S. C. J. Phys. Chem. 1993, 97, 1596–1609. (61) H€anggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 1990, 62, 251–341. (62) Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1989, 90, 3537–3558. (63) Hynes, J. T. The Theory of Reactions in Solution. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. 4, pp171-234. (64) Grote, R. F.; Hynes, J. T. J. Chem. Phys. 1980, 73, 2715–2732. (65) Rey, R.; Hynes, J. T. J. Phys. Chem. 1996, 100, 5611–5615. (66) Straub, J. E.; Borkovec, M.; Berne, B. J. J. Phys. Chem. 1987, 91, 4995–4998. (67) Isbister, D. J.; McQuarrie, D. A. J. Chem. Phys. 1972, 56, 736–738. (68) Douglass, D. C. J. Chem. Phys. 1961, 35, 81–90. (69) Schwarzer, D.; Vikhrenko, D. V.; Vikhrenko, V. S. Chem. Phys. 2004, 300, 197–208. (70) Bursulaya, B. D.; Zichi, D. A.; Kim, H. J. J. Phys. Chem. 1996, 100, 1392–1405. (71) Hynes, J. T. J. Phys. Chem. 1986, 90, 3701–3706. (72) Grote, R. F.; van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 1984, 88, 4676–4680. (73) These could become relevant in situations where free-energy gradients become sufficiently small. (74) Kumar, P. V.; Maroncelli, M. J. Chem. Phys. 2000, 112, 5370– 5381. (75) Horng, M. L.; Gardecki, J. A.; Papazyan, A.; Maroncelli, M. J. Phys. Chem. 1995, 99, 17311–17337. (76) Zichi, D. A.; Ciccotti, G.; Hynes, J. T.; Ferrario, M. J. Phys. Chem. 1989, 93, 6261–6265. (77) Drukker, K. J. Comput. Phys. 1999, 153, 225–272. (78) Tully, J. C. J. Chem. Phys. 1990, 93, 1061–1071. (79) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562–572. (80) Bjerre, A.; Nikitin, E. E. Chem. Phys. Lett. 1967, 1, 179–181. (81) Coker, D. F. Computer Simulation Methods for Nonadiabatic Dynamics in Condensed Systems. In Computer Simulation in Chemical Physics; Allen, M. P., Tildesley, D. J., Eds.; Kluwer Academic Publishers: Norwell, MA, 1993; Vol. 397, pp 315-377. (82) Jasper, A. W.; Truhlar, D. G. J. Chem. Phys. 2005, 122, 044101. (83) Topaler, M. S.; Hack, M. D.; Allison, T. C.; Liu, Y.-P.; Mielke, S. L.; Schwenke, D. W.; Truhlar, D. G. J. Chem. Phys. 1997, 106, 8699– 8709. (84) M€uller, U.; Stock, G. J. Chem. Phys. 1997, 107, 6230–6245. (85) Smith, D. E.; Harris, C. B. J. Chem. Phys. 1990, 92, 1304–1311. (86) Nilsson, L. G.; Padro, J. A. Mol. Phys. 1990, 71, 355–367. (87) Balk, M. W.; Brooks, C. L., III; Adelman, S. A. J. Chem. Phys. 1983, 79, 804–815. (88) Cattaneo, P.; Granucci, G.; Persico, M. J. Phys. Chem. A 1999, 103, 3364–3371. 3734


The Journal of Physical Chemistry A (89) Neither of the previously referenced studies87,88 involving surface hopping with dissipation is affected by the memory effect issues resulting from the hops. In the latter case, this is because the LE has no memory, and in the former case, it is because eq 9 was not used directly, but rather an equivalent treatment involving coupling to extra degrees of freedom was used. (90) From the practical point of view, velocity discontinuities do not prevent the numerical computation of the integral in eq 9. Further, the unphysical aspect of velocity discontinuities is also present in inertial surface-hopping treatments and is unavoidable with the imposition of a energy conservation restriction to a mixed quantum-classical system. (91) Mei, H. S.; Coker, D. F. J. Chem. Phys. 1996, 104, 4755–4767; Misprints: From equations 2.7-2.11 of the article, v(t - Δ) should read v(t - d). (92) Tuckerman, M. E.; Berne, B. J. J. Chem. Phys. 1991, 95, 4389– 4396. (93) Rey, R.; Møller, K. B.; Hynes, J. T. Chem. Rev. 2004, 104, 1015– 1928. (94) Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd ed.; Prentice Hall: New York, 1999. (95) Fukui, K. Acc. Chem. Res. 1981, 14, 363–368. (96) There is a certain partial analogy here with the no- or very low friction lack of stabilization of a high-energy isomerizing species in a product well. (97) We have explored the robustness of this conclusion by considering the case of a time-local LE friction rather than a GLE description, accomplished by increasing the parameter R in the r and θ friction kernels eq 11 by a factor of approximately 25.134 The rate constant then drops slightly but identically for each solvent, from 0.015 fs-1 in Table 3 to 0.013 fs-1. (98) Zgrablic, G.; Voïtchovsky, K.; Kindermann, M.; Haacke, S.; Chergui, M. Biophys. J. 2005, 88, 2779–2788. (99) The FC region is thus defined as {1.29 < r < 1.48 Å, -0.25 < θ < 0.25, -0.13 < z < 0.55 }. (100) For a related effect, see ref 7. (101) Balzer, B.; Hahn, S.; Stock, G. Chem. Phys. Lett. 2003, 379, 351–358. (102) Hahn, S.; Stock, G. J. Chem. Phys. 2002, 116, 1085–1091. (103) In particular, in refs 101 and 102, the coupling and tuning functions of the angular and vibrational coordinates are reversed compared to the present work, no solvent coordinate is introduced, the coordinates’ coupling is to a harmonic “bath” with an Ohmic spectral density, and the dynamics are described via a Markovian Redfield equation. Due to these extensive differences, it is difficult to usefully compare and contrast, for example, the physical origins of the enhanced ground-state trans product found here and in refs 101 and 102. (104) Zgrablic, G.; Ricci, M.; Novello, A. M.; Parmigiani, F. Photochem. Photobiol. 2010, 86, 507–512. (105) A more detailed study of the influence of friction amplitude and time scales may be found in ref 134. (106) Logunov, S. L.; Song, L.; El-Sayed, M. A. J. Phys. Chem. 1996, 100, 18586–18591. (107) Becker, R. S.; Freedman, K. J. Am. Chem. Soc. 1985, 107, 1477– 1485. (108) Zgrablic, G.; Haacke, S.; Chergui, M. J. Phys. Chem. B 2009, 113, 4384–4393. (109) Hamm, P.; Zurek, M.; R€oschinger, T.; Patzelt, H.; Oesterhelt, D.; Zinth, W. Chem. Phys. Lett. 1996, 263, 613–621. (110) Kandori, H.; Sasabe, H. Chem. Phys. Lett. 1993, 216, 126–132. (111) Kennis, J. T. M.; Larsen, D. S.; Ohta, K.; Facciotti, M. T.; Glaeser, R. M.; Fleming, G. R. J. Phys. Chem. B 2002, 106, 6067–6080. (112) Tomasello, G.; Olaso-Gonzalez, G.; Altoe, P.; Stenta, M.; Serrano-Andres, L.; Merchan, M.; Orlandi, G.; Bottoni, A.; Garavelli, M. J. Am. Chem. Soc. 2009, 131, 5172–5186. (113) Valsson, O.; Filippi, C. J. Chem. Theo. Comp. 2010, 6, 1275– 1292. (114) Tannor, D. J. Introduction to Quantum Mechanics: A TimeDependent Perspective; University Science Books: Sausalito, CA, 2007.

ARTICLE

(115) Kosloff, R. Annu. Rev. Phys. Chem. 1994, 45, 145–178. (116) Shalashilin, D. V.; Child, M. S. Chem. Phys. 2004, 304, 103– 120. (117) Beck, M. H.; J€ackle, A.; Worth, G. A.; Meyer, H.-D. Phys. Rep. 2000, 324, 1–105. (118) Raab, A.; Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. J. Chem. Phys. 1999, 110, 936–946. (119) Brill, M. R.; Gatti, F.; Lauvergnat, D.; Meyer, H.-D. Chem. Phys. 2007, 338, 186–199. (120) Burghardt, I.; Spezia, R.; Hynes, J. T. Solvation and Photochemical Funnels: Environmental Effects on Conical Intersections Structure and Dynamics. In Femtochemistry 7; Castleman, A. W., Jr., Kimble, M. L., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; pp 143-153. (121) In addition, less energy will be initially deposited in designated coordinates such as r and θ due to some Franck-Condon energy deposition into the remaining solute internal coordinates. (122) Levine, B. G.; Martínez, T. J. J. Phys. Chem. A 2009, 113, 12815–12824; see Figure 6. (123) Lan, Z.; Domcke, W. Chem. Phys. 2008, 350, 125–138. (124) Yu, X.; Leitner, D. M. Chem. Phys. Lett. 2004, 391, 181–186. (125) Leitner, D. M.; Wolynes, P. G. J. Phys. Chem. A 1997, 101, 541–548. (126) Logan, D. E.; Wolynes, P. G. J. Chem. Pys. 1990, 93, 4994– 5012. (127) Leitner, D. M.; Levine, B.; Quenneville, J.; Martínez, T. J.; Wolynes, P. G. J. Phys. Chem. A 2003, 107, 10706–10716. (128) From a methodological point of view, the introduction of a generalized friction on the r and θ coordinates arising from internal solute degrees of freedom requires attention to the feature that energy125,126,135 or action-angle variables136,137 are a more natural description that the coordinates and momenta in a conventional GLE. (129) Adelman, S. A. Adv. Chem. Phys. 1983, 53, 61–223. (130) Benjamin, I.; Barbara, P. F.; Gertner, B. J.; Hynes, J. T. J. Phys. Chem. 1995, 99, 7557–7567. (131) Ando, K.; Hynes, J. T. J. Phys. Chem. B 1997, 101, 10464– 10478. (132) Keirstead, W. P.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1991, 95, 5256–5267. (133) Honig, B.; Nicholls, A. Science 1995, 268, 1144–1149. (134) Malhado, J. P. Ph.D. thesis, UPMC/ENS, Paris, France, 2009 (135) Sibert, E. L., III; Reinhardt, W. P.; Hynes, J. T. J. Chem. Phys. 1984, 81, 1115–1134. (136) Grote, R. F.; Hynes, J. T. J. Chem. Pys. 1982, 77, 3736–3743. (137) Sibert, E. L., III; Hynes, J. T.; Reinhardt, W. P. J. Chem. Phys. 1984, 81, 1135–1144.

3735


Dynamical Friction Effects on the Photoisomerization of a Model

Recommend Documents