Ab Initio Quantum Study of Nonadiabatic S1–S2 Photodynamics of s

Jul 8, 2013 - The nonadiabatic photodynamics of s-trans-butadiene in its lowest singlet excited states is studied theoretically, using a fully quantal...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCA

Ab Initio Quantum Study of Nonadiabatic S1−S2 Photodynamics of s‑trans-Butadiene A. Komainda,† B. Ostojić,‡ and H. Köppel*,† †

Physikalisch-Chemisches Institut, Universität Heidelberg, Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany Institute of Chemistry, Technology and Metallurgy, University of Belgrade, Studentski trg 14-16, 11 000 Belgrade, Serbia



ABSTRACT: The nonadiabatic photodynamics of s-trans-butadiene in its lowest singlet excited states is studied theoretically, using a fully quantal approach. The coupled 1Bu and 2Ag states are considered in the calculation, representing the lowest dipole-allowed electronic transition, and the dipole-forbidden state with substantial double-excitation character, respectively. Up to six nuclear degrees of freedom, including out-of-plane dihedral angles, are included. The calculation of the underlying potential energy surfaces relies on the CASPT2 method, where widely different CAS spaces have been compared. The ultrafast electronic population decay is confirmed, proceeding on a time scale of 30− 40 fs. Pronounced out-of-plane distortions are obtained for the first time from fully quantal calculations. The complexity of the electronic absorption spectrum increases substantially upon including additional vibrational modes in the calculation. Further computations were performed to facilitate inclusion of the coupling to the ground state in subsequent work. the respective equilibrium geometry.9,17,18 The crossing and near-degeneracy are indicative of a conical intersection (CoIn) between the corresponding potential energy surfaces (PESs) which has also been characterized theoretically.8,9,19−22 In combination with another intersection, between the 2Ag and 1Ag state PESs,8,10 this gives rise to the ultrafast relaxation dynamics of trans-butadiene and other short polyenes. Experimental manifestations are a substantial broadening of the UV excitation spectrum,23,24 absence of fluorescence,25 and cis−trans isomerization26,27 during the relaxation process. There exist relatively few theoretical studies on the excitedstate dynamics and relaxation behavior of trans-butadiene and related systems. The early trajectory-like calculations of Robb and co-workers focused on the ground-state motion following the 2Ag−1Ag nonadiabatic transition.28 In the later quantum wavepacket study of Krawczyk et al., ground-state normal coordinates were used to construct a linear vibronic coupling (LVC) model for the interacting and 1Bu and 2Ag states.22 By virtue of symmetry selection rules only, planar vibrational modes enter the description, and seven modes were employed to compute absorption and resonance Raman spectra as well as the electronic population transfer. The most recent study to date seems to be the ab initio multiple spawning (AIMS) investigation of Levine and Martinez, which included all three electronic states discussed above and identified several types of conical intersections, partly not addressed before in the literature.10 The role of charge transfer intersections was

I. INTRODUCTION In this work we investigate theoretically the photodynamics of the smallest all-trans polyene, s-trans-1,3-butadiene, hereafter referred to in short as trans-butadiene. Polyenes being representatives and structural motifs of many different biomolecules, a vast amount of experimental and theoretical work has been devoted to clarify the nature of their excited electronic states and their photobiological function.1−5 Even for the smallest systems like butadiene many open questions remain, and we recall here some of the key aspects of this complexity. The lowest excited states differ strongly in their electronic structure,6,7 which makes a balanced theoretical treatment difficult to achieve. At least three electronic states are important to describe the photochemistry of short and midsized polyenes, the 1Ag ground state and the 2Ag and 1Bu singlet excited states.8−10 The covalent 2Ag state is of marked doubleexcitation character and optically dark whereas the 1Bu zwitterionic state carries most of the oscillator strength in the ~6 eV excitation energy range (see, for example refs 6, 11, and 12). Their energetic order has been controversial for a long time in the literature, partly due to difficulties to precisely locate the dipole-forbidden 2Ag state experimentally.13−15 For transbutadiene, the issue has been settled only recently by the EOMCC benchmark study of Watson and Chan, which gave best estimates of 6.21 ± 0.02 eV for the vertical excitation energy of the 1Bu state and 6.39 ± 0.07 eV for that of the 2Ag state.16 Although the dipole-allowed state is thus placed vertically below the dipole-forbidden doubly excited state, due to their energetic proximity, they may well cross in or near the Franck− Condon (FC) zone and their ordering thus be reversed during the nuclear motion. The 2Ag state is placed below the state at © 2013 American Chemical Society

Special Issue: Structure and Dynamics: ESDMC, IACS-2013 Received: May 2, 2013 Revised: July 6, 2013 Published: July 8, 2013 8782

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

emphasized and also a comprehensive account of the literature given. We refer to this paper for further references on the relaxation dynamics of trans-butadiene and other small polyenes. In the present paper we want to extend the scope of previous quantum wavepacket investigations on trans-butadiene by using curvilinear (symmetry-adapted internal) coordinates to construct the vibronic Hamiltonian and by including out-of-plane vibrational degrees of freedom in the analysis. This is motivated by the desire to better understand the mechanisms underlying the experimental absorption spectrum and to clarify the nature of the short-time wavepacket motion as it moves out of the FC zone; the latter was predicted10 to involve out-of-plane (torsional) motion, which is important to verify from a fully quantal calculation. The impact of an increasing number of degrees of freedom on the electronic population dynamics and on the excitation spectrum of trans-butadiene will also be elucidated. We further investigate the structural distortions relevant to induce a 2Ag−1Ag state conical intersection and thus provide the basis for three-state wavepacket calculations, to be undertaken in subsequent work. Regarding the general scope of this work, we note that we do not aim at reproducing Rydberg excited states and will omit diffuse functions in the one-particle basis, specified below. This is motivated by our overall goal of revealing general trends in the polyene series of molecules upon increasing chain length (see also note in the Conclusion). It reduces the density of excited states and facilitates to introduce the S2 and S1 shorthand notation for the 1Bu and 2Ag states, respectively. The latter does not follow the energetic ordering in the FC zone center, where the 1Bu state is below the 2Ag state, but rather the ordering at the C2h-constrained minima (see above). It is also motivated by the fact that, after the excursion of the photoexcited system from the FC zone, the 2Ag state is generally below the 1Bu state and is the lowest excited state prior to the internal conversion to the S0 (1Ag) ground state.10,28 This paper is organized as follows. In section II we detail the theoretical framework adopted, namely the underlying coordinates and the form of the vibronic Hamiltonian, the quantum-dynamical treatment (wavepacket propagation scheme) and the electronic-structure methodology. Section III is devoted to the results obtained for the potential energy surfaces and coupling elements, the time-dependent quantities like electronic populations and the spectral intensity distribution. Section IV concludes.

Table 1. Definition of the Symmetry-Adapted Internal Coordinates of trans-Butadiene in-plane S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17

S18 S19 S20 S21 S22 S23 S24

Δτ7,2,3,8 + Δτ1,2,3,4 + Δτ7,2,3,4 + Δτ1,2,3,8 Δτ5,1,2,3 + Δτ6,1,2,7 + Δτ5,1,2,7 + Δτ6,1,2,3 + Δτ2,3,4,9 + Δτ8,3,4,10 + Δτ2,3,4,10 + Δτ8,3,4,9 Δτ5,1,2,3 + Δτ6,1,2,7 + Δτ5,1,2,7 + Δτ6,1,2,3 − Δτ2,3,4,9 − Δτ8,3,4,10 − Δτ2,3,4,10 − Δτ8,3,4,9 Δγ10,9,3,4 + Δγ6,5,2,1 Δγ10,9,3,4 − Δγ6,5,2,1 Δγ8,2,4,3 + Δγ7,3,1,2 Δγ8,2,4,3 − Δγ7,3,1,2

Figure 1. Definition of atom numbering, bond lengths and angles, and torsional angles for trans-butadiene.

the most relevant one. It comprises the coordinates S3 (ag, symmetric stretching of the terminal CC bonds), S13 (ag, symmetric CCC bending), S17 (bu, antisymmetric C CC bending), S19 (au, disrotatory twisting of the CH2 groups around the terminal CC bonds), and S20 (bg, conrotatory twisting of the CH2 groups around the terminal CC bonds). The other coordinates were found to be less important in the cross sections considered. Here we follow these earlier findings but augment the coordinate space by S5 (ag, symmetric stretching of the central CC bond). The energy profiles were represented earlier by sums of elementary transcendental functions, which provide a faithful representation along the cross sections in question. These were used with slight modifications owing to different CAS spaces in the present work (see next subsections). To define the vibronic Hamiltonian, the kinetic energy operator needs to be specified. Because we employ a (quasi-) diabatic electronic basis, the residual derivative couplings in this basis are neglected and the usual kinetic energy operator for noninteracting states is used. Here we employ the familiar G-

II. THEORETICAL METHODOLOGY At its ground-state equilibrium geometry, trans-butadiene has C2h symmetry, and its 24 vibrational modes transform as Γvib = 9a g + 4a u + 3bg + 8bu

out-of-plane

Δr4,10 + Δr1,6 Δr4,9 + Δr1,5 Δr4,3 + Δr1,2 Δr3,8 + Δr2,7 Δr2,3 Δr4,10 − Δr1,6 Δr4,9 − Δr1,5 Δr4,3 − Δr1,2 Δr3,8 − Δr2,7 Δϕ3,4,10 + Δϕ2,1,6 Δϕ3,4,9 + Δϕ2,1,5 Δϕ4,3,8 + Δϕ1,2,7 Δϕ4,3,2 + Δϕ3,2,1 Δϕ3,4,10 − Δϕ2,1,6 Δϕ3,4,9 − Δϕ2,1,5 Δϕ4,3,8 − Δϕ1,2,7 Δϕ4,3,2 − Δϕ3,2,1

(1)

The au vibrational modes reduce the molecular symmetry to C2, whereas the bg and bu modes reduce it to Ci and Cs, respectively. A. Vibronic Hamiltonian. The vibronic Hamiltonian adopted follows the general lines of ref 9. The potential energy matrix is written in a (quasi-)diabatic electronic basis employing symmetry-adapted internal coordinates as listed in Table 1 and depicted in Figure 1. Cuts through the adiabatic PES for all degrees of freedom, using the CASPT2 method, have been computed in ref 9 and a subset has been identified as 8783

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

matrix technique of Wilson29,30 and derive the G-matrix elements by symmetry adaptation from Tabs. III−IV of ref 31, where they are tabulated for the localized internal coordinates as given in Table 1. The results are available upon request. Thus we arrive at the total Hamiltonian adopted in this work: H = TN1 + W(Q)

The solution of the equations of motion requires the frequent computation of the mean fields. The efficiency of the MCTDH method thus demands their fast evaluation and necessitates avoiding the explicit calculation of high-dimensional integrals. We emphasize that the present form of the Hamiltonian represents a sum of elementary functions of the Qi that is exactly the form that makes the application of the MCTDH algorithm efficient. Given the time-dependent wavepacket (5) various timedependent and time-independent quantities can be computed directly or indirectly. Electronic spectra can be obtained directly from a time-dependent treatment as the Fourier transform of the autocorrelation function C(t), assuming a direct transition from the initial to the final states within the framework of Fermi’s golden rule:38,39

(2)

with 2TN = p TGp

(3)

and ⎛Wa(Q ) Wab(Q )⎞ ⎟ W(Q) = ⎜⎜ ⎟ ⎝Wab(Q ) Wb(Q ) ⎠

(4)

P(E) ∝

Here p denotes the vector of momenta conjugate to the symmetry coordinates of Table 1 and 1 stands for the 2 × 2 unit matrix. The meaning of the other quantities is obvious or will be specified below. Note that the matrices 1 and W refer to electronic function space, whereas within TN boldface notation refers to nuclear coordinate space. B. Quantum Dynamical Treatment. The dynamical calculations performed in this work and reported below rely on fully quantal, time-dependent methods, namely wavepacket propagation techniques. The multiconfiguration time-dependent Hartree (MCTDH) method32−37 uses a time development of the wave function expanded in a basis of sets of variationally optimized time-dependent functions φ(κ) called single-particle functions (SPFs). The MCTDH equations of motion are obtained from the Dirac−Frenkel variational principle. By virtue of this optimization, the length of this expansion can be much smaller than in standard integration schemes (so-called MCTDH contraction effect). The efficiency is even enhanced by two important additional features: each of the coordinates used in the integration scheme (called particles) can comprise several physical coordinates Qi. Furthermore, for vibronically coupled systems the wave function is written as a sum of several wave functions which is called multiset formulation,

C(t ) = ⟨Ψ(0)|Ψ(t )⟩ = ⟨0|τ †e−iHt τ |0⟩ = ⟨Ψ(t /2)*|Ψ(t /2)⟩

(5)

The SPFs may then be optimized separately for each electronic state α, and therefore fewer coefficients are needed in the wave function expansion. Both choices are employed in this work and, in combination, lead to an MCTDH contraction35 effect of about 6 orders of magnitude (Table 2). This is of the same order as in other typical applications and makes the quantum dynamical calculations feasible at all. Table 2. Number of Basis Functions for the Primitive as Well as the Time-Dependent (SPF) Basis Used in the MCTDH Calculations for the 1Bu and 2Ag Excited Electronic States of Butadienea

a

modes

primitive basis

SPF basis

(S3, S20) (S13, S19) (S5, S17)

(60, 70) (70, 60) (70, 70)

[35, 35] [30, 30] [30, 30]

(7) (8)

In eq 7, |0⟩ is the vibrational ground state of the initial electronic state (the 1Ag ground state in our case) and τ† is the vector of individual transition matrix elements τα between the initial state and the final electronic states labeled by α. The dependence of the transition dipole matrix elements on the nuclear coordinates is neglected. The ground-state wave function is obtained by a relaxation scheme in imaginary time35,40 using the ground-state PES (computed as described below). This wave function is subsequently “lifted” to the 1Bu excited state. The autocorrelation function C(t) measures the overlap between this time-evolving wave packet and the initial one, and its Fourier transform gives the corresponding spectrum according to eq 6. The scalar product involving the vector τ of transition matrix elements in general implies a summation over various partial spectra, each being proportional to |τα|2 (different final electronic states). However, because one of the electronic transitions is dipole-forbidden in the C2h symmetry group, only the contribution from the 1Bu ← 1Ag transition remains, whereas the other corresponds to the 2Ag ← 1Ag transition and vanishes. Equation 8 is valid here because our Hamiltonian is Hermitian and the initial wavepacket is real41,42 and allows us to reduce the propagation time by a factor of 2. Due to the finite propagation time T of the wavepacket, the Fourier transformation causes artifacts known as the Gibbs phenomenon.43 To reduce this effect, the autocorrelation function is first multiplied by a damping function cos2(πt/ 2T).35,44 Furthermore, to simulate the experimental line broadening, the autocorrelation functions are damped by an additional multiplication with a Gaussian function exp[−(t/ τd)2], where τd is the damping parameter, also called dephasing time below. This multiplication is equivalent to a convolution of the spectrum with a Gaussian with a full width at halfmaximum (fwhm) of 4(ln 2)1/2/τd. The convolution thus simulates the experimental spectrometer resolution used in experiments, plus intrinsic line broadening effects. The two other quantities of interest are the time-evolving (diabatic) electronic populations, Pα(t), and two-dimensional reduced densities ρα(Qi,Qj,t) for the electronic state α. These

∑ Ψα(t )|α⟩ α

(6)

with

ns

Ψ(t ) =

∫ eiEt C(t ) dt

For all modes harmonic primitive basis functions were used. 8784

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

of 0.2 Eh and energies obtained in the RS2C calculations were corrected for this shift. Three series of active spaces involved in the CASSCF procedure on the ground and first two excited states of butadiene were employed. We refer here to orbitals that are symmetric with respect to reflection in the plane of the molecule as σ and those that are antisymmetric upon reflection as π. Cave and Davidson56 first pointed out that the σ−π correlation effects can be important for the description of the 1Bu state. The first active space denoted by CAS(4,4) involved four active orbitals, two valence bonding π molecular orbitals, and two valence antibonding π* molecular orbitals. The second active space denoted by CAS(8,8) involved two valence bonding π orbitals, two valence antibonding π* orbitals, two highest-lying occupied σ orbitals at the equilibrium geometry of the ground electronic state, and two valence antibonding σ* orbitals. The third active space denoted by CAS(10,7) (ten electrons are distributed among seven active orbitals) involved two valence bonding π orbitals, two valence antibonding π* orbitals and three additional σ orbitals that are the highest-lying occupied σ orbitals at the equilibrium geometry of the ground electronic state. It is related to the CAS space used in earlier work.9 The basis set used was the correlation-consistent polarized core−valence triple-ζ basis set (cc-pCVTZ) on carbon57,58 and the correlation-consistent polarized valence triple-ζ basis set (cc-pVTZ) on hydrogen.57 In the case of calculating the vertical energy differences on the equilibrium geometry of the ground electronic state, a modified augmented correlation-consistent polarized core−valence triple-ζ basis set on carbon was used (aug′-cc-pCVTZ) too. The diffuse functions derived from ref 59 corresponding to 3s Rydberg atomic orbitals with the exponents 0.04370, 0.01725, 0.01045, and 0.004125 were centered on carbon and added to the cc-pCVTZ basis set (aug′cc-pCVTZ). The vertical excitation energies (VE) were calculated at our computed equilibrium geometry of the ground electronic state. For the 1Bu and 2Ag states the VE values obtained at the CAS(8,8)-RS2C level of theory were 6.13 and 6.42 eV, respectively. The VE values calculated at the same level of theory and employing the aug′-cc-pCVTZ basis set on carbon and cc-pVTZ basis set on hydrogen were VE(1Bu) = 6.12 eV and VE(2Ag) = 6.44 eV. It is noteworthy that these results are close to each other and to the benchmark values of ref 16. The potential energy surfaces have been calculated along cross sections for the relevant coordinates. They were investigated along internal symmetry coordinates and the reference geometry was the equilibrium geometry of the molecule in its ground electronic state. The set of internal coordinates chosen for the description of the distortions from the reference geometry was the same as used in our earlier paper.9 It is the recommended set of valence coordinates defined in the paper of Pulay et al.60 The symmetry-adapted internal coordinates are defined as suitable linear combinations of these internal coordinates (Table 1; see also Figure 1). The computations have been carried out in the framework of four point groups: C2h, C2, Ci, and Cs (same axis choice as before). Within the framework of the C2, Ci, and Cs point groups the 1Ag, 2Ag, and 1Bu electronic states correlate with the 1A, 2A, 1B; 1Ag, 2Ag, 1Au; and 1A′, 2A′, 3A′ electronic states, respectively. The potential energy curves as a function of the S3, S5, S13, S19, and S20 coordinates have been calculated at the CAS(4,4)-RS2C and CAS(10,7)-RS2C level of theory. For the

quantities are defined as follows, using the wave function given by eq 5: Pα(t ) = ⟨Ψα(t )|Ψα(t )⟩ ρα (Q i ,Q j ,t ) =

(9)

∫ Ψ*α(t ) Ψα(t ) ∏

dQ l (10)

l≠i ,j

In complete analogy, also one-dimensional reduced densities can be obtained and will be presented below. The timedependent wavepacket (WP) calculations have been performed using the Heidelberg MCTDH package.45 C. Electronic Structure Calculations. 1. Ground-State Equilibrium Geometry. As a starting point of the calculations the optimization of the electronic ground state, 1Ag, was performed at the MP2/aug-cc-pVTZ level of theory. The C2h symmetry was imposed during the geometry optimization and the molecule was oriented to lie in the xy-plane so that the zaxis associated with the C2 axis was perpendicular to the molecular plane. Table 3 compares the calculated C−C and C− Table 3. Parameters of the Ground-State Equilibrium Geometry of trans-Butadiene

a

bond length/bond angle

expa

expb

thc

semiexpd

this work

rC1C2 rC2C3 φC1−C2−C3 rC2H7 φC1−C2H7 rC1H5 φC2−C1H5 rC1H6 φC2C1H6

1.348 1.468 124.3 1.107 120.7 1.107 120.7 1.107 120.7

1.337 1.467 123.5 1.089 119.6 1.085 121.5 1.087 121.6

1.3377 1.4548 123.5 1.0846 119.9 1.0823 121.0 1.0799 121.5

1.3376 1.4539 123.62 1.0847 119.91 1.0819 120.97 1.0793 121.47

1.341 1.453 123.6 1.085 119.5 1.083 120.9 1.080 121.4

In ref 46. bIn ref 47. cIn ref 49. dIn ref 48.

H bond lengths and bond angles and contains the experimental,46,47 semiexperimental equilibrium structure,48 and ab initio values based on complete basis set extrapolations of the frozen core CCSD(T) parameters using up to the aug-ccpVQZ basis set with core/valence and Douglas−Kroll relativistic corrections applied.49 Craig et al. determined the equilibrium structure of trans-1,3-butadiene after adjusting the rotational constants obtained from rotational spectroscopy by vibration−rotation constants calculated from the results of quantum chemical calculations.48 This is denoted as “semiexp” in Table 3 and agrees rather well, mostly very well, with our results (last column in the table). 2. Calculation of the Potential Energy Surfaces for the 11Ag, 21Ag, and 11Bu Electronic States. We have computed the potential energy surfaces of the 1Ag, 2Ag, and 1Bu electronic states of by employing the complete active space (CASSCF) technique50,51 followed by a modified version of CASPT2 (complete active space with second-order perturbation theory) developed by Celani and Werner52 and referred to as RS2C. The potential energy curves for the S17 symmetry-adapted internal coordinates were calculated using the CASSCF method followed by a multireference configuration interaction (MRCI) treatment.53−55 For calculations of electronic states we used the CASSCF state-averaging (SA) procedure; i.e., the three states were included with equal weights. To avoid intruder-state problems in the RS2C calculations, we applied an energy shift 8785

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

S5 coordinate we performed additional calculations at the CAS(8,8)-RS2C level. In the case of bu vibrations, the molecule possesses Cs symmetry and the ground and both of the excited states are of A′ symmetry. The calculations for the coupling between the 2Ag and 1Bu electronic states via the bu mode (S17) coordinate have been performed at the CASSCF-MRCI level of theory using the CAS(8,8) orbital space. The CI expansion was generated within the internally contracted method using single and double substitutions (MRCISD) from each reference determinant. In these MRCISD calculations eight electrons were correlated and the effect of higher excitations was taken into account by using the Davidson correction61 (denoted as CAS(8,8)-MRCISD+Q). All electronic structure calculations were carried out using the MOLPRO 2008.162 suite of programs.

III. RESULTS AND DISCUSSION A. Potential Energy Matrix for the S1 and S2 States. As pointed out in the previous section, we obtain the CASPT2

Figure 3. Same as Figure 2, but along S13.

longer bond lengths (similar for the minima of both potentials). Even more evident is the different curvature of these potentials. The CAS(10,7) space yields much steeper potentials, whereas the CAS(4,4)-based calculations lead to very flat and broad potentials. Of course, the different curvature of the potentials can result in a quite different dynamical behavior. Besides the mode S3, also the two results for the S13-mode (Figure 3) show some (minor) differences. The CAS(4,4) results exhibit an intersection at S13 = −20°. The CAS(10,7)-based calculations lead to a steeper potential for the 1Bu state and a small shift of the intersection point to a C−C−C angle closer to the equilibrium value of the electronic ground state. Generally, the CASPT2 method based on the CAS(10,7) active space appears to be most appropriate for extended calculations of PESs for this set of electronic states of trans-butadiene. This finding is supported, as mentioned before, by the work of Davidson and Cave,56 emphasizing the importance of the inclusion of σ−π correlations to achieve a better description of the 1Bu state. The CAS(10,7) active space incorporates this correlation, whereas the CAS(4,4) space contains only π-orbitals in the active space. For the treatment of the vibronic interaction between the 1Bu and 2Ag states we need to construct a model Hamiltonian, i.e., specify the matrix elements of eq 4. There, the diagonal elements Wa and Wb represent the PES of the 2Ag and 1Bu excited states, respectively. They are taken to be the sums of the cross sections through the potential profiles according to eqs 11 and 12.

Figure 2. Potentials along the S3-coordinate for the 1Bu (green line) and 2Ag (red line) excited states of trans-butadiene. The upper panel presents the results of the CAS(10,7) calculations; the lower one present those obtained with CAS(4,4). The full lines always represent the fit results of eqs 13, 19, 25, and 31; the circles denote the ab initio data points.

potentials for the modes S3, S5, S13, S19, and S20 based on either the CAS(4,4) or CAS(10,7) wave function. The results with both active spaces agree qualitatively quite well, although it is noteworthy that there are some differences. Especially for the totally symmetric mode S3 we find very different potential curves (Figure 2). For the smaller CAS(4,4) space the intersection of the 1Bu and the 2Ag state PESs is shifted to

Wa = Va(S3) + Va(S5) + Va(S13) + Va(S17) + Va(S19) + Va(S20)

(11)

Wb = Vb(S3) + Vb(S5) + Vb(S13) + Vb(S17) + Vb(S19) + Vb(S20) 8786

(12)

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

Va(S19) = 10.3414 − 7.0386 cos(S19) + 3.2765 cos(2S19) − 0.3173 cos(3S19)

(17)

Va(S20) = 19.4878 − 21.7024 cos(S20) + 10.4439 cos(2S20) − 1.9681 cos(3S20) Vb(S3) = 62.3849 − 123.9020e−S3 + 67.5354e−2S3

(18) (19)

Vb(S5) = 11.2529 − 6.5670e−S5 − 2.3166e−2S5 + 4.7233e−3S5 − 1.0050e−4S5

(20)

Vb(S13) = 6.0334 − 0.0073S13 + 0.0034S132

(21)

Vb(S17) = 6.6068 + 0.0019S17 2

(22)

Vb(S19) = 5.2402 + 2.0156 cos(S19) − 2.1439 cos(2S19) + 0.9545 cos(3S19)

(23)

Vb(S20) = 6.5056 − 0.0877 cos(S20) − 1.5208 cos(2S20) + 1.1284 cos(3S20)

(24)

Examples for the quality of the fits can be found in Figures 2 and 3. The results of the CAS(4,4) calculations are well fitted by the following functions: Va(S3) = 6.4866 − 7.5928S3 + 18.1944S32 − 12.5221S33 + 4.0783S3 4

(25)

Va(S5) = 14.4882 − 16.9059e−S5 + 11.3050e−2S5 − 2.0968e−3S5 + 0.5242e−4S5

(26)

Va(S13) = 6.4828 + 0.0054S13 + 0.0026S132

(27)

Va(S17) = 6.6557 + 0.0039S17 2

(28)

Va(S19) = 18.9907 − 21.1956 cos(S19) Figure 4. Time-dependent populations of the 1Bu state following a vertical excitation to this state. The upper panel displays the diabatic results obtained with the CAS(10,7) PES, on the middle one those from the CAS(4,4) PES. The lower panel contains the adiabatic populations for the CAS(10,7); they represent the lower adiabatic state, which correlates with the 1Bu state at the FC zone center. In all panels the green line represents the six-mode case and the black one the four-mode case.

+ 11.4403 cos(2S19) − 3.3384 cos(3S19) + 0.5834 cos(4S19)

(29)

Va(S20) = 11.3987 − 8.0316 cos(S20) + 2.9172 cos(2S20) + 0.5348 cos(3S20) − 0.3362 cos(4S20)

(30)

Vb(S3) = 6.0238 − 4.6628S3 + 18.4803S32 − 15.6061S33

The single-mode potentials for the CAS(10,7) results are well described by the following expressions (all energies in electronvolts, bond lengths in bohr, and bond angles in degrees): Va(S3) = 48.2334 − 99.9906e−S3 + 58.0371e−2S3

+ 5.7281S3 4

Vb(S5) = 12.9105 − 14.1163e−S5 + 8.9141e−2S5 − 2.0968e−3S5 + 0.4103e−4S5

(13)

Va(S13) = 6.2726 + 0.0038S13 + 0.0026S13

Va(S17) = 6.6557 + 0.0039S17 2

2

(32)

Vb(S13) = 6.0148 − 0.0022S13 + 0.0031S132

(33)

(14)

Vb(S17) = 6.6068 + 0.0019S17 2

(34)

(15)

Vb(S19) = 10.6899 − 6.6752 cos(S19) + 2.2209 cos(2S19)

Va(S5) = 11.7304 − 6.9994e−S5 − 1.9783e−2S5 + 4.5861e−3S5 − 0.9594e−4S5

(31)

− 0.3392 cos(3S19) + 0.1283 cos(4S19)

(16) 8787

(35)

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

Figure 5. 1D-densities along the S3-coordinate [CAS(10,7)-results]. On the left side of the figure the densities of the 1Bu state are presented (upper panel for t = 10 fs and lower panel for t = 18 fs). The right side shows the corresponding results for the 2Ag state.

Figure 6. 1D-densities along the S3-coordinate [CAS(4,4)-results]. On the left side of the figure the densities of the 1Bu state are presented (upper panel for t = 13 fs and lower panel for t = 20 fs). The right side shows the corresponding results for the 2Ag state.

absorption spectra and also seems to play a role for the location of the S0−S1 intersection. Because all potentials were fitted to the obtained ab initio results, every single-mode potential contains the vertical excitation energy. Accordingly, for the following dynamical calculations it is necessary to subtract the excess energy.

Vb(S20) = 9.0110 − 4.0581 cos(S20) + 0.6323 cos(2S20) + 0.4368 cos(3S20)

(36)

Modes S3 and S13 serve as total-symmetric tuning modes; although displacement along S5 alone does not lead to a CoIn (see Figure 3 of ref 9), it is expected to be important for the 8788

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

Figure 7. 2D-wave packet along the modes S19 and S20. The potentials (obtained with CAS(10,7)) are plotted as dashed lines, and the wave-packets as full lines. The left side of the figure shows the 1Bu state and the right one the 2Ag state.

As previously mentioned, only mode S17 leads to a significant coupling between the excited states. Following earlier work, the off-diagonal element Wab is taken to be a linear function of S17, Wab = λS17

population transfer of the 1Bu state to the dark 2Ag state, proceeding on a time-scale of 30−40 fs. The transfer is rather incomplete, amounting to 50−70%, depending on the details of the calculation. It should be born in mind that the excess energies in question are moderate only because the ground state is not included in the calculation: relative to the 1Bu state, which is the lower state in the FC zone center, it amounts to approximately 1 eV within our Hamiltonian adopted. Figure 4 compares the result with all six modes and the four planar modes in each panel, and the two panels give evidence of the influence of the CAS space on the population dynamics. Not surprisingly, the fluctuations are generally larger with four than with six vibrational modes. This is in line with earlier work showing that a larger nuclear conformational space tends to smooth these quantities.63,64 Nevertheless, given the rather small excess energies, the substantial change when going from four to six modes seems remarkable. It is probably related to the big displacements along the out-of-plane modes that occur; see below. Also, the influence of a bigger CAS space is much stronger in the four-mode than in the six-mode calculation. The reasons for this behavior are not entirely clear at present. We also computed the adiabatic electronic populations (see lowest panel in Figure 4). At time t = 0 the system is located mostly on the lower adiabatic state, which is due to the fact that the 1Bu state is below the 2Ag state at the center of the FC zone (Figures 2 and 3). The calculation reveals that also for longer times the system evolves mostly on the lower adiabatic state. This is in line with other systems with a similar location of the CI with respect to the FC zone center65,66 and can be understood with the aid of phase-space arguments.67,68 It is also consistent with the fact that for longer times the WP evolves to the right of the crossing point in Figure 2 (see next paragraph), where the lower adiabatic state correlates with the 2Ag state.

(37)

The coupling parameter λ is determined by the working equation for the adiabatic potentials V±:38 V± =

Va(S17) + Vb(S17) 2 ±

⎛ Va(S17) + Vb(S17) ⎞2 2 2 ⎜ ⎟ − λ S17 ⎝ ⎠ 2

(38)

The CAS(8,8)-MRCISD+Q calculations for this mode lead to a coupling constant λ = 0.025 eV/deg. To explore the influence of the out-of-plane modes S19 and S20, we consider two cases for our quantum dynamical calculations. On the one hand we use the full coupling matrix as described in eqs 4, 11, and 12. And on the other hand we drop S19 and S20 and limit the Hamiltonian to the four in-plane modes. Below we will present the results of the full six-mode case and the reduced four-mode case. B. Time-Dependent Electronic Populations and Nuclear Densities. Using the model Hamiltonian as described above, we have applied the MCTDH algorithm to solve the time-dependent Schrödinger equation. Some basis set details are collected in Table 2. A primitive basis set size of 8.6 × 1010 leads to 63 000 time-dependent SPFs, amounting to the aforementioned contraction effect. For a propagation time of 500 fs, a CPU time of nearly 170 h on an AMD dual-core Opteron processor (2.6 GHz, 10 GB memory) is required. This leads to the time-dependent population of the optically bright, diabatic 1Bu state, following a vertical transition to this state, as depicted in Figure 4. Figure 4 gives evidence of the ultrafast 8789

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

Figure 9. Same as Figure 8 (without insets), but for the CAS(4,4) PES.

Figure 8. Calculated absorption spectra obtained from the CAS(10,7) PES. In the upper panel the six-mode case is displayed, together with the experimental one23 in the inset, whereas the lower panel shows the absorption spectra of the four-mode case and the reduced S3/S17 spectrum in the inset. In both panels the full lines represent the lowresolution spectra, the dashed lines the high-resolution spectra with a dephasing time as given in the top right of the panels.

excited states in Figure 7, where they are superimposed on WP densities for representative propagation times. Figure 7 demonstrates the large distortions and complex behavior of the WP for both out-of-plane modes. The large distortions are caused by the flat and wide double-mimimum potentials for both modes, leading to the 4-fold minima in the 2D contour plots of the figure. The two left panels display the WP and PES for the 1Bu state, where it is mostly located for these short times (Figure 4); the two right panels are for the 2Ag state according to the then-dominant population of this state. Although for the initial electronic state the geometric distortions and wavepacket complexity are still rather moderate, both increase quite dramatically for later time when the WP is located mostly on the dark 2Ag state. This helps to understand the almost smooth, nonoscillatory behavior of the electronic populations when the out-of-plane modes are included (Figure 4). It also bears on the irregularity of the electronic absorption spectra, to be discussed in the next subsection. C. Vibronic Structure of the Absorption Spectrum. The absorption spectrum of trans-butadiene in the 6 eV energy range represents a major source of experimental information on the photodynamics of the system. The spectra of Vaida et al.23 show a substantial vibronic progression with a dominant spacing of 1400 cm−1. In the jet spectra the 580 cm−1 line width observed in the gas phase is reduced to 490 cm−1, which should be mostly due to homogeneous broadening. This amounts to a

Figures 5−7 give snapshots of the time-dependent wavepacket “accompanying” this population dynamics. For the 1D densities along S3 (Figures 5 and 6) we again document the impact of different CAS spaces on the WP motion. Owing to the different PES along S3, also the WP displacements differ markedly, being much larger for the CAS(4,4) than for the CAS(10,7) approach. In addition, the motion is more pronounced in the 2Ag than in the initial 1Bu state and also displays a splitting of the probability density in this state. The former is due to the stronger shift of the 2Ag equilibrium geometry relative to the 1Ag ground state (irrespective of the CAS space, Figure 2), the latter due to the possibility that the WP can take two different pathways when undergoing the 1Bu− 2Ag nonradiative transition. The coordinate S3 is the most important totally symmetric mode, because of its combined FC and tuning activity. That is, it exhibits a significant shift of the equilibrium geometry in the optically bright state and also modulates the S1−S2 energy gap, leading to a crossing (=cut through a conical intersection) and reversed-state ordering in Figure 2. Two further important modes are the out-of-plane coordinates S19 and S20. The potential energy profiles are given as contour plots for both 8790

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

represents a major step toward a fully microscopic description of the absorption spectrum of trans-butadiene. As before, the importance of the CAS space is also documented for the absorption spectrum, by the results displayed in Figure 9. This shows the same comparison regarding dimensionality and resolution as in Figure 8, but now for the CAS(4,4) potential energy data. The impact of the outof-plane modes is similar as before, but in detail the spectral structures differ markedly from the larger CAS space (Figure 8). This is partly caused by the different potential energy curves along S3 (Figure 2), which is the dominant mode according to the inset in Figure 8 (lower panel). The excitation of S3 is thus weaker in the larger CAS space than in the smaller one (cf. Figure 8 with Figure 9), which improves the agreement with experiment.23 D. Remark on Coupling to the S0 State. Another scope of the present work is to provide the basis for including the coupling to the S0 ground state in future quantum-dynamical calculations. This requires to include the S0−S1 conical intersection in the nuclear coordinate space covered by the quantum treatment. When the present modeling of eqs 11−36 is used for the various diabatic matrix elements; no degeneracy or near-degeneracy of the S0−S1 PES is found (rather, their energetic separation always exceeds ∼2 eV). To provide the basis for a more extended modeling, we have followed the results of refs 10 and 28 and computed 2D S0−S1 PES along the out-of-plane angles S19 and S20, but with other coordinates suitably displaced. In the results presented in Figure 10, S20 is fixed at 60° and S19 is varied. The S3 coordinate is displaced by 0.1 Å and furthermore the antisymmetric stretching mode S8 (Table 1) is activated. The symmetric bending mode S13 in this structure is smaller than at equilibrium, whereas the antisymmetric bending mode is bigger. Additionally the C−C torsional coordinate S18 is increased by 63.6°; see the inset in the figure. This is reminiscent of the transoid structure of ref 10. As inferred from Figure 10, we used the CAS(4,4) and the CAS(10,7) approach in this study. It is seen that indeed the two PES approach each other very closely for the displacement ranges adopted. The CAS(4,4) case essentially reproduces the results of Martinez et al.;10 the CAS(10,7) data provide additional confirmation within the larger CAS space adopted. A similar behavior is also found for other values of S20. However, the crossing in the upper panel must be considered accidental and only confirms the proximity of the S0 and S1 states. Generally, the two PES cannot cross upon varying a single coordinate (noncrossing rule). Note, however, that the noncrossing rule is relaxed in higher-dimensional coordinate space and, generally, two coordinates need to be varied to lead to an energetic degeneracy of two states of the same symmetry. The results show that only two additional coordinates, namely S8 and S18 need to be included to describe the coupling between the S0 and S1 states. The treatment of these two coordinates will be the scope of our future investigation.

Figure 10. Potential along the S19 coordinate for fixed displacements of S3, S5, S8, S13, S17, S18, and S20 as given in the figure. In the upper panel the CAS(10,7) treatment is used and in the lower one CAS(4,4). In both cases the blue line stands for the 1Ag ground state and the red line represents the 2Ag state.

lifetime of 10 fs, in the same range as inferred above from the electronic populations.18,27 In Figures 8 and 9 we present the absorption spectra obtained theoretically at various levels of treatment. In Figure 8 the six-mode (upper panel) and four-mode (lower-panel) spectra are compared as resulting from the CAS(10,7) scheme. Lower-resolution (full lines) and higher-resolution (dashed lines) results are compared in each panel. They amount to dephasing times as given there in the top right. The lowerresolution spectra may be compared to the experimental recording23 given in the inset to the upper panel. They are seen to reproduce it very satisfactorily. This holds, in particular, for the six-mode calculation where a phenomenological broadening with a dephasing time of 30 fs has been applied. The four-mode calculation, on the other hand, requires a dephasing time of 15 fs to achieve similar agreement with experiment. The reason for the different phenomenological broadening becomes evident from the high-resolution spectra (dashed lines), which exhibit a much more complex vibronic structure in the six-mode than in the four-mode calculation. This is necessarily due to the inclusion of the out-of-plane modes S19 and S20, which are discarded in the four-mode treatment. Apparently, these modes are strongly excited in the 1Bu−1Ag transition (see dashed lines in the upper panels of Figures 8 and 9); though not resolved, they represent a major cause for its broad absorption features. An earlier theoretical study, considering only in-plane modes, obtained agreement not only similar with experiment but also based on a short phenomenological dephasing time of 15 fs.22 Apparently, the inclusion of out-of-plane modes substantially reduces the need for phenomenological broadening and thus

IV. CONCLUSIONS In this work we have presented an ab initio quantum dynamical study of the nonadiabatic photodynamics of the smallest alltrans polyene. Accurate CASPT2 potential energy data and multidimensional MCTDH wavepacket calculations have been combined to arrive at an extended, fully quantal description of the system’s motion after the dipole-allowed transition from the 1Ag ground state to the 1Bu excited state. The system is seen to undergo an ultrafast internal conversion process to the doubly 8791

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

(13) Doering, J. P.; McDiarmid, R. J. Chem. Phys. 1980, 73, 3617− 3624. (14) Doering, J. P.; McDiarmid, R. J. Chem. Phys. 1981, 75, 2477− 2478. (15) Chadwick, R. R.; Gerrity, D. P.; Hudson, B. S. Chem. Phys. Lett. 1985, 115, 24−28. (16) Watson, M. A.; Chan, G. K.-L. J. Chem. Theory Comput. 2012, 8, 4013−4018. (17) Aoyagi, M.; Osamura, Y.; Iwata, S. J. Chem. Phys. 1985, 83, 1140−1148. (18) Ito, M.; Ohmine, I. J. Chem. Phys. 1997, 106, 3159−3173. (19) Graham, R. L.; Freed, K. F. J. Chem. Phys. 1992, 96, 1304−1316. (20) Dallos, M.; Lischka, H. Theor. Chem. Acc. 2004, 112, 16−26. (21) Cave, R. J.; Davidson, E. R. Chem. Phys. Lett. 1988, 148, 190− 196. (22) Krawczyk, R. P.; Malsch, K.; Hohlneicher, G.; Gillen, R. C.; Domcke, W. Chem. Phys. Lett. 2000, 320, 535−541. (23) Leopold, D. G.; Pendley, R. D.; Roebber, J. L.; Hemley, R. J.; Vaida, V. J. Chem. Phys. 1984, 81, 4218−4229. (24) Vaida, V. Acc. Chem. Res. 1986, 19, 114−120. (25) Zerbetto, F.; Zgierski, M. Z. J. Chem. Phys. 1990, 93, 1235− 1245. (26) Dou, Y.; Torralva, B. R.; Allen, R. E. J. Phys. Chem. A 2003, 107, 8817−8824. (27) Fuß, W.; Schmid, W. E.; Trushin, S. A. Chem. Phys. Lett. 2001, 342, 91−98. (28) Olivucci, M.; Ragazos, I.; Bernardi, F.; Robb, M. A. J. Am. Chem. Soc. 1993, 115, 3710−3721. (29) Wilson, E. B. J. Chem. Phys. 1939, 7, 1047−1052. (30) Wilson, E. B.; Decius, J.; Cross, P. C. Molecular Vibrations - The Theory of Infrared and Raman Vibrational Spectra; Dover Publications: New York, 1980. (31) Frederick, J. H.; Woywood, C. J. Chem. Phys. 1999, 111, 7255− 7271. (32) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. Chem. Phys. Lett. 1990, 165, 73−78. (33) Manthe, U.; Meyer, H.-D.; Cederbaum, L. S. J. Chem. Phys. 1992, 97, 3199−3213. (34) Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. J. Chem. Phys. 1996, 105, 4412−4426. (35) Beck, M.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. Phys. Rep. 2000, 324, 1−105. (36) Meyer, H.-D.; Worth, G. A. Theor. Chem. Acc. 2003, 109, 251− 267. (37) Meyer, H.-D.; Gatti, F.; Worth, G. A. Multidimensional Quantum Dynamics: MCTDH Theory and Applications; Wiley-VCH: Weinheim, 2009. (38) Köppel, H.; Domcke, W.; Cederbaum, L. S. Adv. Chem. Phys. 1984, 57, 59−246. (39) Schinke, R. Photodissociation Dynamics; Cambridge University Press: Cambridge, U.K., 1991. (40) Kosloff, R.; Tal-Ezer, H. Chem. Phys. Lett. 1986, 127, 223−230. (41) Manthe, U.; Meyer, H.-D.; Cederbaum, L. S. J. Chem. Phys. 1992, 97, 9062−9071. (42) Engel, V. Chem. Phys. Lett. 1992, 189, 76−78. (43) Jerri, A. J. The Gibbs Phenomenon in Fourier Analysis, Splines and Wavelet Approximations; Kluwer: New York, 1998. (44) Raab, A.; Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. J. Chem. Phys. 1999, 110, 936−946. (45) Worth, G. A.; Beck, M. H.; Jäckle, A.; Meyer, H.-D. The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). See http://mctdh.uni-hd.de/. (46) Kveseth, K.; Seip, R.; Kohl, D. A. Acta Chem. Scand. Ser. A. 1980, 34, 31−42. (47) Caminati, W.; Grassi, G.; Bauder, A. Chem. Phys. Lett. 1988, 148, 13−16. (48) Craig, N. C.; Groner, P.; McKean, D. C. J. Phys. Chem. A 2006, 110, 7461−7469.

excited (optically dark) 2Ag state. This is caused by a low-lying conical intersection being reached, for example, along the totally symmetric C−C stretching coordinate S3, which causes a (partial) population transfer from the bright to the dark state. Although the conical intersection has been established before in the literature, this is the first fully quantal treatment where also out-of-plane vibrational modes have been included in the modeling. The out-of-plane modes, in particular, have been demonstrated to lead to a highly irregular computed vibronic structure that represents a major cause of the diffuseness of the absorption spectrum. The experimental spectrum could be well reproduced, and less phenomenological broadening proved necessary than before. This is considered an essential step toward a fully microscopic description and understanding of the system. For a full treatment of the relaxation dynamics of transbutadiene also the interaction with the S0 ground state needs to be included. In earlier such work the 1Bu−2Ag nonadiabatic coupling was often ignored (i.e., only the two Ag states were treated). In the present line of work we aim at a fully quantal three-state study comprising the S0−S1−S2 set of states. For this purpose the range of coordinates is to be further increased as documented by Figure 10. A suitable reduced-dimensionality treatment is to be developed and to be combined with efficient diabatization techniques like the concept of regularized diabatic states. In this way we hope to develop models and tools that help to understand also bigger members of the polyene series of molecules.



AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged. B. O. acknowledges the financial support of the Ministry of Education, Science and Technological Development of the Republic of Serbia (Contract No. 172001). We are indebted to A. Dreuw for stimulating discussions.



REFERENCES

(1) Levine, B. G.; Martinez, T. Annu. Rev. Phys. Chem. 2007, 58, 613−634. (2) van der Horst, M. A.; Hellingwerf, K. J. Acc. Chem. Res. 2004, 37, 13−20. (3) Dugave, C.; Demange, L. Chem. Rev. 2003, 103, 2475−2532. (4) Bernardi, F.; Olivucci, M.; Robb, M. A. Chem. Soc. Rev. 1996, 25, 321−328. (5) Orlandi, G.; Zerbetto, F.; Zgierski, M. Z. Chem. Rev. 1991, 91, 867−891. (6) Starcke, J. H.; Wormit, M.; Schirmer, J.; Dreuw, A. Chem. Phys. 2006, 329, 39−49. (7) Boggio-Pasqua, M.; Bearpark, M. J.; Klene, M.; Robb, M. A. J. Chem. Phys. 2004, 120, 7849−7860. (8) Serrano-Andrés, A.; Merchán, M.; Nebot-Gil, I.; Lindh, R.; Roos, B. O. J. Chem. Phys. 1993, 98, 3151−3162. (9) Ostojic, B.; Domcke, W. Chem. Phys. 2001, 269, 1−10. (10) Levine, B. G.; Martinez, T. J. Phys. Chem. A 2009, 113, 12815− 12824. (11) Nakayama, K.; Nakano, H.; Hirao, K. Int. J. Quantum Chem. 1998, 66, 157−175. (12) Cave, R. J.; Zhang, F.; Maitra, N. T.; Burke, K. Chem. Phys. Lett. 2004, 389, 39−42. 8792

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793

The Journal of Physical Chemistry A

Article

(49) Feller, D.; Peterson, K. A. J. Chem. Phys. 2007, 126, 114105−1− 114105−13. (50) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053− 5063. (51) Werner, H.-J.; Knowles, P. J. Chem. Phys. Lett. 1985, 115, 259− 267. (52) Celani, P.; Werner, H.-J. J. Chem. Phys. 2000, 112, 5546−5557. (53) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1988, 89, 5803− 5814. (54) Werner, H.-J.; Knowles, P. J. Chem. Phys. Lett. 1988, 145, 514− 522. (55) Werner, H.-J.; Knowles, P. J. Theor. Chim. Acta 1992, 84, 95− 103. (56) Cave, R. J.; Davidson, E. R. J. Phys. Chem. A 1987, 91, 4481− 4490. (57) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007−1023. (58) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1995, 103, 4572− 4585. (59) Dunning, T. H.; Hay, P. J. In Methods of Electronic Structure Theory; Schaefer, H. F., Ed.; Plenum Press: New York, 1977. (60) Pulay, P.; Fogarasi, G.; Pongar, G.; Boggs, J.; Vargha, A. J. Am. Chem. Soc. 1983, 105, 7037−7047. (61) Langhoff, S. R.; Davidson, E. R. Int. J. Quantum Chem. 1974, 8, 61−72. (62) Werner, H.-J.; et al. MOLPRO, version 2008.1, a package of ab initio programs. See http://www.molpro.net. (63) Manthe, U.; Köppel, H. J. Chem. Phys. 1990, 93, 345−356. (64) Schneider, R.; Domcke, W.; Köppel, H. J. Chem. Phys. 1990, 92, 1045−1061. (65) Lévêque, C.; Komainda, A.; Taïeb, R.; Köppel, H. J. Chem. Phys. 2013, 138, 044320−1−044320−11. (66) Müller, H.; Köppel, H. Chem. Phys. 1994, 183, 107−116. (67) Manthe, U.; Köppel, H. J. Chem. Phys. 1990, 93, 345−356. (68) Manthe, U.; Köppel, H. J. Chem. Phys. 1990, 93, 1658−1669.

8793

dx.doi.org/10.1021/jp404340m | J. Phys. Chem. A 2013, 117, 8782−8793