A variational approach for calculating Franck-Condon factors including mode-mode anharmonic coupling

We have implemented our new procedure for computing Franck-Condon factors utilizing vibrational configuration interaction based on a vibrational self-consistent field reference. Both Duschinsky rotations and anharmonic three-mode coupling are taken into account. Simulations of the first ionization band of ClO2 and C4H4O furan using up to quadruple excitations in treating anharmonicity are reported and analyzed. A developer version of the MIDASCPP code was employed to obtain the required anharmonic vibrational integrals and transition frequencies. © 2006 American Institute of Physics. DOI: 10.1063/1.2360944


I. INTRODUCTION
Ab initio simulations of UV and photoelectron spectra can play a key role in the assignment and interpretation of high-resolution experimental measurements.Within the Born-Oppenheimer approximation the dominant contribution to transition intensities is given by the Franck-Condon factors ͑FCFs͒. 1,2Although there are other contributions, usually quite small, given by Herzberg-Teller factors, 3 possible non-Born-Oppenheimer effects, 4 and rotational contributions we focus in this paper on the theoretical evaluation of vibronic FCFs.
In the literature there is a wide range of methods to calculate FCFs, but nearly all of them are based on the harmonic approximation for the potential energy surface ͑PES͒ of both the initial and final electronic states.][9][10][11][12][13][14] Just last year, for instance, a new efficient procedure for calculating Duschinsky rotations 15 for large molecules was presented. 16ccounting for the effect of anharmonicity in the initial and final electronic states PES dramatically increases the computational cost, and in addition, usually requires a more complex formulation.Recently, several new procedures to include anharmonicity have been formulated, [17][18][19][20] allowing the accurate reproduction of experimental spectra for small molecules.][22] We have recently formulated and implemented an alternative approach [23][24][25] including Duschinsky rotations and anharmonicity ͑including nondiagonal mode-mode coupling͒.It has been successfully applied to ClO 2 , ethylene, 24 and furan, 25 even though the latter molecule has 21 vibrational degrees of freedom.In ethylene and furan, calculations were feasible, thanks to an algorithm to truncate the vibrational basis set that reduces drastically the computational cost without significant loss of accuracy.Our method introduces mode-mode anharmonicity through vibrational perturbation theory using the harmonic oscillator zeroth-order model. 23,24his perturbation theory treatment is the least computationally demanding option, but will often diverge when the anharmonicity is large as in van der Waals species, molecular clusters, 26 molecules with multiple minima separated by low barriers, 24 and, in general, chemical systems with large amplitude motions.The well-known alternative that can be used to circumvent the divergence problem is the variational method. 26ibrational self-consistent field ͑VSCF͒ theory is the lowest level variational treatment that takes into account mode-mode anharmonicity. 27In the VSCF approach each mode is governed by a mean field potential obtained by averaging over the coupling with all other modes.This neglects the instantaneous coupling, or correlation, between modes.As in electronic structure theory there are a number of ways to include vibrational correlation.Perhaps the most straightforward, generally effective, method is vibrational configuration interaction [27][28][29][30][31] ͑VCI͒ which is the primary focus of this paper.The work presented here represents the first at-tempt to evaluate FCFs for a system containing more than four atoms with mode-mode coupling included by means of VCI.
In Sec.II a brief general review of the theory is provided.Then, Sec.III describes the more technical aspects of the calculations.In Sec.IV the results for ClO 2 and furan are analyzed in detail.Finally, our key conclusions and future plans are given in Sec.V.

A. Simultaneous equations for determination of Franck-Condon integrals
In the calculation of vibronic Franck-Condon integrals we seek to calculate the overlap between vibrational wave functions obtained for different electronic states.Accordingly we consider the calculation of the vibrational eigenstates of two electronic states, Here T ˆg, V ˆg, v g g , and E v g g are the vibrational kinetic energy operator, the vibrational potential energy function, wave function, and energy of the ground electronic state, respectively, and T ˆe, V ˆe, v e e , and E v e e are their counterparts for the electronic excited state.We consider the vibronic case ͑i.e., the effect of the rotations are neglected͒ and use the simple kinetic energy operator T =−ប 2 ͚ k 3N−6 ‫ץ‬ 2 / ‫ץ‬Q k 2 in terms of the normal coordinates Q k .The additional terms in Watson's Hamiltonian give small, though not always insignificant, 32 contributions to the energy, and their effect on the FCFs through the wave functions is anticipated to be small.
Our new procedure for calculating FCFs is based on the fact that the vibrational wave function for the initial electronic state may be expanded in terms of the complete set of vibrational eigenfunctions of the final electronic state.If the system is originally in the ground electronic state, for example, it is then straightforward to derive the following set of homogenous linear simultaneous equations 23,24 for the Franck-Condon integrals, Here ␦ e v e is the Kronecker delta.In this derivation the vibrational kinetic operators we use for the two electronic states exactly cancel one another. 23Using the Duschinsky rotation Q g = K + JQ e , relating the normal coordinates of the electronic excited state ͑Q e ͒ to those of the ground state ͑Q g ͒ it is easy to verify this cancellation.Indeed the total nuclear kinetic energy operator, expressed in terms of momenta conjugate to external Cartesian coordinates, is also independent of electronic state.For any v g , the normalized eigenvector

2
. Of course, the set of simultaneous equations to be solved is infinite.Thus, in practice, the set of all v e must be truncated to a finite dimension M ͑see below͒.
Using again the Duschinsky rotation, one can write the potential energy difference V ˆg − V ˆe solely in terms of excited electronic state normal coordinates. 23,24Then all quantities in Eq. ͑3͒ except for E v g g and the Franck-Condon integrals become properties of the excited electronic state.The quantities in Eq. ͑3͒ are easily evaluated at the harmonic oscillator level of approximation.However, when anharmonicity must be taken into account the calculations become more complex and expensive.Vibrational perturbation theory, using the normal coordinate harmonic oscillator description as zeroth-order model, allows the effect of anharmonicity to be introduced at relatively low computational cost. 23,24Unfortunately, the perturbation treatment will often diverge. 26We present here an alternative approach where the effect of anharmonicity is included by means of VCI based on VSCF functions for the individual modes.Such an approach is more elaborate but, in principle, can account for anharmonicity to a high degree of accuracy for bound vibrational states.

B. Truncation of the basis set
The number of vibrational states increases exponentially with the number of normal modes.For that reason the solution of Eq. ͑3͒ rapidly becomes unfeasible as the molecular size increases, even when only low energy vibrational states are considered.One of the most critical steps in our method is the truncation of the vibrational basis set in the face of competing requirements.On one hand, for an accurate determination of the Franck-Condon integrals all vibrational basis functions that contribute significantly to the Franck-Condon intensity need to be taken into account.On the other hand, efficient computation requires limiting the calculation to a relatively small number of vibrational excited states.Our algorithm includes only the necessary vibrational excited states.It involves an iterative buildup of the basis that increases the range of vibrational quantum numbers while, simultaneously, removing unnecessary functions. 23,24All vibrational states that, in a particular iteration, lead to FCFs smaller than a certain threshold are excluded in basis set augmentations for the following iterations.This algorithm drastically reduces the growth of the basis set and was vital, for instance, in calculating mode-mode anharmonicity contributions to FCFs for the photoelectron spectrum of furan. 25s far as we know, except for our own previous simulations, 24,25 no FCF calculations have been reported that include the effect of anharmonic mode-mode coupling for molecules containing more than four atoms.
The iterative approach described here is essentially the same as that used in our previous work. 24,25The principal difference in the present work is that another engine is now used for calculating the necessary energies and integrals for Eq.͑3͒ as will now be described.

C. Vibrational configuration interaction wave functions, energies, and determination of anharmonic
Š e e ͦV ˆg − V ˆeͦ v e e ‹ integrals Our description of anharmonicity in this paper is based on VCI using the VSCF reference state.The VSCF reference state, in turn, is taken to be a direct product of one-mode functions.Applying the variational principle to this ansatz provides the working equations for calculating the optimal one-mode functions called modals. 27The VSCF treatment accounts for anharmonicities in an average field sense.To provide a more accurate treatment we employ the VCI approach which has been documented in a number of works.Our particular version is described in Ref. 33 and implemented in MIDASCPP. 34Thus, we begin with a VSCF calculation for the vibrational ground state.Subsequently a VCI is performed where we find a number of requested roots.The definition of the VCI calculation requires a definition of the VCI space.This is done in terms of the number of modals per mode and an excitation level.Thus, for each mode we choose the energetically lowest VSCF modals up to a certain number as the basis.The number used for each mode, as dictated by the iterative procedure in Sec.II B, is the highest quantum number for that mode in the particular set of manymode basis states.The excitation level is the maximum number of modes that are allowed to be simultaneously excited out of the VSCF reference state.Thus VCI͓3͔ includes up to three mode couplings in the wave function, and can also be denoted VCISDT following a nomenclature that is standard for electronic wave functions.
Our implementation of VCI is based on a direct CI approach wherein the VCI Hamiltonian is never explicitly constructed.Instead the relevant VCI states are found using it-erative methods that require only the transformation of vectors with the VCI Hamiltonian matrix.The program takes as input a set of target states defined in terms of a set of quantum numbers.This set of quantum numbers is similar to describing states in the harmonic approximation in terms of the harmonic quantum numbers.The difference is that the quantum numbers refer to the VSCF modals and, furthermore, that the states are allowed to have other components.Thus, given an initial target vector the program iterates using Davidson's algorithm to find the true VCI eigenstate with this target vector as the largest component.The true VCI vectors incorporate anharmonicity in the approximation defined by the VCI wave function level.
Once the VCI states have been found the integrals required for Eq.͑3͒ must be calculated.In the context of a direct VCI it is relatively straightforward to replace the vibrational Hamiltonian operator with the V ˆg − V ˆe operator to facilitate calculation of the ͗ e e ͉V ˆg − V ˆe͉ e e ͘ matrix elements.
Thus, for each converged VCI state we perform one direct VCI transformation with the V ˆg − V ˆe operator, and the requested integrals are subsequently calculated as dot products of the resulting vector with the converged VCI vectors.

III. COMPUTATIONAL DETAILS
In order to test our method for calculating FCFs using VCI we computed the ClO 2 The new variational values are compared with our previously validated results obtained by means of vibrational perturbation theory. 23,25The same geometries as in the perturbation study were utilized for the neutral and cation in all calculations.Whereas the geometry of the neutral molecule was taken to be the experimental one, [35][36][37] for the cation the geometry was obtained from an iterative Franck-Condon analysis 17 ͑IFCA͒ in which the FCFs were determined by second-order vibrational perturbation theory ͑PT2͒.
For ClO 2 the harmonic and anharmonic force constants were derived from the PESs calculated by Peterson and Werner for the neutral 38 and cationic 39 ground states by means of a complete active space self-consistent-field/ multireference configuration interaction ͑CASSCF/MRCI͒ treatment with a cc-pVQZ basis set. 23For furan, ab initio quadratic, cubic, and quartic vibrational force constants, expressed in terms of normal modes, were obtained at the ͑U͒B3LYP 6-311+ + G͑d , p͒ level using the Romberg method 40 to reduce numerical error. 25All ab initio electronic structure calculations were performed with the GAUSSIAN 03 suite of programs. 41heoretical FCFs and vibrational energies were used to determine the relative peak intensities and positions in the He I PE spectrum except that the first peak was shifted to the experimental value, which occurs at the adiabatic ionization energy ͑AIE͒ ͑10.345 eV for ClO 2 as measured by Flesch et al. 42 and 8.883 eV for furan as measured by Derrick et al. 43 ͒.In order to simulate the experimental spectra we used Gaussian functions with a full width at half maximum of 30 meV for ClO 2 and 15 meV for furan.
In the FCF calculations we perform preliminary tests on the threshold ͑10 −n ͒ FCF value that determines which vibrational states are excluded ͑see Sec.II B͒ in order to find what is adequate for the desired level of accuracy.The value of n was chosen so that the height of the most intense peak obtained with the thresholds 10 −n and 10 −͑n+1͒ differs by less than 0.2%.For ClO 2 this led to n = 4 whereas for furan n =6.
The VSCF modals were obtained by expanding in terms of the set of harmonic oscillator eigenfunctions for that mode.All harmonic oscillator functions are included up to vibrational quantum number of 20.We have checked that the change in the energy of the ground and excited vibrational states needed for evaluation of the FCFs is negligible when the maximum quantum number is reduced to 15.Our vibrational configurations are constructed from these modals using a maximum excitation level from the input and the modal space determined by the algorithm described in Sec.II B. The maximum error in the vibrational energies obtained with the above thresholds was 9 ϫ 10 −4 % for ClO 2 and 10 −5 % for furan.

A. Chlorine dioxide
In Fig. 1 we present both our harmonic and anharmonic VCISDT simulations of the first band in the ClO 2 He I PE spectrum.All three modes, corresponding to calculated harmonic frequencies of 1012 cm −1 ͑ 1 ͒, 511 cm −1 ͑ 2 ͒, and 1409 cm −1 ͑ 3 ͒, were included in the vibrational basis set of each simulation.Here 1 and 2 are the totally symmetric ͑a 1 ͒ stretch and bend, respectively, while 3 is the asymmetric ͑b 2 ͒ stretch.
The most important FCFs and vibronic energy levels used to generate the VCISDT simulation are shown in Table I, whereas our assignment of peaks is given in Table II.The results indicate that the largest FCFs are associated with the two a 1 modals.The progression with the most intense peaks is ͑ 1 00͒ = ͑000͒ , ͑100͒ , ͑200͒ , ͑300͒ ,..., although the progression ͑ 1 00͒ = ͑010͒ , ͑110͒ , ͑210͒ , ͑310͒ ,... is also clearly observable.Apart from the two dominant progressions, all other transitions listed in Table II have an intensity less than 1% of the major ͑100͒ transition.States with an odd number of quanta in the b 2 mode have a zero Franck-Condon integral.Nonetheless, transitions where b 2 is involved become the dominant ones in the three highest energy peaks presented in Table II.
The data presented in Fig. 1, as well as Tables II and III, evidence the important difference between harmonic and VCISDT transition energies, relative peak intensities, and peak assignments.As expected, the harmonic vibrational transition frequencies are larger than their VCISDT counterparts, particularly at the high energy end of the spectrum where the difference can be as large as 5%.The difference between the harmonic and VCISDT transition energies follows the same pattern for both progressions described above.The largest FCF associated with each peak is in boldface.

154114-4
On the contrary, the difference in relative intensities exhibits a unique pattern for each progression.For the ͑ 1 10͒ progression the maximum difference in relative intensities is only 3.1% and it is associated with the ͑010͒ peak which, in the case of VCISDT, has a relative intensity of 12.5%.On the other hand, the relative intensity differences in the ͑ 1 00͒ progression are quite a bit larger.The maximum is 13.2%, which occurs for the ͑300͒ peak which has a relative intensity ͑versus the main peak͒ of 20.7%.For the ͑400͒ peak of this series, the difference in relative intensities ͑8.6%͒ is larger than the relative intensity ͑4.2%͒ itself.
For the main peaks of the ClO 2 band the harmonic and VCISDT simulations lead to equivalent assignments.However, there are a couple of differences to be noted.One is that the harmonic simulation includes more transitions that are weak, but still greater than the threshold of 10 −4 .For instance, the 11.061 eV VCISDT peak is generated solely by the ͑302͒ transition, whereas the cor-responding harmonic peak has contributions from the ͑600͒ and ͑520͒ transitions as well.A second difference is that, in the harmonic simulations, the ͑ 1 00͒ and ͑ 1 10͒ progressions also make the dominant contribution to the high energy peaks.
In Fig. 2 we depict the VCISDT and PT2 simulations of the first band in the ClO 2 PE spectrum while Table IV contains the transition energies, relative intensities, and assignments for the PT2 calculation.The results show very good agreement between the VCISDT and PT2 simulations.This validates the VCI implementation of our FCF computational method and also indicates that, for this chemical system, PT2 is an excellent procedure to include the effect of anharmonicity.In spite of the excellent overall agreement some interesting minor discrepancies can be observed in the assignments.The most evident is that the number of PT2 transitions that are weak, yet have intensity greater than the threshold of 10 −4 , is even larger than in the harmonic simu- lation.Maybe more important is the presence of negative FCFs in the PT2 assignments.Although the negative FCFs are used to obtain the relative intensities presented in Table IV, we did not use them in the PT2 simulation since they have no physical meaning.A detailed analysis shows that the FCFs which are negative in the PT2 simulation can change their value radically depending on the treatment of anharmonicity.For instance, the harmonic FCF value for the ͑020͒ transition ͑7.4ϫ 10 −3 ͒ is reduced to 2.3ϫ 10 −4 at the VCISDT level and becomes negative at the PT2 level ͑−2.4ϫ 10 −4 ͒.
Finally, we performed a VCISD simulation with the goal of determining the effect of the triple excitations.As can be seen in Fig. 3, the differences between the VCISD and VCISDT simulations are negligible for the ͑ 1 00͒ progression ͑Ͻ0.2% ͒, and very small for the ͑ 1 10͒ progression ͑Ͻ1.3% ͒.

B. Furan
Our previous PT2 simulations of the C 4 H 4 O + X ˜2A 2 ← C 4 H 4 O X ˜1A 1 band in the photoelectron spectrum 25 reproduced the observed spectrum quite accurately 43 and led to the conclusion that a basis set containing excitations of just three a 1 modes ͑ 4 , 6 , 8 ͒ was sufficient to accurately reproduce the experimental spectrum.Excitations of a fourth a 1 mode ͑ 5 ͒ were needed to obtain the anharmonic frequencies but not the Franck-Condon factors themselves.Taking into account these previous PT2 results, we included only the modals corresponding to 4 , 6 , and 8 in the FCF basis set for our first series of variational simulations.Although the resulting VCISDT modal space, i.e., VCISDT-3, involved excitation of only these three modals, all 21 modes were included in optimizing the VSCF state used as reference for the VCI calculation.The resulting modals are optimal for describing anharmonicity in the ground state and presumably constitute a good basis for the VCI calculation.In the VCI calculation some modes are simply frozen in the ground state modal and the anharmonic interaction with the other modes is included within this restriction.A more complete account of the vibrational correlation would include VCI for all modes.Such a procedure is possible but computationally expensive.First of all each step in optimizing the VCI wave function becomes significantly more expensive.In addition, the convergence of high-lying states necessary in our iterative procedure becomes difficult requiring many iterations and/or specialized methods.In this first work, therefore, we restrict our study to a few modes based on previous studies and later verify that adding more modes does not greatly affect the results.In  The largest FCF associated with each peak in boldface.Only the FCF larger than 10 −4 units has been used to label the peaks.
come dominant for the high energy peaks.For instance, the 4 = 1 peak of the ͑ 4 01͒ progression overlaps the 4 =0 peak of the ͑ 4 20͒ progression with the former FCF being larger.On the other hand, for the higher energy members of these overlapping progressions we find that ͑201͒ is weaker than ͑120͒, ͑301͒ is weaker than ͑220͒, etc.
Although the harmonic-3 spectrum reproduces all the main features of the VCISDT-3 spectrum, the difference in relative intensities can be as large as 15%.This occurs for the second peak of the ͑ 4 00͒ progression which has a relative intensity of 73.8% compared to the main peak.The most important cases, as seen in Fig. 4, are the second and third peaks of the ͑ 4 00͒ progression and the first and second peaks of the ͑ 4 10͒ progression.For the other 24 peaks the difference between the harmonic-3 and VCISDT-3 relative intensities is always smaller than 2.9%, although for the high frequency peaks that error can be larger than the relative intensity versus the main peak.
For C 4 H 4 O + the harmonic vibrational transition frequencies are again larger than their VCISDT counterparts.However, the difference is always smaller than 2.1%.There is also a very good agreement between the harmonic-3 and VCISDT-3 assignments.Indeed, for the first 15 peaks in Tables V and VI, which include all peaks with a relative intensity larger than 10% of the ͑000͒ transition, the assignments are identical.In contrast with the scenario presented for ClO 2 , the VCISDT simulation results in more weak transitions contributing to the high energy peaks than in the harmonic simulation.
For test purposes we also compared the literature PT2 simulation with the equivalent VCISDT simulation including excitations of only the three key a 1 modes in the PT2 and VSCF basis sets.The agreement between both simulations is virtually exact with the maximum difference in relative intensities being smaller than 0.7%.As mentioned before, excitations of the 5 a 1 mode do not affect the Franck-Condon integral results but they are needed for accurate determination of the vibrational frequencies.Next, we performed a VCISDT-4͑ 5 ͒ calculation by adding excitations of the 5 modal to the three modal FCF basis set and VCISDT-3 space.Once more, however, excitations in all 21 modes were included in the VSCF calculation.In Table VII we present transition energies, relative intensities, and assignments for the VCISDT-4͑ 5 ͒ calculation.The comparison of this simulation with the equivalent PT2 treatment 25 again shows excellent agreement.For instance, the difference between the VCISDT-4͑ 5 ͒ and PT2-4͑ 5 ͒ relative intensities is always smaller than 1.2%.These results once more validate our implementation in the VCI context and they also confirm that PT2 is an adequate procedure to evaluate the FCFs of the were considered in the reference state VSCF calculation for the VCISDT-4͑ 5 ͒ simulation, the good agreement with the PT2-4͑ 5 ͒ simulation provides further evidence for the minor role of the other 17 modes in determining the vibrational transition frequencies.To study the effect of excitations of these other modes on the FCFs they would also have to be included in the FCF and VCI spaces.
There is a very good agreement between the VCISDT-3 and VCISDT-4͑ 5 ͒ simulations.The maximum differences in the relative intensities and in the vibrational transition frequencies are 1.7% and 0.8%, respectively.Thus, the VSCF reference state calculation sufficiently accounts for coupling of 5 with 4 , 6 , and 8 as far as evaluation of vibrational transition frequencies is concerned.On the other hand, as previously observed for the PT2 calculations, 25 it is unnecessary to include excitations of 5 in the FCF basis set in order to obtain accurate Franck-Condon integrals.As confirmation of the good agreement between the VCISDT-3 and VCISDT-4͑ 5 ͒ simulations we may consider the peak assignments shown in Tables V and VII.The largest FCF associated with each peak agrees in both simulations.Although the VCISDT-4͑ 5 ͒ simulation includes a few states with excitations of 5 , these states play a minor role.
We also carried out four modal calculations obtained by adding excitations of either 1 , 2 , 3 , or 7 to VCISDT-3.For 1 , 2 , and 7 the effect is small with the differences in relative intensities being always smaller than 1.4%.How- The largest FCF associated with each peak is in boldface.Only the FCF larger than 10 −4 units has been used to label the peaks.
ever, for 3 there are differences as large as 6.0% and 4.3% in the ͑100͒ and ͑200͒ peaks, which have VCISDT-4͑ 5 ͒ relative intensities ͑versus the main peak͒ of 79.8% and 35.8%, respectively.In fact this non-negligible role of 3 excitations can also be seen in harmonic and PT2 simulations.However, the differences between the harmonic-3 ͑PT2-3͒ and harmonic-4͑ 3 ͒ ͓PT2-4͑ 3 ͔͒ results become negligible when the IFCA procedure is applied. 25In Table VIII we present transition energies, relative intensities, and assignments for the VCISDT-4͑ 3 ͒ calculation whereas Fig. 5 compares the VCISDT-3 and VCISDT-4͑ 3 ͒ simulations.
The comparison between the assignments shown in Tables V and VIII permits us to rationalize the effect of 3 .The 3 and 4 vibrations have a very similar fundamental frequency, i.e., 184 vs 176 meV.Consequently, the ͑0 4 00͒ progression overlaps the ͑ 3 000͒ progression leading to a significant increase in the intensity of the peaks assigned to this progression.Although with less magnitude, the same effect also oc-curs because of the overlapping of the ͑0 4 10͒ and ͑ 3 010͒ progressions.Finally, we performed a VCISDT-5 simulation including the excitations of 3 , 4 , 5 , 6 , and 8 .The maximum difference between the VCISDT-4͑ 3 ͒ and VCISDT-5 is 3.3% in the relative intensities and 0.7% in the vibrational transition frequencies.Contrary to the VCISDT simulation of the first band in the ClO 2 PE spectrum, triple excitations play an essential role in the first band of the furan PE spectrum ͑see Fig. 6͒.Obviously, the VCISD-4͑ 3 ͒ calculation cannot simulate the ͑0111͒ triple excitation peak which has a relative intensity of 15.8%.In addition, however, the omission of the triple excitations prevents the correct description of single and double excitation peaks.For instance, the relative intensity of ͑0001͒ in the VCISDT-4͑ 3 ͒ simulation is 27.0%, while its VCISD-4͑ 3 ͒ counterpart is only 15.3%.Considering the important differences between VCISD-4͑ 3 ͒ and VCISDT -͑ 3 ͒ one might think that still higher-order excitations are shows that both simulations agree very well.The differences in relative intensity are always smaller than 2.0%.Not by chance, the largest differences are in the ͑0 4 11͒ triple excitation progression, since these peaks overlap the ͑1 4 11͒ quadruple excitation progression.

V. CONCLUSIONS
Our procedure for calculating vibronic Franck-Condon factors has been generalized so that vibrational configuration interaction ͑VCI͒ based on a vibrational self-consistent field ͑VSCF͒ reference can be employed.The required anharmonicity integrals and anharmonic vibrational transition frequencies are computed using a developer version of   PE spectra as a test of the new method.Both anharmonic mode-mode coupling and Duschinsky rotations effects are taken into account.This procedure is the first to be employed for incorporating mode-mode anharmonicity in chemical systems with more than four atoms as far as we know.The algorithm to select the minimal sufficient basis set for a given accuracy is a vital feature of our approach.Inclusion of anharmonicity effects in the PE spectra is essential.For the two bands investigated in this paper the differences between VCI and harmonic relative intensities are larger than 12%.Nonetheless, in both cases the harmonic simulation describes correctly the most important features and assignments of the VCI spectra.
The VCISDT simulations of the ClO bands presented in this paper show excellent agreement with their literature PT2 counterparts. 23,25This accord between VCISDT and PT2 not only validates our implementation of the VCI method but also proves that, for these two bands, PT2 is an excellent way to treat the effect of anharmonicity on FCFs.Although the PT2 procedure is computationally far cheaper than VCI, it can diverge when the anharmonicity is large.The results presented here suggest that when the differences between harmonic and PT2 relative intensities are smaller than 20%, PT2 is likely to be reliable.Of course, more comparisons will be necessary to validate this result.The effect of triple and higher-order excitations in VCI depends on the band.Whereas for ClO 2 the effect of triple excitations is very small, in furan the triple excitations generate the largest FCF associated with some important peaks.Nevertheless, for furan the effect of selected quadruple excitations was found to be very small.
For furan, with 21 normal modes, the selection of a minimal, but reliable, basis set is essential to account for mode-mode anharmonicity.Using the IFCA procedure the harmonic, PT2, and VCI simulations agree that a basis set containing excitations of just three a 1 modes is sufficient to reproduce the experimental relative intensities.However, the IFCA procedure hides the non-negligible role of the excitations of 3 , which lead to differences in the relative intensity as large as 6%.The excitations of another mode, i.e., 5 , are required for accurate determination of the vibration transition frequencies, although they have no role in the evaluation of FCFs.Including these effects through the reference state VSCF calculation seems to be sufficient to take care of the deficiency in this case.
In the near future we will undertake VCI simulations of PE spectra for systems that are highly anharmonic and cannot be successfully treated by the PT2 version of our method.Furthermore, we plan to extend our code so as to calculate the vibrational contribution to two-photon absorption and Herzberg-Teller contributions to one-photon spectra.

FIG. 1 .
FIG. 1. Simulation of the first band in the ClO 2 He I PE spectrum.The dashed and solid lines represent the harmonic and the VCISDT anharmonic spectrum, respectively.

Fig. 4 FIG. 4 .
FIG. 4. Simulation of the first band in the furan PE spectrum.The dashed and solid lines represent the harmonic-3 and VCISDT-3 spectra, respectively.

FIG. 5 .
FIG. 5. Simulation of the first band in the furan PE spectrum.The solid and dashed lines represent the VCISDT-3 and VCISDT-4͑ 3 ͒ spectra, respectively.
, and Christiansen J. Chem.Phys.125, 154114 ͑2006͒ Downloaded 02 Dec 2010 to 84.88.138.106.Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions determines the entire set of FCFs given by the square of the Franck-Condon integrals, i.e., F v g e = S v g e

TABLE I .
VCISDT frequencies and FCFs for the ClO 2 + X ˜1A 1← ClO 2 X ˜2B 1 PE band.The FCFs are in units so that the FCF vector is normalized to unity ͓see the discussion following Eq.͑3͔͒.

TABLE II .
VCISDT frequencies, relative intensities, and assignments of peaks in the ClO 2 + X ˜1A 1 ← ClO 2 X ˜2B 1 PE band with a relative intensity greater than 0.1% of the ͑100͒ peak.

TABLE III .
Harmonic frequencies, relative intensities, and assignments of peaks in the ClO 2 + X ˜1A 1← ClO 2 X ˜2B 1 PE band with a relative intensity greater than 0.1% of the ͑100͒ peak.
a The largest FCF associated with each peak is in boldface.FIG. 2. Simulation of the first band in the ClO 2 He I PE spectrum.The dashed and solid lines represent the PT2 and VCISDT anharmonic spectra, respectively.

TABLE IV .
PT2 frequencies, relative intensities, and assignments of peaks in the ClO 2 b Negative FCFs ͑see text͒ are in italics.FIG. 3. Simulation of the first band in the ClO 2 He I PE spectrum.The dashed and solid lines represent the VCISD and VCISDT anharmonic spectra, respectively.

TABLE V .
VCISDT-3 frequencies, relative intensities, and assignments of peaks in the C 4 H 4 O + X ˜2A 2

TABLE VI .
Harmonic-3 frequencies, relative intensities, and assignments of peaks in the C 4 H 4 O + X ˜2A 2 ← C 4 H 4 O X ˜1A 1 PE band with a relative intensity greater than 0.5% of the ͑000͒ peak.

TABLE VII .
VCISDT-4͑ 5 ͒ frequencies, relative intensities, and assignments of peaks in the C 4 H 4 O + X ˜2A 2 ← C 4 H 4 O X ˜1A 1 PE band with a relative intensity greater than 0.5% of the ͑000͒ peak.The largest FCF associated with each peak is in boldface.Only the FCF larger than 10 −4 units has been used to label the peaks.
a a The largest FCF associated with each peak is in boldface.Only the FCF larger than 10 −4 units has been used to label the peaks.Downloaded 02 Dec 2010 to 84.88.138.106.Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissionsrequired.However, a comparison of VCISDTQ-4͑ 3 ͒ transition energies, relative peak intensities, and peak assignments in Table IX with the VCISDT-4͑ 3 ͒ data in Table VIII

TABLE IX .
VCISDT-4͑ 3 ͒ frequencies, relative intensities, and assignments of peaks in the C 4 H 4 O + X ˜2A 2 ← C 4 H 4 O X ˜1A 1 PE band with a relative intensity greater than 0.5% of the ͑000͒ peak.The largest FCF associated with each peak is in boldface.Only the FCF larger than 10 −4 units has been used to label the peaks.