Nonadiabatic Dynamics of Uracil: Population Split among Different

ARTICLE pubs.acs.org/JPCA

Nonadiabatic Dynamics of Uracil: Population Split among Different Decay Mechanisms Dana Nachtigallov
a,*,† Ad
elia J. A. Aquino,‡ Jaroslaw J. Szymczak,‡,^ Mario Barbatti,*,§ Pavel Hobza,† and Hans Lischka*,‡ †

Institute of Organic Chemistry and Biochemistry, Flemingovo n
am. 2, 166 10 Praha 6, Czech Republic Institute of Theoretical Chemistry, University of Vienna, Waehringerstrasse 17, A 1090 Vienna, Austria § Max-Planck-Institut f€ur Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 M€ulheim, Germany ‡

bS Supporting Information ABSTRACT: Nonadiabatic dynamics simulations performed at the state-averaged CASSCF method are reported for uracil. Supporting calculations on stationary points and minima on the crossing seams have been performed at the MR-CISD and CASPT2 levels. The dominant mechanism is characterized by relaxation into the S2 minimum of ππ* character followed by the relaxation to the S1 minimum of nπ* character. This mechanism contributes to the slower relaxation with a decay constant larger than 1.5 ps, in good agreement with the long time constants experimentally observed. A minor fraction of trajectories decay to the ground state with a time constant of about 0.7 ps, which should be compared to the experimentally observed short constant. The major part of trajectories decaying with this time constant follows the ππ* channel and hops to the ground state via an ethylenic conical intersection. A contribution of the relaxation proceeding via a ring-opening conical intersection was also observed. The existence of these two latter channels together with a reduced long time constant is responsible for a significantly shorter lifetime of uracil compared to that of thymine.

1. INTRODUCTION Upon UV excitation all five naturally occurring nucleobases relax into the ground state through internal conversion mechanisms on time scales not exceeding a few picoseconds,1
6 in contrast to some of their analogues not found in nucleic acids.1,7 This ultrafast decay serves as a protection mechanism, which prevents photochemical damage. This remarkable photostability of the nucleobases has been investigated intensively by means of experimental and theoretical approaches. Ultrafast pump
probe techniques1
4,6
8 have been applied in order to obtain detailed information on the specific time scales of the deactivation processes. For the purine bases, adenine and guanine,1,2,4,6,8,9 experimental data indicate a single homogeneous mechanism with a faster decay compared to pyrimidine bases. The experimental studies performed on the pyrimidine bases, thymine, uracil, and cytosine, reveal more complicated photodynamics with generally longer lifetimes; the longest one is for thymine.1,2,4 For uracil, using gas phase time-resolved photoelectron spectra obtained with a pump wavelength of 250 nm, three time constants were determined: a decay of less than 0.05 ps, followed by intermediary (0.53 ps) and longer (2.4 ps) decay constants.4 Using the same initial energy, mass-selected femtosecond-resolved pump
probe ionization measurements were fitted with two exponentials, giving short and long decay times of 0.13 and 1.05 ps, respectively.1 Kang et al. used a monoexponential fit in the analysis of femtosecond-resolved pump
probe resonant ionization spectra with 267 nm pump wavelength. They reported an exponential decay of 2.4 ps.2 All these experimental results agree on a faster decay than that of thymine,1,2,4 for which r 2011 American Chemical Society

lifetimes of 6.4 and 5.1 ps have been obtained. A dark state with a lifetime of several nanoseconds found in time-resolved photoelectron spectra of uracil and thymine was explained by trapping in the nπ* state.10,11 Many theoretical studies have been performed to obtain a molecular picture explaining the observed ultrafast decay mechanisms. Reaction paths and conical intersections between energy surfaces of different electronic states have been computed for the natural nucleobases,12
24 nucleobase models,25
27 isomers,17,28 and substituted species.7,20,23,26 Special emphasis must be attributed to these conical intersections and their energetic accessibility since they are responsible for the ultrafast nonradiative decay occurring in the nucleobases. Concentrating the discussion on uracil, the energetically lowest conical intersection found so far on the S2/S1 crossing seam has a boat conformation.19,29,30 Two main conical intersections for S1 f S0 deactivation were identified for uracil. The first one is an ethylenic conical intersection formed by twisting of the C5
C6 bond, which promotes a crossing between the ππ* and closed shell (cs) states. The second conical intersection, found by Lan et al. using the semiempirical OM2 approach, is formed by out-ofplane distortion of the oxygen atom connected to the C4 atom, causing a crossing between the nπ* and the closed shell states.31 On the basis of these investigations, different mechanisms have been proposed to explain the photophysical deactivation to the Received: February 10, 2011 Revised: April 15, 2011 Published: May 06, 2011 5247

dx.doi.org/10.1021/jp201327w | J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A ground state after excitation to the first ππ* (S2) state. These mechanisms can be summarized as follows: (i) a direct path relaxing along the ππ* state to the ethylenic conical intersection;19,20,32 (ii) an indirect path with initial conversion to the S1 minimum (nπ*) and subsequent motion to one of the conical intersections;19,33 and (iii) an early trapping in the S2 minimum (ππ*),30 crossing to S1 followed by either the direct path or the indirect path described above. Hudock et al. also showed that the initial fast time constant does not correspond to a change in electronic state, but corresponds to an increase in the ionization energy and change in the electronic character of the ionized state.30 Conventional quantum-chemical calculations have provided detailed information about the several reaction mechanisms. However, a final word on the actual deactivation path cannot be given in that way. Photodynamic investigations constitute a logical next step to build a more precise picture of the internal conversion of nucleobases. For that purpose, nonadiabatic dynamics simulations have been performed for all five natural nucleobases using ab initio,30,34
38 semiempirical,31,39
41 and density functional40,42,43 methods. The relaxation mechanisms of all five nucleobases have been overviewed recently in ref 38. Multiple spawning calculations at the complete active space self-consistent-field (CASSCF) level, restricted to the initial phase of the uracil photodynamics, predicted that an S2 trapping should occur similar to thymine.30 In contrast thereto, semiempirical surface hopping dynamics did not show S2 trapping at all. It was predicted that the indirect path with decay at the conical intersection characterized by strong out-of-plane displacement of oxygen atom should dominate the dynamics.31 This work aims at a significantly extended nonadiabatic dynamics simulation of the photophysics of uracil, covering a longer temporal range of deactivation through the different deactivation pathways than in previous investigations. It will be shown that the above-mentioned S2 trapping occurs for a major fraction of trajectories, but that also important contributions of the direct path are observed. The present simulations have additionally revealed a third mechanism not found in previous studies: deactivation via ring-opening conical intersections.

2. COMPUTATIONAL DETAILS On-the-fly ab initio dynamics simulations were performed using Tully’s surface hopping approach44 by solving Newton’s equations for nuclear motion in time steps of 0.5 fs, using the velocity-Verlet algorithm.45 The time-dependent Schroedinger equation (TDSE) was integrated with the fifth-order Butcher algorithm46 with a time step of 0.01 fs using linear interpolation between the time steps of the classical dynamics. The adiabatic representation was used for computing the surface hopping probabilities for nonadiabatic transitions according the fewest switches algorithm.44 Analytic energy gradients necessary for the integration of Newton’s equations and nonadiabatic coupling vectors required for the integration of TDSE were computed by the procedure described in refs 47
51. Decoherence correction with a decay parameter of 0.1 hartree was included according to the prescription of ref 52. Dynamics simulations were carried out at the CASSCF level, with the active space composed of ten electrons in eight orbitals (CASSCF(10,8)). This orbital space contains one n, four π, and three π* orbitals at the ground-state geometry. State averaging was performed over three states (SA-3), and the 6-31G* basis set53 was used.

ARTICLE

The initial conditions were taken from the calculated absorption spectrum that was simulated within a semiclassical approach based on a harmonic oscillator Wigner distribution for the ground state.54 The Wigner distribution composed of 1500 different geometries was projected on two excited states. For the spectrum simulation, a Lorentzian line shape with phenomenological broadening of 0.05 eV was chosen. The absorption spectrum is characterized by a single band with the maximum at 7.00 eV. In the present calculations only the first absorption peak was simulated. The simulation of the second peak would require more states to be included in the state-averaging procedure. To perform dynamics simulation close to experimental pump
probe conditions, the trajectories were selected from that part of the absorption spectra, which corresponds to an excitation (pump) wavelength of 267 nm, as used in time-resolved spectroscopy.1,2 The vapor spectrum of uracil is centered at 244 nm;55 thus the pump energy was located 0.44 eV below the absorption maximum. The same shift was used to restrict the initial conditions for the dynamics. Thus, 59 trajectories were generated with excitation energy of 6.56 ( 0.25 eV and additional 31 trajectories were generated with excitation of 7.06 ( 0.25 eV, i.e., slightly above the computed band maximum in order to study the dependence of the dynamics on the pump energy.56 All trajectories were initiated in the bright S2 state and were propagated for 1.2 ps. For a subset of 47 trajectories, the simulation was continued to 2.5 ps. Ring puckering was characterized by adopting the Cremer
Pople (CP) parameters57 and the Boeyens classification scheme.58 The CP parameters quantify the degree (Q) and the kind (θ and φ) of ring puckering. In particular, Q = 0 Å indicates planar structures and θ and φ parameters can be used to classify the ring conformation as boat (B), envelope (E), chair (C), half-chair (H), screw boat (S), and twisted boat (T). Stationary points and minima on the crossing seam were optimized at the CASSCF level with the same active space as described above and at the multireference configuration interaction level with single and double excitations (MR-CISD). In the CI procedure the orbitals with natural occupations lower than 0.1 and higher than 0.9 were moved to virtual and doubly occupied spaces, respectively, which resulted in the reference space composed by a complete active space with six electrons in five orbitals (MR-CISD(6,5)). Twelve orbitals were kept frozen, and generalized interacting space restrictions were adopted.59 Single-point energies were corrected for higher-order excitation effects using the Davidson correction (þQ).48,60,61 Additional static calculations were performed with complete active space self-consistentfield second-order perturbation theory in its multistate version (MS-CASPT2)62 and the resolution-of-identity approximate coupled cluster to the second-order (RI-CC2)63
65 method, employing the TZVPP basis sets.66 Dynamics simulations, MRCI, MS-CASPT2, and RI-CC2 calculations were carried out by means of the NEWTONX,67,68 COLUMBUS,69
71 MOLCAS,72 and TURBOMOLE73 suites of programs, respectively.

3. RESULTS AND DISCUSSION 3.1. Potential Energy Surface (PES) of Uracil. Previously reported theoretical studies of the excited states of uracil74
77 show that at the ground-state geometry the first singlet excited state (S1) is of nπ* character followed by a bright S2 (ππ*) state. The latter state dominates the first band of the UV absorption spectrum.54,55,76 The vertical excitation energies and the energies 5248

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A

ARTICLE

Table 1. Vertical Excitation Energies (in eV) and Energies of the S1 and S2 Minima (in eV) Obtained with Several Methods geometry min S0

min S1

min S2

CAS(10,8)a

MS-CASPT2b

CC2c

MR-CISDd

previous results

S1(nπ*)

4.84

5.05

4.83 (0.00)

5.08 (5.20)

5.15,e 4.83f

S2(ππ*)

6.92

5.78

5.41 (0.19)

6.50 (5.91)

5.51,e 5.20,f 5.1g

state S0(cs)

S0(cs)

1.08

1.04

1.25 (0.89)

S1(nπ*)

3.80

4.41

4.11 (4.44)

S2(ππ*)

5.96

5.95

6.32 (5.58)

S0(cs)

1.75

1.65

1.55 (1.41)

0.85f

S1(nπ*) S2(πOπ*)

4.03 5.66

4.65 5.71

4.20 (4.46) 5.52 (5.40)

4.29f 4.56f

a SA-3-CASSCF(10,8)/6-31G*. b MS-3-CASPT2/SA-3-CASSCF(10,8)/6-31G* using CASSCF geometries. c RI-CC2/TZVPP; values of oscillator strengths calculated at the CC2 method are given in parentheses. d MR-CISD(6,5)/SA-3-CASSCF(10,8)/6-31G* using the geometry optimized at the same level; results with Pople corrections are given in parentheses. e SA-3-CASSCF(12,9)-MS-CASPT2.33 f SA-5-CASSCF(8,7)-MS-CASPT2.30 g Maximum of vapor spectrum (refs 55 and 80).

of two lowest excited-state minima localized on the PES (S1(nπ*) and S2(πOπ*)) obtained in this work are presented in Table 1 and are compared with previous results. The numbering scheme and the main orbitals involved in the excitations are displayed in Figure 1. In the vertical excitation region, the first excited state is due to the excitation from the n orbital localized at the O8 atom. The nπ* excitation energies obtained with different methods agree within the range of 0.4 eV. The energy of the second excited state, which is of ππ* character, is more sensitive to the method used. At the CASPT2 and MR-CISDþQ levels, the excitation energies are 0.4 and 0.5 eV higher than the value obtained with the CC2 method, which is in best agreement with experiment. The result obtained with the CASSCF(10,8) method used in the dynamics studies is 1.5 eV too high. This disagreement reflects a wellknown deficiency in the treatment of the ionic ππ* state at the CASSCF level.30,78 Nevertheless, as will be shown later, this method is comparable with higher level methods when exploring the features of the potential energy surface. Selected structural parameters of excited state S1 and S2 minima obtained at the CASSCF(10,8) and MR-CISD(6,5) levels are collected in Table 2. Cartesian coordinates and total energies are provided in the Supporting Information. Note that the geometries obtained at the MRCI level confirm very well the CASSCF result (see Table 2). In the S1 minimum at the CASSCF level, the C4
O8 and C5
C6 bonds are elongated by 0.16 and 0.06 Å, respectively, and the C4
C5 bond is shortened by 0.1 Å with respect to the ground-state minimum structure. These results agree well with previously reported results obtained at the MRCI19 and OM231 levels. The S2 minimum shows similar geometry changes with reference to the ground state as computed for the S1 minimum. In addition, the structure is slightly pyramidalized at the C6 atom. The same changes of geometric parameters in the S2 minimum were obtained also by Hudock et al.30 at the MS-CASPT2 level. Interestingly, it was reported by several authors19,20,31,33 that the unconstrained geometry optimization leads directly to the conical intersection. Energies and characterization of wave functions of the minima on the crossing seams (MXSs) between the S1 and S0 states and between the S2 and S1 states optimized at the CASSCF(10,8) and MR-CISD(6,5) levels are given in Table 3. Structures are shown in Figure 2. Cartesian coordinates and total energies are provided in the Supporting Information. Due to the split of the energies at

Figure 1. Valence bond structure of uracil and molecular orbitals involved in the excitations.

Table 2. Selected Geometric Parameters for Minima Obtained at the CASSCF(10,8)/6-31G*/MR-CISD(6,5)/6-31G* Levels (in Å) S0a S1a S2b a

C4
C5

C5
C6

C4
O8

1.463/1.455

1.345/1.348

1.200/1.241

1.364/1.365

1.408/1.405

1.361/1.365

1.369/1.351

1.486/1.483

1.350/1.365

Planar. b C4C5C6H dihedral angle 139/145.

the MS-CASPT2 level between S1 and S0 for MXS(S1/S0) and S2 and S1 states for MXS(S2/S1) structures optimized at the CASSCF level, the average between the two energy values is given in Table 3. Five MXSs were localized on the S1/S0 crossing seam, corresponding to three qualitatively different structures: two ethylenic structures with C6 or C5 and C6 atoms displaced from the ring plane resulting in ππ*/cs crossing (6E and 6S5), two structures with the O8 atom moved out of plane with boat (3,6B) and planar (oop-O) conformations resulting in the nπ*/cs crossing, and a N3
C4 ring-opening structure. In the ringopening structure the N3 atom is displaced from the ring plane. 5249

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A

ARTICLE

Table 3. Energies (in eV) Relative to the Ground State Minimum and Wave Function Character of the Minima on the S1/S0 and S2/S1 Crossing Seams Calculated with Several Methods

MXS(S1/S0) ring opening 6

S5 oop-O

wave function

CASSCF(10,8)a

MS-CASPT2b,e

MR-CISDc,e

σ(n
π)π*/cs

4.43

4.76

5.25 (4.88)

ππ*/cs nπ*/cs

4.55 5.03

4.28 5.62

4.32 (4.08) 5.33 (5.40)

3,6

nπ*/(ππ*þcs)

5.32

5.91

5.40 (5.62)

6

ππ*/cs

5.82

5.36

5.54 (5.20)

B

E

wave function

CASSCF(10,8)a

MS-CASPT2b,e

MR-CISDc,e

ππ*/σ(n
π)π*/cs

5.25, 4.89d

5.84, 5.01d

5.97, 5.69d (5.57, 5.49)d

MXS(S2/S1) ring opening 3,6

nπ*/ππ*

5.54

5.26

5.25 (5.00)

3

(πþn)π*/(πþn)π*

5.88

6.75

6.21 (6.37)

5

(πþn)π*/(πþn)π*

6.24

6.78

6.40 (6.65)

B

T1 S6

a

SA-3-CASSCF(10,8)/6-31G*. b MS-3-CASPT2/SA-3-CASSCF(10,8)/6-31G* using CASSCF geometries. c MR-CISD(6,5)/SA-3-CASSCF(10,8)/ 6-31G* using geometry optimized at the same level; results with Pople corrections (þQ) are given in parentheses. d S0 energies are given in italics. e Due to the splitting of energies at the MS-CASPT2 and MR-CISDþQ levels, the average value is given.

Due to the distortion of the N3 atom from planarity, the orbital with π character localized mainly on this atom makes an antibonding combination with the nonbonding orbital of the O8 atom forming an antibonding orbital of σ character, in the text labeled further on as σ(n
π) (Figure 1). The C5-puckered conical intersection characterized in previous studies19,31,32 corresponds to our 6S5 MXS. In addition, the conical intersection described by Lan et al.31 characterized by out-of-plane displacement of the O8 atom corresponds to the oop-O conical intersection localized in the present study. The 6S5 conical intersection lies 4.55 eV above the ground state at the CASSCF level. When optimized at the MRCI level, it has E5 conformation and it is stabilized to 4.32 eV, or 4.08 eV if þQ corrections are added. At CASPT2, this conical intersection occurs at about 4.28 eV. The ring-opening conical intersection is the lowest at the CASSCF level, but it is strongly destabilized at the CASPT2 and MRCISþQ levels, becoming the second lowest. As we shall discuss, only a minor fraction of trajectories deactivates at this intersection at the present excitation energies, which means that its too low energy at the CASSCF level does not interfere with the simulations. Similar conclusions can be drawn about the oop-O and 3,6B conical intersections, which are about 0.3
0.4 eV lower at the CASSCF than at the MRCISDþQ. Additionally, these two conical intersections are already too high to be expected to contribute significantly to the dynamics of uracil excited in the first band. Concerning the S2/S1 crossing seam, three of the four structures found in the current study are characterized by strong puckering of the ring. The lowest-energy MXS according to the CASPT2 and MR-CISD calculations is the 3,6B structure, which has been previously characterized also at the MR-CI19,29 and MSCASPT230 levels. In the structures of the energetically higher 3T1 and 6S5 conical intersections, the O8 atom is also displaced from the molecular plane in addition to ring puckering. The lowest S2/ S1 conical intersection at the OM2 level is planar.31 In addition to the puckered structures, an MXS characterized by a broken C4
N3 bond was located at low energy. It corresponds to the crossing of σ(n
π)π*/ππ*/cs states with the MRCI energies of

Figure 2. Structures of (a) S1/S0 and (b) S2/S1 MXSs optimized at the SA-3-CASSCF(10,8)/6-31G* level.

all three states almost degenerate, which means that this structure should be near a section of the crossing seam involving three states, as discussed by Matsika.79 The structures localized on the S2/S1 crossing seam can be grouped in a way similar to that discussed above for S1/S0 intersections. The two most stable structures, ring opening and 3,6 B, are separated from higher-lying structures by 0.3, 0.9, and 0.8 eV at the CASSCF, CASPT2, and MR-CISDþQ levels. The energy of the ring-opening structure is destabilized by 0.6 and 0.3 eV at the CASPT2 and MR-CISDþQ levels, respectively, with respect to the CASSCF level. Consequences of the reversed ordering of the ring-opening and 3,6B conical intersections on the dynamical behavior will be discussed later in the text. Optimized geometries and relative energies of the stationary points and minima on the crossing seams given above show that, except for the overestimation of the excitation energy to the S2(ππ*) state at the vertical region, the critical points on the potential energy surface calculated at the CASSCF(10,8) method are in a good agreement with results obtained with more reliable CASPT2 and MR-CISD(6,5) methods, at least for uracil excited 5250

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A

ARTICLE

Table 4. Experimental and Simulated Time Constantsa and Corresponding Fraction of Trajectories Deactivating via Different Mechanismsb λpump (nm) window τ1 (ps)

τ2 (ps)

c

f2 (%) c

τ3 (ps) c

f3 (%)

250

H

(1.5 (2.4 )

55

267

L

(0.13d) 0.74

23

>1.8 (1.1,d 2.4e)

77

a

Short (τ1), intermediary (τ2), and long (τ3) decay constants. b Experimental results are given in parentheses. c Reference 4. d Reference 1. e Obtained by monoexponential fitting, ref 2.

Figure 3. Potential energy curves along the linear interpolation paths connecting Franck
Condon region, S2 minimum, S2/S1 3,6B MXS, and S1/S0 6S5 MXS. The energies are computed at the SA-3-CASSCF(10,8)/6-31G* level (a) and MR-CISD(6,5)/SA-3-CASSCF(10,8,)/ 6-31G* level (b) with geometries optimized at the same level.

Figure 4. Time evolution of the average population of S0, S1, and S2 states running in the L (6.56 ( 0.25 eV) and H (7.06 ( 0.25 eV) energy windows.

into the first absorption band. The only significant difference is the energetic ordering of ring-opening conical intersections with respect to other conical intersections. Inclusion of dynamic electronic correlation destabilizes these intersections on both crossing seams. To confirm the reliability of the CASSCF method for the dynamics simulation studies, we have also calculated potential energy curves along the linearly interpolated paths connecting S0 min f S2 min f S2/S1 MXS(6,3B) f S1/S0 MXS(6S5) points. The interpolation curves were calculated using the CASSCF(10,8) method (Figure 3a) and MRCI-SD(6,5) method (Figure 3b) with excited state S2 minimum and conical intersections optimized at the same level. In agreement with previously reported interpolation curves connecting S0 min f S2

min f S2/S1 MXS(6,3B), the inclusion of the electron correlation decreases the excitation energy of S2 state in the vertical region and shifts slightly the S2/S1 crossing to lower distance values.30 The main features of the PES of all three states are, however, similarly described for all three methods. This finding justifies the use of the CASSCF method for the dynamics simulations. The small barrier on the path between S2 minimum and S2/S1 crossing obtained with the CASSCF method is only an upper limit since the path was not optimized. In addition, due to the large excess of the energy given to the system to reach the state due to the ππ* transition in the Franck
Condon region, this small barrier can be easily overcome and should not influence the results of dynamics studies. 3.2. Dynamics Results. The time evolution of the population in S0, S1, and S2 states for trajectories started in both energy windows [6.56 ( 0.25 eV (low, L) and 7.06 ( 0.25 eV (high, H)] is shown in Figure 4. The relaxation mechanisms and lifetimes were analyzed for both energy windows separately. In Table 4 the values of the experimental and simulated time constants, labeled as τ1 (short), τ2 (intermediary), and τ3 (long), are reported. The experimental value of τ1 is highly sensitive to data deconvolution and influenced by factors such as coherent processes and instrumental response. Since the present dynamics simulations do not take these factors into account, we do not intend to interpret the dynamics during this time period. As will be discussed below, both direct and indirect decay mechanisms were found to be responsible for the excited-state relaxation of uracil. Therefore, dynamics results were interpreted using τ2 and τ3 decay constants. To obtain the time constant τ2 for populating S0, the quantity 1
F(S0) [F(S0) is the fraction of trajectories in the S0 state] is fitted with the function f ðtÞ ¼ ð1
f2 Þ þ f2 expð
t=τ2 Þ

ð1Þ

where the fitting parameter f2 gives the fraction of the population that decays with time constant τ2. For the case of the L window, values f2 = 0.23 and τ2 = 0.74 ps were obtained (Table 4). This means that 23% of trajectories returned to S0 with a 0.74 ps time constant. To estimate a lower limit for the decay to the ground state of the remaining 77% of trajectories, the fraction of trajectories in S2 was fitted by the function f ðtÞ ¼ A expð
t=τS2 Þ

ð2Þ

which gives τS2 = 1.8 ps for the L window. Since the decay to S0 will occur after this step, then τ3 > τS2, which gives the estimate of the lower limit of τ3. The fraction of trajectories decaying with this constant is given as f3 = 1
f2, i.e., f3 = 77%. In the H window 45% of trajectories decay with τ2 = 0.65 ps and 55% of trajectories with τ3 > 1.5 ps. 5251

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A

ARTICLE

Figure 5. Time evolution of the C5C6 bond for (a) the whole set of trajectories and (b) for trajectories trapped in the S1 state.

As one can see from the slow decay of the S2 population (Figure 4), the reason for the long time constant is the trapping in the S2 state. This can also be seen in the time evolution of the average C5
C6 bond distance. Among the parameters which show the largest distortion from the ground-state geometry (see Table 2), the C5
C6 bond distance is especially well suited to explore details of the time evolution because the differences between optimized values of S0, S1, and S2 minima are well pronounced. The time evolution of the C5
C6 bond distance averaged over the whole set of trajectories is shown in Figure 5a. In the beginning of the dynamics, the C5
C6 bond is quickly stretched as a consequence of the excitation and reaches the value of the S2 minimum within 100 fs. After that, it oscillates around this value for the whole simulation time, which is a consequence of the dominance of the S2-trapping mechanism. When the average value is restricted to the subset of trajectories that decay to the S1 state within the simulation time, the average value is reduced to the shorter values of the S1 minimum (Figure 5b). In Figure 6b, the distribution of the structures occurring at the S2/S1 hopping is shown. For both windows, the deactivation occurs mainly at 6S5, 3,6B, and 3T1 conformations. The distribution of Cremer
Pople parameters of structures at S1 f S0 hoppings (Figure 6a) shows that the majority of the trajectories decay at the structures characterized by out-of-plane distortion of the C6 atom and rotation around the corresponding C6
C5 and C6
N1 bonds. This means that the deactivation occurs preferentially in the ethylenic region of the crossing seam, with a large variety of conformations, centered mainly at 3,6B and 6 S5 conformations. Note that the participation of 3,6B conformations in the deactivation may be overestimated in our simulations due to the relatively low energy of the 3,6B conical intersection at the CASSCF level, as already discussed. Among the trajectories decaying with the time constant τ2, conical intersection structures characterized by ring opening at the N3
C4 bond were also found. They appear more frequently among the trajectories starting in the H window. Three typical situations were observed and are illustrated in Figure 7 for selected trajectories. In Figure 7, each column presents one different trajectory. The upper panels show the evolution of the adiabatic states near the hopping time. The current state at each step is indicated by circles. The bottom panels illustrate the evolution of the weight of main configuration state functions contributing to the wave function at each time. In the first situation (Figure 7A), the ππ* state switches to the nπ* state when hopping from S2 to S1. After that, the dynamics will continue to the S1 minimum before returning to the ground state. Trajectories behaving like this one correspond to the

Figure 6. Cremer
Pople analysis of structures at (a) S1 f S0 and (b) S2 f S1 hoppings. Structures found for H and L energy windows are given by green filled circles and blue open squares, respectively. Optimized structures (MXS) are identified by red crosses.

indirect path discussed in the Introduction. Among the trajectories contributing to the τ3 decay time, a vast majority follow this mechanism. Thus we can expect that most of the S2-trapped trajectories evolve into this indirect pathway. However, longer simulation times would be necessary to confirm this hypothesis conclusively. For the fraction of trajectories contributing to τ2, a different behavior was observed. Once the S1 state is reached, it quickly relaxes into the ground state. The time spent in the S1 state does not exceed 100 fs. Two mechanisms are responsible for this fast decay. The first mechanism is illustrated in Figure 7B. After the change from the S2 to the S1 state, the initial ππ* character experiences a mixing of closed shell contributions. Shortly after the deactivation to the S1 state, the molecule finds its way to the ground state via an ethylenic conical intersection between the ππ* and cs states. This means that the diabatic ππ* character is maintained during the entire deactivation process and no trapping in S1(nπ*) occurs. This kind of trajectory corresponds to the direct pathway referred to in the Introduction. 35% of trajectories in the H window and 21% of trajectories in the L window belong to this case. The second type of the fast mechanisms is illustrated in Figure 7C. With the hopping from S2 to S1, the nπ* configuration mixes into the wave function. During this time the N3
C4 bond breaks. Due to the out-of-plane distortion of the N3 atom, the π orbital localized on this atom makes an antibonding combination with the nonbonding orbital localized on the O8 atom, resulting in orbital labeled as σ(n
π) (see Figure 1). The S1 state can be thus characterized as σ(n
π)π*. In addition, the S0 state destabilizes and the system proceeds to the S1/S0 crossing until it finally relaxes to the S0 state. At this time, the configuration with the doubly occupied σ(n
π) orbital mixes into the final wave function. This mechanism, which has not been previously described, occurs for about 10% of trajectories in the H window and 2.5% of trajectories in the L window. Due to the destabilization of ring-opening conical intersections with inclusion of higher-level correlation effects, the fraction of trajectories decaying with this mechanism is probably overestimated at the CASSCF level. 3.3. Photophysics of Uracil in Comparison with Thymine. Between the two pyrimidine nucleobases thymine and uracil, the 5252

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A

ARTICLE

Figure 7. Analysis of three trajectories following different deactivation mechanisms. (A) Indirect path ππ* f nπ* f S0 (crossing to S0 is not shown). (B) Direct path ππ* f S0. (C) Ring-opening path. Top panels: adiabatic energies. The current state is indicated by circles. Bottom panels: diabatic character of the current state given by the weight of configuration. Hopping events are indicated by arrows.

Figure 8. Scheme of the reaction mechanisms of uracil.

former possesses significantly longer excited-state lifetimes.1,2,4 The experimentally observed excited-state decay of thymine was interpreted in terms of only one relaxation decay constant. The experimental data observed for uracil were interpreted either using two decay time constants, intermediary (τ2) of 0.5 ps4 and

long (τ3) of 1.11,2 and 2.4 ps,4 or using only one decay constant (τ3) of 2.4 ps.2 Results of dynamics simulations reported in the present study suggest that the photophysics of uracil is more complex compared to thymine and the use of two decay constants, one in the intermediary range and another in the long time range, is the most appropriate to fit the experimental data. It has been proposed for thymine that, after the excitation into the S2(ππ*) state, the system quickly relaxes from the Franck
Condon region into the shallow S2(πOπ*) minimum, where it remains for a few picoseconds.30,37 For uracil, a similar occurrence of two types of ππ* states in S2 is found. The dominant mechanism, which contributes to the long time constant observed in the experiment, is the same for thymine and uracil (see Figure 8a) and corresponds to the trapping in the S2 state, in agreement with results of the short dynamics simulations of Hudock et al.30 This is in contradiction to the results of semiempirical dynamics simulations for thymine31,41 and uracil,31 according to which the ground state is significantly populated in less than 1 ps. The simulation time of the present photodynamics calculations is sufficiently long to allow a detailed study of the more complex behavior of uracil compared to thymine. Static calculations on thymine and uracil performed at the same level of calculation show quite similar results concerning the location of conical intersections and deactivation pathways.16,37 The major difference found for the relaxation mechanisms for the two bases mentioned above comes from the direct ππ* f gs (ground state) decay channel (Figure 8b). The ππ* f gs mechanism was not observed in previous dynamics of thymine.37 This difference is caused very likely by a kinematic effect of CH3 group of thymine, as already suggested by Hudock et al.30 On the other hand, the nonoccurrence of the direct path in thymine could also indicate a bias of the CASSCF surface toward the indirect 5253

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A path, given that an intermediary time constant τ2 has also been reported for this nucleobase.4 In addition, ring-opening conical intersections open a new decay channel characterized as ππ* f σ(n
π)π* f gs, contributing to the intermediary time constant (Figure 8b). Although only a small part of trajectories follows this path for the present excitation energies, it may have important consequences for the deactivation process of uracil excited at higher energies. While the deactivation at the puckered structures should enhance the photostability of uracil, the deactivation at the ring-opening conical intersections leads almost certainly to new photochemical products. Matsika19 and Mercier et al.33 predicted the coexistence of two channels, in particular indirect ππ* f nπ* f gs and direct ππ* f gs, in agreement with the present study. According to the analysis of Mercier et al., the direct ππ* f gs channel requires a time scale of 100 fs.33 Our simulations show that the time spent in the direct channel is not so short, in spite of the large CASSCF initial excitation energy. In fact, the results indicate that it is rather responsible for the intermediary time constant.4

4. CONCLUSIONS The photophysics of uracil has been investigated by means of nonadiabatic dynamics simulations performed at the CASSCF level including three electronic states. To justify the use of the CASSCF method for nonadiabatic dynamics simulations, representative points on the potential energy surfaces were computed at the MRCI and MS-CASPT2 levels. The results of optimizations, relative energies of the critical points, and the reaction paths connecting them obtained at the MRCI level confirm very well the CASSCF results. In addition, comparison of the results obtained at the CASSCF and CASPT2 levels shows that the former method adequately describes the character of the potential energy surface and provides a reliable method to be used in dynamics simulations studies of uracil excited in the first absorption band. The present dynamics simulations show that several different pathways and not a single decay mechanism are responsible for the relaxation of uracil. The dominant mechanism (55 and 77% of trajectories, started in the high (H) and low (L) energy windows, respectively) contributes to the relaxation with a time constant larger than 1.5 ps and is related to the experimentally observed time constant of 2.4 ps.1,2,4 It is characterized by the indirect ππ* f nπ* f gs relaxation. Following this path, shortly after excitation in the Franck
Condon region, the system relaxes to the S2 minimum, where it gets trapped for some time before it decays to the S1 state of nπ* character. This mechanism is the same as that found for thymine. Excited-state decay through a direct ππ* f gs channel at an ethylenic conical intersection contributes to the intermediary time constant of about 0.7 ps and should correspond to the experimentally observed decay constant of 0.5 ps. In addition, a new ππ* f σ(n
π)π* f gs channel through a ring-opening conical intersection has been observed for a minor fraction of trajectories. This channel, which also contributes to the intermediary time constant, should lead to new photochemical products rather than back to ground-state uracil. The existence of these two channels not observed in thymine together with a reduced long time constant in uracil is found to be responsible for the significantly shorter lifetime of uracil compared to thymine.

ARTICLE

’ ASSOCIATED CONTENT

bS

Supporting Information. Cartesian coordinates of excited state minima and conical intersections optimized at the SA3-CASSCF(10,8)/6-31G* level. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected] (D.N.); barbatti@ kofo.mpg.de (M.B.); [email protected] (H.L.). Present Address ^

Department of Chemistry, University of Basel, Klingelbergstrasse 80, 4056 Basel, Switzerland.

’ ACKNOWLEDGMENT This work has been supported by the Austrian Science Fund within the framework of the Special Research Program F41 (ViCoM) and by the German Research Foundation, Priority Program SPP 1315, project No. GE 1676/1-1. This work was part of the research project Z40550506 of the Institute of Organic Chemistry and Biochemistry of the Academy of Sciences of the Czech Republic (AS CR) and supported by a grant of the Ministry of Education of the Czech Republic (Center for Biomolecules and Complex Molecular Systems, LC512). Support from the Praemium Academiae of the AS CR, awarded to P.H. in 2007, and computer time at the Vienna Scientific Cluster (Project Nos. 70019 and 70151) are gratefully acknowledged. ’ REFERENCES (1) Canuel, C.; Mons, M.; Piuzzi, F.; Tardivel, B.; Dimicoli, I.; Elhanine, M. J. Chem. Phys. 2005, 122, 074316. (2) Kang, H.; Lee, K. T.; Jung, B.; Ko, Y. J.; Kim, S. K. J. Am. Chem. Soc. 2002, 124, 12958. (3) Crespo-Hernandez, C. E.; Cohen, B.; Hare, P. M.; Kohler, B. Chem. Rev. 2004, 104, 1977. (4) Ullrich, S.; Schultz, T.; Zgierski, M. Z.; Stollow, A. Phys. Chem. Chem. Phys. 2004, 6, 2796. (5) Saigusa, H. J. Photochem. Photobiol. C: Photochem. Rev. 2006, 7, 197. (6) Ullrich, S.; Schultz, T.; Zgierski, M. Z.; Stolow, A. J. Am. Chem. Soc. 2004, 126, 2262. (7) Blancafort, L.; Cohen, B.; Hare, P. M.; Kohler, B.; Robb, M. A. J. Phys. Chem. A 2005, 109, 4431. (8) Kang, H.; Jung, B.; Kim, S. K. J. Chem. Phys. 2003, 118, 6717. (9) Chin, C. H.; Mebel, A. M.; Kim, G. S.; Baek, K. Y.; Hayashi, M.; Liang, K. K.; Lin, S. H. Chem. Phys. Lett. 2007, 445, 361. (10) He, Y. G.; Wu, C. Y.; Kong, W. J. Phys. Chem. A 2003, 107, 5145. (11) He, Y. G.; Wu, C. Y.; Kong, W. J. Phys. Chem. A 2004, 108, 943. (12) Chen, H.; Lis, S. J. Phys. Chem. A 2005, 109, 8443. (13) Perun, S.; Sobolewski, A. L.; Domcke, W. Chem. Phys. 2005, 313, 107. (14) Perun, S.; Sobolewski, A. L.; Domcke, W. J. Am. Chem. Soc. 2005, 127, 6257. (15) Perun, S.; Sobolewski, A. L.; Domcke, W. J. Phys. Chem. A 2006, 110, 13238. (16) Zechmann, G.; Barbatti, M. J. Phys. Chem. A 2008, 112, 8273. (17) Serrano-Andres, L.; Merchan, M.; Borin, A. C. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 8691. (18) Ismail, N.; Blancafort, L.; Olivucci, M.; Kohler, B.; Robb, M. A. J. Am. Chem. Soc. 2002, 124, 6818. (19) Matsika, S. J. Phys. Chem. A 2004, 108, 7584. 5254

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255

The Journal of Physical Chemistry A (20) Zgierski, M. Z.; Patchkovskii, S.; Fujiwara, T.; Lim, E. C. J. Phys. Chem. A 2005, 109, 9384. (21) Marian, C. M. J. Chem. Phys. 2005, 122, 104314. (22) Marian, C. M. J. Phys. Chem. A 2007, 111, 1545. (23) Gustavsson, T.; Banyasz, A.; Lazzarotto, E.; Markovitsi, D.; Scalmani, G.; Frisch, M. J.; Barone, V.; Improta, R. J. Am. Chem. Soc. 2006, 128, 607. (24) Chung, W. C.; Lan, Z. G.; Ohtsuki, Y.; Shimakura, N.; Domcke, W.; Fujimura, Y. Phys. Chem. Chem. Phys. 2007, 9, 2075. (25) Barbatti, M.; Lischka, H. J. Phys. Chem. A 2007, 111, 2852. (26) Nielsen, S. B.; Solling, T. I. ChemPhysChem 2005, 6, 1276. (27) Nachtigallova, D.; Zeleny, T.; Ruckenbauer, M.; Muller, T.; Barbatti, M.; Hobza, P.; Lischka, H. J. Am. Chem. Soc. 2010, 132, 8261. (28) Perun, S.; Sobolewski, A. L.; Domcke, W. Mol. Phys. 2006, 104, 1113. (29) Epifanovsky, E.; Kowalski, K.; Fan, P. D.; Valiev, M.; Matsika, S.; Krylov, A. I. J. Phys. Chem. A 2008, 112, 9983. (30) Hudock, H. R.; Levine, B. G.; Thompson, A. L.; Satzger, H.; Townsend, D.; Gador, N.; Ullrich, S.; Stolow, A.; Martinez, T. J. J. Phys. Chem. A. 2007, 111, 8500. (31) Lan, Z. G.; Fabiano, E.; Thiel, W. J. Phys. Chem. B 2009, 113, 3548. (32) Merchan, M.; Gonzalez-Luque, R.; Climent, T.; SerranoAndres, L.; Rodriuguez, E.; Reguero, M.; Pelaez, D. J. Phys. Chem. B 2006, 110, 26471. (33) Mercier, Y.; Santoro, F.; Reguero, M.; Improta, R. J. Phys. Chem. B 2008, 112, 10769. (34) Hudock, H. R.; Martinez, T. J. ChemPhysChem 2008, 9, 2486. (35) Barbatti, M.; Lischka, H. J. Am. Chem. Soc. 2008, 130, 6831. (36) Groenhof, G.; Schafer, L. V.; Boggio-Pasqua, M.; Goette, M.; Grubmuller, H.; Robb, M. A. J. Am. Chem. Soc. 2007, 129, 6812. (37) Szymczak, J. J.; Barbatti, M.; Hoo, J. T. S.; Adkins, J. A.; Windus, T. L.; Nachtigallova, D.; Lischka, H. J. Phys. Chem. A 2009, 113, 12686. (38) Barbatti, M.; Aquino, A. J. A.; Szymczak, J. J.; Nachtigallova, D.; Hobza, P.; Lischka, H. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 21453. (39) Fabiano, E.; Thiel, W. J. Phys. Chem. A 2008, 112, 6859. (40) Lei, Y. B.; Yuan, S. A.; Dou, Y. S.; Wang, Y. B.; Wen, Z. Y. J. Phys. Chem. A 2008, 112, 8497. (41) Alexandrova, A. N.; Tully, J. C.; Granucci, G. J. Phys. Chem. B 2010, 114, 12116. (42) Langer, H.; Doltsinis, N. L.; Marx, D. ChemPhysChem 2005, 6, 1734. (43) Improta, R.; Barone, V.; Lami, A.; Santoro, F. J. Phys. Chem. B 2009, 113, 14491. (44) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (45) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637. (46) Butcher, J. J. Assoc. Comput. Mach. 1965, 12, 124. (47) Shepard, R.; Lischka, H.; Szalay, P. G.; Kovar, T.; Ernzerhof, M. J. Chem. Phys. 1992, 96, 2085. (48) Shepard, R. The Analytic Gradient Method for Configuration Interaction Wave Functions. In Modern Electronic Structure Theory; Yarkony, D. R., Ed.; World Scientific: Singapore, 1995; Vol. 1, p 345. (49) Lischka, H.; Dallos, M.; Shepard, R. Mol. Phys. 2002, 100, 1647. (50) Lischka, H.; Dallos, M.; Szalay, P. G.; Yarkony, D. R.; Shepard, R. J. Chem. Phys. 2004, 120, 7322. (51) Dallos, M.; Lischka, H.; Shepard, R.; Yarkony, D. R.; Szalay, P. G. J. Chem. Phys. 2004, 120, 7330. (52) Granucci, G.; Persico, M. J. Chem. Phys. 2007, 126, 134114. (53) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (54) Barbatti, M.; Aquino, A. J. A.; Lischka, H. Phys. Chem. Chem. Phys. 2010, 12, 4959. (55) Clark, L. B.; Reschel, G. G.; Tinoco, I. J. J. Phys. Chem. 1965, 69, 3615. (56) Nachtigallova, D.; Barbatti, M.; Szymczak, J. J.; Hobza, P.; Lischka, H. Chem. Phys. Lett. 2010, 497, 129. (57) Cremer, D.; Pople, J. A. J. Am. Chem. Soc. 1975, 97, 1354.

ARTICLE

(58) Boeyens, J. C. A. J. Chem. Crystallogr. 1978, 8, 317. (59) Bunge, A. J. Chem. Phys. 1970, 53, 20. (60) Langhoff, S. R.; Davidson, E. R. Int. J. Quantum Chem. 1974, 8, 61. (61) Bruna, P. J.; Peyerimhoff, S. D.; Buenker, R. J. Chem. Phys. Lett. 1980, 72, 278. (62) Finley, J.; Malmqvist, P. A.; Roos, B. O.; SerranoAndres, L. Chem. Phys. Lett. 1998, 288, 299. (63) Hattig, C.; Weigend, F. J. Chem. Phys. 2000, 113, 5154. (64) Hattig, C.; Kohn, A. J. Chem. Phys. 2002, 117, 6939. (65) Christiansen, O.; Koch, H.; Jorgensen, P. Chem. Phys. Lett. 1995, 243, 409. (66) Schafer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (67) Barbatti, M.; Granucci, G.; Lischka, H.; Ruckenbauer, M.; Persico, M. NEWTON-X: a package for Newtonian dynamics close to the crossing seam, version 0.14b; 2007; www.univie.ac.at/newtonx. (68) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; Eckert-Maksic, M.; Lischka, H. J. Photochem. Photobiol., A: Chem. 2007, 190, 228. (69) Lischka, H.; Shepard, R.; Brown, F. B.; Shavitt, I. Int. J. Quantum Chem. 1981, S15, 91. (70) Lischka, H.; Shepard, R.; Pitzer, R. M.; Shavitt, I.; Dallos, M.; Muller, T.; Szalay, P. G.; Seth, M.; Kedziora, G. S.; Yabushita, S.; Zhang, Z. Y. Phys. Chem. Chem. Phys. 2001, 3, 664. (71) Lischka, H.; Shepard, R.; Shavitt, I.; Pitzer, R. M.; Dallos, M.; Mueller, T.; Szalay, P. G.; Brown, F. B.; Ahlrichs, R.; Boehm, J. G.; Chang, A.; Comeau, D. C.; Gdanitz, R.; Dachsel, H.; Ehrhardt, C.; Ernzerhof, M.; Hoechtl, P.; Irle, S.; Kedziora, G.; Kovar, T.; Parasuk, V.; Pepper, M. J. M.; Scharf, P.; Schiffer, H.; Schindler, M.; Schueler, M.; Seth, M.; Stahlberg, E. A.; Zhao, J.-G.; Yabushita, S.; Zhang, C. L.; Barbatti, M.; Matsika, S.; Schuurmann, M.; Yarkony, D. R.; Brozell, S. R.; Beck, E. V.; Blaudeau, J.-P.; Ruckenbauer, M.; Sellner, B.; Plasser, F.; Szymczak, J. J. COLUMBUS, an ab initio electronic structure program, release 5.9.2; 2008; www.univie.ac.at/columbus. (72) Karlstrom, G.; Lindh, R.; Malmqvist, P. A.; Roos, B. O.; Ryde, U.; Veryazov, V.; Widmark, P. O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L. Comput. Mater. Sci. 2003, 28, 222. (73) Ahlrichs, R.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Chem. Phys. Lett. 1989, 162, 165. (74) Lorentzon, J.; Fulscher, M. P.; Roos, B. O. J. Am. Chem. Soc. 1995, 117, 9265. (75) Marian, C. M.; Schneider, F.; Kleinschmidt, M.; Tatchen, J. Eur. Phys. J. D 2002, 20, 357. (76) Petke, J. D.; Maggiora, G. M.; Christoffersen, R. E. J. Phys. Chem. 1992, 96, 6992. (77) Neiss, C.; Saalfrank, P.; Parac, M.; Grimme, S. J. Phys. Chem. A 2003, 107, 140. (78) Muller, T.; Dallos, M.; Lischka, H. J. Chem. Phys. 1999, 110, 7176. (79) Matsika, S. J. Phys. Chem. A 2005, 109, 7538. (80) Clark, L. B.; Tinoco, I. J. Am. Chem. Soc. 1965, 87, 11.

5255

dx.doi.org/10.1021/jp201327w |J. Phys. Chem. A 2011, 115, 5247–5255