First principles simulation of the UV absorption

JOURNAL OF CHEMICAL PHYSICS

VOLUME 121, NUMBER 5

1 AUGUST 2004

First principles simulation of the UV absorption spectrum of ethylene using the vertical Franck-Condon approach Anirban Hazra Department of Chemistry, Princeton University, Princeton, New Jersey 08544 and Guelph-Waterloo Centre for Graduate Work in Chemistry, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1

Hannah H. Chang Department of Chemistry, Princeton University, Princeton, New Jersey 08544

Marcel Nooijena) Guelph-Waterloo Centre for Graduate Work in Chemistry, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1

共Received 29 April 2004; accepted 12 May 2004兲 A new method which we refer to as vertical Franck-Condon is proposed to calculate electronic absorption spectra of polyatomic molecules. In accord with the short-time picture of spectroscopy, the excited-state potential energy surface is expanded at the ground-state equilibrium geometry and the focus of the approach is more on the overall shape of the spectrum and the positions of the band maxima, rather than the precise position of the 0-0 lines. The Born-Oppenheimer approximation and the separability of the excited-state potential energy surface along the excited-state normal mode coordinates are assumed. However, the potential surface is not necessarily approximated as harmonic oscillator potentials along the individual normal modes. Instead, depending upon the nature of the potential surface along a particular normal mode, it is treated either in the harmonic approximation or the full one-dimensional potential is considered along this mode. The vertical Franck-Condon approach is applicable therefore even in cases where the excited state potential energy surface is highly anharmonic and the conventional harmonic Franck-Condon approach is inadequate. As an application of the method, the ultraviolet spectrum of ethylene between 6.2 eV 共50 000 cm⫺1兲 and 8.7 eV 共70 000 cm⫺1兲 is simulated, using the Similarity Transformed Equation of Motion Coupled-Cluster method to describe the required features of the potential energy surfaces. The spectrum is shown to be a result of sharp doublet structures stemming from the ␲ →3s 共Rydberg兲 state superimposed on top of a broad band resulting from the ␲→␲* 共valence兲 state. For the Rydberg state, the symmetric CvC stretch and the torsion mode contribute to the spectrum, while the broad valence band results from excitation into the CvC stretch, CH2 scissors, and the torsion mode. For both states, the potential along the torsion mode is highly anharmonic and the full treatment of the potential along this mode in the vertical Franck-Condon method is required. © 2004 American Institute of Physics. 关DOI: 10.1063/1.1768173兴

I. INTRODUCTION

tween the two approaches, we will refer to the former as adiabatic FC and the latter as vertical FC. The suggestion of a vertical FC model, additionally assuming the same force constants for ground- and excited-states, has been made many years ago by Domcke et al.7 based on a moment analysis of the intensity distribution in the spectrum.8 The assumption of linearly displaced harmonic potentials to describe excited-state potentials is also widely used to simulate resonance Raman spectra following work by Tannor and Heller.9 In this work we will consider this alternative FC approach but including more sophisticated approximations for the excited-state potential energy surface than displaced harmonic oscillators. Of the two approaches, the vertical FC approximation is expected to give a better description of the broad features 共position of intensity maximum, overall width, and envelope兲 of the spectrum than the adiabatic FC approach. This can be understood based on Heller’s time-dependent perspective on molecular spectroscopy.9–11 Let 兩 ␹ 0 典 be the ground-state vi-

Ever since its enunciation in the 1920s, the FranckCondon 共FC兲 approximation1–3 has played an important role in the field of electronic spectroscopy. It is used both for interpretation of experimental spectra and in ab initio simulations. The procedure that is frequently used to calculate electronic absorption spectra involves optimizing the ground- and excited-state geometries, making a harmonic approximation 共second-order Taylor expansion兲 of the groundand excited-state potential energy surfaces at their respective minimum energy geometries and evaluating the FC overlap factors between the ground- and excited-states to construct the spectrum 共for example Refs. 4 – 6兲. Instead of optimizing the excited state geometry, one might expand the excited state potential energy surface in a Taylor series around the equilibrium geometry of the ground state. To distinguish bea兲

Electronic mail: [email protected]

0021-9606/2004/121(5)/2125/12/$22.00

2125

© 2004 American Institute of Physics

Downloaded 07 Nov 2006 to 134.99.82.26. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

2126

J. Chem. Phys., Vol. 121, No. 5, 1 August 2004

Hazra, Chang, and Nooijen

brational wave function and ␮ the electronic transition moment between the ground- and excited-electronic states. Then writing 兩 ␾ 0 典 as 兩 ␾ 0典 ⫽ ␮ 兩 ␹ 0典 ,

共1兲

the absorption spectrum is given by

⑀ 共 ␻ 兲 ⫽C ␻

冕

T

⫺T

e i 共 ␻ ⫹E 0 兲 t 具 ␾ 0 兩 ␾ 共 t 兲典 dt,

共2兲

where C is a constant, ␻ is the frequency of incident radiation, E 0 is the energy of 兩 ␹ 0 典 , and 兩 ␾ (t) 典 satisfies the timedependent Schro¨dinger equation iប

⳵兩␾共 t 兲典 ⫽H ex兩 ␾ 共 t 兲典 ⳵t

共兩 ␾ 共 0 兲典 ⬅ 兩 ␾ 0 典 ),

共3兲

where H ex is the vibrational Hamiltonian for the excited state.10 The broad features of the spectrum are obtained by integrating Eq. 共2兲 over a short time T. By propagating the wave packet 兩 ␾ (t) 典 for longer times T more details in the structure of the spectrum emerge, but the envelope of the spectrum remains essentially unchanged.10 Since, in practice we do not use the full excited-state Hamiltonian H ex, the approximation that we use for it is critical in determining the spectrum. If we approximate the Hamiltonian by writing the potential as a Taylor series around the equilibrium geometry of the ground state such as in the vertical FC approach, then the difference between the effective potential and exact potential will be least in the region where the wave function 兩 ␾ (t) 典 has its largest value at a short time. This is because at time equal to zero, 兩 ␾ (t) 典 is a Gaussian function with its center at the equilibrium geometry of the ground state. Thus the broad features of the spectrum are expected to be obtained more accurately from the vertical FC approach than from an adiabatic FC calculation. On the other hand, specific low-lying vibrational lines in the spectrum, in particular the 0-0 transition may be better described using adiabatic FC as here the potential around the excited state equilibrium geometry is most critical. Though physically appealing for calculating spectra, the vertical FC approach in the harmonic approximation runs into difficulties if one or more of the excited-state frequencies is imaginary at the ground-state equilibrium geometry. Along a normal mode with an imaginary frequency, a harmonic approximation to the potential energy surface then gives a downward parabola corresponding to a potential that tends to negative infinity. A Hamiltonian with such a potential is clearly unphysical and does not have real eigenfunctions. It is a poor approximation to the actual potential corresponding to a bound state and cannot be used to calculate a spectrum in a meaningful way. In reality, the second derivative of such a potential surface must change sign with change in geometry and the surface is therefore highly anharmonic. In molecules, negative force constants most often occur for nonsymmetric modes indicating that the equilibrium geometry of the excited state has a lower symmetry than the ground state. The low-lying excited states of ethylene are a typical example, where the molecule breaks the D 2h molecu-

lar symmetry of the ground state along the torsion coordinate12 to achieve a stable minimum in these excited states. In the present paper, we present a method to calculate electronic spectra using the vertical FC approximation, which can be used even in cases when the excited-state potential has negative force constants at the ground-state equilibrium geometry. In this method, the separability of the excited-state potential along its normal mode coordinates, defined using the force constant matrix at the ground-state geometry, is still assumed. For those normal modes that correspond to real frequencies, the harmonic approximation is made for the potential while for the modes that correspond to imaginary frequencies, the complete one-dimensional potential in the directions of those normal modes are calculated to construct the excited-state vibrational Hamiltonian. While calculating the one-dimensional potential as a function of a particular 共imaginary frequency兲 normal mode coordinate, the remaining coordinates are frozen at their values corresponding to the ground-state geometry. This is consistent with the real frequency modes, where the ground-state geometry is taken as the reference point in the second-order Taylor expansion while making the harmonic approximation for these modes. This approach is motivated by the shorttime picture of spectroscopy and aims to provide an accurate representation of the potential energy surface in the Condon region. The vertical FC approach discussed here thus retains some of the simplicity of the harmonic approximation, while introducing more sophistication when this approximation is inadequate. Using the above form of the excited-state Hamiltonian, FC overlap integrals are calculated. The procedure for calculating FC overlaps is not separable in the normal modes because the ground- and excited-state normal modes are not parallel in general. The ground- and excited-state normal mode coordinates are related by a linear transformation, the so called Duschinsky rotation.13–15 Since the excited-state force constant matrix is defined at the groundstate geometry, the eigenvectors that describe overall translation and rotation are identical in both states. Therefore axis-switching effects,14,16 which make the transformation between the ground- and excited-state vibrational normal mode coordinates nonlinear, do not arise in the present approach. Because of the Duschinsky rotation, calculating FC overlaps involves evaluating multidimensional overlap integrals over functions that are products of partly harmonic, and partly completely general one-dimensional functions. In our implementation of the vertical FC method presented here, this is done in an efficient manner as will be described below. The vertical FC method is applied to calculate the electronic absorption spectrum of ethylene between 6.2 and 8.7 eV which corresponds to the first broad band in the overall spectrum. This energy range includes the ␲ →3s 共Rydberg兲 and ␲→␲* 共valence兲 states and has been examined by several researchers, e.g. 共Refs. 4, 12, and 17–31兲. Various ab initio simulations of the spectrum have been carried out, but have not been fully successful in reproducing the spectrum.4,30–32 Mebel, Chen, and Lin 共MCL兲 have calculated the spectrum by optimizing excited-state geometries,4 that is, using the adiabatic FC approach. Their calculations



showed the critical role of dynamical correlation. At the complete active space self consistent field 共CASSCF兲 level with the 6-311共2⫹兲G* basis set, which includes primarily nondynamic correlation, the calculated geometry of the Rydberg state is planar. However, at the equation of motion coupled-clustered singles and doubles 共EOM-CCSD兲 level of theory, the Rydberg state breaks symmetry along the torsion mode and the twisting angle after freezing the other geometric parameters at the CASSCF optimized values is found to be about 20°. The subsequent FC calculations by MCL used a combination of the CASSCF and EOM-CCSD data to simulate the spectrum. Still the adiabatic FC simulation of the Rydberg state spectrum was not fully satisfactory as MCL note that the energy gaps between the major peaks in the calculated spectrum correspond to the excitation of the torsion ( ␯ 4 ) vibrational mode and are evenly spaced, while in the experiment the vibrational progression shows a series of doublet structures. From this study it is clear that the harmonic potential along the torsion coordinate provides a poor representation of the actual potential energy surface. For the valence state, MCL have used the D 2d geometry in which the torsion angle is equal to 90°. In this case the equilibrium structure in the valence excited state is displaced enormously compared to the ground-state geometry. For such cases the adiabatic Franck-Condon approximation is dubious and it seems more promising to describe the broad valence band starting from a short-time dynamics perspective, as the envelope of the spectrum is the quantity of interest, not so much the individual lines. Baeck and Martinez31 have used the ab initio multiple spawning method 共AIMS兲,33,34 using a quasidiabatic scheme to follow the adiabatic states and have calculated the spectrum from the time-correlation function using a small number of trajectories that covered about 50 fs of dynamics. The vibrational dynamics was calculated on potential energy surfaces that were generated on the fly at the EOM-CCSD/631⫹G level. Nonadiabatic dynamics plays a minor role in this part of the spectrum. The authors point out that their calculated Rydberg and valence state intensities are weighted equally while comparing the total spectrum to experiment, but in principle, that each spectrum should be weighted by its oscillator strength. They further note that the major difference between the calculated and observed spectra is the peak doubling attributed to the Rydberg state, which is seen in the experimental spectrum but is absent in their calculated spectrum. For the valence state, Baeck and Martinez31 and before Ben-Nun and Martinez30 concluded in their work that the main vibrational progression is due to the symmetric stretch. Moreover, the time-correlation function in their work rapidly decays and does not show any recurrences provided that the trajectories are not constrained to be planar and sample the torsion motion. It appears that the torsion mode plays an important role in suppressing the recurrences in the timecorrelation function, resulting in a broad featureless spectrum, which is indeed the main characteristic of the valence band. In this work we calculate the absorption spectrum of ethylene from first principles using the vertical FC method and good agreement is obtained with experiment. The elec-

UV absorption spectrum of ethylene

2127

tronic structure calculations use the similarity transformed equation of motion coupled-cluster singles and doubles 共STEOM-CCSD兲 method.35–38 In both the Rydberg and the valence state, all modes are treated in the harmonic approximation except for the torsion mode where we use the full ab initio potential. Since the torsion mode is the same in the ground- and excited-states due to symmetry, we can analyze the contribution to the spectrum from the torsion and the other modes separately and this will be used to interpret the spectrum. This paper is organized as follows: Sec. II describes the theory of the vertical FC method for a general excited-state potential energy surface, which may have negative force constants at the ground-state geometry. Section III discusses the results from the calculation of the absorption spectrum of ethylene using the vertical FC method and in Sec. IV we summarize our conclusions. II. THEORY A. Vibrational Hamiltonian for the excited electronic state and overlap integrals for calculating spectra

The calculation of the positions and intensities of lines in the electronic absorption spectrum involves obtaining the vibrational energies and the overlaps between the vibrational eigenfunctions of the ground and excited electronic states. For this we need the vibrational Hamiltonians of the two states. First we define the coordinates. Let us suppose that A and Ã are the ground and excited electronic states in question where the ⬃ symbol is associated with the excited state. The ground-state normal mode coordinates are denoted by Q i (i ⫽1, N兲. For the excited state, the normal mode coordinates are defined for a model surface obtained by making a Taylor expansion using the vertical excitation energy, gradient, and Hessian of the excited state surface at the ground-state geometry. Suppose that f is the Hessian matrix and g is the gradient vector in the ground-state normal mode coordinates. Let U be the unitary matrix containing eigenvectors of f as columns. The excited-state normal mode coordinates denoted ˜ (i⫽1, N兲 are related to the ground-state normal coorby Q i dinates by N

˜ ⫽ Q j

兺 S ji Q i ⫹d j ,

共4兲

i⫽1

where 共5兲

S⫽UT and d⫽UTf⫺1 g⫽ 共 UTf U兲 ⫺1 UTg.

共6兲

In Eq. 共4兲, S is an N⫻N Duschinsky matrix that describes the rotation between the ground- and excited-state normal modes and d is a vector of length N that describes the displacement of the excited-state equilibrium geometry, which is the origin of the excited-state normal coordinate system. From Eq. 共6兲, we see that the eigenvalues of f which are elements of the diagonal matrix (UTf U) and the gradient in the terms of the eigenvectors of f, i.e., UTg determine the 13–15


2128



displacement d of the excited-state normal coordinates with respect to the ground-state ones. If the Hessian matrix f has negative eigenvalues then the model surface with negative force constants will not have a well defined equilibrium geometry. In this case, we set the components of the vector UTg corresponding to the negative eigenvalues of f to zero to define d in Eq. 共6兲 and thus the excited-state normal coordinates 关Eq. 共4兲兴. We assume that the potential energy functions of the states are separable in their respective normal mode coordinates. In the usual harmonic FC method, the ground- and excited-state potentials are approximated as quadratic 共or harmonic兲 potentials along each normal mode. In the present method, however, if the excited-state potential along certain normal modes is poorly represented by the harmonic approximation, then along those modes, FC factors are calculated using the full one-dimensional potential. The groundstate potential on the other hand is assumed to be harmonic as before. For the sake of distinction, we will refer to the normal modes along which the excited-state potential is approximated to be harmonic as harmonic normal modes and those along which the potential is anharmonic as anharmonic nor˜h mal modes. We denote the harmonic normal modes as Q ˜ a . A integer subscript labels and the anharmonic ones as Q ˜ a corresponds to the the particular mode. So 共for instance兲 Q i ith mode which happens to be anharmonic. Let us assume there are M anharmonic modes out of a total of N modes. Since the excited-state Hamiltonian is assumed to be separable in the normal modes, we partition this Hamiltonian into two parts—one which depends on the anharmonic normal modes and the other on the harmonic ones. The former is a sum of one-dimensional Hamiltonians where the potential is a function of the particular anharmonic normal mode coordinate while the latter is simply a sum of one-dimensional harmonic oscillator Hamiltonians. M

˜⫽ H

N

˜h 兺兵 KE共 Q˜ ai 兲 ⫹V i共 Q˜ ai 兲其 ⫹ j⫽M兺⫹1 H HO j 共 Q j 兲,

i⫽1

共7兲

KE represents the kinetic energy operator and V i is the excited-state potential energy as function of the normal co˜ a and is obtained by restricting all the other norordinate Q i mal coordinates to have values that are fixed at the groundis the one-dimensional harmonic state geometry. H HO j ˜h. oscillator Hamiltonian in the coordinate Q j The eigenfunctions of this multidimensional but sepa˜ have the following product form: rable Hamiltonian H ˜ ⫽兩␾ ˜ a 1共 Q ˜ a兲 ␾ ˜ a 2共 Q ˜ a兲¯␾ ˜ aM共 Q ˜ a 兲 ˜␹ h M ⫹1 共 Q ˜h 兲 ⌿ 1 2 M M ⫹1 1 2 M M ⫹1 h ⫹2 ˜ h h ˜h ⫻ ˜␹ MM⫹2 共 Q M ⫹2 兲 ¯ ˜␹ NN 共 Q N 兲典,

共8兲

˜ a i (Q ˜ a ) for i⫽1, M and ˜␹ h j (Q ˜ h ) for j⫽M ⫹1, N are where ␾ i j j i eigenfunctions of one-dimensional Hamiltonians in anharmonic and harmonic normal modes, respectively. That is, for i⫽1, M a

a

a

˜ i共 Q ˜ a 兲典 ˜ a 兲 ⫹V 共 Q ˜ a 兲其兩␾ ˜ i共 Q ˜ a 兲典 ⫽˜⑀ i 兩 ␾ 兵 KE共 Q i i i i i i i i a i ⫽0,1,2,...,

共9兲

and for j⫽M ⫹1, N h

h

h

˜h ˜ j ˜h ˜ j ˜ j ˜h 兵 H HO j 共 Q j 兲其兩␹ j 共 Q j 兲典 ⫽ ⑀ j 兩␹ j 共 Q j 兲典 ˜⑀ ai i

h j ⫽0,1,2,...,

共10兲

˜⑀ hj j

and are the corresponding eigenvalues. For where convenience we introduce a shorter notation for the eigenfunctions a

˜ i共 Q ˜ a 兲典, 兩ã i 典 ⬅ 兩 ␾ i i 共11兲

h ˜h 兩˜h i 典 ⬅ 兩 ˜␹ j j 共 Q j 兲典.

˜ of Eq. 共8兲, we see that it can be For the product function ⌿ fully specified by a list of the individual one-dimensional ˜ as eigenfunctions in the product. We therefore denote ⌿ ˜ ⬅ 兩ã ,a ˜ ˜ ˜ ˜ M ,h ⌿ 1 ˜ 2 ,...,a M ⫹1 ,h M ⫹2 ,...,h N 典 ˜典. ˜ 2 ,...,a ˜ M ,h ⬅ 兩ã 1 ,a

共12兲

The absorption spectrum is calculated from the overlaps of the above eigenfunctions with the ground-state vibrational eigenfunction. The vibrational Hamiltonian of the ground state is a sum of one-dimensional harmonic oscillator Hamiltonians in the ground-state normal mode coordinates. Using similar notation to that discussed above for the excited state, the ground-state vibrational eigenfunction ⌿ is written as h

h

h

⌿⫽ 兩 ␹ 1 1 共 Q h1 兲 ␹ 2 2 共 Q h2 兲 ¯ ␹ NN 共 Q Nh 兲典 ⬅ 兩 h 1 ,h 2 ,...,h N 典 ⬅ 兩 h典 .

共13兲

For low temperatures, the harmonic oscillator eigenfunctions along all the modes correspond to the ground vibrational state. In our present simulations we use this condition, although this is not a fundamental limitation for the formalism. The ground state is then ⌿⫽ 兩 0,0,...,0 典 ⬅ 兩 0典 .

共14兲

The FC overlap integrals that we are interested in calculating are ˜ 典 ⫽ 具 0兩ã ,a ˜典, ˜ M ,h 具⌿兩⌿ 1 ˜ 2 ,...,a

共15兲

for different vibrational excitations in the excited electronic state corresponding to different integer values of ã 1 , ˜ . Note that ⌿ is a product ã 2 ,...,a ˜ M and ˜h M ⫹1 , ˜h M ⫹2 ,...,h N of functions of the ground-state normal mode coordinates, ˜ is a product of functions of the excited-state normal while ⌿ mode coordinates. The excited-state normal mode coordi˜ are linear combinations of the ground-state normal nates Q i mode coordinates Q i as indicated by Eq. 共4兲. The integrals in Eq. 共15兲 are thus N-dimensional integrals in the variables Q i ˜ (i⫽1, N). or Q i B. Evaluation of FC overlap integrals and calculation of the spectrum

If the one-dimensional excited-state potentials are harmonic along all the normal modes, then the multidimensional FC overlaps of Eq. 共15兲 can be efficiently calculated using a recursion relation 共for example, Ref. 39兲. In the present more general case the excited-state potential is not harmonic along all the normal modes. To calculate the overlaps efficiently,




we express the eigenfunctions along each anharmonic mode in a basis of harmonic oscillator eigenfunctions along that mode. We can thereby express the multidimensional excitedstate vibrational eigenfunction in an intermediate harmonic basis with which the overlap with the ground state can be calculated efficiently using a recursion relation. This provides an efficient mechanism to deal with Duschinsky rotation effects, even in the case of anharmonic potentials. n ã Let ˜␩ i i (Q i ) for n i ⫽0,1,2,... be harmonic oscillator eigenfunctions in the anharmonic normal mode coordinate ˜ a . The n in the superscript n is associated with a particular Q i i harmonic oscillator eigenfunction in the ith anharmonic normal mode coordinate. Then using notation developed earlier n ã 兩ñ i 典 ⬅ 兩 ˜ ␩ i i共 Q i 兲典.

共16兲

The eigenfunctions in the anharmonic modes can be written in terms of harmonic oscillator eigenfunctions 兩ã i 典 ⫽

兩ñ i 典具ñ i 兩ã i 典 , 兺 ñ

i⫽1, M .

共17兲

i

The multidimensional FC overlap integral of Eq. 共15兲 then becomes ˜ 典⫽ 具⌿兩⌿

兺

ñ1 ,n ˜ 2 ,...,n ˜M

˜ 典具ñ 兩ã 典 ˜ 2 ,...,n ˜ M ,h 具 0兩ñ 1 ,n 1 1

⫻具ñ 2 兩ã 2 典 ... 具ñ M 兩ã M 典 .

共18兲

˜ 典 is a product of harmonic os˜ 2 ,...,n ˜ M ,h The function 兩ñ 1 ,n cillator functions and for different integer values of ñ 1 , ˜ forms an intermediate ñ 2 ,...,n ˜ M and ˜h M ⫹1 , ˜h M ⫹2 ,...,h N harmonic basis in which the excited-state vibrational eigen˜ 典 is expressed. Equation 共18兲 gives a prescription function 兩 ⌿ ˜ 典 . We will now for calculating the FC overlap integrals 具 ⌿ 兩 ⌿ discuss how each piece on the right-hand side of Eq. 共18兲 is calculated and how they are all put together to calculate the overlap integrals and the spectrum. To begin with, we have the ground-state energy, frequencies, and normal modes and the excited-state vertical excitation energy, gradient, and Hessian or force constant matrix for the excited state evaluated at the ground-state equilibrium geometry. These are obtained from the electronic structure calculation. Using the excited-state Hessian and gradient, the excited-state normal mode coordinates are defined in the way described in Sec. II A. If a normal coordinate corresponds to a negative or very small force constant, we treat the normal mode as an anharmonic mode. To be able to define a set of harmonic eigenfunctions which constitute an intermediate basis in the anharmonic normal mode coordinates, we replace the negative 共or small兲 force constant by a positive number. This number is taken to be the ground-state force constant corresponding to the normal mode that has the largest overlap with the anharmonic mode. This particular choice of intermediate force constant is not necessary in principle, but leads to a relatively small number of multidimensional overlap integrals in the intermediate basis. We also set the component of the gradient along this normal mode direction to zero. The component of the gradient is already zero if the

2129

normal mode is not totally symmetric and it breaks the ground-state molecular symmetry in the excited state along that direction. We thus obtain a multidimensional harmonic oscillator potential in the excited-state normal mode coordinates. The multidimensional harmonic overlap integrals between the eigenfunctions of this harmonic potential and the ˜ 典 in Eq. 共18兲 ˜ 2 ,...,n ˜ M ,h ground vibrational state, i.e., 具 0兩ñ 1 ,n can now be calculated efficiently using a recursion relation.39 Next, we calculate the one-dimensional overlaps between the harmonic oscillator eigenfunctions and the actual eigenfunctions in the anharmonic normal mode coordinates, i.e., 具ñ 1 兩ã 1 典 , 具ñ 2 兩ã 2 典 ,..., 具ñ M 兩ã M 典 in Eq. 共18兲. The actual eigenfunctions are solutions of the Schro¨dinger equation in Eq. 共9兲. The potential energy in the Hamiltonian is a slice or scan of the multidimensional potential energy surface along the particular anharmonic normal mode. A scan of the potential energy surface is obtained by calculating the electronic energy of the excited state keeping all normal coordinates fixed at their ground-state value, except for the scanned coordinate. The potential energy function can have a general shape but will be bound, unless it involves dissociation. The latter case is not considered here. For a bound potential, we solve the one-dimensional Schro¨dinger equation by representing the Hamiltonian in a finite basis. The basis set that we use are Legendre polynomials.40 We take the range of the polynomials 共which is scaled as necessary兲 to be large enough so that we get a near zero value of the calculated eigenfunction at large geometrical distortions. The Legendre functions have the advantage that it is numerically stable to calculate their values within their domain using a recursion relation. Hermite polynomials are not very suitable for this purpose, as they rapidly take on very large values. Nevertheless, an issue with Legendre polynomials is that they do not satisfy the boundary conditions associated with the Hamiltonian. As a result, we have to be careful with the evaluation of the matrix representation of the kinetic energy operator. The Legendre basis representation of the usual form of the kinetic energy operator, i.e., ⫺d 2 /dx 2 is in general not a Hermitian matrix. Consider this kinetic energy operator sandwiched between two basis functions b 1 (x) and b 2 (x). The off-diagonal element of the matrix representation is K 1,2⫽⫺ ⫹

冕冕

x2

x1

dx b * 1 共x兲

x2

x1

dx

d 2b 2共 x 兲 dx 2

冋

⫽ ⫺b * 1 共x兲

db 1* 共 x 兲 db 2 共 x 兲 . dx dx

db 2 共 x 兲 dx

册

x2 x1

共19兲

* ⫽K 2,1 The kinetic energy matrix is Hermitian, that is K 1,2 only if

冋

⫺b 1* 共 x 兲

db 2 共 x 兲 dx

册

x2

⫽0.

共20兲

x1

This is not the case for our basis of Legendre polynomials, which do not satisfy the boundary conditions of the true eigenfunctions of the Hamiltonian. For this reason, we use an alternative expression for the kinetic energy matrix element,


2130


K 1,2⫽

冕

x2

x1

dx

db 1* 共 x 兲 db 2 共 x 兲 , dx dx

共21兲

which is explicitly Hermitian, and which is equivalent to the usual definition if the basis functions would satisfy the boundary conditions in Eq. 共20兲. We note that similar issues regarding the hermiticity of the kinetic energy operator may arise in the discrete variable representation, which is a particular kind of finite basis representation. A discussion is provided in Ref. 41. The integrals that involve the Legendre polynomials and the arbitrary potential are evaluated using a numerical integration scheme based on Legendre points and weights,40 scaled to the appropriate interval. The eigenfunctions in the Legendre basis, consisting typically of a 100 or so basis functions are obtained using diagonalization of the matrix representation of the Hamiltonian. Finally the overlap integrals involving the eigenfunctions of the general 1-D Hamiltonian and the intermediate harmonic oscillator functions are evaluated using again numerical integration. Since the intermediate basis typically does not exceed ten basis functions, numerical stability issues do not arise here. We have checked that our procedure yields proper eigenvalues and overlap integrals by comparing to the results from the program LEVEL Ref. 42, which uses a fully numerical representation of the eigenfunctions. We have now discussed how to obtain the different pieces on the right hand side of Eq. 共18兲 from which the FC overlap integrals, given by Eq. 共15兲, can be calculated. Calculating the overlap integrals for vibrational states, one by ˜ M , ˜h M ⫹1 , one, i.e., for a particular set of ã 1 , ã 2 ,...,a ˜h ˜ M ⫹2 ,...,h N , is not efficient. Rather, Eq. 共18兲 can be viewed as a transformation of the N-dimensional harmonic overlap integrals ˜典 ˜ 2 ,...,n ˜ M ,h 具 0兩ñ 1 ,n

共22兲

to the actual FC overlap integrals ˜典 ˜ 2 ,...,a ˜ M ,h 具 0兩ã 1 ,a

共23兲

through the transformation matrices

具 ñ1 兩ã1 典 , 具 ñ2 兩ã2 典 ,..., 具 ñM兩ãM典 .

共24兲

The efficient way of computer implementation of this transformation is to do it for one mode at a time and to transform all the harmonic eigenfunctions for a particular mode simultaneously to all the actual eigenfunctions for that mode. Let us consider the example of a particular normal mode, say mode 2 which we transform. We assume that mode 1 is already transformed to the actual eigenfunctions but none of the higher modes are transformed. If, for mode 2, there are B harmonic eigenfunctions, i.e., ñ 2 ⫽0,1,2,..., B⫺1 and A actual eigenfunctions, i.e., ã 2 ⫽0,1,2,..., A⫺1, then the transformation matrix is a B⫻A overlap matrix between the onedimensional harmonic oscillator eigenfunctions and the actual eigenfunctions. The transformation for mode 2 is then the multiplication of a vector of length B

冏冏


˜典 ˜ 2 ⫽0, ñ 3 ,...,n ˜ M ,h 具 0兩ã 1 ,n ˜典 ˜ 2 ⫽1, ñ 3 ,...,n ˜ M ,h 具 0兩ã 1 ,n

冏冏

共25兲

] ˜典 ˜ ˜ ñ 3 ,...,n ˜ M ,h 0 兩 a ,n ⫽B⫺1, 具 1 2

with the B⫻A matrix 具 ñ2 兩ã2 典 to give a vector of length A ˜典 ˜ 2 ⫽0, ñ 3 ,...,n ˜ M ,h 具 0兩ã 1 ,a ˜典 ˜ 2 ⫽1, ñ 3 ,...,n ˜ M ,h 具 0兩ã 1 ,a

] ˜典 ˜ 2 ⫽A⫺1, ñ 3 ,...,n ˜ M ,h 具 0兩ã 1 ,a

.

共26兲

From the calculated FC overlap integrals, the intensities in the spectrum can be calculated. The energy of the lines can be calculated keeping in mind that the energy is additively separable in the excited-state normal mode coordinates because the Hamiltonian 关Eq. 共7兲兴 is separable in these coordinates.

III. ELECTRONIC ABSORPTION SPECTRUM OF ETHYLENE A. Computational details

The electronic absorption spectrum of ethylene was calculated using the vertical FC approach. The input parameters for this approach are obtained from an electronic calculation and are the following. 共1兲 The ground-state energy, frequencies, and normal modes after optimizing the geometry. 共2兲 The excited-state energy, gradient, and force constant matrix at the ground-state equilibrium geometry. 共3兲 For every excited-state normal mode corresponding to an imaginary frequency, a one-dimensional scan of the excited-state potential energy surface obtained by changing only that particular normal mode coordinate and keeping the others fixed at the ground-state geometry. We obtain the quantities enumerated above using CCSD for the ground-state and the STEOM-CCSD35–38 method for the excited states. In the STEOM-CCSD approach, the excited states are calculated by diagonalizing a doubly transformed Hamiltonian over singly excited configurations. The computational cost of calculating the transformed Hamiltonian is in practice about twice the cost of a ground-state CCSD calculation, while the subsequent cost of calculating all of the excited states is essentially negligible. In the STEOM approach the calculation of a large set of excited states is highly efficient, and the gradients and force constant matrices calculated at the ground-state geometry for all excited states of interest are evaluated using finite difference techniques. As for the basis set, we noticed that the higher lying excited states of ethylene are particularly sensitive to the inclusion of diffuse basis functions. However, numerical instabilities occur if diffuse functions are included on all atoms in ethylene. The problems arise from the fact that significant overlap exists between diffuse functions, and the near-linear dependency between them results in convergence difficulties in the electronic structure calculation, in particular in the numerical evaluation of the force constant matrix. In order to




TABLE I. Exponents in the augmented PBS⫹diffuse basis set used for the ghost atom. s 0.015 0.005 0.001 667 0.000 555

p 0.015 0.005 0.001 667

TABLE III. Vibrational normal modes and vibrational frequencies of ethylene.

d 0.035 0.016 67 0.004 00

work around this problem, we created two ground-state models for ethylene. The Ghost model contains a ghost atom of large mass between the two carbons. A polarized basis set 共PBS兲 from Sadlej43 was used on the C and H atoms while a special set of diffuse functions containing four s, three p, and three d functions evenly tempered was used on the ghost atom. The exponents of this basis set are listed in Table I. The second model, referred to as the C 2 v model, has PBS basis functions on all atoms except for one of the carbons, for which a PBS⫹diffuse basis containing seven s, six p, and five d functions was applied.44 We found that both models achieve numerical stability in the electronic structure calculations and that the results obtained from all calculations including spectra simulations do not differ significantly between the two models. Thus for the purpose of the present discussion, the distinction between the two models can be treated as merely a technical detail. In what follows we will present results from the Ghost model which is more appealing for interpretative purposes as it retains D 2h symmetry. In the STEOM calculation, a total of 5 occupied orbitals and 41 virtual orbitals were included in the active space. In terms of the Hartree-Fock eigenvalues, all occupied orbitals above ⫺23 eV and all virtual orbitals below 7.5 eV were included in the active space. B. Results and discussion

1. Electronic structure data

The ethylene molecule has D 2h symmetry in the ground state. To assign the symmetry labels for the electronic states and normal modes, the z axis is taken along the CvC bond while the x axis is perpendicular to the molecular plane. The optimized ground-state geometry is given in Table II, along with experimental results. The 12 vibrational normal modes and their frequencies have been provided in Table III, also juxtaposed with experimental measurements. Vertical excitation energies have been reported in the literature for a variety of electronic structure methods. A small group of vertical energies has been selected to compare to our STEOM-CCSD results as shown in Table IV. It can be seen that there is general agreement across methods. The results from

Mode

Symmetry label

Description

Calculated frequency 共cm⫺1兲

Experimental frequency 共cm⫺1兲a

1 2 3 4 5 6 7 8 9 10 11 12

ag ag ag au b 3g b 3g b 3u b 2g b 2u b 2u b 1u b 1u

CH2 stretch CvC stretch CH2 scissors Torsion CH2 stretch CH2 rock CH2 wag CH2 wag CH2 stretch CH2 rock CH2 stretch CH2 scissors

3131 1675 1355 1034 3215 1225 987 954 3235 825 3107 1446

3022 1625 1344 1026 3083 1222 949 940 3105 826 2989 1444

a

Reference 56.

STEOM-CCSD are comparable to those obtained from EOM-CCSD and multireference configuration interaction MRD-CI. The CASPT2 value for the valence ␲→␲* ( 1 B 1u ) state is significantly higher, however, than that from other methods. This is likely due to unaccounted for Rydberg/ valence mixing, a problem in early CASPT2 calculations, that has since been resolved.45 Note that the STEOM-CCSD value of 7.78 eV for this 1 B 1u state is close to the best theoretical estimate of 7.7 eV by Mu¨ller, Dallos, and Lischka who performed extensive multireference CI singles and doubles 共MR-CISD兲 and multireference averaged quadratic coupled cluster 共MR-AQCC兲 calculations46 and with the MR-CISD value of 7.8 eV calculated by Ben-Nun and Martinez.47 It is also in agreement with the vibrationally corrected value of 7.8 eV calculated by Davidson and Jarzecki.48 In Table V, our calculated excited-state energies, oscillator strengths, spatial extent, i.e., 具 r 2 典 values, and assignments concerning valence or Rydberg character of the states is presented. In addition to the valence ␲→␲* transition, we also found another dipole forbidden valence ␴→␲* state of B 1g symmetry at 8.49 eV, which has not been obtained by other theoretical studies 共Table IV兲. It is possible that this state was missed by other calculations due to computational issues. For TABLE IV. Vertical excitation energies of ethylene in eV. ˜ ) CASPT2 EOM-CCSD EOM-CCSD(T State STEOM-CCSD 1 1 1 1 1 1 1

TABLE II. Ground-state geometry of ethylene: calculated and experimental results. Parameter r(CvC) r(CuH) ⬔共H-CvC兲 a

Reference 55.

Current calculation 1.348Å 1.098Å 121.3°

Experimenta 1.339Å 1.086Å 121.2°

2131

1 1 1 1 1

B 3u B 1u B 1g B 2g Ag B 1g B 3u B 3u B 2u Au B 3u B 1u

7.21 7.78 7.91 7.97 8.27 8.49 8.71 8.93 8.95 8.96 9.11 9.24

a

b

c

MRD-CId

7.28 7.98 7.94 7.99 8.45

7.10 7.74 7.76 7.80 8.28

7.17 8.40 7.85 7.95 8.40

7.13 7.96 7.86 7.89 8.21

8.79 9.08 9.26 9.02

8.61 8.90 9.08 8.85

8.66 9.03 9.18 8.94

8.73 9.31 8.99 9.04

9.30

9.13

9.31

a

Reference 57. Reference 57. c Reference 58. d Reference 25. b


2132



TABLE V. Assignment, spatial extent, and oscillator strength of calculated excited states in ethylene. R⬅Rydberg, V⬅valence.

State 1 1 1 1 1 1 1 1 1 1 1 1

B 3u B 1u B 1g B 2g Ag B 1g B 3u B 3u B 2u Au B 3u B 1u

STEOM-CCSD energy 共eV兲 7.21 7.78 7.91 7.97 8.27 8.49 8.71 8.93 8.95 8.96 9.11 9.24

具 r 2典

Assignment

共a.u.兲

Oscillator strength 共a.u.兲

R 关 b 3u ( ␲ )→a g (3s) 兴 V 关 b 3u ( ␲ )→b 2g ( ␲ * ) 兴 R 关 b 3u ( ␲ )→b 2u (3p y ) 兴 R 关 b 3u ( ␲ )→b 1u (3p z ) 兴 R 关 b 3u ( ␲ )→b 3u (3p x ) 兴 V 关 b 3g ( ␴ )→b 2g ( ␲ * ) 兴 R 关 b 3u ( ␲ )→a g 兴 R 关 b 3u ( ␲ )→a g 兴 R 关 b 3u ( ␲ )→b 1g 兴 R 关 b 3u ( ␲ )→b 3g 兴 R 关 b 3u ( ␲ )→a g 兴 R 关 b 3u ( ␲ )→b 2g 兴

121.3 98.7 135.4 142.3 154.9 89.8 172.6 186.9 190.3 201.3 298.1 278.7

8.00⫻10⫺2 3.67⫻10⫺1 0 0 0 0 2.41⫻10⫺3 5.41⫻10⫺2 3.01⫻10⫺2 0 1.45⫻10⫺2 2.38⫻10⫺2

example, the excitation energies of valence states are typically overestimated in EOM-CCSD 共Ref. 49兲 while in CASPT2 and MR-CI the active space may not have been chosen appropriately for this particular state. The STEOM method, on the other hand, automatically places the valence state at a presumably reasonable energy by implicitly including triple excitations, which are known to lower the valence state energies.38 2. Absorption spectra

In the energy range under study, i.e., 6.2– 8.7 eV, there are only two electronic states with nonzero oscillator strength, namely, the ␲ →3s ( 1 B 3u ) and ␲→␲* ( 1 B 1u ) states. In the literature these states are referred to as R and V states, respectively, and the ground electronic state is called N.12 The absorption spectrum is simulated for the R and V excited states. The R and V states are found to have one imaginary frequency along the direction of the torsion normal mode ( ␯ 4 ), and this mode is subjected to the special anharmonic treatment in the vertical FC approach. The scan of the potential energy curves along the Cartesian torsion normal coordinate for the R and V states with all other normal coordinates frozen to their value at the groundstate geometry is shown in Fig. 1. It should be realized that the torsion in the Cartesian coordinate system is not a pure rotation around the CvC bond for larger than infinitesimal distortion of the molecule. In fact, in Cartesian coordinates, the torsion coordinate corresponds to a displacement of the hydrogen atoms out of the molecular plane in a way such that one pair of trans-hydrogen atoms have equal displacement and the other pair have the same magnitude of displacement, but in the opposite direction. For large diplacements the potential energy curve corresponds to the simultaneous breaking of the four CH bonds, although this regime is not reached in the scan. The potential is not periodic, however, as it would be in the case of a rotational coordinate. We chose this coordinate to define the scan, rather than the torsional angle, as this is the normal mode in the ground state and the kinetic energy takes a simple form in terms of the orthogonal normal mode coordinates. The geometrical parameters of the molecule corresponding to the

FIG. 1. Energies of the R 共Rydberg; full lines兲 and V 共valence; dashed lines兲 excited states of ethylene as a function of the torsion normal mode coordinate. All other normal coordinates are frozen to their ground-state geometry.

minimum energy point on the R and V scans are shown in Table VI along with their ground-state values. As the molecule is distorted in the torsion coordinate, the CvC distance remains virtually unchanged but the torsion angle primarily changes and this is accompanied by changes in the C—H distance and H—CvC angle. These latter changes are small for the Rydberg state, but rather substantial for the valence state. The minimum in the potential energy curve for the valence state is determined by the CH stretching motion rather than the 90° minimum that would be expected for a variation along the torsion angle. The torsion mode in ethylene has a special property which makes the interpretation of the spectrum easier; the torsion mode is the only normal mode in the a u irreducible representation of the D 2h point group. Thus, as long as force constant matrices are calculated at the ground-state molecular geometry with D 2h symmetry, the torsion normal mode cannot mix with any other normal mode, and Duschinsky rotation does not involve the torsion mode. Within the vertical FC approximation the calculated spectrum is therefore identical to that obtained by convoluting the spectrum for only the torsion mode with the spectrum for all other modes. We will use this fact to interpret the spectrum. From the vertical FC calculation we get a line spectrum, where the height of a line corresponds to the intensity. Each of the lines is then broadened by a scaled Lorentzian function of a width of 0.02 eV 共161 cm⫺1兲 to obtain our final representation of the spectrum. The calculated spectrum for the ␲ →3s Rydberg state or R state is shown in the bottom panel of Fig. 2 while the corresponding experimental spectrum is shown in the top TABLE VI. Geometry of ethylene corresponding to minimum energy on the one-dimensional scan of the R and V state potential energy surfaces as shown in Fig. 1. Parameter

N state

R state

V state

r(CvC) r(CuH) ⬔共H—CvC兲 ⬔共Torsion兲

1.348Å 1.098Å 121.3° 0°

1.348Å 1.119Å 120.6° 25.9°

1.346Å 1.179Å 118.9° 49.4°




2133

FIG. 2. Absorption spectrum for the R ( 1 B 3u ) state of ethylene. The spectrum calculated using the vertical FC method 共bottom panel兲 is compared to the experimental spectrum 共Refs. 17 and 21兲共top panel兲.

panel. The calculated spectrum is shifted by about 1000 cm⫺1 with respect to the experimental one. It is seen that all features in the spectrum are rather well reproduced although the intensity ratio of the peaks that make up the repeated doublet structure is not quite correct. For the sake of interpretation of the spectrum, two further calculations were performed. In one case, the spectrum was calculated with only harmonic normal modes included in the simulation of the spectrum, while in the other case it was calculated with only the torsion mode or anharmonic mode included. These spectra are shown in Fig. 3. It is seen from the top figure that only the CvC stretch ( ␯ 2 ) contributes to the purely harmonic spectrum. The doubling of the peaks in the full spectrum comes from excitation in the torsion ( ␯ 4 ) mode. In addition the higher excitation in the torsion mode with four vibrational quanta is also clearly observed, and is identified also in the full spectrum of the Rydberg state. Note that, along the torsion mode, by symmetry only the vibrations with even number of quanta contribute to the spectrum. The frequencies of the individual lines are compared to the experimental values in Table VII. In the fourth column 共Z兲 of the table, for ease of comparison, the calculated frequencies shifted overall by 1005 cm⫺1 so that the 0-0 frequency from the calculation coincides with experiment, is shown. The last column 共Z-X兲 in the table is indicative of the error in the vibrational energy calculation. It is seen that the error grows steadily as we approach higher lying vibrational levels. The error in the symmetric modes is likely due to the neglect of anharmonicities. Regarding the vibrational structure due to the torsion mode, it is likely that the decoupling approxima-

FIG. 3. Absorption spectrum for the R ( 1 B 3u ) state of ethylene calculated with only Harmonic normal modes 共top panel兲 and only the torsion mode 共bottom panel兲.

tion of the normal modes gives rise to some of the error. Moreover it is possible that the Hamiltonian would provide a better representation of the full potential if the scan were made along the torsion angle, rather than along the Cartesian

TABLE VII. Experimental and calculated frequencies of lines in the spectrum of the Rydberg state.

a

Assignment ␯ 1␯ 2␯ 3␯ 4

Experiment X 共cm⫺1兲a

Calculation Y 共cm⫺1兲

Y⫹1005 Z 共cm⫺1兲

Z-X 共cm⫺1兲

0000 0002 0004 0100 0102 0104 0200 0202 0204 0300 0302 0304 0400 0402 0404

57 336 57 833 58 420 58 706 59 199 59 905 60 055 60 513 61 226 61 379 61 839 62 547 62 689 63 152 63 879

56 331 56 939 57 878 58 080 58 688 59 627 59 829 60 437 61 376 61 577 62 185 63 124 63 326 63 934 64 873

57 336 57 944 58 883 59 085 59 693 60 632 60 834 61 442 62 381 62 582 63 190 64 129 64 331 64 939 65 878

0 111 463 379 494 727 779 929 1155 1203 1351 1582 1642 1787 1999

Reference 17.


2134



FIG. 4. Absorption spectrum for the V ( 1 B 1u ) state of ethylene calculated using the vertical FC method.

normal mode. This may also contribute to the error in the relative intensities of the doublet pattern. Turning to the ␲→␲* valence or V state, we point out that there has been a long-standing controversy over assignment of the vibrational structure of this V state.21,27 Wilkinson and Mulliken first assigned the vibrational progression to the CvC stretch,17 although they provided also a detailed early discussion of the torsion mode, and these authors were well aware that it could not be represented by a harmonic potential. McDiarmid and Charney reassigned the structure to the torsion mode.18 Their analysis of the spectrum appeared consistent with earlier suggestions50 that the geometry of ethylene in the valence state is twisted by 90° along the CvC bond, so that the molecule has D 2d symmetry. This twisted equilibrium geometry was supported by theoretical calculations22,51 and continued to be the basis for the assignment of the spectrum52 for a long time. Quite recently BenNun and Martinez 共BM兲33,47 have argued that the D 2d geometry is actually not a minimum, but rather a saddle point. Furthermore BM assigned the vibrational progressions in the V state to the CvC stretch.31 Our calculated spectrum for the V state is shown in Fig. 4. The overall broad shape of the valence spectrum is well reproduced. The vibrational structure in the spectrum involves the CvC stretch ( ␯ 2 ), CH2 scissors ( ␯ 3 ), and the torsional ( ␯ 4 ) motion. This is more clearly seen by bisecting the spectrum into contributions from only harmonic modes and only the anharmonic torsion mode, as shown in Fig. 5. From the top figure and the assignments it is seen that both the CvC stretch and the CH2 scissors are highly activated in the transition to the valence state, and the width of the valence band is captured more or less by the progressions in the symmetric modes. However, this simulated harmonic spectrum involving primarily the symmetric modes has far more detailed structure than the experimental valence band. From the contributions in the torsion mode shown in the bottom panel of the figure it becomes clear why. Each of the many symmetric bands is convoluted with about nine clearly visible transitions in the torsion mode, and this gives rise to an overall very dense band for the valence state, in which little structure can be discerned, except at low energy. At the onset

FIG. 5. Absorption spectrum for the V ( 1 B 1u ) state of ethylene calculated with only Harmonic normal modes 共top panel兲 and only the torsion mode 共bottom panel兲.

of the spectrum near 55 000 cm⫺1, we can identify the transitions into the torsion mode, that are not accompanied by additional excitation of the symmetric modes. Using our potential the early visible features in the spectrum correspond to 8 and 10 quanta in the torsion mode, while the most intense feature in the torsion-only spectrum involves 14 quanta. This is indicative of the huge displacement that occurs along this coordinate, and it provides a justification for the short-time picture that underlies our vertical FC framework. Next we superimpose the two spectra of the valence and Rydberg state. The total intensity of each spectrum is determined by the transition dipole moment and no adjustment or rescaling of the intensities is made while combining the spectra. As for the Lorentzian line broadening function, we noticed that if the width is taken as 0.05 eV for the V state instead of the 0.02 eV broadening we have employed thus far, we get a somewhat better agreement with the experimental spectrum, as the valence spectrum shows less structure. Using a bigger line broadening for the valence state spectrum corresponds to a shorter lifetime of the state as compared to the Rydberg state, which is in agreement with its higher oscillator strength. In addition, the large displacement along the torsion mode in the valence state would eventually result in coupling to the other normal modes, thus leading to further congestion of the spectrum if the full potential energy sur-



FIG. 6. Absorption spectrum of ethylene in the energy range 6.2 eV 共50 000 cm⫺1兲 to 8.7 eV 共70 000 cm⫺1兲. The spectrum calculated using the vertical FC method 共bottom panel兲 is compared to the experimental spectrum 共Refs. 21 and 59兲共top panel兲.

face would be taken into account. This is mimicked to some extent by a broader line width. The total simulated spectrum along with the experimental spectrum is shown in Fig. 6. We see that the simulated spectrum has the prominent doublet structures riding on top of a broad continuum. The features near the onset of the spectrum before the first doublet as well as the long tail of the spectrum is well reproduced. There is an inconsistency of the relative intensities of the two peaks in a doublet, but apart from that, the agreement of the calculated spectrum with experiment is quite good. We note that no adjustments have been made of the vertical line positions, and this provides support for the vertical excitation energy of 7.78 eV for the valence state as calculated at the STEOM-CCSD level. The vertical excitation energy of the Rydberg state is likely somewhat higher than the 7.21 eV calculated using STEOMCCSD and the EOM-CCSD value of 7.28 eV may be more accurate as our current value for the 0-0 transition for the Rydberg state based on STEOM, is lower than experiment by about 1000 cm⫺1 共Table VII兲. IV. CONCLUDING REMARKS

A method for calculating electronic absorption spectrum of polyatomic molecules is presented which is in the spirit of the familiar FC method, but goes beyond the harmonic ap-


2135

proximation for the potential energy surfaces. We refer to this new method as the vertical Franck-Condon approach. The procedure for efficient computer implementation is discussed. The method is very efficient for calculating spectra as compared to approaches based on vibronic model Hamiltonians 共for example, Refs. 44, 53, and 54兲, which are sometimes used to calculate spectra when surfaces are highly anharmonic. Of course, unlike the vibronic model Hamiltonian approaches, the present approach is based on the BornOppenheimer approximation and does not capture nonadiabatic effects. Moreover the potentials underlying the vertical FC model are approximated to be separable and this is another limitation of the present approach that is overcome to some extent in vibronic models. The vibronic models can provide anharmonic model potentials that can be used in model type vertical FC approaches, and this would alleviate the computationally expensive ab initio evaluation of the one-dimensional potential energy surface along anharmonic normal modes. Work along these lines is in progress. The first principles simulation of the electronic absorption spectrum of ethylene in the 6.2 eV 共50 000 cm⫺1兲 to 8.7 eV 共70 000 cm⫺1兲 range using the STEOM-CCSD method in conjunction with the vertical FC approach developed in this paper is in nice agreement with experiment. Excitation to the ␲ →3s Rydberg and ␲→␲* valence states determines the spectrum in this frequency range. The Rydberg state is responsible for the four sets of prominent doublets that ride on top of the broad valence band. The vibrational progressions in the Rydberg state are found to be due to the CvC stretch and torsion modes, with the torsion imprinting the prominent doublet structure on the symmetric progressions. For the valence transition, the CvC stretch, CH2 scissors, and the torsion modes are all found to be important in generating the featureless continuum, that constitutes the valence band. The torsion mode effectively smoothens the symmetric progressions and is also responsible for the onset of the spectrum at 55 000 cm⫺1. ACKNOWLEDGMENTS

This work was supported by the NSERC Discovery Grants program. We would like to thank Professor R. J. Le Roy for his assistance in using the program LEVEL,42 which we used to test a part of our code. J. Franck, Trans. Faraday Soc. 21, 536 共1925兲. E. Condon, Phys. Rev. 28, 1182 共1926兲. 3 E. U. Condon, Phys. Rev. 32, 858 共1928兲. 4 A. M. Mebel, Y. T. Chen, and S. H. Lin, J. Chem. Phys. 105, 9007 共1996兲. 5 B. O. Roos, P. A. Malmqvist, V. Molina, L. Serrano-Andre´s, and M. Merchań, J. Chem. Phys. 116, 7526 共2002兲. 6 P. Calaminici, A. M. Ko¨ster, and D. R. Salahub, J. Chem. Phys. 118, 4913 共2003兲. 7 W. Domcke, L. S. Cederbaum, H. Ko¨ppel, and W. Von Niessen, Mol. Phys. 34, 1759 共1977兲. 8 L. S. Cederbaum and W. Domcke, Adv. Chem. Phys. 36, 205 共1977兲. 9 D. J. Tannor and E. J. Heller, J. Chem. Phys. 77, 202 共1982兲. 10 E. J. Heller, Acc. Chem. Res. 14, 368 共1981兲. 11 E. J. Heller, J. Chem. Phys. 62, 1544 共1975兲. 12 A. J. Merer and R. S. Mulliken, Chem. Rev. 共Washington, D.C.兲 69, 639 共1969兲. 13 F. Duschinsky, Acta Physicochim. URSS 7, 551 共1937兲. 14 I. Ozkan, J. Mol. Spectrosc. 139, 147 共1990兲. 1 2


2136


G. M. Sando and K. G. Spears, J. Phys. Chem. A 105, 5326 共2001兲. J. T. Hougen and J. K. G. Watson, Can. J. Phys. 43, 298 共1965兲. 17 P. G. Wilkinson and R. S. Mulliken, J. Chem. Phys. 23, 1895 共1955兲. 18 R. McDiarmid and E. Charney, J. Chem. Phys. 47, 1517 共1967兲. 19 R. J. Buenker, S. D. Peyerimhoff, and H. L. Hsu, Chem. Phys. Lett. 11, 65 共1971兲. 20 P. D. Foo and K. K. Innes, J. Chem. Phys. 60, 4582 共1974兲. 21 M. B. Robin, Higher Excited States of Polyatomic Molecules 共Academic, New York, 1975兲, Vol. 2, pp. 2–22. 22 R. J. Buenker and S. D. Peyerimhoff, Chem. Phys. 9, 75 共1975兲. 23 R. McDiarmid, J. Phys. Chem. 84, 64 共1980兲. 24 D. G. Wilden and J. Comer, J. Phys. B 13, 1009 共1980兲. 25 C. Petrongolo, R. J. Buenker, and S. D. Peyerimhoff, J. Chem. Phys. 76, 3655 共1982兲. 26 A. Gedanken, N. A. Kuebler, and M. B. Robin, J. Chem. Phys. 76, 46 共1982兲. 27 M. B. Robin, Higher Excited States of Polyatomic Molecules 共Academic, New York, 1985兲, Vol. 3, pp. 213–227. 28 B. A. Williams and T. A. Cool, J. Chem. Phys. 94, 6358 共1991兲. 29 J. Ryu and B. S. Hudson, Chem. Phys. Lett. 245, 448 共1995兲. 30 M. Ben-Nun and T. J. Martinez, J. Phys. Chem. A 103, 10517 共1999兲. 31 K. K. Baeck and T. J. Martinez, Chem. Phys. Lett. 375, 299 共2003兲. 32 A. M. Mebel, M. Hayashi, and S. H. Lin, Chem. Phys. Lett. 274, 281 共1997兲. 33 M. Ben-Nun, J. Quenneville, and T. J. Martinez, J. Phys. Chem. A 104, 5161 共2000兲. 34 M. Ben-Nun and T. J. Martinez, Adv. Chem. Phys. 121, 439 共2002兲. 35 M. Nooijen and R. J. Bartlett, J. Chem. Phys. 106, 6441 共1997兲. 36 M. Nooijen and R. J. Bartlett, J. Chem. Phys. 107, 6812 共1997兲. 37 M. Nooijen and R. J. Bartlett, J. Chem. Phys. 106, 6449 共1997兲. 38 M. Nooijen, Spectrochim. Acta 55, 539 共1999兲. 39 A. Hazra and M. Nooijen, Int. J. Quantum Chem. 95, 643 共2003兲. 40 V. I. Krylov, Approximate Calculation of Integrals 共Macmillan, New York, 1962兲, 共translated by A. H. Stroud兲.

Hazra, Chang, and Nooijen H. Wei and T. Carrington, J. Chem. Phys. 101, 1343 共1994兲. R. J. Le Roy, LEVEL 7.5: A Computer Program for Solving the Radial Schro¨dinger Equation for Bound and Quasibound Levels, University of Waterloo Chemical Physics Research Report CP-655 共2002兲. The source code and manual for this program may be obtained from the ‘‘Computer Programs’’ link at http://leroy.uwaterloo.ca 43 A. J. Sadlej, Theor. Chim. Acta. 79, 123 共1991兲. 44 M. Nooijen, Int. J. Quantum Chem. 95, 768 共2003兲. 45 B. O. Roos, K. Andersson, M. P. Fulscher, P. A. Malmqvist, L. SerranoAndre´s, K. Pierloot, and M. Merchań, Adv. Chem. Phys. 93, 219 共1996兲. 46 T. Mu¨ller, M. Dallos, and H. Lischka, J. Chem. Phys. 110, 7176 共1999兲. 47 M. Ben-Nun and T. J. Martinez, Chem. Phys. 259, 237 共2000兲. 48 E. R. Davidson and A. A. Jarzecki, Chem. Phys. Lett. 285, 155 共1998兲. 49 J. E. Del Bene, J. D. Watts, and R. J. Bartlett, J. Chem. Phys. 106, 6051 共1997兲. 50 A. D. Walsh, J. Chem. Soc. 1953, 2325. 51 K. B. Wiberg, C. M. Hadad, J. B. Foresman, and W. A. Chupka, J. Phys. Chem. 96, 10756 共1992兲. 52 W. Siebrand, F. Zerbetto, and M. Z. Zgierski, Chem. Phys. Lett. 174, 119 共1990兲. 53 H. Ko¨ppel, W. Domcke, L. S. Cederbaum, and W. von Niessen, J. Chem. Phys. 69, 4252 共1978兲. 54 H. Ko¨ppel, W. Domcke, and L. S. Cederbaum, Adv. Chem. Phys. 57, 59 共1984兲. 55 G. Herzberg, Molecular Spectra and Molecular Structure (III) 共Van Nostrand, Princeton, NJ, 1966兲. 56 R. J. Sension and B. S. Hudson, J. Chem. Phys. 90, 1377 共1989兲. 57 J. D. Watts, S. R. Gwaltney, and R. J. Bartlett, J. Chem. Phys. 105, 6979 共1996兲. 58 L. Serrano-Andre´s, M. Merchań, I. Nebot-Gil, R. Lindh, and B. O. Roos, J. Chem. Phys. 98, 3151 共1993兲. 59 J. Geiger and K. Wittmaack, Z. Naturforsch. 20A, 628 共1965兲.

15

41

16

42
