Photoelectron Spectrum and Dynamics of the Uracil Cation - The

Article pubs.acs.org/JPCA

Photoelectron Spectrum and Dynamics of the Uracil Cation Mariana Assmann,*,† Horst Köppel,‡ and Spiridoula Matsika† †

Department of Chemistry, Temple University, Philadelphia, Pennsylvania 19122, United States Theoretische Chemie, Universität Heidelberg, Im Neuenheimer Feld 229, 69120 Heidelberg, Germany

‡

S Supporting Information *

ABSTRACT: The photoelectron spectrum of uracil and the molecular dynamics of its radical cation are investigated using the multiconfigurational time-dependent Hartree (MCTDH) method. For this aim, the vibronic coupling model Hamiltonian is used including up to ten important a′ modes. Moreover, to account for coupling through conical intersections between states of different symmetry in the system, coupling constants of two a″ modes are taken into account. The parameters used in the model are obtained by fitting to ab inito data obtained with extensive EOM-IP-CCSD calculations. The first four cationic states were investigated, which are either of A″ (hole in a π orbital) or A′ (hole in a nO orbital) symmetry. The results of the wavepacket propagations were used to calculate the corresponding photoelectron spectrum and compare to the experimental spectrum. The MCTDH simulations reproduce the experimental spectrum well. The dynamics starting from the D2 and D3 ionic states show a fast relaxation to the cationic ground state often involving direct D2−D0 or D3−D1 transitions.

■

INTRODUCTION The formation and subsequent dynamics of cations play an important role in ionization processes. Radiation damage in DNA/RNA creates radical cations of the bases which are the initial reactive intermediates of one-electron oxidation.1,2 Radical-cation states are also important in interpreting pump− probe experiments where ionization is used as the probe, such as in time-resolved photoelectron (PE) spectroscopy, or strongfield dissociative ionization. In order to interpret the results of the experiments, we need to understand the dynamics of the cations.3 Detailed studies of the dynamics in radical cations involving several ionic states are limited. Recent examples are studies on the propanal4 and ethylene cations5 using the multiple spawning method. Furthermore, radical cations, and particularly PE spectra, have been investigated using the multiconfiguration time-dependent Hartree (MCTDH) method for several systems, such as benzene,6,7 allene,8 cyclobutadiene,9 butatriene10 and smaller heterocycles.11 The MCTDH method has proven to be particularly efficient for treating nonadiabatic dynamics of highdimensional systems like those mentioned above and even bigger molecules.12−14 Furthermore, wavepacket propagations can be used to calculate the PE spectrum through the corresponding MCTDH program.15 Our motivation for this work comes from previous studies where theory was combined with experiments on uracil using strong-field dissociative ionization.16,17 Strong-field ionization can populate several ionic states,17,18 and the dynamics following this population determine the dissociation patterns. An assumption that we have made in our work so far is that relaxation from the higher ionic states to the ground state is very fast and precedes fragmentation. Thus, fragmentation occurs on the ground state. Following this assumption we calculated the © 2015 American Chemical Society

possible fragmentation pathways in the uracil cation in order to obtain dissociation energies.19 The assumption of fast relaxation before fragmentation is a reasonable one, but its generality should be tested. In this work we are studying the relaxation of the uracil cation starting from different ionic states in order to check the validity of fast decay. We have some information about critical features of the ionic states; specifically we have previously located two- and three-state conical intersections (CIs) among ionic states in uracil. We have found that minima on the two-state CI seam and the three-state CI seam exist about 0.5 and 1 eV above the D0 minimum, respectively.20 These CIs likely modulate the dynamics in the uracil cation, however, the dynamics behind the processes is still unknown. Both CIs are planar and, thus, retain the Cs symmetry. The first four cationic states of uracil have alternative A″ and A′ symmetries: D0 (A″), D1 (A′), D2 (A″) and D3 (A′). From our dynamics we can calculate the PE spectrum using the autocorrelation function and compare to the experimental one. The first PE spectrum of uracil was reported by Padva et al.21 The bands of the experimental spectrum have been assigned to the four mentioned states in several previous studies with the help of different theoretical methods, such as SCF-MO and HAM/3,22 CNDO/s23 as well as configuration interaction.24 The aim of this work is to unravel the relaxation dynamics and their time scales in the uracil radical cation using the MCTDH method taking into account the first four cationic states. Other spin states than doublets have not been considered in this work. We are using the vibronic coupling Hamiltonian25−27 consisting of up to 10 a′ modes. In order to account for coupling between Received: December 8, 2014 Revised: January 6, 2015 Published: January 7, 2015 866

DOI: 10.1021/jp512221x J. Phys. Chem. A 2015, 119, 866−875

Article

The Journal of Physical Chemistry A

Here, γ(α) ij are the second derivatives of the adiabatic PESs and μ(αβ) are the second order nonadiabatic couplings between states ij α and β. All parameters, κ, λ, γ and μ correspond to the values at q0. In this work all μ are neglected and γ are only considered for i = j. Historically, these terms have been neglected for states of different symmetry. Furthermore, fitting properly these terms would require multidimensional grids, which would increase the computational effort dramatically. Finally, as our results will show, these neglected terms are not necessary for an adequate representation of the spectra. In our case, q0 corresponds to the geometry of the S0 minimum of uracil and Eα correspond to the ionization potentials. In most cases, the expansion of the potentials could be cut after including the linear terms and the single mode and single state γ(α) ii . In a few a′ modes, anharmonicity had to be accounted for by replacing the potential term in eq 2 with a Morse potential of the form

neighboring states of different symmetry, two selected a″ modes are added. From the wavepacket propagations the PE spectrum is calculated and compared to experimental results as a reference for the accurateness of our models. In the following section, we will give an overview of the methods and computational details used. Afterward, we present the results of the fitting of the ab initio data, the PE spectra and the populations dynamics. Finally, the last section summarizes our work.

■

METHODOLOGY Vibronic Coupling Hamiltonian and fitting of the PESs. The vibronic coupling Hamiltonian25−27 is a diabatic expression of the potential energy surfaces (PESs) in terms of the massweighted vibrational normal modes q around a reference geometry q0. The diabatic and adiabatic descriptions of the PESs coincide at q0. The electronic Hamiltonian comprises the zeroth order Hamiltonian H(0) and a Taylor series of coupling matrices W(m): (0)

H=H

+W

(1)

+W

(2)

+W

(3)

+ ...

V i(α) = d0,(αi )[exp(ai(α)(qi − q0,(αi ))) − 1]2 + e0,(αi )

(1)

Additionally, for fitting the a″ modes, quartic terms of the form

The elements of the zeroth order Hamiltonian are the sum of the kinetic and potential energy terms of the vibrational modes qi expressed as harmonic oscillators for each electronic state α (in atomic units): (0) Hαα = Eα + Ti + V i(α) = Eα +

∑ i

⎞ ωi ⎛ ∂ 2 2⎟ ⎜ q + i ⎟ 2 ⎜⎝ ∂qi2 ⎠

(4) W αα =

=

∑

κi(α)qi

i (1) W αβ =

∑ λi(αβ)qi i

(2)

∂H ϕ ∂qi α

ϕα

λiαβ =

∂H ϕ ∂qi β

ϕα

(4)

(2) W αβ

1 2

1 = 2

∑ γij(α)qiqj i,j

∑ i,j

(10)

κi(α) ≠ 0, Γi ⊃ ΓA

(11)

λi(αβ) ≠ 0, Γα ⊗ Γi ⊗ Γβ ⊃ ΓA

(12)

(13)

Since we only considered γ(αα) ii , the multiplication of all the modes and states involved in that case always yields the totally symmetric representation. All PESs and coupling parameters are obtained by fitting the model to the ab initio data with the help of a set of programs in the MCTDH suite (VCHAM) which perform a least-squares fit. MCTDH method. For the quantum dynamics simulations the MCTDH method was used, which is an efficient tool to treat high dimensionality systems. Details of the MCTDH algorithm have been described comprehensively elsewhere (refs 12, 13) and, therefore, will only be outlined briefly here. The multistate ansatz for the wave function in the MCTDH approach is written as

(5)

(6)

ns

ψ=

(7)

∑ ψ (α)|α⟩ α=1

μij(αβ)qiqj

i

γij(αβ) ≠ 0, Γα ⊗ Γβ ⊗ Γi ⊗ Γj ⊃ ΓA

with ϕα being the diabatic electronic wave functions. The truncation of the coupling matrices after the linear terms gives rise to the linear vibronic coupling (LVC) model.25,28 The terms of the second order coupling matrices, which lead to the quadratic coupling model (QVC), are written in a similar way as (2) W αα =

∑ ki(α)qi4

with Γα and Γβ being the irreducible representations of the state, Γi the symmetry of the mode and ΓA the totally symmetric representation. For uracil that possesses Cs symmetry this means that only for the a′ modes the on-diagonal κ(α) i are nonzero. The coupling elements γ(αβ) are nonzero when multiplying the state i symmetries yields the symmetry of the mode. For the second order γ(αβ) it is ij

(3)

The linear intrastate couplings κ(α) i are related to the gradients of the adiabatic PESs with respect to the nuclear coordinates, and λ(αβ) are the linear couplings between the states: i κ iα =

1 24

were found necessary to be added in the PES description in order to better describe the double wells. If symmetry is present in a molecule, certain terms of the coupling matrices vanish due to symmetry reasons. The following restrictions apply for the first order parameters

where ωi is the vibrational frequency for mode i which is taken to be the same for all states. Eα is the energy of the electronic state α at q0. The off-diagonal terms of the zeroth order Hamiltonian are equal to 0. The first order coupling matrix elements are expressed as (1) W αα

(9)

(14)

where the wave function on each state, ψ(α), is expanded as a linear combination of Hartree products:

(8) 867


Article

The Journal of Physical Chemistry A nαf

n1α

(α)

ψ (q1 , ..., qf , t ) =

∑ j1α = 1

...

multireference character in the ionic states. It provides very good agreement with experimental data, which can be seen in the following section. The same method was used for single point calculations to scan the cationic PESs of all vibrational modes i around q0. Four cationic states were taken into account.

f

∑ A jα ...j (t ) ∏ φ j(k ,α)(qk , t )

j fα = 1

α 1

α f

k=1

α k

(15)

■

where q1,...,qf are the nuclear coordinates for the f degrees of freedom. Every coordinate qk can be a composite coordinate of one or more system coordinates. However, here we do not contract any coordinates. Aj1αα···jfα are the time-dependent expansion coefficients and φ(k,α) are the nk time-dependent jαk basis functions for the degree of freedom k, which are the socalled single particle functions (SPFs). Here, the SPFs are represented in a time-independent discrete variable representation (DVR),29 more precisely the harmonic oscillator DVR12 where 32 grid points have been used for all DOFs. The number of SPFs was adapted to every wavepacket propagation and a full list of the SPF numbers can be found in the Supporting Information. The equations of motion to describe the evolution of the wavepacket are then derived for the expansion coefficients and the SPFs employing the Dirac−Frenkel variational principle. A detailed description of the derivation of these equations can be found in ref 12. Dynamics simulations were run starting from the first four cationic states, each for 120 fs. The initial wavepackets were obtained from relaxation of a guess wavepacket in the S0 potentials of the included modes. The initial wavepacket was transferred to the cationic states via instantaneous vertical excitation. For relevant propagations the adiabatic state populations were calculated using the adpop program that is implemented in the MCTDH program package; in the cases with six and more a′ modes Monte Carlo integration was employed. Calculation of the PE Spectra. The PE spectra Pα(E) describing the excitation to the state α were calculated from the respective wavepacket dynamics by Fourier transform of the autocorrelation function Cα(t): P α(E) ∼

1 Re π

∫0

∞

C α(t )eiEt dt

RESULTS AND DISCUSSION Fitting of the Potential Energy Surfaces. PESs were calculated for the first four cationic states, which are described by configurations with an electron hole in orbitals πCC (A″), nO,2 (A′), πNO (A″), and nO,1 (A′) for D0, D1, D2, and D3, respectively. Figure 1 shows the orbitals involved in these states. The first

Figure 1. Hartree−Fock orbitals participating in the EOM-IP-CCSD description of states.

vertical IP at the EOM-IP-CCSD level of theory is 9.42 eV compared to the experimental 9.59 eV.21 The IPs to the higher states D1, D2, and D3 are 10.11 (10.11), 10.48 (10.56), and 11.08 (11.16) eV, respectively, with the experimental values given in parentheses. The errors between the experimental and theoretical IPs are between 0.00 and 0.17 eV. Initial crude scans along all normal modes were used to identify the modes that may be contributing to the spectra and dynamics. The modes that were determined to be important and were included in our calculations are shown in Figure 2, and the labeling of atoms is shown in Figure 3. The most important a′ modes are the ones corresponding to carbonyl stretch (denoted ν25 and ν26), CN stretch plus NCH/CNH bending motions (ν21, ν20, ν18), and the C5C6 stretch (ν24). Further, the NCO bend (ν3) and another CN stretch + bending motion (ν19), the ring breathing deformation (ν11) and another ring distortion (ν7) were considered. These latter modes were not as important as the first six modes but they also contribute to the spectrum as will be discussed below. In total, up to 10 a′ modes have been included in the calculations. Table 1 summarizes which fitting model was used and the coordinate range of the fit for each mode. The ab initio data of modes ν3 and ν18−ν26 were fitted using the second order model (QVC), i.e. via κ, λ and γ. For ν7 and ν11, the linear model was used. For ν24−ν26 Morse potentials were needed instead of the harmonic oscillator potentials in eq 2 to account for the anharmonicity of the PESs. The fitted parameters can be found in Table 2 for ν3−ν26 and in Table 3 for ν24−ν26. The fitted PESs are plotted together with the ab initio data in Figure 4. The a″ modes provide the coupling between the A′ and A″ ionic states. Identifying the a″ modes was more difficult because the corresponding coupling constants are weak in this case. To investigate their influence, scans along the a″ modes were performed at different displaced geometries of the a′ modes, that is at points that exhibit degeneracies of the corresponding a′ mode PESs. Cuts at q21 = −4 were used. This point corresponds to a three-state CI of D0, D1, and D2 (see the three curves approaching the three-state CI in Figure 4g. The range on the x axis has been extended beyond the point used in the fitting in order to show this CI). The most significant couplings were found to be along ν10 and ν12 at q21 = −4. The importance of the modes is determined by the coupling strength that is reflected in

(16)

The autocorrelation function is given through α

α

C α(t ) = ⟨ψgs|Ô e−i Ht Ô |ψgs⟩

(17)

with ψgs being the S0 vibrational ground state, i.e., initial wavepacket. Ô α is the excitation operator to the state α. The broadening of the experimental spectrum was simulated by multiplication with a damping function exp(−t/τ) to convolute the spectral lines. Here, we chose the damping time τ to be 14.1 fs. Furthermore, a weight function cos (πt/2T) was used to reduce problems arising from the finite propagation time T in the Fourier transformation. Quantum Chemistry Calculations. The optimization of the S 0 minimum of uracil and the calculation of the corresponding harmonic vibrational frequencies was done using B3LYP/aug-cc-pVTZ to have the best agreement with experimental data according to ref 30. The calculations were performed with the Gaussian03 program package.31 The cationic states were obtained with the equation-of-motion coupledcluster method with single and double excitations for the calculation of ionization potentials32 (EOM-IP-CCSD) and the basis set 6-311 + G* using the Q-chem program.33 This method calculates the ionization potentials using as reference the neutral molecule which is a well behaved closed shell system, and thus, avoids any problems due to the high density of states or 868


Article


Figure 2. Important vibrational normal modes and corresponding vibrational frequencies.

the shape of the PESs. The bigger the coupling, the more repulsion can be seen between the single PESs. It should be noted that if scans starting at a different point are used then the couplings for these two modes change somewhat, even though our model dictates that they should be universal. This is a weakness of the model, but as it will be seen later, the spectra are not sensitive to the a″ modes, so we expect this to have a small effect in our results. The PESs of the a″ modes were fitted using the quartic potentials for the diagonal terms and the linear model for the couplings. The fitted coupling and potential parameters can be found in Table 3. The PESs are shown in Figure 4k,l. The PES of D0 flattens and has a double well while the PESs of the other states are steeper. All fittings are in excellent agreement with the ab initio data in the ranges considered. Photoelectron Spectrum. Analysis of Reduced Dimensionality Spectra. Figure 5 shows the PE spectra calculated with an increasing number of a′ modes included. Spectra in which the a″ modes are not considered yet will be referred to as PE1 spectra. The most important modes that are responsible for the main vibrational progressions of the D1−D3 bands are ν21, ν25 and ν26. The corresponding three-dimensional (3D) spectrum is depicted in Figure 5a. To determine which bands are produced by which modes, single mode spectra for these three modes were calculated as well. These spectra are depicted in Figure 6. The contribution from the modes to the different bands correlates with the gradient of the respective modes at q0 which is given by κ. The contributions are what one would expect based on the character of each state. The vibronic structure of the D1 state is dominated by mode ν25 which has the largest gradient at q0 in D1. D1 corresponds to a hole in nO,2, which is the n orbital at O8. The stretch of C4O8 is described by ν25, therefore it is reasonable that the vibrational progression of the D1 band is determined by this mode. A similar relation applies to D3 being the hole in nO,1

Figure 3. Atom labeling in uracil.

Table 1. Summary of the Models Used for Each Mode mode

range

model A′

ν21 ν25 ν26 ν24 ν18 ν20 ν3 ν19 ν11 ν7

[−3:3] [−4:2] [−2:4] [−3:3.5] [−2:2.25] [−3:3] [−4:4] [−3:3] [−3:2] [−2:3] A″ [−2:2] [−2:2]

ν10 ν12

QVC QVC + Morse QVC + Morse QVC + Morse QVC QVC QVC QVC LVC LVC LVC + quartic LVC + quartic

a (α,β) Table 2. Fitted Parameters κ(α) , and γ(α) i , λi ii in eV for States D0−D3 (α = 0−3) and Modes ν3, ν7, ν11, ν18−ν21 (i = 3, 7, 11, 18−21)

mode

κ(0)

κ(1)

κ(2)

κ(3)

λ(02)

ν3 ν7 ν11 ν18 ν19 ν20 ν21

0.041 39 −0.053 67 0.053 57 −0.022 03 0.074 72 −0.121 47 −0.094 68

−0.026 88 −0.005 03 −0.024 56 0.090 74 −0.005 82 0.053 16 0.044 54

−0.058 53 0.007 75 0.036 97 0.027 48 0.008 89 0.112 33 0.145 39

0.001 32 −0.045 81 −0.021 30 −0.040 54 −0.021 36 0.007 47 0.000 50

0.021 41 0.020 82 −0.035 38 −0.037 63 −0.020 49

λ(13)

0.018 45 0.080 77 0.058 34 0.072 84

γ(0)

γ(1)

γ(2)

γ(3)

0.003 83

−0.004 26

−0.003 66

−0.009 47

0.019 38 0.015 90 0.014 89 0.009 70

0.006 94 0.013 48 0.008 28 0.000 96

−0.002 94 0.009 01 0.001 83 −0.001 14

0.007 52 0.004 97 0.005 46 0.011 08

When a value is missing from the table, it is either because it was not included in the fitting or if it was included the fitting procedure returned a value of zero. a

869


Article

The Journal of Physical Chemistry A Table 3. Fitted Parameters for States D0−D3 (α = 0−3) and Modes ν24−ν26, ν10, and ν12 (i = 24, 26, 10, 12)a 2 (α) (α) (α) (α) Morse potentials V(α) i = d0,i [exp(ai (qi − q0,i )) − 1] + e0,i

mode (state) ν25 (D0) (D1) (D2) (D3) ν26 (D0) (D1) (D2) (D3) ν24 (D0) (D1) (D2) (D3) mode ν10 ν12 a

d0

a

4.802 70 −0.136 75 74.159 95 −0.030 64 90.769 28 −0.033 74 20.560 79 −0.080 44 22.928 02 0.074 38 18.274 40 0.079 11 9.468 94 0.086 53 65.096 78 0.036 60 41.897 04 0.047 19 38.371 22 0.052 31 39.256 91 0.052 86 37.978 47 0.054 31 (α) 4 quartic potentials V(α) i = (1/(24))ki qi (0)

k

0.033 17 0.029 79

coupl. const.

q0

e0

λ(02)

λ(13)

0.028 83 −1.344 68 −0.299 23 0.388 41 −0.320 69 −0.017 11 0.376 35 1.663 12 0.814 40 0.374 88 0.148 59 −0.181 52

−0.000 07 −0.120 82 −0.009 16 −0.020 71 −0.012 74 −0.000 03 −0.010 37 −0.256 39 −0.064 31 −0.015 05 −0.002 44 −0.003 66

0.001 14

0.126 06

0.130 35

0.142 72

−0.018 32

coupl. const.

(1)

(2)

k

k

0.011 57 0.014 88

0.015 34 0.016 71

λ

(01)

0.046 33 0.035 40

λ(12) 0.031 48 0.036 07

(α) (α,β) (α) All κ(α) , and e(α) and q(α) i , d0,i , λi 0,i are given in eV. a 0,i are dimensionless.

which is the n orbital at O7 and, thus, related to the stretch of C3O7, i.e., ν26. As Figure 6 shows the D3 band structure is caused by mode ν26. Although when only this mode is included three distinct peaks appear, in the combined 3D spectrum a main peak appears, while the others broaden and overlap. The electronic wave function of D2 is described by a hole in the π orbital located at the CN bonds and the oxygen atoms. Hence, ν21 that includes CN stretch and NCH/CNH bending motions has the most significant contribution on the D2 band with a κ of 0.14539 eV. These three modes are used initially to produce the 3D spectrum. The next set of modes that should contribute to the spectrum based on their κ values are modes, ν24, ν18, and ν20, followed by ν3, ν19, ν11, and ν7. The corresponding six-dimensional (6D) and 10dimensional (10D) PE spectra can be found in Figure 5b,c. From these figures it can be seen that the influence of the additional modes is less important for the D1 and D3 bands but it is quite important for D0 and less so for D2. When ν24, ν18, and ν20 are included (Figure 5b), the shoulder of the D0 band becomes a distinct peak with an intensity comparable to that of the first peak. A third peak also appears at 10.1 eV. Since ν24 is the C5C6 stretching mode it is obvious that the vibrational structure of the D0 πCC state should be influenced by ν24. Furthermore, a decrease in the relative height of the main peak of the D0 and also D2 bands with respect to the other bands is also seen. In the D2 band the second peak becomes broader merging somewhat with the first peak. The inclusion of the modes ν3, ν19, ν11, and ν7 in the 10D spectrum affects mainly the D0 and D2 bands. Both bands partially lose their clear shape (Figure 5c) as more modes are included, with the D0 band being much more affected. This is because these modes pose a certain displacement on their D0 PESs, and thus, have significant κ(0) values (0.04139, 0.07472, i 0.05357, and −0.05367 eV) and κ(0) i /ωi ratios. Adding a″ modes changes dramatically the dynamics (as will be shown in the section Dynamics), but it has much smaller effects on the spectrum. Figure 7 shows the 6D spectra with the presence of λ10 only (PE2) or both λ10 and λ12 (PE3) in direct comparison to PE1. The effect is very small when including only λ10 (what can be seen in the similarity of the black and the blue

spectrum) except for the D1 band in which the relative intensities change by a small amount. Adding λ12 visibly alters the D1 vibrational progressions and slightly the one of D2. Comparison of the Spectra with Experimental Data. The experimental PE spectrum taken from ref 21 in the region between 9 and 12 eV is shown in Figure 8a. The spectrum shows a broad band between 9 and 10 eV with a peak at 9.6 eV corresponding to D0, then there is a structured band with three peaks, which has been assigned to D1. The next band associated with D2 is peaked at 10.56 eV and shows no structure again, and finally the D3 band shows a maximum at 11.16 eV and a shoulder on the right side of the spectrum. In order to compare the experimental spectrum to our calculations we choose the spectrum with the higher number of degrees of freedom. It is shown in Figure 7 that the a″ modes do not have an important effect on the spectrum. On the other hand, all 10 a′ modes are important for D0, so we will use the 10D PE1 spectrum for our comparisons. One discrepancy between theory and experiment that we would like to take into account is that we predict the spacing between the D1 and D2 electronic energies to be somewhat smaller than the experimental spacing. In the experimental spectrum the gap is 0.45 eV, whereas the calculated states are only 0.37 eV apart. In order to take into account this effect we have shifted the spectra for each state according to the experimental spacing between them rather than the theoretical one. The resulting spectrum for 10D PE1 is shown in Figure 8b while Figure 8a shows the experimental and shifted theoretical spectra overlaid. It can be seen that the agreement with the experimental spectrum is quite good. More detail comparisons follow. The main difference between the theoretical and experimental spectra is in the shape of D0. Even though the experimental D0 band shows no structure, our calculated spectrum shows a structure for all bands. D0 is actually the hardest band to describe. When only three degrees of freedom are included, in 3D PE1, a single narrow peak with a small shoulder appears. In that spectrum it is the highest peak in the calculated spectrum while in the experimental one it is much broader and not as high. Adding a′ modes causes the band to broaden with several peaks contributing to it. Also, as the number of modes increases the 870


Article


Figure 4. Fitted adiabatic PESs of the a′ modes (a−j) and the a″ modes (k,l). The ab initio points are indicated as dots. In part g, the plot range of the PESs of ν21 were extended beyond the fitting range to allow for showing the three-state CI. For the a″ modes in parts k and l, qi = 0 corresponds to the point of q21 = −4.

relative intensity of the vibrational peaks changes, and the band approaches more the experimental one. It appears that including even more a′ modes could lead to the peak overlapping even more and eventually become one broad peak as seen experimentally. In general, many low frequency modes seem to contribute to D0. D1 shows three peaks in all of our calculated spectra similarly to what is seen in the experimental one. The fact that only one vibrational mode (ν25) contributes mostly to the progression in

this band makes its description much easier than D0. The theoretical spectrum contains the three distinct peaks of the experimental D1 band. The separation between the first two peaks is around 1300 cm−1 in PE1 (1500 cm−1 in 6D PE2/PE3), while the experimental one is reported to be 1200 cm−1. The separation between the second and third peak is calculated to be around 1300 cm−1 in PE1 (1500 cm−1 in 6D PE2/PE3), while the experimental one is 1450 cm−1. The values of these splittings are smaller than the vibrational frequency of mode ν25 (1761 871


Article


Figure 5. PE1 spectra: (a) 3D including modes ν21, ν25, and ν26; (b) 6D including modes ν24, ν18, and ν20, in addition to the 3D ones; (c) 10D including modes ν3, ν7, ν11, and ν19 in addition to 6D. Figure 8. (a) Experimental spectrum (modified with permission from ref 21 Copyright 1974 Elsevier) and (b) 10D PE1 spectrum shifted to the experimental values.

are used. Both D0 and D2 are A″ states involving orbitals in the ring and they are affected a lot by all the low frequency modes. On the other hand the D1 and D3 states are A′ states and they are mainly affected by the carbonyl stretches producing easier and more clean spectra. The D3 band in our calculations shows very good agreement with experiment. A strong peak followed by a shoulder shifted about 1700 cm−1 is seen in both our calculations and the experiment. The qualitative shape is actually already reproduced at the 3D PE1 spectrum, and there is no need for this band to include any more than 3 modes. The value of the separation here is very similar to the vibrational frequency for mode ν26 which is the mode primarily responsible for the shoulder. Dynamics. The diabatic and adiabatic state populations corresponding to the PE1 and PE3 spectra are shown in Figure 9 and Figure 10, respectively, for selected 3D, 6D, and 10D wavepacket propagations. When only three degrees of freedom are used in the 3D PE1 case (Figure 9a, Figure 10a), and particularly no a″ modes are used, no population transfer between A′ and A″ states can occur, so only states with the same symmetry as the initial state are populated. Therefore, the diabatic state populations of the dynamics starting in D1 are constant. When starting in D2 there is some population transfer directly to D0, bypassing D1. Similarly when the propagation is initiated in D3 there is population decay from D3 to D1 but the population on D0 remains zero. When transforming from the diabatic to the adiabatic populations the picture changes somewhat since the states may change symmetry when a crossing occurs. For example, when D0 crosses with D1 the character of the two states will switch and therefore the symmetry. The diabatic populations describing dynamics from D1 and D3 do not change significantly when adding modes ν24, ν18 and ν20 (i.e., going from 3D to 6D), but they change more for D2 where a lot more population transfer to D0 occurs. The remaining additional modes do not contribute to a greater population transfer, hence the population dynamics of the higher dimensional wavepacket propagations do not differ significantly from the 6D case. In the adiabatic picture the populations are affected more when adding more a′ modes. In all cases the dynamics

Figure 6. Single mode spectra for (a) ν21, (b) ν25, and (c) ν26.

Figure 7. Overlay of the 6D PE1, PE2, and PE3 spectra.

cm−1 according to our calculations), which is primarily responsible for the splitting. This indicates that the other modes play some role in the progression. D 2 shows no structure but some asymmetry in the experimental spectrum. In our spectrum, the asymmetry is reproduced although it is more pronounced. The two peaks are distinct when only three degrees of freedom are used in 3D PE1 while they are fused with each other as more degrees of freedom 872


Article


Figure 9. Diabatic state populations corresponding to the dynamics of selected (a) 3D, (b−d) 6D, and (e) 10D systems. Shown are the dynamics with the initial wavepacket in D1 (left column), D2 (middle), and D3 (right column).

Figure 10. Adiabatic state populations corresponding to the dynamics of selected (a) 3D, (b−d) 6D, (e) 10D systems. Shown are the dynamics with the initial wavepacket in D1 (left column), D2 (middle), and D3 (right column). For the 6D and 10D cases, the black lines in the adiabatic state populations are the sums of all single populations to estimate the error due to Monte Carlo integration.

simulations. In general, the presence of more modes allows for a

show an increased adiabatic population transfer to D1 in the 6D dynamics with respect to the corresponding 3D dynamics. Moreover, D2 and D3 get depopulated much faster and to a greater extend, and D1 obtains more population in the same

greater distribution of the wavepacket and an easier population transfer. 873


Article

The Journal of Physical Chemistry A Notes

The a″ modes have a more significant effect on the dynamics since they allow transitions between the A′ and A″ states. The comparison of the PE2/PE3 state populations with those of the PE1 simulations reveals the change in the dynamics. The most dramatic effect appears on dynamics starting from D1 where decay to the ground state can now occur. D2 and D3 also get depopulated faster with population moving to all states. In view of the effect of a″ modes, we expect that the 6D PE3 case represents better the dynamics compared to 10D PE1, contrary to what was seen in the PE spectra.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge the support from the Department of Energy under Award Number DE-FG02-08ER15983. S.M. acknowledges support from the Humboldt Foundation through an Experienced Researcher Fellowship during her visit to Germany. Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged (H.K.).

■

■

SUMMARY AND CONCLUSIONS The PE spectrum of uracil and the dynamics of the ionic states produced after ionization to the first four cationic states have been studied with molecular quantum dynamics. The MCTDH method in combination with the vibronic coupling model were used. The PE spectrum resulting from this work is in good agreement with the experimental spectrum, validating the approach. The theoretical analysis also offers a detailed understanding of which modes contribute to the bands and their structure. We have seen that 10 a′ modes were needed for a good reproduction of the spectrum while the a″ modes do not contribute much. On the other hand, the a″ modes are very important for the dynamics since they allow for coupling between neighboring states, which have different symmetries in this case. The results also show that direct nonadiabatic transitions from D2 to D0 and from D3 to D1 contribute significantly to the dynamics, indicating that the three-state conical intersections found previously in uracil cation play a crucial role. The presence of the low energy three-state CI in combination with the small coupling between D0−D1 and D1−D2 are possibly responsible for the big contribution from the direct D2-D0 transitions. The uracil cation, therefore, gives a strong example of the importance of three-state CIs in cases of high density of electronic states. The present study cannot address fragmentation of the radical cation directly since we are using a Hamiltonian appropriate for describing the bound region only. The results however can be used to make indirect arguments. Specifically, the results show that radiationless decay from higher ionic states to D0 occurs very fast, so fragmentation will have to occur in less than 100 fs in order to compete with the radiationless decay. This may be possible to occur in cases where the ionic PES is directly dissociative, but most likely not possible in other cases. In a future publication we will present studies of the dynamics of uracil cation using surface hopping dynamics. These studies will serve as a comparison of the quantum dynamics with limited degrees of freedom to semiclassical dynamics where all degrees of freedom are included. Furthermore, surface hopping studies could potentially address fragmentation, although even in that case one should use the appropriate underlying ab initio method for that to occur.

■

ADDITIONAL NOTE This paper was originally submitted for the “David R. Yarkony Festschrift”, published as the December 26, 2014, issue of J. Phys. Chem. A (Vol. 118, No. 51).

■

(1) Becker, D.; Sevilla, M. D. In The Chemical Consequences of Radiation Damage in DNA; Lett, J. T., Sinclair, W. K., Eds.; Academic Press: San Diego, CA, 1993; Vol. 17, pp 121−180. (2) Steenken, S. Purine Bases, Nucleosides, and Nucleotides: Aqueous Solution Redox Chemistry and Transformation Reactions of their Radical Cations and e− and OH Adducts. Chem. Rev. 1989, 89, 503−520. (3) Allison, T. K.; Tao, H.; Glover, W. J.; Wright, T. W.; Stooke, A. M.; Khurmi, C.; Tilborg, J. v.; Liu, Y.; Falcone, R. W.; Martínez, T. J.; Belkacem, A. Ultrafast internal conversion in ethylene. II. Mechanisms and pathways for quenching and hydrogen elimination. J. Chem. Phys. 2012, 136, 124317. (4) Tao, H.; Shen, L.; Kim, M. H.; Suits, A. G.; Martínez, T. J. Conformationally selective photodissociation dynamics of propanal cation. J. Chem. Phys. 2011, 134, 054313. (5) Joalland, B.; Mori, T.; Martínez, T. J.; Suits, A. G. Photochemical Dynamics of Ethylene Cation C2H+4 . J. Phys. Chem. Lett. 2014, 5, 1467− 1471. (6) Köppel, H.; Döscher, M.; Bâldea, I.; Meyer, H.-D.; Szalay, P. G. Multistate vibronic interactions in the benzene radical cation. II. Quantum dynamical simulations. J. Chem. Phys. 2002, 117, 2657−2671. (7) Worth, G.; Carley, R.; Fielding, H. Using photoelectron spectroscopy to unravel the excited-state dynamics of benzene. Chem. Phys. 2007, 338, 220−227. (8) Mahapatra, S.; Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S.; Köppel, H. The Ã 2E/B̃ 2B2 Photoelectron Bands of Allene beyond the Linear Coupling Scheme: An ab Initio Dynamical Study Including All Fifteen Vibrational Modes. J. Phys. Chem. A 2001, 105, 5567−5576. (9) Saddique, S.; Worth, G. Applying the vibronic coupling model Hamiltonian to the photoelectron spectrum of cyclobutadiene. Chem. Phys. 2006, 329, 99−108. (10) Cattarius, C.; Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. All mode dynamics at the conical intersection of an octa-atomic molecule: Multi-configuration time-dependent Hartree (MCTDH) investigation on the butatriene cation. J. Chem. Phys. 2001, 115, 2088−2100. (11) Trofimov, A. B.; Köppel, H.; Schirmer, J. Vibronic structure of the valence π-photoelectron bands in furan, pyrrole, and thiophene. J. Chem. Phys. 1998, 109, 1025−1040. (12) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The multiconfiguration time-dependent Hartree method: A highly efficient algorithm for propagating wavepackets. Phys. Rep. 2000, 324, 1−105. (13) Multidimensional Quantum Dynamics: MCTDH Theory and Applications; Meyer, H.-D., Gatti, F., Worth, G. A., Eds.; VCH: Weinheim, Germany, 2009. (14) Worth, G. A.; Meyer, H.-D.; Köppel, H.; Cederbaum, L. S.; Burghardt, I. Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics. Int. Rev. Phys. Chem. 2008, 27, 569−606. (15) Worth, G. A.; Beck, M. H.; Jäckle, A.; Burghardt, I.; Lasorne, B.; Meyer, H.-D. The MCTDH Package, Development Version 10.0, 2011.

ASSOCIATED CONTENT

S Supporting Information *

Table giving details of the dynamics simulations setup. This material is available free of charge via the Internet at http://pubs. acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*(M.A.) E-mail: [email protected]. 874


Article

The Journal of Physical Chemistry A (16) Matsika, S.; Spanner, M.; Kotur, M.; Weinacht, T. C. Ultrafast Relaxation Dynamics of Uracil Probed via Strong Field Dissociative Ionization. J. Phys. Chem. A 2013, 117, 12796−12801. (17) Kotur, M.; Zhou, C.; Matsika, S.; Patchkovskii, S.; Spanner, M.; Weinacht, T. C. Neutral-Ionic State Correlations in Strong-Field Molecular Ionization. Phys. Rev. Lett. 2012, 109, 203007. (18) Spanner, M.; Patchkovskii, S.; Zhou, C.; Matsika, S.; Kotur, M.; Weinacht, T. C. Dyson norms in XUV and strong-field ionization of polyatomics: Cytosine and uracil. Phys. Rev. A 2012, 86, 053406. (19) Zhou, C.; Matsika, S.; Kotur, M.; Weinacht, T. C. Fragmentation Pathways in the Uracil Radical Cation. J. Phys. Chem. A 2012, 116, 9217−9227. (20) Matsika, S. Two- and three-state conical intersections in the uracil cation. Chem. Phys. 2008, 349, 356−362. (21) Padva, A.; LeBreton, P. R.; Dinerstein, R. J.; Ridyard, J. N. A. Uv photoelectron studies of biological pyrimidines: The electronic structure of uracil. Biochem. Biophys. Res. Commun. 1974, 60, 1262− 1268. (22) Urano, S.; Yang, X.; LeBreton, P. R. UV photoelectron and quantum mechanical characterization of DNA and RNA bases: valence electronic structures of adenine, 1,9-dimethyl-guanine, 1-methylcytosine, thymine and uracil. J. Mol. Struct. 1989, 214, 315−328. (23) Dougherty, D.; Wittel, K.; Meeks, J.; McGlynn, S. P. Photoelectron spectroscopy of carbonyls. Ureas, uracils, and thymine. J. Am. Chem. Soc. 1976, 98, 3815−3820. (24) O’Donnell, T. J.; LeBreton, P. R.; Petke, J. D.; Shipman, L. L. Ab initio quantum mechanical characterization of the low-lying cation doublet states of uracil. Interpretation of UV and x-ray photoelectron spectra. J. Phys. Chem. 1980, 84, 1975−1982. (25) Köppel, H.; Domcke, W.; Cederbaum, L. S. Multimode Molecular Dynamics Beyond the Born-Oppenheimer Approximation. Adv. Chem. Phys. 1984, 57, 59. (26) Domcke, W.; Yarkony, D. R.; Köppel, H. Conical Intersections: Electronic Structure, Dynamics & Spectroscopy; World Scientific: Singapore, 2004. (27) Domcke, W.; Yarkony, D. R.; Köppel, H. Conical Intersections: Theory, Computation and Experiment; World Scientific: Singapore, 2011. (28) Köppel, H.; Domcke, W.; Cederbaum, L. S. The Multi-Mode Vibronic-Coupling Approach, In Conical Intersections: Electronic Structure, Dynamics and Spectroscopy; Domcke, W., Yarkony, D. R., Köppel, H., Eds.; World Scientific: Singapore, 2004; pp 323−367. (29) Light, J. C.; Hamilton, I. P.; Lill, J. V. Generalized discrete variable approximation in quantum mechanics. J. Chem. Phys. 1985, 82, 1400− 1409. (30) Puzzarini, C.; Biczysko, M.; Barone, V. Accurate Anharmonic Vibrational Frequencies for Uracil: The Performance of Composite Schemes and Hybrid CC/DFT Model. J. Chem. Theory Comput. 2011, 7, 3702−3710. (31) Frisch, M. J. et al. Gaussian 03, Revision D.01, Gaussian, Inc.: Wallingford, CT, 2004. (32) Krylov, A. I. Equation-of-Motion Coupled-Cluster Methods for Open-Shell and Electronically Excited Species: The Hitchhiker’s Guide to Fock Space. Annu. Rev. Phys. Chem. 2008, 59, 433−462. (33) Shao, Y.; et al. Advances in methods and algorithms in a modern quantum chemistry program package. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191.

875


Photoelectron Spectrum and Dynamics of the Uracil Cation - The

Recommend Documents