Applying the vibronic coupling model Hamiltonian to the photoelectron spectrum of cyclobutadiene

Chemical Physics 329 (2006) 99–108 www.elsevier.com/locate/chemphys

Applying the vibronic coupling model Hamiltonian to the photoelectron spectrum of cyclobutadiene S. Saddique, G.A. Worth

*

School of Chemistry, University of Birmingham, Birmingham B15 2TT, UK Received 11 May 2006; accepted 28 June 2006 Available online 12 July 2006

Abstract The vibronic coupling model Hamiltonian has been set up to provide global potential energy surfaces for the lowest three singlet states of the neutral molecule and the doubly-degenerate ground state of the radical cation. Both systems are dominated by non-adiabatic coupling between the states, and the model is able to describe the anharmonic adiabatic surfaces using simple coupled diabatic states. The model is fitted to electronic structure calculations at the CASSCF level. To test the model, the photoelectron spectrum is calculated and compared to experiment. It is found that while the simple model is able to reproduce qualitatively the known topology of the potential energy surfaces, it is unable to accurately reproduce the photoelectron spectrum. It is concluded that this is due to inaccuracies in the electronic structure calculations.  2006 Elsevier B.V. All rights reserved. Keywords: Vibronic coupling; Wavepacket dynamics; Cyclobutadiene photoelectron spectrum

1. Introduction Cyclobutadiene, C4H4, is a classic case in which vibronic coupling plays a crucial role in its structure and stability. At first glance it would be expected that this molecule will have a square planar, D4h, structure with an aromatic psystem. However, as there are four electrons in the p-system, the Hu¨ckel 4n + 2 rule expects this to be unstable. Calculations [1–4] and experiment [5] find that in fact the equilibrium geometry is rectangular, D2h. This has been explained by vibronic coupling between the ground state ~ 1 A1g at D4h [1], ~ 1 B1g and A and the first-excited state, X which induces symmetry breaking and leads to a stable structure with lower symmetry. Vibronic coupling also plays a role in the structure of the cyclobutadienyl radical cation, C4 Hþ 4 . The ground state at D4h is now 2Eg, and this structure is again subject to symmetry breaking. This time the coupling is classified as E  b Jahn–Teller, and the branching space for the symme*

Corresponding author. E-mail address: [email protected] (G.A. Worth).

0301-0104/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2006.06.040

try lowering is provided by vibrational modes with b1g and b2g symmetry [2,6]. In this paper, we use the vibronic coupling model Hamiltonian [7] to set up potential energy surfaces for both the neutral molecule and radical cation. To test the model Hamiltonians, quantum dynamics calculations are then used to calculate the photoelectron spectrum, which can be compared to the experimental one. Vibronic coupling is due to the coupling between electronic and nuclear motion, and gives rise to pathways for radiationless transitions between electronic states [8–10]. A common result is the existence of a conical intersection between adiabatic potential energy surfaces [11,12]. Systems in which vibronic coupling plays a role are inherently multi-mode problems. At least two degrees of freedom must be present, and in the majority of interesting photochemical systems many nuclear modes are coupled. The vibronic coupling model provides a simple, yet powerful, framework for describing multi-dimensional potential energy surfaces of molecules where these non-adiabatic effects are important. A simple Taylor expansion of the diabatic potential surfaces and couplings is made around a

100

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

suitable point, taking symmetry into account. This approach, which uses the fact that diabatic states are smooth, is able to reproduce the complicated topology seen in the adiabatic picture often using a surprisingly small number of terms. Developed by Cederbaum and co-workers [13,7], the model has been used to great effect particularly in modelling photoelectron spectra such as butatriene [14,13], benzene [15], allene [16] and small heterocycles [17]. The study of the butatriene photoelectron spectrum can claim to be the first instance where the importance of a conical intersection was shown in a non-Jahn–Teller molecule. The model has also been used to describe other systems where a conical intersection is present, an example being the absorption spectrum of pyrazine [18]. Vibronic coupling problems are also inherently dynamical problems, and require quantum dynamics calculations for a good treatment. This is a significant problem as standard quantum dynamics methods, such as grid-based wavepacket propagation [19–21] are severely restricted with respect to the size of system that can be treated. Usually only 3–4 nuclear modes can be treated. For this reason the multi-configuration time-dependent Hartree (MCTDH) method plays an important role [22– 24]. Also originating in the Cederbaum group, this is an extension to wavepacket propagation methods able to treat many modes. For example, we have been able to treat the dynamics of pyrazine after excitation to the S2 state as it passes through a conical intersection to the S1 state including all 24-vibrational degrees of freedom [18]. The Hamiltonian for this calculation was provided by the vibronic coupling model. Thus in combination the model and the MCTDH method allow a system passing through a conical intersection to be investigated in detail. 2. Methods 2.1. The vibronic coupling model The Hamiltonian for a set of N coupled states can be written as an N · N matrix, H. Assuming that a diabatic basis is being used, this can be written H ¼ Hð0Þ þ Wð0Þ þ Wð1Þ þ    ;

ð1Þ

where H(0) is a diagonal matrix containing the kinetic energy operators and a potential. For a bound-state problem harmonic oscillators provide a suitable zero-order picture, ! X x a o2 ð0Þ 2 H ii ¼ þ Qa ; ð2Þ 2 oQ2a a where the sum is over the nuclear coordinates, Qa, which are dimensionless (mass-frequency scaled) normal coordinates. The subsequent matrices include the effects of electronic excitation and vibronic coupling as a Taylor expansion around a point Q0, the zero point of the normal modes. Thus W(0) has matrix elements

ð0Þ

W ij ¼ hUi ðQ0 ÞjH el jUj ðQ0 Þi;

ð3Þ

where Hel is the standard clamped nucleus electronic Hamiltonian and Ui the diabatic electronic functions. This diabatic basis is defined by the condition that it transforms the singular non-adiabatic operator in the adiabatic Hamiltonian into a local potential-type operator [8,25]. The ‘‘global gauge’’ must also be fixed by making the diabatic basis equivalent to the adiabatic functions at a point. Choosing this point as Q0, the matrix W(0) simply contains the adiabatic excitation energies on the diagonal. The linear coupling matrix elements can be written X oH el ð1Þ W ij ¼ hUi ðQ0 Þj jUj ðQ0 ÞiQa : ð4Þ oQa a Thus the diagonal elements are the forces, and the off-diagonal elements the non-adiabatic coupling at Q0. Higher order terms may be added, but experience has shown that the major physics of many systems can be described by this simple linear model. This is especially the case when a conical intersection is present. Symmetry can now be used to further simplify and add a descriptive element to the model. An integral is zero unless the integrand is totally symmetrical with respect to the symmetry of the molecule. Thus an element of ð1Þ the linear coupling matrix, W ij;a , is zero unless the product of the symmetries of the two states and the vibrational mode contains the totally symmetric irreducible representation: Ci  Ca  Cj  Ag :

ð5Þ

Thus terms linear on the diagonal of the linear coupling matrix involve only symmetric vibrations, and, if the two states have different symmetries, only terms involving a particular non-symmetric irrep appear on the off-diagonal. The expansion coefficients in the coupling matrices are the parameters describing the potential energy surfaces. These can be calculated by evaluating derivatives of the surfaces at Q0. A more satisfactory approach is to choose the parameters so that the model optimally fits the calculated surfaces. As in previous applications on butatriene and allene [14,26], this was done by calculating points along the various normal modes and optimising the parameters. A set of programs, called VCHAM, have been developed to do this [27]. A database of points and corresponding energies is set up. The Hamiltonian parameters are then calculated optimising the least-squares fit function F ¼

X

 2 wi V calc  V mod ; i i

ð6Þ

i

where V calc and V mod are the ab initio calculated and model i i potential energies at point i. The latter are the eigenvalues of the diabatic Hamiltonian matrix at the point. The weight function   ð7Þ wi ¼ exp 2V calc i

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

101

ensures that the fit is best in the more important low energy regions. The initial guess used the neutral molecule frequencies along with the state energies and gradients at Q 0.

Further details of the method are in recent reviews [32,24].

2.2. Calculating the photoelectron spectrum

3.1. Neutral cyclobutadiene

A spectrum can be calculated from the Fourier transform of the autocorrelation function Z 1 IðxÞ  x dtCðtÞeixt ; ð8Þ

Previous studies show that the equilibrium geometry of cyclobutadiene is a rectangular, D2h, structure with a square planar, D4h, transition state on the reaction path for the interconversion between the two equivalent forms. The equilibrium structure was calculated and characterised by a normal mode analysis using RHF, MP2 and CASSCF methods with a 6-31G* basis set. The active space taken for the CAS calculations used four electrons distributed in the four p-orbitals. The orbitals are shown schematically in Fig. 1. All calculations used the GAUSSIAN program [33]. The frequencies are shown in Table 1 and compared to experimental values. Symmetry labels correspond to the molecule in the standard orientation with the molecule in the yz-plane and the z-axis bisecting the shorter pair of bonds. The RHF values are clearly inadequate, with errors of over 300 cm1 in some modes. Even the low frequency modes have significant errors. In general, MP2 provides reasonable values. Excluding the four high frequency modes, errors are less than 100 cm1. The exception is 1b3g with an experimental value of 723 cm1 compared to a calculated value of 858.7 cm1. The CASSCF values sit between the MP2 and RHF in terms of quality. To set up the vibronic model Hamiltonian, the energy of neutral cyclobutadiene was minimised in a square planar, D4h, geometry using the CASSCF method. The molecule is defined with the C2 axes bisecting the bonds, and C 02 bisecting the atoms. If this definition is reversed, the labels b1 and b2 swap over. The convention chosen here matches that used by Nakamura et al. [1]. This geometry is taken as Q0, the point around which the expansion is made. The frequencies and normal modes at this geometry were then calculated. These are the nuclear coordinates to be used in the Hamiltonian. The D4h minimum energy geometry is a transition state between two rectangular,

1

where CðtÞ ¼ hWð0ÞjWðtÞi ¼

      t  t W : W 2 2

ð9Þ

The initial wavefunction, W(0), is formed by a vertical excitation, i.e. the ground-state neutral wavefunction is placed on a potential energy surface in the manifold of states describing the radical cation. This is the Condon approximation, and is equivalent to assuming that the nuclei do not move as an electron is ejected. The nuclear wavefunction at time t, W(t), is then obtained by solving the time-dependent Schro¨dinger equation, allowing the initial wavefunction to evolve over the set of surfaces. Note that the right-hand side of Eq. (9) allows one to obtain an autocorrelation function twice the length of the propagation time [28,29]. The solution of the time-dependent Schro¨dinger equation was done using the Heidelberg MCTDH package [30]. This implements the highly efficient MCTDH wavepacket propagation method [23,24] in a convenient computer program. The method takes as a wavefunction ansatz X ð1Þ ðf Þ WðQ1 . . . Qf ; tÞ ¼ Aj1 ...jf uj1 ðQ1 ; tÞ . . . ujf ðQf ; tÞ; ð10Þ j1 ...jf

writing the function as a sum of products of low-dimensional basis functions. Equations of motion for the basis functions, called single-particle functions, and the expansion coefficients are then set up using a variational method. The power comes from the fact that both the coefficients and the basis functions are time-dependent and evolve so that the MCTDH wavefunction optimally reproduces the numerically exact solution of the Schro¨dinger equation. Flexibility is also available in the ansatz as the singleparticle functions may be functions of more than one-dimension. In the example below, two-dimensional functions were used, with two system degrees of freedom ‘‘coupled together as a single-particle’’. As a result the number of coefficients required if n basis functions are used for f degrees of freedom increases as nf/2 rather than nf. A huge saving of effort. The functions themselves are described by a set of onedimensional discrete variable representation (DVR) functions [31]. Here, a harmonic oscillator DVR was used with frequencies matching the zero-order oscillator frequencies of the problem.

3. Results

Fig. 1. The four p molecular orbitals, /1  /4 of cyclobutadiene defining the active space and the three configurations |1i, |2i and |3i which can be used to describe the lowest three singlet states.

102

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

Table 1 Frequencies of neutral cyclobutadiene calculated at the D2h minimum energy geometry compared to experimental results Mode

MP2

CAS(4, 4)

RHF

Experimenta

1au 1b2g 1b3u 2au 1b2u 1b1g 1b3g 1ag 1b1u 2ag 2b3g 2b2u 3ag 2b1u 3b3g 3b2u 3b1u 4ag

472.0 490.4 553.9 709.9 749.8 777.2 858.7 991.5 1080.1 1149.5 1201.3 1291.7 1600.7 1615.0 3263.8 3278.0 3298.1 3307.6

535.1 432.4 524.1 734.4 852.6 698.2 966.3 999.3 1130.8 1215.0 1307.4 1410.5 1513.3 1639.5 3402.5 3417.3 3430.7 3445.1

607.2 708.2 682.8 1012.0 809.7 995.9 925.0 1034.5 1162.2 1219.3 1302.8 1404.9 1782.6 1806.9 3398.5 3415.3 3438.1 3449.0

– 531 576 – 721 – 723 989 1028 1059 – 1245 1678 1526 3093 3107 3124 3140

All values in cm1. All calculations use a 6-31G* basis set. a Ref. [36], except 1b1u Ref. [37].

D2h minima and a single imaginary frequency was found along a mode with b1g symmetry. The modes and frequencies are listed in Table 2. The transition mode is m4. A correlation with the D2h is also given. The lowest energy electron configuration in the p-orbital space is a22u e2g , which can give rise to 1B1g, 1A1g, 1B2g and 3 A2g. The energies of the lowest three singlet electronic states were calculated using a state-averaged CASSCF ~ 1 B1g and the first and wavefunction. The ground-state is X ~ 1 A1g and B ~ 1 B2g . Using the configusecond excited states A rations shown schematically in Fig. 1, the dominant config-

Table 2 Definition of the modes used in the vibronic coupling model of cyclobutadiene Mode

Symmetry

x (cm1)

D2h equivalent

x (cm1)

m1 m2 m3 m4a m5 m6 m7 m8a m8b m9 m10 m11 m12a m12b m13a m13b m14a m14b

1a1g 2a1g 1a2g 2b1g 1b1g 1b2g 2b2g 1eg 1eg 1a2u 1b1u 2b1u 1eu 1eu 2eu 2eu 3eu 3eu

0.165 0.429 0.164 0.176 0.144 0.126 0.424 0.049 0.049 0.062 0.062 0.074 0.124 0.124 0.181 0.181 0.426 0.426

2ag 4ag 2b3g 3ag 1ag 1b3g 3b3g 1b2g 1b1g 1b3u 1au 2au 1b2u 1b1u 2b2u 2b1u 3b2u 3b1u

0.151 0.427 0.162 0.188 0.124 0.120 0.422 0.054 0.087 0.065 0.066 0.091 0.106 0.140 0.175 0.203 0.424 0.425

Column 2 uses the symmetry defined at D4h, and column 4 is the corresponding mode at D2h (see Table 1). a Transition mode with imaginary frequency.

Table 3 Energies of lowest lying states of neutral cyclobutadiene as calculated using a state-averaged CASSCF with a 6-31G* basis set and an active space of four electrons in four orbitals State

Dominant configurationsa

E (CAS) (eV)

E (EOM-CC) [4] (eV)

~ 1 B1g X ~ 1 Ag A ~ 1 B2g B

|1i  |3i |1i + |3i |2i

0.000 2.230 3.840

0.000 1.814 2.137

a

Defined in Fig. 1.

urations are listed in Table 3. Due to the degeneracy of the eg orbitals, other choices of molecular orbitals are possible (see e.g. [1]) and this would affect the description of the configurations. The ones chosen here correlate directly with the orbitals at D2h. The state energies are also given in Table 3. The CAS calculations give vertical excitation energies of EA~ ¼ 2:23 eV and EB~ ¼ 3:84 eV. As a comparison, the results from high level calculations using equation-ofmotion coupled cluster techniques give values of 1.81 and 2.14 eV [4]. The vibronic model Hamiltonian, Eq. (1), can now be set up for neutral cyclobutadiene. The zero-order Hamiltonian is the set of harmonic oscillators defined in Eq. (2) with the frequencies in Table 2. For the transition mode the calculated imaginary frequency was taken as real. The matrix of state energies is 0 1 0 EX~ 0 B C W ð0Þ ¼ @ 0 EA~ 0 A ð11Þ 0

0

EB~

with the energies as listed in Table 3. The linear coupling matrix W(1) can only have on-diagonal elements due to totally symmetric modes. As the point of expansion is at a minimum for all these modes, these terms are all zero for the ground state. They are also found to be negligible for the two excited states being considered. After a consideration of the symmetries of the states and ~ states are coupled by modes ~ and A normal modes, the X ~ and B ~ states by modes with a2g with b1g symmetry, the X ~ and B ~ states by modes with b2g symmetry. symmetry and A Using the mode numbering scheme in Table 2, the nonzero linear coupling matrix elements are thus X ð01Þ ð1Þ ka Q a ; ð12Þ W 12 ¼ a¼4;5

W W

ð1Þ 13 ð1Þ 23

ð02Þ

¼ k3 Q 3 ; X ð12Þ ¼ ka Q a :

ð13Þ ð14Þ

a¼6;7

W(1) is, like all the coupling matrices, symmetric. A final consideration must be made for the quadratic and bilinear terms. Neglecting the ungerade and doubly degenerate modes, which do not appear in the linear model and so play a minor role in the physics, useful information can be obtained by correlating the modes at D4h with those

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

at D2h (see Table 2). Modes m2, m3, m6 and m7 correlate directly with a single modes without change in frequency. Thus they have negligible second-order parameters. In contrast, on going to D2h modes m1, m4 and m5 mix together to become the first three ag modes. This Duschinsky rotation is an important feature in describing molecules as they change geometry and must be included. The important terms in the second-order matrix are X X ð0Þ ð2Þ 2 W 11 ¼ cð0Þ cab Qa Qb ; ð15Þ a Qa þ a¼1;4;5

W

ð2Þ 33

¼

X

a¼4;b¼5 2 cð1Þ a Qa

þ

a¼1;4;5

¼

X

X

ð1Þ

ð16Þ

cab Qa Qb ;

a¼4;b¼5 2 cð2Þ a Qa ;

ð17Þ

ð01Þ

ð18Þ

a¼1;4;5 ð2Þ

W 12 ¼

X

lab Qa Qb :

a¼1;b¼4;5

All other parameters are zero by symmetry or found to be negligible. The parameters were calculated by the fitting procedure described above. At 10 points along various cuts, both along and between modes where required, the energies of the three states were calculated using state-averaged CASSCF calculations, giving a database of energies at 101 points to which the model Hamiltonian was fitted. The final values are given in Table 4. The adiabatic surfaces from the model and the ab initio points are shown along the three most important modes, m1, m4 and m6. The modes themselves are shown in Fig. 2. In all cases the fit is seen to be very good. Along mode m4 ~ (kð01Þ ) leads ~ and A the strong vibronic coupling between X 4 to the double minima on the ground state. This is the symmetry breaking mode, with the minima at D2h and the transition state at D4h. Along mode m6, weak coupling between Table 4 Parameters for the vibronic coupling model Hamiltonian of neutral cyclobutadiene k

c

ð02Þ

¼ 0:000

ð01Þ

¼ 0:509

k3

m4

k4

(0.583)

~ and B ~ (kð12Þ ) leads to the widening of the firststates A 6 excited adiabatic state potential well, but the coupling is not strong enough to break the symmetry. Along m1 an error is present due to a significant anharmonicity of the ab initio potential surfaces that cannot be fitted by the model potentials (see Fig. 3). Fig. 4 shows two-dimensional cuts through the lowest adiabatic surface as contour plots to demonstrate the effect of the second-order terms. A separate fit including only linear terms in the model was also made. This results in ð01Þ slightly different values of k4 ¼ 0:583 eV and ð01Þ k5 ¼ 0:057 eV to compensate for the missing terms. (a) and (c) are in the space of modes m4 and m6: the modes that provide the strong vibronic coupling. The double minima are clear. (b) and (d) show the surface in the space of the ~ In both cases, the sec~ and A. b1g modes which couple X ond-order terms make the wells in (c) and (d) wider and deeper. The quality of the fit can be measured by the root-meansquare displacement (RMSD) between the CASSCF energies and the values from the model at the points. For the linear model, the value is 0.099 eV. For the second-order model, this value is 0.069 eV. These errors are small. The critical points on the ground-state surface are the minima positions and the barrier height. From the CASSCF calculations, the D2h minima are at (Q1, Q2, Q4, Q5) = (0.215, 0.015, ±3.147, ±0.414) with an energy of 0.382 eV below the D4h barrier. From the second-order model, the minima

l

ð0Þ

m1

m3

Fig. 2. The most important vibrational modes of cyclobutadiene: (a) m1(a1g), (b) m4(b1g) and (c) m6(b2g).

ð01Þ

c1 ¼ 0:011

l1;4 ¼ 0:018

7.0

ð1Þ c1

ð01Þ l1;5

6.0

¼ 0:014

¼ 0:023

ð1Þ

c4 ¼ 0:049 ð1Þ

c4;5 ¼ 0:049

energy [eV ]

W

ð2Þ 22

103

a

b

c

-4 -2 0 2 4 Q4

-4 -2 0 2 4 Q6

5.0 4.0 3.0 2.0

ð2Þ

c4;5 ¼ 0:042 m5

ð01Þ

k5

¼ 0:031

(0.057)

ð0Þ

c5 ¼ 0:021 ð2Þ

k6

ð12Þ

¼ 0:250

m7

ð12Þ k7

¼ 0:079

0.0

ð1Þ

c5 ¼ 0:023 c5 ¼ 0:021

m6

1.0

The expansion is around the D4h minimum energy geometry. Values for ð01Þ ð01Þ k4 and k5 in brackets are for the pure linear model. All values in eV.

-4 -2 0 2 Q1

4

Fig. 3. Cuts through potential energy surfaces for the lowest three singlet ~ 1 Ag and B ~ 1 B1g , A ~ 1 B2g of neutral cyclobutadiene along the D4h states X vibrational modes (a) m1(1ag), (b) m4(2b1g) and (c) m6(1b2g). Points are energies calculated using SA-CAS(4,4)/6-31G*, lines are from the vibronic coupling model Hamiltonian.

104

a

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

b3

4

2 1 Q5 0 -1 -2 -3

2 Q6 0 -2 -4 -8 -6 -4 -2 0 2 4 6 8 Q4

-8 -6 -4 -2 0 2 4 6 8 Q4

To decide which of the b-modes go on the diagonal, and which off-diagonal in the linear coupling matrix, symmetry labels using the D2h point group, the largest non-abelian subset of the full D4h, are used. Under this reduction of symmetry, b1g becomes ag while b2g becomes b3g. Thus the linear coupling matrix has the elements X X ð1Þ ja Qa þ ja Qa ; ð19Þ W 11 ¼ a¼1;2

c

d

W

3 2 1 Q5 0 -1 -2 -3

4 2

Q6 0 -2 -4 -8 -6 -4 -2 0 2 4 6 8 Q4

ð1Þ 22

¼

X

a¼4;5

ja Qa 

a¼1;2 ð1Þ

W 12 ¼

X

X

ð20Þ

ja Qa ;

a¼4;5

ð21Þ

ka Qa

a¼6;7

-8 -6 -4 -2 0 2 4 6 8 Q4

Fig. 4. Contour plots of the potential energy surface for the lowest ~ 1 B1g , of neutral cyclobutadiene in the space of (b), adiabatic singlet state, X (d) the b1g modes, m4 and m5, and (a), (c) the strongest vibronically coupled modes, m4 and m6. Surfaces are from the vibronic coupling model Hamiltonian fitted to SA-CAS(4,4)/6-31G* data: (a) and (b) include only linear coupling, (c) and (d) including second-order terms.

~ in the neu~ and A with b1g modes, those that couple the X ~ tral molecule, on-diagonal, and b2g, those that couple the A ~ and B in the neutral molecule, off-diagonal. The totally symmetrical modes also appear on the diagonal, with the same coupling constant for both states in contrast to the b1g Jahn–Teller active modes which have different signs. The non-negligible second-order terms, taking account of symmetry, are X ð2Þ ca Q2a þ c1;4 Q1 Q4 þ c4;5 Q4 Q5 ; ð22Þ W 11 ¼ a¼2;4;5;7

are at (Q1, Q2, Q4, Q5) = (0.496, 0.000, ±3.052, ±0.551) with a barrier height of 0.181 eV. In the linear model, the minima are at (Q1, Q2, Q4, Q5) = (0.100, 0.072, ±2.724, ±0.325) with the same barrier height of 0.181 eV, very similar to the second-order model minima. The b2g modes are both at zero in all cases. In the model the sign is wrong on the symmetric mode Q1. This is probably due to the anharmonicity of the potential along these modes which is hard to fit with a quadratic potential. The double well minimum is in the correct place along the b1g modes, although the barrier is a bit flat. 3.2. Cyclobutadiene radical cation The vibronic coupling model Hamiltonian for the radical cation was obtained as for the neutral molecule. The expansion is around the neutral molecule D4h equilibrium geometry, and the same normal modes were used as coordinates. Thus the zero-order Hamiltonian is the same set of harmonic oscillators. The linear, quadratic, and bilinear parameters were again obtained by fitting to points calculated using CASSCF calculations. The same active space was used, now with three electrons distributed in the four orbitals. At D4h, the state is 2Eg, with a single electron in the doubly degenerate orbital. This degeneracy can be lifted along modes with b1g and b2g symmetry, in a E  b Jahn–Teller interaction. The first excited state is calculated to lie 4.55 eV higher in energy. The pseudo-Jahn–Teller coupling that takes place between the doubly degenerate state and this singly degenerate state has been shown to be important in a full analysis of the properties of the radical cation [6], but this is a weak effect.

W

ð2Þ 22

¼

X

ca Q2a þ c1;4 Q1 Q4 þ c4;5 Q4 Q5 :

ð23Þ

a¼2;4;5;7

The parameters for the model after fitting to a suitable database of points are listed in Table 5. Parameters both for the purely linear model (fitted without second-order terms) and the full second-order model are shown. The RMSD of the linear and second-order models are 0.075 and 0.066 eV, respectively. Fig. 5 shows cuts through the potential energy surface demonstrating how well the model fits the ab initio points along the two main Jahn–Teller active modes, m4 and m6. The much stronger coupling along the m4 modes means that this dominates, and as the contour plot shows the double minima occur along this mode either side of Q0. At the origin, (Q4, Q6) = (0, 0), is the conical intersection point with the cone formed by the upper adiabatic surface (see Table 5 Parameters for the vibronic coupling model Hamiltonian of the cyclobutadiene radical cation Mode

j

m1 m2 m4 m5 m6 m7

0.038 0.173 0.301 0.030

k

j

k

0.065 0.147 0.275 0.024 0.115 0.035

c 0.033 0.011 0.029

0.115 0.039

0.039

c1,4 =  0.025 c4,5 =  0.033 Ecat = 7.170 The expansion is around the D4h neutral molecule minimum energy geometry. The second and third columns are for the purely linear model. Columns 4–6 are for the second-order model. All values in eV.

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

a

b

-4 -2 0 2 4 Q4

-4 -2 0 2 4 Q6

12.0

energy [eV]

11.0 10.0 9.0 8.0 7.0

c

4 2

Q6

0 -2 -4 -8 -6 -4 -2 0 2 4 6 8 Q4

~ 2 Eg ground-state Fig. 5. Cuts through potential energy surfaces for the X of the cyclobutadiene radical cation along the D4h vibrational modes (a) m4(2b1g) and (b) m6(1b2g). Points are energies calculated using SACAS(4,4)/6-31G*, lines are from the vibronic coupling model Hamiltonian. (c) Contour plot of the lower adiabatic surface from the model Hamiltonian, in the space of m4 and m6.

Fig. 1 Ref. [8] for a 3D representation of an E  b Jahn– Teller surface). The D2h minima in the second-order model lie at (Q1, Q2, Q4, Q5) = (0.238, 0.370, ±1.795, ±0.481) with a Jahn–Teller stabilisation energy of 0.279 eV. In the linear model, this is (Q1, Q2, Q4, Q5) = (0.023, 0.403, ±1.710, ±0.205) with a Jahn–Teller stabilisation energy of 0.295 eV. These values should be compared to the CASSCF minimum geometry, which lies at (Q1, Q2, Q4, Q5) = (0.350, 0.067, ±1.736, ±0.234) with a Jahn–Teller stabilisation energy of 0.327 eV. The seemingly large difference on Q2 is related to the high frequency of this mode. A large change in coordinate does not correspond to a large change in geometry. 3.3. The photoelectron spectrum The photoelectron spectrum was calculated from wavepacket dynamics simulations using the MCTDH method. The first step was to obtain the ground-state nuclear wavefunction for the neutral molecule. This was done by applying energy relaxation (propagation in imaginary time) [24,34] to a guess wavepacket using the vibronic coupling model Hamiltonian for the neutral molecule. The initial

105

guess was taken as the ground-state harmonic oscillator eigenfunction of the zero-order Hamiltonian. While the functions in the diabatic states have density in both wells, the lowest eigenfunction created by adding the parts is localised in one well, centred at (Q1, Q2, Q4, Q5, Q6, Q7) = (0.91, 0.18, 4.33, 0.85, 0.00, 0.00). The next step is to perform the dynamics of the radical cation after its formation. A vertical excitation, which corresponds to the instantaneous removal of an electron with no relaxation of the nuclear framework, was used to form the initial wavepacket. This involves taking the lowest energy eigenfunction of the neutral molecule and placing it in one of the diabatic states of the cation. Note that the eigenfunction formed by energy relaxation is in a diabatic representation, and all components need to be simultaneously excited. The photoelectron spectrum was then obtained from the subsequent dynamics, again solved using the MCTDH program, from the Fourier transform of the autocorrelation function. Three calculations were made. One used the full secondorder Hamiltonians and included six modes (m1, m2, m4, m5, m6, m7), which comprise the set of modes with significant linear coupling in both the neutral molecule and radical cation. Pairs of modes were combined together to form particles, thus effectively reducing the problem to a threemode calculation. The number of basis functions used are listed in Table 6. Note that a total of 25.9 · 106 primitive basis functions are being used for each state. This makes the calculation impossible for standard grid-based wavepacket methods, but the calculations are fairly straightforward for MCTDH – the largest calculation made needed 45 MB of memory and took less than half an hour on a standard linux workstation. The second calculation used the same model Hamiltonians, but only included the two most important Jahn–Teller active modes m4 and m6. The third included all six modes, but used the purely linear model Hamiltonian. These calculations help to show the importance of various terms in the model. The photoelectron spectra are shown in Fig. 6. The wavepacket has been propagated for 150 fs, producing an autocorrelation function of 300 fs. To allow for the finite length of the propagation, the autocorrelation functions has been multiplied by an exponential damping function with a time constant of s = 100 fs. The spectra have also Table 6 Number of basis functions used in the MCTDH wavepacket propagation calculations used to simulate the photoelectron spectrum of cyclobutadiene mode

(ni, nj)

N0, N1, N2

N1, N2

m4, m6 m1, m5 m7, m2

50, 18 10, 12 20, 12

4, 4, 4 4, 4, 4 4, 4, 4

8, 8 8, 8 8, 8

ni are the no. of primitive basis functions used to describe each mode. Ni are the number of single-particle functions used for the wavepacket on each state. Column 3 refers to the neutral molecule, column 4 the radical cation.

106

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

Intensity

a

8

8.5

9

9.5

10

9.5

10

9.5

10

Intensity

b

8

8.5

9

Intensity

c

8

8.5

9 Energy [ev]

Fig. 6. The photoelectron spectra of cyclobutadiene obtained from wavepacket calculations using (a) second-order six-mode model, (b) linear six-mode model (c) second-order two-mode model.

been shifted so that the first peak is close to the experimental peak at 8.16 eV [5]. This involved the addition of 1.0 eV to the six-mode model spectra and 0.5 eV to the two-mode model spectrum. The three spectra differ significantly, and show how sensitive the spectrum is to the model and parameters used. The six-mode second-order spectrum, (Fig. 6a), is a broad, structured band with 11 peaks approximately 0.18 eV apart. The most intense peak is the third. The two-mode model spectrum (Fig. 6c) has far less structure, which shows that all six modes are excited in the ionisation, and the higher peaks are very weak. Only seven peaks are really clear, and the second peak is the strongest. The linear model spectrum in Fig. 6b is simpler still. Only three peaks are clearly visible, and the strongest peak is the first. Thus without the second-order terms little excitation of the modes takes place. The second-order terms play only a marginal role in the radical cation model. They, however, make a significant difference for the neutral molecule surfaces as can be seen by the positions of the minima described above. A brief analysis of the dynamics can be made by calculating the expectation values of the coordinates as a function of time. Some results from the six-mode second-

order model are shown in Fig. 7. The a1g and b1g modes m1, m4 and m5 are all seen to oscillate with a large amplitude: they are highly excited by the removal of the electron. Interestingly, the strongly coupled mode m4 is seen to remain trapped in the initial well as it barely moves to positive values. It also is not damped in the same way as m1 and m5. The off-diagonal b2g mode m6 has a very different behaviour. The oscillations grow as time progresses and energy flows into this mode. As can be seen, though, the excitation of this mode is small. The diabatic state populations are shown in Fig. 8. There is an oscillatory transfer of population between the two components of the degenerate state which has a period around 20 fs. Comparing this with the expectation values in Fig. 7 this corresponds to the frequencies of the three strongly excited modes, which all have similar periods. Comparing the results to the published experimental spectra [5,35] the best fit is Fig. 6c, the two-mode model. The experimental spectrum is much less well resolved, but approximately seven peaks can be distinguished under the envelope, and the second peak is the most intense. The discrepancy of the more complete model with respect to this is due to the exact location of the neutral molecule minima compared to the cation minima. These are very similar in the linear model, hence little excitation occurs, but quite different along all the modes in the second-order six-mode model. A more significant difference is the spacing between the peaks. The experimental gap [5] is 0.08 eV (645 cm1), under half that found in the model. The width of the experimental spectrum is accordingly narrower. This provides a mystery. The simulated spectrum of Ref. [5] does not report the frequencies of their model, but the gap appears to be about 0.13 eV, also wider than 0.08 eV. There are two sources of error present in the calculations presented here. The first is the simple nature of the reduced dimensionality model used: harmonic diabatic states and couplings only up to second-order. The second is the level of electronic structure theory used to calculate the adiabatic potential energy surfaces. The modes not included in the model all couple weakly and will play only minor roles in the spectrum. Furthermore, the quality of the fit along the included modes means that the frequencies from the model must be very similar to those of the original. Thus the discrepancy between the experimental and calculated spectrum must be due not to the model but the underlying ab initio potential energy surfaces. Looking at Table 1 there are a few vibrational modes close to this frequency. Ignoring the ungerade modes which cannot be directly excited by this transition, m8a, m8b and m6 are possible candidates. The pair degenerate at D4h, however, also cannot be excited directly, leaving m6. The calculations at the level used here firmly indicate that its Jahn–Teller coupling is not strong enough for its progression to dominate over those from the b1g and a1g modes. This may be a failure of the method. A clue might be that the calculated frequency for this mode at D2h has a large

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

a

1

b

0.5 0

-1

Q4

Q1

-0.5

-1.5 -2 -2.5

C

0

20

40

60

80 100 Time [fs]

120

140

160

d

1.5

Q6

Q5

0.5

20

40

60

80 100 Time [fs]

120

140

160

120

140

160

0.004 0.002 0 -0.002 -0.004 -0.006

-0.5 -1 -1.5

0

0.01 0.008 0.006

1

0

1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5

107

0

20

40

60

80

100

120

140

160

-0.008 -0.01

0

20

Time [fs]

40

60

80

100

Time [fs]

~ state from the neutral molecule. Taken from a Fig. 7. Expectation values of the coordinates for the cyclobutadiene radical cation after formation in the X wavepacket propagation calculation using a six-mode second-order vibronic coupling model Hamiltonian (a) m1(1a1g), (b) m4(2b1g), (c) m5(1b1g) and (d) m6(1b2g).

1

course be possible, but might lead to a physically incorrect model. Further calculations at a higher level are required to resolve this problem.

0.8 0.6

4. Conclusions

0.4 0.2 0

0

20

40

60

80 100 time [fs]

120

140

160

~ state of the Fig. 8. Diabatic state populations for the components of the X cyclobutadiene radical cation after formation from the neutral molecule. Taken from a wavepacket propagation calculation using a six-mode second-order vibronic coupling model Hamiltonian.

error, even at the MP2 level. High level calculations (CCSD(T)/cc-pVTZ) [2] also indicate that the rectangular structure is not the global minimum, but a puckered rhomboid–rhomboid due to distortion along b2g (b1g in [6]) followed by pseudo-Jahn–Teller coupling puckering the ring. This indicates that indeed the coupling between the states along m6 is being underestimated. A final possibility is that the coupling between the modes is also not correctly described and the frequency of m1 may be lowered to match that seen. Given the choices, it is not sensible to manually optimise parameters to fit the spectrum. This would of

The calculation of potential energy functions for multidimensional systems is a serious bottleneck in the study of molecular dynamics. The most satisfactory approach is to fit analytic functions to the results of high level ab initio electronic structure calculations. If such a function is able to reproduce experimental data then we can be sure that the physical interpretation from the simulations is correct. The vibronic coupling model Hamiltonian provides a suitable functional form for the description of systems in which vibronic coupling is important. Experience shows that often a very simple structure in the diabatic picture used in the model is able to describe very anharmonic potential surfaces in the adiabatic picture. A second bottleneck is the need for an accurate simulation method able to treat multi-dimensional systems. Nonadiabatic effects such as conical intersections are now widely accepted as being important in a range of photophysical and photochemical processes. These require quantum dynamical simulation for a complete treatment. Few methods can be applied to systems of the size interesting to chemists. The MCTDH method is at present the best

108

S. Saddique, G.A. Worth / Chemical Physics 329 (2006) 99–108

suited for accurate quantum calculations on large systems. This is true especially in combination with the vibronic coupling model Hamiltonian. In the system studied here, cyclobutadiene, the model Hamiltonian has been parameterised so that it fits the adiabatic potential energy surfaces at the CASSCF/6-31G* level for both the neutral molecule and the radical cation. It is able to reproduce the topological features of these surfaces in an analytic form suitable for dynamics simulations. To test the accuracy of the model, it was compared to the measured photoelectron spectrum. Unfortunately, the dynamics calculations do not agree with experiment in a satisfactory way. The ‘‘correct’’ answer could be obtained by suitable adjustment of the parameters, or even by using a reduced dimensionality (the result from a two-dimensional model resembles the experimental spectrum closer than the one from a sixdimensional model). There is, however, a number of ways the adjustment could be done, and each would give a different interpretation of the spectrum. For example to get the gap between spectral peaks right a change could be made to the frequency of a strongly-coupled mode, or else to the coupling of a mode with the correct frequency. More accurate calculations are required to distinguish which is the right way. In conclusion, the combined power of the model Hamiltonian with an automated fitting procedure together with the MCTDH method able to provide a test of the model, provides a set of tools suitable for the study of photoinduced processes. Acknowledgement This paper is dedicated to Prof. Lenz Cederbaum on the occasion of his 60th birthday. With thanks for his continued support and friendship and in appreciation of his scientific insight and endless enthusiasm. References [1] K. Nakamura, Y. Osamura, S. Iwata, Chem. Phys. 136 (1989) 67. ˇ a´rsky, Chem. Phys. Lett. [2] M. Roeselova´, T. Bally, P. Jungwirth, P. C 234 (1995) 3295. [3] A. Metropolous, Y.-N. Chiu, J. Mol. Struct. (Theochem.) 365 (1996) 119. [4] S.V. Levchenko, A. Krylov, J. Chem. Phys. 120 (2004) 175. [5] D.W. Kohn, P. Chen, J. Am. Chem. Soc. 115 (1993) 2844. [6] I.B. Bersuker, Chem. Rev. 101 (2001) 1067. [7] H. Ko¨ppel, W. Domcke, L.S. Cederbaum, Adv. Chem. Phys. 57 (1984) 59. [8] G.A. Worth, L.S. Cederbaum, Ann. Rev. Phys. Chem. 55 (2004) 127. [9] F. Bernardi, M. Olivucci, M.A. Robb, Chem. Soc. Rev. 25 (1996) 321. [10] M. Klessinger, Angew. Chem. 107 (1995) 597.

[11] W. Domcke, D.R. Yarkony, H. Ko¨ppel (Eds.), Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, Singapore, 2004. [12] D.R. Yarkony, Rev. Mod. Phys. 68 (1996) 985. [13] L.S. Cederbaum, W. Domcke, H. Ko¨ppel, W. von Niessen, Chem. Phys. 26 (1977) 169. [14] C. Cattarius, G.A. Worth, H.-D. Meyer, L.S. Cederbaum, J. Chem. Phys. 115 (2001) 2088. [15] M. Do¨scher, H. Ko¨ppel, Chem. Phys. 225 (1997) 93. [16] S. Mahapatra, G.A. Worth, H.-D. Meyer, L.S. Cederbaum, H. Ko¨nppel, J. Phys. Chem. A 105 (2001) 5567. [17] A.B. Trofimov, H. Ko¨ppel, J. Schirmer, J. Chem. Phys. 109 (1998) 1025. [18] A. Raab, G. Worth, H.-D. Meyer, L.S. Cederbaum, J. Chem. Phys. 110 (1999) 936. [19] K.C. Kulander (Ed.), Time-dependent Methods for Quantum Dynamics, Elsevier, Amsterdam, 1991. [20] R. Kosloff, J. Phys. Chem. 92 (1988) 2087. [21] C. Cerjan (Ed.), Numerical Grid Methods and their Application to Schro¨dinger’s Equation, Kluwer Academic Publishers., Dordrecht, 1993. [22] H.-D. Meyer, U. Manthe, L.S. Cederbaum, Chem. Phys. Lett. 165 (1990) 73. [23] U. Manthe, H.-D. Meyer, L.S. Cederbaum, J. Chem. Phys. 97 (1992) 3199. [24] M.H. Beck, A. Ja¨ckle, G.A. Worth, H.-D. Meyer, Phys. Rep. 324 (2000) 1. [25] M. Baer, Chem. Phys. 259 (2000) 123. [26] A. Markmann, G.A. Worth, L.S. Cederbaum, J. Chem. Phys. 122 (2005) 144320. [27] C. Cattarius, A. Markmann, G.A. Worth, The VCHAM program, see http://www.pci.uni-heidelberg.de/tc/usr/mctdh/ (2005). [28] V. Engel, Chem. Phys. Lett. 189 (1992) 76. [29] U. Manthe, H.-D. Meyer, L.S. Cederbaum, J. Chem. Phys. 97 (1992) 9062. [30] G.A. Worth, M.H. Beck, A. Ja¨ckle, H.-D. Meyer, The MCTDH Package, Version 8.3, see http://www.pci.uni-heidelberg.de/tc/usr/ mctdh/, 2004. [31] J.C. Light, I.P. Hamilton, J.V. Lill, J. Chem. Phys. 82 (1985) 1400. [32] H.-D. Meyer, G.A. Worth, Theor. Chem. Acc. 109 (2003) 251, feature article. [33] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery, Jr., K.N. Kudin, J.C. Burant, J.M. Millam, R.E. Stratmann, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, S. Iyengar, G.A. Petersson, M. Ehara, K. Toyota, H. Nakatsuji, C. Adamo, J. Jaramillo, R. Cammi, C. Pomelli, J. Ochterski, P.Y. Ayala, K. Morokuma, P. Salvador, J.J. Dannenberg, 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. Gomperts, 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, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, GAUSSIAN 03, 2003. [34] R. Kosloff, H. Tal-Ezer, Chem. Phys. Lett. 127 (1986) 223. [35] J. Kreile, N. Mu¨nzel, A. Schweig, H. Specht, Chem. Phys. Lett. 124 (1986) 140. [36] B. Arnold, J. Michl, J. Phys. Chem. 97 (1993) 13348. [37] G. Meier, Angew. Chem. 100 (1988) 317.