PCCP View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

PAPER

Cite this: Phys. Chem. Chem. Phys., 2014, 16, 15957

View Journal | View Issue

The role of the low-lying dark np* states in the photophysics of pyrazine: a quantum dynamics study a Matthieu Sala,*a Benjamin Lasorne,b Fabien Gattia and Ste ´phane Gue ´rin

The excited state dynamics of pyrazine has attracted considerable attention in the last three decades. It has long been recognized that after UV excitation, the dynamics of the molecule is impacted by strong non-adiabatic effects due to the existence of a conical intersection between the B2u(pp*) and B3u(np*) electronic states. However, a recent study based on trajectory surface hopping dynamics simulations suggested the participation of the Au(np*) and B2g(np*) low-lying dark electronic states in the ultrafast radiationless decay of the molecule after excitation to the B2u(pp*) state. The purpose of this work was to pursue the investigation of the role of the Au(np*) and B2g(np*) states in the photophysics of pyrazine. A linear vibronic coupling model hamiltonian including the four lowest excited electronic states and the Received 18th May 2014, Accepted 16th June 2014

sixteen most relevant vibrational degrees of freedom was constructed using high level XMCQDPT2

DOI: 10.1039/c4cp02165g

performed and used to simulate the absorption spectrum and the electronic state population dynamics

www.rsc.org/pccp

of the system. Our results show that the Au(np*) state plays an important role in the photophysics of pyrazine.

electronic structure calculations. Wavepacket propagations using the MCTDH method were then

1 Introduction Non-adiabatic processes are ubiquitous in the photophysics and photochemistry of molecular systems.1,2 These processes are associated with a breakdown of the Born–Oppenheimer approximation and involve nuclear dynamics on several coupled potential energy surfaces (PESs). In the last few decades, the concept of conical intersection (CIs), which denotes the situation where two PESs become degenerate in a subspace of the molecular geometries, has appeared as a key concept in the study of the excited state dynamics of polyatomic molecules. In particular, conical intersections are often responsible for a very fast (from femtoseconds to picoseconds) non-radiative relaxation of ¨ppel in ref. 3 excited electronic states (see e.g. the chapter of Ko and references therein). The excited state structure and dynamics of pyrazine has attracted considerable attention in the last three decades. It has long been recognized that after UV excitation, the dynamics of the molecule is impacted by strong non-adiabatic effects. In particular, the broad band observed4–10 in the B2u(pp*) state

a

Laboratoire Interdisciplinaire Carnot de Bourgogne UMR 6303 CNRS, Universite´ de Bourgogne, BP 47870, F-21078 Dijon, France. E-mail: [email protected] b CTMM, Institut Charles Gerhardt UMR 5253 CNRS, CC 15001, Universite´ Montpellier 2, F-34095 Montpellier, France

This journal is © the Owner Societies 2014

region of the spectrum is an indication of a fast relaxation process. Recent time-resolved photoelectron spectroscopy studies11,12 gave more insight into this relaxation process. Specifically, an ultrafast (B20 fs) non-radiative decay from the B2u(pp*) state to the B3u(np*) state has been observed. On the theoretical side, extensive ab initio electronic structure calculations of the excited state energies and properties have been reported13–20 using various levels of theory. In addition, a number of quantum dynamics simulations of the absorption spectrum using multi-mode treatments of the vibronic coupling between the B3u(np*) and B2u(pp*) states have been reported. Most of these studies have been performed within the framework of the vibronic coupling model ¨ppel, Domcke and Cederbaum.21 Early developed by Ko 22–26 studies reported three- and four-mode models. These calculations gave confirmation that the broad band observed in the absorption spectrum is the signature of the existence of a CI between the B3u(np*) and B2u(pp*) states. In addition, it has been shown that a four-dimensional model allowed for a qualitatively correct simulation of the UV absorption spectrum.26 Simulations of real-time transient absorption,27–29 two-photon ionization,30 resonance Raman and fluorescence spectra31 have also been reported. Further improvements of the model by inclusion of more vibrational modes32 have culminated with the work of Worth et al.33–36 who reported accurate quantum dynamics simulations of the absorption

Phys. Chem. Chem. Phys., 2014, 16, 15957--15967 | 15957

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

Paper

spectrum using a model including the twenty-four vibrational modes of the molecule. Rigorous quantum dynamics calculations on such large systems have been made possible by the development of the multi configuration time-dependent Hartree (MCTDH) method.37–42 These benchmark results have then been used to test various approximate methods for the simulation of non-adiabatic dynamics of molecular systems.43–55 Recently, an on-the-fly trajectory surface hopping study based on time-dependent density functional theory (TDDFT) electronic structure calculations54 has suggested that the two low-lying Au(np*) and B2g(np*) dark states were significantly populated after excitation to the B2u(pp*) state. However the ability of the TDDFT method to accurately describe excited state PESs in the vicinity of CIs can be questioned. Therefore, further investigations of the role of the low-lying dark electronic states in the photophysics of pyrazine, based on multi-reference electronic structure calculations and accurate wavepacket propagation techniques, are of great interest for a better understanding of the non-radiative deactivation mechanism of this important system. In this paper, as a first step toward this goal, we derive a vibronic coupling model Hamiltonian including the four lowest lying excited electronic states, namely the B3u(np*), Au(np*), B2u(pp*) and B2g(np*) states, and the sixteen most relevant vibrational modes. This model is then used to simulate the low resolution UV absorption spectrum and the electronic state population dynamics after excitation to the B2u(pp*) state. In contrast to the results of ref. 54, our results suggest that the B2g(np*) state plays a negligible role in the non-adiabatic dynamics of pyrazine. However, we predict that the Au(np*) is significantly populated through a low-lying B2u(pp*)/Au(np*) CI. This decay channel is found to compete with the well established decay through the B2u(pp*)/B3u(np*) CI. In addition, we show the existence of a strong vibronic coupling between the Au(np*) and B3u(np*) states which influences the topography of the S1 adiabatic PES. The rest of the paper is organized as follows. The details of the electronic structure calculations performed in this work and the model Hamiltonians used for the dynamics calculations are presented in Section 2. Our results are presented and discussed in section 3 and section 4 concludes the paper.

2 Construction of the model 2.1

The diabatic model Hamiltonian

In the usual adiabatic representation, the derivative coupling terms have a singular behavior at CIs.21,56 For this reason, the quantum description of the non-adiabatic dynamics of molecular systems is most often performed in a so-called quasidiabatic representation, where the potential energy matrix elements are smooth functions of the nuclear coordinates. ¨ppel et al., the Within the vibronic coupling model of Ko diabatic Hamiltonian for several coupled electronic states is written as the sum of a reference Hamiltonian H0(Q) and a

15958 | Phys. Chem. Chem. Phys., 2014, 16, 15957--15967

PCCP

potential energy matrix W(Q), i.e H(Q) = H0(Q) + W(Q). The Hamiltonian is conveniently expressed in terms of dimensionless normal coordinates,21,26 gathered in the vector Q. The reference Hamiltonian H0(Q) is usually the ground electronic state Hamiltonian in the harmonic approximation  X oi  @ 2 2 H 0 ðQÞ ¼  2 þ Qi I; (1) 2 @Qi i where the oi are the harmonic vibrational frequencies and I is the n  n identity matrix, n being the number of electronic states included in the model. The W(Q) matrix expresses the changes in the excited state potential energy with respect to the ground state as a Taylor expansion around the ground state equilibrium geometry X ðnÞ X ðnÞ W nn ðQÞ ¼ En þ ki Qi þ gij Qi Qj þ    i

W nn0 ðQÞ ¼

X i

ðnn0 Þ

li

i;j

Qi þ

X

ðnn0 Þ

mij

(2) Qi Qj þ    ;

i;j

and where n a n 0 , En are the vertical excitation energies, k(n) i g(n) are, respectively, the linear and quadratic intrastate coupling ij constants for the nth electronic state and l(nn0) and m(nn0) are, i ij respectively, the linear and quadratic interstate coupling constants between the nth and n 0 th electronic states. In the case of highly symmetric molecules, many terms of these last two equations vanish. When only linear terms are retained (linear vibronic coupling model) the non-vanishing terms fulfill the following condition: Gn#GQ#Gn 0 *GA

(3)

where Gn and Gn 0 refer to the electronic state symmetry, GQ to the normal mode symmetry and GA is the totally symmetric irreducible representation of the symmetry point group of the molecule. Pyrazine is planar with D2h symmetry in its ground state equilibrium geometry. Its twenty-four normal modes can be classified as Gvib = 5Ag + 1B1g + 2B2g + 4B3g + 2Au + 4B1u + 4B2u + 2B3u. (4) It follows immediately from eqn (3) that only totally symmetric modes give rise to non-vanishing k(n) constants. In addition, i two given states of different symmetries are coupled by modes of a unique symmetry. In this work, we consider the four lowest excited electronic states, namely the B3u(np*), Au(np*), B2u(pp*) and B2g(np*) states, corresponding to n = 1, 2, 3 and 4 respectively. It also follows from eqn (3) that the B3g modes give rise to non-vanishing l12 i constants, the unique B1g mode gives rise to a non-vanishing l13 i constant, the B2g modes give rise to non-vanishing l23 constants and the B1u, B2u and Au i 24 modes give rise to non-vanishing l14 and l34 constants, i , li i respectively. Therefore, in a first-order description, the two B3u modes can be neglected, a priori. We construct several linear vibronic coupling model Hamiltonians including up to four electronic states and sixteen vibrational modes (see Section 2.3 for more details). Following ref. 26, quadratic diagonal terms

This journal is © the Owner Societies 2014

View Article Online

PCCP

Paper

2 g(n) i Qi are added for the non-totally symmetric modes, i.e. the modes for which the diagonal linear terms vanish by symmetry. The diabatic potential matrix elements then take the form X ðnÞ X ðnÞ ki Qi þ gj Qj2 W nn ðQÞ ¼ En þ

j

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

i

W nn0 ðQÞ ¼

X

ðnn0 Þ

lk

(5) Qi ;

k

where i runs over the totally symmetric modes, j runs over the nontotally symmetric modes and k runs over the modes of appropriate symmetry, depending on n and n0 , as explained above. The linear intrastate coupling constants k(n) i are obtained as the numerical derivative of the energy of the nth adiabatic electronic state with respect to Qi at the ground state equilibrium geometry  @Enad  ðnÞ ki ¼ : (6) @Qi Q¼0 Any displacement along a normal coordinate associated with a non-totally symmetric mode leads to a coupling of two electronic states only. Let us consider the n10a mode of B1g symmetry as an example. Displacement along the corresponding Q10a coordinate leads to a coupling of the B3u(np*) and B2u(pp*) electronic states. We obtain the following diabatic potential matrix 0 1 0 0 E1 0 B C B 0 E 0 0 C 2 B C o10a C Q10a2 I þ B W eff ðQ10a Þ ¼ B C 2 B 0 C 0 0 E 3 @ A 0 0

ð1Þ

g10a Q10a2

B B B 0 B þB B B l10a Q10a B @ 0

0

0

E4 l10a Q10a

0

g10a Q10a2

0

0

0

g10a Q10a2

0

0

0

g10a Q10a2

ð3Þ

ð4Þ

C C C C C: C C C A (7)

The eigenvalues of this matrix are the corresponding adiabatic potentials ad;ð1;3Þ ðQ10a Þ Veff

ad;ð2;4Þ ðQ10a Þ Veff

 1 ¼ o10a Q10a2 þ E1 þ E3 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ðDE þ Dg10a Q10a2 Þ2 þ 4l10a 2 Q10a2

(8)

o10a ð2;4Þ ¼ Q10a2 þ E2 þ g10a Q10a2 ; 2 g(3) 10a

g(1) 10a.

g(n) 10a

where DE = E3  E1 and Dg10a =  The l10a and constants with n = 1, 2, 3, 4 are obtained as a least-square fit of the expressions of eqn (8) to the computed ab initio data points. 2.2

Ab initio electronic structure calculations

All the electronic structure calculations presented in this work are performed with the aug-cc-pVDZ basis set of Dunning.57 The ground state geometry optimization and normal mode

This journal is © the Owner Societies 2014

CC CN CH NCC CNC NCH

This work

MP226

CASPT259

Exp.7

1.393 1.339 1.083 122.3 115.4 117.0

1.402 1.346 1.085 122.3 115.3 116.8

1.396 1.341 1.081 122.4 115.3 116.9

1.403 1.339 1.115 122.2 115.6 113.9

calculations are performed at the second-order Møller–Plesset (MP2) level of theory using the Gaussian 03 program package.58 The resulting equilibrium geometry and harmonic frequencies, compared with previous calculations and experimental values, are presented in Tables 1 and 2 respectively. Excited-state electronic structure calculations are performed using the extended multi-configuration quasi-degenerate second-order perturbation theory (XMCQDPT2) method60 using the Firefly QC package61 which is partially based on the GAMESS (US) source code.62 For the underlying state-average complete active space (SA-CASSCF) wavefunction, an active space of ten electrons in eight orbitals, including the full p orbital subset and the two nitrogen lone-pair orbitals is used. The orbitals are averaged over the five lowest CASSCF states. In addition, for the sake of comparison, multi-reference configuration interaction with single and double excitations (MRCISD) vertical excitation energies including the Davidson correction Table 2 Harmonic vibrational frequencies (in cm1) obtained at the MP2/ aug-cc-pVDZ level

1

0 ð2Þ

Table 1 Equilibrium geometry obtained in this work at the MP2/aug-ccpVDZ level compared with previous calculations and experimental values. Bond lengths are given in Å and bond angles in degrees

Symmetry

This work

MP232

MP235

Exp.7

n6a n1 n9a n8a n2

Ag Ag Ag Ag Ag

593 1017 1242 1605 3226

597 1027 1264 1633 3280

594 1024 1261 1629 3276

596 1015 1230 1582 3055

n10a

B1g

936

914

914

919

n4 n5

B2g B2g

734 942

761 913

758 915

756 983

n6b n3 n8b n7b

B3g B3g B3g B3g

700 1352 1553 3205

711 1384 1592 3254

707 1379 1586 3250

704 1346 1525 3040

n16a n17a

Au Au

337 966

343 900

343 903

341 960

n12 n18a n19a n13

B1u B1u B1u B1u

1022 1148 1486 3206

1032 1166 1456 3253

1029 1163 1516 3250

1021 1136 1484 3012

n18b n14 n19b n20b

B2u B2u B2u B2u

1079 1364 1440 3221

1092 1369 1456 3277

1089 1367 1453 3272

1063 1149 1416 3063

n16b n11

B3u B3u

417 797

426 781

423 781

420 785

Phys. Chem. Chem. Phys., 2014, 16, 15957--15967 | 15959

View Article Online

Paper

PCCP

Table 3 Vertical excitation energies computed using the XMCQDPT2 method and a SA5-CASSCF reference wavefunction compared with previous theoretical results and experimental data

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

a

XMCQDPT2 MRCISD+Qa CASPT220 CASPT264 CASPT215 CASPT217 CC265 CC365 EOM-CCSD(T)16 SAC-CI19 TDDFT/B3LYP54 Exp.10 a

This work.

b

B3u(np*)

Au(np*)

B2u(pp*)

B2g(np*)

3.93 4.35 3.86 4.02 3.85 3.83 4.26 4.24 3.83 4.25 3.96 3.83b

4.45 5.26 4.52 4.75 4.63 4.36 4.95 5.05 4.81 5.24 4.6 —

4.79 4.99 4.81 4.80 4.76 4.79 5.13 5.02 4.64 4.84 5.46 4.81c(4.69b)

5.38 5.77 5.48 5.56 — 5.50 5.92 5.74 5.56 6.04 6.3 —

0–0 transition. c Band maximum.

(denoted +Q) and using the same reference SA-CASSCF wavefunction are computed using the MOLPRO program package.63 Our results for the vertical excitation energies of the four lowest lying singlet excited electronic states, in comparison with previous calculations performed using a variety of methods and experimental data, are presented in Table 3. The 0–0 transition of the B3u(np*) state has been consistently measured at 3.83 eV in previous experimental studies.4–7,9,10 The 0–0 transition being the highest intensity transition, the equilibrium geometry of the B3u(np*) state is expected to be rather close to that of the ground state and the vertical excitation energy of the B3u(np*) state should not be much higher than 3.83 eV. Most previous theoretical work, including various CASPT2 calculations as well as TDDFT and EOM-CCSD(T) calculations, predicted a B3u(np*) state vertical excitation energy between 3.83 eV and 4.02 eV, in good agreement with our value of 3.93 eV. In contrast, the CC2, CC3, SAC-CI and MRCI methods yielded significantly higher vertical excitation energies of 4.26, 4.24, 4.25 and 4.35 eV respectively. Similarly, the 0–0 transition and absorption band maximum of the B2u(pp*) state has been consistently measured at 4.69 and 4.81 eV, respectively.4,5,10 Again, our vertical excitation energy of 4.79 eV is in good agreement with previous CASPT2 and SAC-CI results. In contrast, the EOM-CCSD(T) value of ref. 16 appears slightly underestimated while the CC2, CC3 values of ref. 65 and our MRCI value appear slightly overestimated. Finally, we note that the TDDFT value of ref. 54 is significantly overestimated. As seen in Table 3 the vertical excitation energy of the dark Au(np*) state is a matter of debate. Walker and Palmer proposed an experimental value of 5.0 eV based on nearthreshold energy loss spectra.9 However, other studies quote a value of 4.72 eV.19,20 From the theoretical point of view, the Au(np*) state vertical excitation energy appears more sensitive to the methodology used for the calculation. Our XMCQDPT2 and previous CASPT2 calculations predicted the Au(np*) state as the second singlet excited state, i.e. below the B2u(pp*) state, with vertical excitation energies ranging from 4.36 to 4.75 eV. Previous CC2 and TDDFT calculations also predicted the Au(np*) state to lie below the B2u(pp*) state. CC3 calculations,

15960 | Phys. Chem. Chem. Phys., 2014, 16, 15957--15967

in contrast, predicted the Au(np*) and B2u(pp*) states to be essentially degenerate whereas the SAC-CI, EOM-CCSD(T) and MRCI calculations predicted the Au(np*) state to lie above the B2u(pp*) state. A similar uncertainty exists for the position of the B2g(np*) state. Okuzawa et al.8 reported a value of 5.19 eV based on UV-IR double resonance dip spectroscopy measurements. However, values of 5.5 eV 7,19 and 6.10 eV 5,20 have also been quoted in the literature. From the theoretical point of view, as for the Au(np*) state, rather large variations exist in the vertical excitation energies calculated using different methods. The CASPT2 and EOM-CCSD(T) values reported in Table 3 are in good agreement with each other, ranging from 5.48 to 5.56 eV. The values obtained from other methods are typically higher, ranging from 5.74 eV to 6.3 eV. Our XMCQDPT2 value of 5.38 eV appears slightly underestimated with respect to previous calculations. The interested reader is referred to the work of Weber and Reimers17,18 for further extensive comparisons of computational methods for the excitation energies of the low-lying singlet and triplet states of pyrazine. 2.3

Details of the model diabatic Hamiltonians

The values of the parameters entering the definition of the vibronic coupling model Hamiltonians used for the quantum dynamics calculations reported in Section 3, obtained as explained in Section 2.1 using our XMCQDPT2 ab initio data are reported below. Table 4 presents the values of the linear intrastate coupling constants k(n) i . These values represent the gradients of the excited states at the Franck–Condon (FC) geometry. We also report the dimensionless quantities k(n) i /oi, which give a better measure of the relative importance of the different totally symmetric modes in the non-adiabatic decay dynamics of the molecule after photoexcitation.66,67 In this work, we are particularly interested in the decay mechanism after excitation to the B2u(pp*) state. An inspection of Table 4 shows that the B2u(pp*) state has large gradients along the n6a and n1 modes. In addition, one can note the very large gradient of the Au(np*) state along the n8a mode. In contrast, the high frequency n2 mode is expected to play a marginal role in the dynamics of the molecule. The relative magnitudes of the interstate coupling constants, presented in Table 5, give further insight into the decay mechanism of the molecule. Our calculations show a strong coupling of the bright B2u(pp*) state with the B3u(np*) and B2g(np*) states, essentially mediated by the n10a and n16a modes, respectively. The coupling of the B2u(pp*) state with the Au(np*) Table 4 Linear intrastate coupling constant k(n) values (in eV) obtained i in this work. The values between parentheses are the dimensionless quantities k(n) i /oi

B3u(np*)

Au(np*)

B2u(pp*)

B2g(np*)

k6a 0.081(1.103) 0.168(2.283) 0.128(1.739) 0.184(2.500) k1 0.038(0.304) 0.083(0.659) 0.183(1.452) 0.117(0.926) 0.165(1.073) k9a 0.117(0.762) 0.071(0.459) 0.045(0.295) 0.172(0.862) k8a 0.087(0.436) 0.465(2.338) 0.026(0.132) 0.022(0.054) 0.060(0.150) 0.018(0.044) 0.030(0.074) k2

This journal is © the Owner Societies 2014

View Article Online

PCCP

Paper

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

Table 5 Quadratic intrastate g(n) and linear interstate coupling constant i l(nn0) values (in eV) obtained in this work. The values between parentheses i are the dimensionless quantities l(nn0) /oi i

Mode

Symm.

l

g(1)

g(2)

g(3)

g(4)

n10a

B1g

0.195(1.684)

0.012

0.048

0.012

0.013

n4 n5

B2g B2g

0.060(0.663) 0.053(0.456)

0.030 0.014

0.031 0.026

0.031 0.026

0.027 0.009

n6b n3 n8b n7b

B3g B3g B3g B3g

o105 0.065(0.385) 0.219(1.140) 0.020(0.050)

0.013 0.006 0.012 0.003

0.013 0.006 0.012 0.003

0.005 0.001 0.007 0.004

0.006 0.004 0.043 0.003

n16a n17a

Au Au

0.112(2.686) 0.018(0.146)

0.013 0.016

0.013 0.041

0.008 0.012

0.008 0.012

n12 n18a n19a n13

B1u B1u B1u B1u

0.207(1.629) 0.090(0.632) 0.094(0.511) o105

0.006 0.006 0.006 0.004

0.022 0.002 0.010 0.004

0.006 0.005 0.002 0.004

0.006 0.006 0.006 0.004

n18b n14 n19b n20b

B2u B2u B2u B2u

0.044(0.332) 0.044(0.263) 0.072(0.404) o105

0.001 0.019 0.013 0.005

0.003 0.021 0.006 0.003

0.002 0.020 0.015 0.004

0.003 0.021 0.006 0.003

in the two-state model, the three-state model includes the n4 and n5 modes of B2g symmetry and the n3 and n8b modes of B3g symmetry. Therefore, a total of nine modes are included in the three-state model. When the model is extended to include the B2g(np*) state, the Au, B1u and B2u modes need to be considered. The four-state model considered in this work includes the n16a mode of Au symmetry, the n12, n18a and n19a modes of B1u symmetry and the n18b, n14 and n19b modes of B2u symmetry, resulting in a total of sixteen modes.

3 Time-dependent nuclear quantum dynamics simulations 3.1

The multi-configuration time-dependent Hartree method

¨dinger equation is solved The nuclear time-dependent Schro using the multi-configuration MCTDH method37–40 as implemented in the Heidelberg MCTDH package.42 The wavefunction, in the so-called multi-set formalism68 which is used in this work, reads jCi ¼

state, through the n4 and n5 modes, appears comparatively weaker. In addition, we predict strong couplings between the B3u(np*) and Au(np*) states, essentially through the n8b mode, and between the B3u(np*) and B2g(np*) states, essentially through the n12 mode. In order to get some insight into the role of the Au(np*) and B2g(np*) states in the non-adiabatic dynamics of pyrazine, three different models, differing in the number of electronic states included, are considered in the quantum dynamics calculations presented in Section 3: a two-state model, including only the B3u(np*) and B2u(pp*) states, a three-state model including the Au(np*) state and a four-state model including the B2g(np*) state. In the two-state model, as explained in Section 2.1, only the Ag and B1g modes give rise to non-vanishing first-order coupling constants. Therefore, the two-state model considered in this work includes the four most important totally symmetric modes (n6a, n1, n9a and n8a) and the n10a mode. Similar two-state models have been previously used26,32,59 to simulate the excited state dynamics of pyrazine. When the Au(np*) state is included in the model, the B2g and B3g modes need to be considered as they give rise to nonvanishing coupling constants. Besides the five modes included

Table 6

na X

CðaÞ ðQ; tÞjai;

(9)

a¼1

where {|ai} denotes the set of electronic states. The nuclear wavefunction C(a)(Q,t) in the electronic state |ai is expanded on the basis of time-dependent functions j(qk,t) called singleparticle functions (SPFs) CðaÞ ðQ; tÞ ¼ CðaÞ ðq; tÞ ðaÞ

¼

n1 X

ðaÞ



ðaÞ

j1

np X

ðaÞ

Aj1 ;;jp

p Y k¼1

ðaÞ

jp

ðk;aÞ

(10)

j ðaÞ ðqk ; tÞ: jk

The SPFs can be multi-dimensional, i.e. the coordinate qk may be a collective one qk = (Qa,Qb,. . .). The SPFs are expanded in a primitive basis, built as a direct product of one-dimensional basis functions for each degree of freedom ðkÞ

jjk ðqk ; tÞ ¼

Nk X

ðkÞ

akjk wkk ðqk Þ:

(11)

k¼1

The equations of motion for the expansion coefficients and for the SPFs are derived variationally, ensuring an optimal description of the wavepacket during its evolution.

Number of SPF and primitive basis functions used in the calculations

Model

Combinations of modes

Numbers of SPFs

Numbers of grid points

Two-state

(n6a, n10a), (n1, n9a, n8a)

[4,34,14], [4,12,8], [4,12,8], [4,10,6]

(32,40), (16,14,24)

Three-state

(n6a, n10a), (n1, n4), (n9a, n3, n8b), (n8a, n5)

[4,30,34,14], [4,12,15,8], [4,13,14,8], [4,20,24,8]

(32,40), (16,20), (14,10,14), (24,10)

Four-state

(n6a, n10a, n16a), (n1, n4, n19b), (n9a, n3, n8b), (n8a, n5, n12), (n18a, n18b, n19a, n14)

[4,36,42,14,6], [4,13,14,8,5], [4,24,28,8,5], [4,18,19,8,5], [4,7,7,4,4]

(32,40,20), (16,20,10), (14,10,14), (24,10,11), (7,7,7,8)

This journal is © the Owner Societies 2014

Phys. Chem. Chem. Phys., 2014, 16, 15957--15967 | 15961

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

Paper

PCCP

In this work, for the representation of the Hamiltonian and the wave function, a Hermite polynomial DVR scheme69 was used for all the degrees of freedom. The number of SPF and primitive basis functions used in the calculations are listed in Table 6. Test calculations with both larger primitive and SPF bases have been performed. The results (populations and spectra) of these calculations were found to be almost identical to those performed with the bases of Table 6, indicating a satisfactory convergence of the results presented in this paper. 3.2

Simulation of the UV absorption spectrum

In this section, we report the results of our simulations of the low resolution UV absorption spectrum in the region between 3.8 and 5.7 eV, using the model Hamiltonians presented in Section 2.3. In this energy range, the spectrum is dominated by two distinct bands attributed to the B3u(np*) and B2u(pp*) states. The absorption spectrum is computed as the Fourier transform of the autocorrelation function a(t) given by a(t) = hC(Q,0)|C(Q,t)i,

(12)

where C(Q,0) is the initial wavefunction. To account for the homogeneous broadening of the experimental spectrum, the autocorrelation function is pre-multiplied by a damping function f (t) = et/t, where the damping time t is a free parameter. This is equivalent to convoluting the spectrum by a Lorentzian function of full width at half maximum of 2/t. In addition, to avoid problems arising from the finite propagation time T (the so-called ‘‘Gibbs phenomenon’’), the autocorrelation function was further pre-multiplied by a filter function g(t) = cos2(pt/2T)Y(t  T), where Y denotes the Heaviside step function. Photoexcitation to the B3u(np*) and B2u(pp*) states are both allowed processes whereas photoexcitation to the Au(np*) and B2g(np*) states are symmetry-forbidden processes at the D2h ground state equilibrium geometry. Two wavepacket propagations are performed, starting on the B3u(np*) and B2u(pp*) states, respectively, from which two partial absorption spectra are computed. In each case, the initial wavefunction is obtained as a vertical electronic excitation, i.e. it corresponds to the ground vibronic state wavefunction projected in each excited electronic state. The full absorption spectrum is then obtained as the oscillator strength-weighted sum of the two partial absorption spectra. Based on extensive test calculations, an adjustment of the Au(np*) state vertical excitation energy was found to be necessary to obtain a reasonable agreement between simulated and experimental spectra. A value of 4.69 eV was found to provide the most satisfactory description of the absorption spectrum. This value is somewhat higher than our ab initio value of 4.45 eV. Nevertheless, the magnitude of this adjustment appears reasonable given the large deviations existing between vertical excitation energies of the Au(np*) state computed using different methods, as seen in Table 3. The experimental estimations of 0.006 and 0.1 reported in ref. 7 for the oscillator

15962 | Phys. Chem. Chem. Phys., 2014, 16, 15957--15967

strengths of the B3u(np*) and B2u(pp*) states, respectively, are used. Each partial absorption spectrum is computed from a propagation of 120 fs. Our simulated absorption spectra, using the three different models presented in Section 2.3, compared with the experimental spectrum4 are presented in Fig. 1. In each case, the whole computed spectrum is blueshifted so that the position of the B3u(np*) 0–0 peak matches the experimental one at approximately 324 nm. The spectrum computed with the two-state model is presented in Fig. 1(b). In this case, damping times t of 400 fs and 40 fs are used for the B3u(np*) and B2u(pp*) components, respectively. The whole spectrum is blue-shifted by 0.04 eV. The two-state model is seen to satisfactorily reproduce the shape of the B2u(pp*) band. However this band appears slightly red-shifted with respect to the experimental one (the experimental spectrum is presented in Fig. 1(a)). In contrast, the agreement for the B3u(np*) band is less satisfactory. In particular, the relative intensities of the two most intense peaks at approximately 324 nm and 319 nm, corresponding to the 0–0 and n16a transitions, respectively, disagree with experiment. The spectrum computed with the three-state model is presented in Fig. 1(c). Damping times t of 400 fs and 100 fs are used for the B3u(np*) and B2u(pp*) components, respectively. The fact that a larger damping time parameter (leading to a less important broadening added to the spectrum) is needed for the

Fig. 1 Comparison of (a) the experimental UV absorption spectrum taken from ref. 4 with the spectra computed with (b) the two-state model, (c) the three-state model and (d) the four-state model.

This journal is © the Owner Societies 2014

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

PCCP

three-state model indicates that the density of vibronic states in the region of the B2u(pp*) state is significantly higher than in the two-state model. The spectrum is blue-shifted by 0.13 eV to approximately match the experimental B3u(np*) 0–0 peak. The B2u(pp*) band position is in better agreement with respect to experiment than in the 5D two-state model. The inclusion of the Au(np*) state in the simulations also has a significant impact on the shape of the spectrum. The B2u(pp*) band shape is slightly improved, especially in the 245–255 nm range, with respect to the two-state model, and is in almost perfect agreement with experiment. The relative intensity of the 0–0 and n16a peaks of the B3u(np*) band is also better reproduced by the three-state model. Nevertheless, the agreement with experiment remains more qualitative than for the B2u(pp*) band. The spectrum computed with the four-state model, presented in Fig. 1(d), is almost identical to the spectrum obtained from the three-state model. Damping times t of 400 fs and 100 fs are used for the B3u(np*) and B2u(pp*) components, respectively. The spectrum is blue-shifted by 0.17 eV to approximately match the experimental B3u(np*) 0–0 peak. This is slightly more than the blue-shift needed for the three-state model. In addition, the B2u(pp*) band now appears slightly blue-shifted with respect to its experimental counterpart. This is the sign of the influence of the B2g(np*) state, through the vibronic coupling, on the topography of the PESs of the lower states. However, as mentioned in Section 2.2, our value of 5.38 eV for the vertical excitation energy of the B2g(np*) state appears underestimated with respect to previous calculations. As a higher vertical excitation energy would imply a weaker impact of the vibronic coupling between the B2g(np*) state and the lower states, this slight influence might be an artifact of our model. The results presented in this section allow us, by taking into account the overall energy shift and the adjustment of the Au(np*) state vertical excitation energy applied in the three-state model, to propose refined estimations of 4.06, 4.82 and 4.92 eV for the vertical excitation energies of the B3u(np*), Au(np*) and B2u(pp*) states, respectively. 3.3

Electronic state populations and decay mechanism

In a recent experimental study based on time-resolved photoelectron spectroscopy,12 the radiationless decay of the molecule after excitation to the B2u(pp*) state has been directly observed using sub-20 fs pump and probe pulses. In this section, we analyze the wavepacket dynamics triggered by a 14 fs pump pulse resonant with the electronic transition from the ground to the B2u(pp*) state. The pulse is characterized by a sinesquared envelope, a peak amplitude of 0.01 a.u. and a photon energy of 4.7 eV. We stress that, contrarily to the calculations presented in the previous section, in this section the ground electronic state is explicitly included in our models. Nevertheless, for clarity, we continue to refer to them as two-, three- and four-state models. The dynamics of the system is analyzed in terms of diabatic and adiabatic electronic state populations. In this section, the diabatic electronic states are labeled by Greek letters whereas

This journal is © the Owner Societies 2014

Paper

the adiabatic electronic states are labeled by Latin letters. The diabatic population for the state |ai is simply the norm of the corresponding wavefunction component (see eqn (9)) Pda = 8Cd,(a)8.

(13)

The adiabatic populations are more difficult to obtain. We first introduce the transformation matrix U(Q) that diagonalizes the diabatic potential matrix W(Q) of eqn (5) U†(Q)W(Q)U(Q) = V(Q),

(14)

where V(Q) is the diagonal matrix of the adiabatic potentials. This transformation matrix relates the adiabatic wavefunction to the diabatic one Cad(Q) = U(Q)Cd(Q).

(15)

Pad a

of the state |ai is obtained as the The adiabatic population ˆad expectation value of the corresponding projection operator P a d ˆad d Pad a = hC |Pa |C i,

(16)

with ad P^a ¼

X

y jbiUba Uag hgj:

(17)

b;g

Unfortunately, the adiabatic projection operator does not have the so-called product form.37–40 Therefore, the integrals of eqn (16) have to be evaluated on the full primitive grid, which can become too much time consuming for systems of more than four or five degrees of freedom. For larger systems, as in the present work, a Monte Carlo integration scheme, implemented in the Heidelberg MCTDH package, can be used. In addition, to reduce the size of the grid, the so-called ‘‘quick algorithm’’ (see the MCTDH html documentation42) is used: all the grid points for which the product of the one-dimensional reduced densities is lower than a given threshold are ignored. The onedimensional reduced density along the coordinate Qi is the density integrated over all the degrees of freedom except Qi Di(Qi) = hC|CiQi.

(18)

To distinguish them from the diabatic states, the adiabatic states are noted Sn with n = 1, 2, 3, 4. We first present the diabatic and adiabatic electronic state populations obtained from the two-state model shown in Fig. 2(a) and (b), respectively. The electronic state populations were analyzed in detail previously using similar two-state models.35,36 The population transferred to the B2u(pp*) state by the laser pulse reaches a maximum of approximately 0.6 at 12 fs and quickly decays to the B3u(np*) state. After 50 fs, the B2u(pp*) state population drops to 0.1. A large recurrence occurs at 90 fs, as noted previously.35,36 Since the B2u(pp*) state is the third excited state at the FC geometry, the two adiabatic excited states in the two-state model are labeled S1 and S3. The procedure explained above used to compute the adiabatic electronic state populations is necessarily approximate. To get some insight into the accuracy of this procedure, two calculations, using different standard

Phys. Chem. Chem. Phys., 2014, 16, 15957--15967 | 15963

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

Paper

Fig. 2 Electronic state populations computed with the two-state model. (a) Populations of the diabatic B3u(np*) (blue) and B2u(pp*) (green) states and (b) populations of the adiabatic S1 (black) and S3 (cyan) states using a Monte Carlo integration scheme combined with the so-called ‘‘quick algorithm’’ (see text for details). The full lines present the adiabatic populations computed using a standard deviation of 0.01 for the Monte Carlo integration and a threshold of 1.0  106 for the quick algorithm. Adiabatic populations computed using values of 0.001 and 1.0  108 for the same quantities are presented as dashed lines.

deviations for the Monte-Carlo integration and thresholds for the quick algorithm, are performed. These quantities are set to 0.01 and 1.0  106, respectively, in the first calculation and to 0.001 and 1.0  108, respectively, in the second calculation. The resulting adiabatic state populations, presented in Fig. 2(b), show some discontinuities. These discontinuities reflect the limited accuracy of the Monte Carlo integration, and are particularly important in the S1 population after 40 fs. This corresponds to the vibrationally hot wavepacket propagating on the S1 PES after the initial decay through the S2/S1 CI. Nevertheless, except for this oscillatory part of the S1 population curve, the results obtained from the two calculations are almost identical. This indicates that the Monte Carlo integration combined with the quick algorithm provides a sufficient accuracy for our purpose. As noted previously,35,36 the decay of the adiabatic population is faster than in the diabatic case, and the final population of the upper state is lower. The diabatic and adiabatic electronic state populations obtained with the three-state model are shown as full lines in Fig. 3(a) and (b), respectively. The diabatic populations obtained with the four-state model are also shown as dashed lines in Fig. 3(a). The comparison of the diabatic populations for the three- and four-state models confirms that the inclusion of the B2g(np*) state in the simulations has a minor effect on the non-adiabatic decay dynamics of the molecule after excitation to the B2u(pp*) state. Only a very small amount of population is transferred to the B2g(np*) state, with a maximum of less than 0.02 at 11 fs. In addition, the populations of the B3u(np*), Au(np*) and B2u(pp*) states are similar in the three- and fourstate models. The population of the B2u(pp*) state reaches a maximum of nearly 0.6 at 11 fs, and then quickly decays to almost zero at 50 fs. A recurrence is then seen at 95 fs, similarly to the two-state model simulation. Between 0 and 20 fs, both the B3u(np*) and Au(np*) state populations rise quickly and

15964 | Phys. Chem. Chem. Phys., 2014, 16, 15957--15967

PCCP

Fig. 3 Electronic state populations computed with the three- and fourstate models. (a) Populations of the diabatic B3u(np*) (blue), Au(np*) (red), B2u(pp*) (green) and B2g(np*) (magenta) states with the three-state model (full lines) and four-state model (dashed lines). The three-state model includes only the B3u(np*), Au(np*) and B2u(pp*) states. (b) Populations of the adiabatic S1 (black), S2 (orange) and S3 (cyan) states computed using a standard deviation of 0.01 for the Monte Carlo integration and a threshold of 1.0  106 for the quick algorithm.

reach approximately 0.15 at 20 fs. Then, between 20 fs and 40 fs, the B3u(np*) state population continues to rise and reaches a value of 0.4 at 40 fs. In contrast, during the same time interval, the Au(np*) state population slightly drops to 0.13. After 40 fs the Au(np*) state population increases again while the B3u(np*) state population starts to drop and the population oscillates between these two states until the end of the simulation. For the calculation of the adiabatic state populations presented in Fig. 3(b), a standard deviation of 0.01 for the Monte Carlo integration and a threshold of 1.0  106 for the quick algorithm are used. Overall, this figure shows a fast and sequential decay from the S3 state to the S1 state via the S2 state. At 40 fs the S1 population reaches a value of approximately 0.6 and remains essentially constant until the end of the simulation. In addition, a residual population of approximately 0.05 is seen to remain on the S2 state. One-dimensional cuts of the diabatic PESs along the Q6a and Q8a totally symmetric coordinates are presented in Fig. 4(a) and

Fig. 4 One-dimensional cuts of the diabatic PESs of the B3u(np*) (blue), Au(np*) (red), B2u(pp*) (green) and B2g(np*) (magenta) states along the (a) Q6a and (b) Q8a totally symmetric coordinates. The full lines represent the potential obtained from our model Hamiltonian whereas the squares are the ab initio computed points.

This journal is © the Owner Societies 2014

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

PCCP

Fig. 5 Populations of the diabatic B3u(np*) (blue), Au(np*) (red) and B2u(pp*) (green) states computed with the three-state model where the l3 and l8b coupling constants, mediating the vibronic coupling between the B3u(np*) and Au(np*) states, are set to zero.

(b) respectively. Both figures show B2u(pp*)/Au(np*) CI points close to the FC geometry. Upon motion along the Q6a coordinate, the B2u(pp*) state then crosses the B3u(np*) state. This B2u(pp*)/B3u(np*) CI is responsible for the decay of the molecule in the two-state model. The B2u(pp*)/Au(np*) CI appears much closer to the FC geometry than the B2u(pp*)/B3u(np*) CI. However, as seen in Table 5, the l4 and l5 coupling constants mediating the B2u(pp*)/Au(np*) coupling are smaller than the l10a constant mediating the B2u(pp*)/B3u(np*) coupling. This explains the competition between population transfer to the B3u(np*) and Au(np*) in the first 40 fs seen in Fig. 3(a). Fig. 4(b) shows the existence of both B2u(pp*)/Au(np*) and Au(np*)/ B3u(np*) low-lying CI points. The Au(np*)/B3u(np*) CI explains the oscillations of population between the B3u(np*) and Au(np*) seen after 40 fs in Fig. 3(a). To confirm this interpretation, we perform a calculation using the three-state model where the l3 and l8b coupling constants are set to zero. Therefore in this model, the B3u(np*) and Au(np*) states are no longer vibronically coupled. The resulting diabatic electronic state populations are shown in Fig. 5. In the first 40 fs the diabatic populations are similar to that obtained in the full three-state model (Fig. 3(a)). At 45 fs, the populations of the B3u(np*) and Au(np*) states are approximately 0.4 and 0.2, respectively. This confirms the competition between a decay to the Au(np*) state via the B2u(pp*)/Au(np*) CI and a decay to the B3u(np*) state via the B2u(pp*)/B3u(np*) CI during the first 40 fs after excitation to the B2u(pp*) state. After 40 fs, the populations of the B3u(np*) and Au(np*) states remain roughly constant, except for the recurrence between the B3u(np*) and B2u(pp*) states seen at 90 fs. Therefore the oscillations between the B3u(np*) and Au(np*) states in the full three-state model calculation (Fig. 3(a)) are due to the vibronic coupling between these states.

4 Conclusion In the last two decades, the pyrazine molecule has been considered as a benchmark system for the study of the radiationless decay of molecules at CIs. The ultrafast decay observed experimentally11,12 has been interpreted within a two-state model including the vibronically coupled B3u(np*) and B2u(pp*) states.

This journal is © the Owner Societies 2014

Paper

However, a recent study based on trajectory surface-hopping calculations54 has suggested the participation of the two lowlying dark Au(np*) and B2g(np*) states in the photophysics of the molecule. In this paper, we have performed quantum dynamics simulations of the excited state non-adiabatic decay of the molecule. For this purpose, a linear vibronic coupling model Hamiltonian augmented with diagonal quadratic terms for the non-totally symmetric modes, including the four lowest lying excited electronic states and up to sixteen vibrational modes, has been constructed. In contrast to the results presented in ref. 54, our calculations predict a negligible influence of the B2g(np*) states on the dynamics of the molecule. However, our calculations suggest that the decay to the Au(np*) state, through a low-lying B2u(pp*)/Au(np*) CI, efficiently competes with the well established decay to the B3u(np*) state through the B2u(pp*)/B3u(np*) CI. In addition, we have shown the existence of a strong vibronic coupling between the Au(np*) and B3u(np*) states which influences the topography of the S1 adiabatic PES. Therefore, the present work provides a significant new insight into the non-adiabatic dynamics of pyrazine after excitation to the B2u(pp*) state. Nevertheless, although our model provides an excellent agreement with experiment in the simulation of the B2u(pp*) band of the UV absorption spectrum, the B3u(np*) band is only qualitatively reproduced. We hope that an extension of the model presented in this paper to a fully quadratic vibronic coupling model will provide a more accurate description of the B3u(np*) band of the UV spectrum. This will be the subject of future work.

Acknowledgements We gratefully acknowledge Prof. Hans-Dieter Meyer for numerous discussions. We acknowledge support from the CoConicS Project (ANR-13-BS08-0013-03). M. S. and S. G. acknowledge support from the Conseil Regional de Bourgogne. Part of the calculations were performed using HPC resources from DSICCUB (Universit de Bourgogne).

References ¨ppel, Conical Intersec1 W. Domcke, D. Yarkony and H. Ko tions: Electronic Structure, Dynamics & Spectroscopy, World Scientific, 2004. ¨ppel, Conical Intersections: 2 W. Domcke, D. Yarkony and H. Ko Theory, Computation and Experiment, World Scientific, 2011. 3 F. Gatti, Molecular Quantum Dynamics, Springer, Heidelberg, 2014. 4 I. Yamazaki, T. Murao, T. Yamanaka and K. Yoshihara, Faraday Discuss. Chem. Soc., 1983, 75, 395. 5 A. Bolovinos, P. Tsekeris, J. Philips, E. Pantos and G. Anditsopoulos, J. Mol. Spectrosc., 1984, 103, 240. 6 J. Kommandeur, W. A. Majewski, W. L. Meerts and D. W. Pratt, Annu. Rev. Phys. Chem., 1987, 38, 433.

Phys. Chem. Chem. Phys., 2014, 16, 15957--15967 | 15965

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

Paper

7 K. K. Innes, I. G. Ross and W. R. Moomaw, J. Mol. Spectrosc., 1988, 132, 492. 8 Y. Okuzawa, M. Fufii and M. Ito, Chem. Phys. Lett., 1990, 171, 341. 9 I. C. Walker and M. H. Palmer, Chem. Phys., 1991, 153, 169. 10 M. Stener, P. Decleva, D. M. P. Holland and D. A. Shaw, J. Phys. B: At., Mol. Opt. Phys., 2011, 44, 075203. 11 V. Stert, P. Farmanara and W. Radloff, J. Chem. Phys., 2000, 112, 4460. 12 Y.-I. Suzuki, T. Fuji, T. Horio and T. Suzuki, J. Chem. Phys., 2010, 132, 174302. 13 W. R. Wadt, W. A. G. III and T. H. Dunning Jr., J. Chem. Phys., 1976, 65, 438. 14 M. P. Flscher, K. Andersson and B. O. Roos, J. Chem. Phys., 1992, 96, 9204. 15 M. P. Flscher and B. O. Roos, Theor. Chim. Acta, 1994, 87, 403. 16 J. E. Delbene, J. D. Watts and R. J. J. Bartlett, J. Chem. Phys., 1997, 106, 6051. 17 P. Weber and J. R. Reimers, J. Phys. Chem. A, 1999, 103, 9821. 18 P. Weber and J. R. Reimers, J. Phys. Chem. A, 1999, 103, 9830. 19 Y. Li, J. Wan and X. Xu, J. Comput. Chem., 2007, 28, 1658. ´sz and A. Vibo ´k, Theor. Chem. 20 C. Woywod, A. Papp, G. J. Hala Acc., 2010, 125, 521. ¨ppel, W. Domcke and L. S. Cederbaum, Adv. Chem. 21 H. Ko Phys., 1984, 54, 59. 22 R. Schneider and W. Domcke, Chem. Phys. Lett., 1988, 150, 235. 23 R. Schneider and W. Domcke, Chem. Phys. Lett., 1989, 159, 61. ¨ppel, J. Chem. Phys., 1990, 93, 1658. 24 U. Manthe and H. Ko 25 R. Seidner, G. Stock, A. L. Sobolewski and W. Domcke, J. Chem. Phys., 1991, 96, 5298. 26 C. Woywod, W. Domcke, A. L. Sobolewski and H.-J. Werner, J. Chem. Phys., 1994, 100, 1400. 27 G. Stock, R. Schneider and W. Domcke, J. Chem. Phys., 1989, 90, 7184. 28 G. Stock and W. Domcke, J. Opt. Soc. Am. B, 1990, 7, 1970. 29 G. Stock and W. Domcke, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 3032. 30 M. Seel and W. Domcke, J. Chem. Phys., 1991, 95, 7806. 31 G. Stock and W. Domcke, J. Chem. Phys., 1990, 93, 5496. 32 G. Stock, C. Woywod, W. Domcke, T. Swinney and B. H. Hudson, J. Chem. Phys., 1995, 103, 6851. 33 G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1996, 105, 4412. 34 G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1998, 109, 3518. 35 G. A. Worth, H.-D. Meyer and L. S. Cederbaum, Chem. Phys. Lett., 1999, 299, 451. 36 A. Raab, G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1999, 110, 936. 37 H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73.

15966 | Phys. Chem. Chem. Phys., 2014, 16, 15957--15967

PCCP

38 U. Manthe, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199. ¨ckle, G. A. Worth and H.-D. Meyer, Phys. 39 M. H. Beck, A. Ja Rep., 2000, 324, 1. 40 H.-D. Meyer and G. A. Worth, Theor. Chem. Acc., 2003, 109, 251. 41 H.-D. Meyer, F. Gatti and G. A. Worth, Multidimensional Quantum Dynamics: MCTDH Theory and Applications, WileyVCH, Weinheim, 2009. ¨ckle and H.-D. Meyer, The 42 G. A. Worth, M. H. Beck, A. Ja MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). See http://mctdh.uni-hd.de. 43 M. Thoss, W. H. Miller and G. Stock, J. Chem. Phys., 2000, 112, 10282. 44 M. Santer, U. Manthe and G. Stock, J. Chem. Phys., 2001, 114, 2001. 45 K. Ando and M. Santer, J. Chem. Phys., 2003, 118, 10399. 46 R. Borrelli and A. Peluso, J. Chem. Phys., 2003, 119, 8437. 47 D. V. Shalashilin and M. S. Child, J. Chem. Phys., 2004, 121, 3563. 48 P. Puzari, R. S. Swathi, B. Sarkar and S. Adhikari, J. Chem. Phys., 2005, 123, 134317. 49 P. S. Christopher, M. Shapiro and P. Brumer, J. Chem. Phys., 2006, 124, 184107. 50 P. Puzari, B. Sarkar and S. Adhikari, J. Chem. Phys., 2006, 125, 194316. 51 X. Chen and V. S. Batista, J. Chem. Phys., 2006, 125, 124313. 52 I. Burghardt, K. Giri and G. A. Worth, J. Chem. Phys., 2008, 129, 174104. 53 C. Lasser and T. Swart, J. Chem. Phys., 2008, 129, 034302. 54 U. Werner, R. Mitric´, T. Suzuki and V. Bonacic´-Koutecky´, Chem. Phys., 2008, 349, 319. 55 D. V. Shalashilin, J. Chem. Phys., 2010, 132, 244111. 56 D. R. Yarkony, Chem. Rev., 2012, 112, 481. 57 R. A. Kendall, T. H. Dunning Jr. and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796. 58 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian, Inc., Wallingford, CT, 2004.

This journal is © the Owner Societies 2014

View Article Online

Published on 18 June 2014. Downloaded by State University of New York at Stony Brook on 23/10/2014 20:24:57.

PCCP

59 T. Shiozaki, C. Woywod and H.-J. Werner, Phys. Chem. Chem. Phys., 2013, 15, 262. 60 A. A. Granovsky, J. Chem. Phys., 2011, 134, 214113. 61 A. A. Granovsky, Firefly version 7.1.G, http://classic.chem. msu.su/gran/firefly/index.html. 62 M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis and J. A. Montgomery, J. Comput. Chem., 1993, 14, 1347. 63 H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, ¨tz, P. Celani, T. Korona, R. Lindh, M. Schu A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, ¨ppl, Y. Liu, A. W. Lloyd, R. A. Mata, G. Jansen, C. Ko

This journal is © the Owner Societies 2014

Paper

64 65 66 67 68 69

A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, ¨ger, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflu R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson and M. Wang, MOLPRO, version 2012.1, a package of ab initio programs, 2012. M. R. Silva-Junior, M. Schreiber, S. P. A. Sauer and W. Thiel, J. Chem. Phys., 2010, 133, 174318. M. Schreiber, M. R. Silva-Junior, S. P. A. Sauer and W. Thiel, J. Chem. Phys., 2008, 128, 134110. F. Wang, S. P. Neville, R. Wang and G. A. Worth, J. Phys. Chem. A, 2013, 117, 7298. S. P. Neville and G. A. Worth, J. Chem. Phys., 2014, 140, 034317. ¨ppel, L. S. Cederbaum and G. A. Worth, H.-D. Meyer, H. Ko I. Burghardt, Int. Rev. Phys. Chem., 2008, 27, 569. J. C. Light, I. P. Hamilton and J. V. Lill, J. Chem. Phys., 1985, 82, 1400.

Phys. Chem. Chem. Phys., 2014, 16, 15957--15967 | 15967