- Email: [email protected]
A time-dependent density functional study of vibrationally resolved excitation, emission, and ionization spectra of the S1 state of phenol
Chemical Physics Letters 420 (2006) 459–464 www.elsevier.com/locate/cplett
A time-dependent density functional study of vibrationally resolved excitation, emission, and ionization spectra of the S1 state of phenol Mattijs de Groot, Wybren Jan Buma
*
Van’t Hoff Institute for Molecular Sciences, Faculty of Science, University of Amsterdam, Nieuwe Achtergracht 166, 127-129, 1018 WV Amsterdam, The Netherlands Received 2 November 2005; in final form 15 December 2005 Available online 2 February 2006
Abstract Franck–Condon simulations of excitation, dispersed emission and ionization spectra, using geometries and force fields obtained with Time Dependent Density Functional Theory, are reported for the lowest excited singlet state of phenol. The quality of the simulations is much higher than that previously obtained with ab initio complete active space self-consistent field (CAS-SCF) calculations, demonstrating the large predictive power of these calculations and their usefulness for the interpretation of experimental, vibrationally resolved spectra. Specifically, the shortening of the C–O bond length in the excited state and the activity of the 6a mode, both missing in the CAS-SCF results, are modeled accurately. Ó 2006 Elsevier B.V. All rights reserved.
1. Introduction One of the great challenges in the current design of new photoresponsive materials is to understand and thereby predict the changes in the electronic properties after electronic excitation that are ultimately responsible for their photophysical and photochemical behavior at a molecular level. These changes can be probed effectively by vibrationally resolved electronic spectroscopy. The changes in the vibrational frequencies are a direct probe for the force constants and thus the bond strengths, while transition intensities are a direct measure of the geometry changes that occur upon excitation and thus contain information about the geometry of the probed excited state. Franck–Condon (FC) simulations can be used to extract this information quantitatively from vibrationally resolved spectra. Ideally, a high-quality FC simulation makes the assignment of the peaks in a spectrum a trivial task and – in as far that minor differences still are present – allows *
Corresponding author. Fax: +31 20 525 6456/6422. E-mail address: [email protected] (W.J. Buma).
0009-2614/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.12.095
one to come even closer to the true geometry changes by using experimental peak intensities [1]. The necessary prerequisites for such simulations are reliable geometries and force fields for the ground and excited states. It has become clear in recent years that for the majority of ground state problems, calculations based on density functional theory are the method of choice. Calculations of the same quality for excited states have proven much more difficult. Up to recently, the most widely used methods for the calculation of optimized excited state geometries and force fields of medium to large size molecules were based on single excitation configuration interaction (CIS) and complete active space self-consistent field (CAS-SCF) methods. On the whole, it is found that the applicability of CIS is limited. CAS-SCF calculations, on the other hand, have been successfully used in various studies, but cannot be employed for larger molecular systems. Moreover, in some cases they fail to give even a good qualitative description of the excited state geometry. One of the notorious cases in this respect is phenol, where FC simulations based on the bare results CASSCF calculations hardly bear any resemblance to the
460
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
experimental spectra [2]. By clever use of ‘chemical intuition’ and additional experimental data available from rotationally resolved excitation spectra, the authors of this study nevertheless succeed to fit the data to an impressive degree, but the fact that the excited state force field is apparently described more accurately than the geometry is surprising and the success of this approach might be called somewhat fortuitous. More distressing, however, is the conclusion that in such cases the FC calculation loses its predictive power and that an easy assignment of peaks becomes impossible. This particularly holds for excitation spectra because here one needs excited state frequencies, which are always calculated less accurately than the ground state frequencies needed for emission spectra. The methods based on Coupled-Cluster (CC) theory [3], such as CC2 and CCSD for the ground state and their Equation-Of-Motion (EOM) analogs for excited states [4], are in principle capable of producing very accurate results. With the development of an analytical gradient technique, EOM–CCSD presently allows for geometry optimization of excited states [5], but because an analytical second derivative technique is not available yet, calculations of excitedstate force fields still represent a challenging computational task. One consequently has to conclude that these methods cannot yet be applied routinely to larger molecules. We have recently started an extensive program in which we study experimentally the photophysics and photochemistry of several chromophores of the Photoactive Yellow Protein (PYP), amongst else by performing high-resolution spectroscopy on chromophores seeded in molecular beams. During these studies, it rapidly became clear that reliable simulations of the vibronic structure are essential to extract the maximum amount of information from the measured spectra. Dierksen and Grimme [6] have shown that with the advent of analytical excited state gradients for TDDFT developed by Furche and Ahlrichs [7] the accurate calculation of geometries and force fields with TD-DFT has become possible and can in fact yield very high quality simulations of vibrational structure in spectra. Since a large amount of experimental and theoretical data is available for phenol, we have employed this molecule as a benchmark system to assess the performance of TD-DFT for interpreting the vibrationally resolved spectra of the molecular systems of interest to us. To this purpose, we have calculated equilibrium geometries and harmonic force fields of the molecule in the S0, S1, and D0 states and present here simulations of the vibrationally resolved S1 S0 fluorescence excitation, S1 ! S0 dispersed emission, and D0 S1 photoelectron spectra. The results demonstrate that an impressive agreement between theory and experiment can be obtained on the basis of the bare theoretical results, and validate the predictive capabilities of the method. 2. Computational details Density functional calculations were performed using the TURBOMOLE 5.7 suite of programs [8–11]. The calcula-
tions used a Gaussian AO basis set of triple-f quality augmented with a double set of polarization functions (defTZVPP [12]), the hybrid B3-LYP functional [13,14], a large grid for numerical quadrature (‘grid m3’ option), and convergence of the ground state energy and density matrix to at least 108 a.u.. The geometry of the molecule in ground and excited states was optimized using analytical DFT and TDDFT [7] gradients, respectively. At the equilibrium geometries, harmonic force fields were calculated using numerical (two-side) differentiation of the gradients for both ground and excited states. Vibronic transition intensities were calculated on the basis of the Condon approximation in which the electronic transition dipole moment is assumed to be independent of the nuclear geometry. Within this approximation, the calculation of the vibronic spectrum of the transition to a single electronic state is reduced to the calculation of the appropriate Franck–Condon factors. These factors were determined employing an in-house developed program that is based on the recursion relations derived by Doktorov et al. [1], which have been applied by numerous other groups as well [2,15,16]. Even though it is found that calculations based on geometries and force fields of good quality can reproduce experimental spectra quite reasonably, there generally remain still differences between the two. If we assume that the primary cause for these differences is that the calculation of the equilibrium geometry of the excited state is not accurate enough, i.e., we assume that the S0 equilibrium geometry and the S0 and S1 force fields are ‘correct’, we can reconstruct the excited state equilibrium geometry from the experimental intensities of the mi ða0 Þ01 transitions using the procedures described by Doktorov et al. [1]. It is obvious that this assumption is not completely valid, but it is justified because the ground state calculations are usually the most accurate ones and the influence of changes in the force fields on the spectrum is small compared to the influence of changes in the equilibrium geometries. In the following, we will employ the labeling of the normal modes as used by Kleinermanns and co-workers [2], which is based on the labeling of Bist et al. [17]. This labeling scheme is, however, somewhat misleading because it suggests that the modes of phenol can easily and uniquely be described in the basis of the benzene modes. Grafton and Wheeler [18] have shown that for most modes an expansion into benzene modes is indeed possible, but because the majority of the modes needs to be described by combinations of two or more benzene modes, the labeling becomes somewhat arbitrary, in particular when different phenol modes have the same benzene mode as the largest contributing mode. To obtain a better guide to the eye in comparing simulated and experimental spectra, we have not employed the frequency scaling factor that is usually employed for B3-LYP calculations (0.9614 [19]), but the average of the scaling factors derived for each of the totally symmetric modes present in the pertaining spectrum.
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
461
3. Results and discussion 3.1. S1 00 ! S0 emission spectrum Dispersed emission and excitation spectra both contain information on the geometry differences between ground and excited state. For several reasons dispersed emission spectra, and in particular that obtained after excitation of the electronic origin can, in general, be simulated with larger accuracy than the excitation spectrum. Firstly, peak intensities in the excitation spectrum are not only determined by transition probabilities between initial and final states, but also by the efficiency of (radiationless) decay processes in the excited state. Because this efficiency can vary between vibrational levels, this can influence the peak intensities in an excitation spectrum, but not the relative peak intensities in a dispersed emission spectrum. Secondly, in recording excitation spectra with acceptable signal-to-noise ratios for the weaker transitions, one is easily lured into increasing the intensity of the excitation light. This can cause saturation of the stronger transitions in the excitation spectrum, and thus an underestimation of their transition probability. Saturation of a particular S0 ! S1 transition, on the other hand, does not affect the relative intensities in the emission spectrum. Thirdly, the errors in the calculated ground state frequencies are usually smaller than the errors in the calculated excited state frequencies. Fig. 1 displays the experimentally obtained S1 00 ! S1 emission spectrum [20] (Fig. 1a) and the simulation based on the present calculations (Fig. 1b). While simulations based on ‘raw’ CAS-SCF calculations [2] do not even reproduce qualitatively the experimental spectrum, it is clear that the present calculations lead to an excellent agreement between experiment and theory, thereby giving evidence for the high quality of the equilibrium geometries and force fields calculated with (TD-)DFT. In particular, we notice that aspects where the CAS-SCF calculations failed – the shortening of the C–O bond length in the excited state and the activity of the 6a mode – are now predicted correctly. Also all other major features in the spectrum are simulated with great accuracy. The only exception is the transition to the (16b) [2] level, which is absent in the simulation. Roth et al. explain the intensity in this transition with the assumption of a Fermi resonance with the nearby fundamental of mode 1. The harmonic model used for the present simulation does not take this anharmonic effect into account, and is therefore not capable of reproducing it. The intensity of the transition to mode 12 is somewhat underestimated, but the correlation between the calculated and measured peaks (and therefore the assignment of the experimental peaks) is evident. An ‘experimental’ excited-state geometry can be obtained from the emission spectrum when the experimentally observed intensities of the transitions to the fundamentals of modes 6a, 12, 1, 18a and 7a are used as input for a Franck–Condon calculation and a back-transformation is
Fig. 1. (a) Experimental S1 00 ! S0 emission spectrum from [21] (figure taken from [2]). (b) Franck–Condon simulation of the spectrum (calculated frequencies have been scaled with a factor of 0.985). (c). Simulated spectrum for the S1 geometry derived from the experimental intensities in the emission spectrum. Simulated spectra have been constructed by convolution of stick spectra with a Gaussian line profile with a width of 6 cm1.
done [1]. Fig. 1c shows the S1 00 ! S0 emission spectrum that is obtained for the fitted geometry; the geometrical parameters for the calculated S0 and S1 states as well as the fitted S1 state are reported in Table 1. From this Table, it can be concluded that the adjustments in the S1 geometry for the fit are small compared with the differences between the calculated S0 and S1 geometrical parameters. The same Table reports a comparison between calculated and experimentally determined rotational constants. The differences between the two are similar for ground and excited states, the largest one being only 1.1%. These results thus show a quasi quantitative agreement between the calculated geometry changes upon excitation and those determined experimentally in rotationally resolved fluorescence excitation spectroscopy studies [21,22]. The ring expands and is elongated along the b-axis, while the C–O bond length is shortened as shown in Fig. 2.
462
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
Table 1 Geometry of phenol in S0, S1 and D0 state S0 a ˚ Bond lengths (A) O–H 0.962 O–C1 1.367 1.393 C1–C2 C1–C6 1.393 C2–C3 1.391 1.388 C6–C5 C3–C4 1.390 C5–C4 1.393 C2–H2 1.085 1.082 C6–H6 C3–H3 1.083 C5–H5 1.083 C4–H4 1.081
S1
S1 fit
D0
0.966 1.348 1.421 1.414 1.421 1.422 1.415 1.413 1.082 1.079 1.079 1.079 1.082
0.966 1.349 1.428 1.419 1.424 1.426 1.420 1.418 1.082 1.080 1.080 1.080 1.082
0.971 1.308 1.428 1.432 1.365 1.366 1.419 1.413 1.083 1.081 1.081 1.081 1.082
Angles (°) C1–O–H C2–C1–O C6–C1–O C1–C2–C3 C1–C6–C5 C2–C3–C4 C6–C5–C4 C2–C1–C6 C3–C4–C5 C1–C2–H2 C1–C6–H6 C2–C3–H3 C6–C5–H5 C3–C4–H4
109.6 120.0 116.3 117.6 117.8 119.1 119.0 123.8 122.7 120.2 119.3 120.7 120.7 118.5
109.5 119.8 116.2 117.7 117.8 118.5 118.4 124.0 123.7 120.2 119.0 121.4 121.3 118.0
113.8 122.6 116.0 118.6 118.9 120.1 119.8 121.4 121.2 119.7 118.3 120.4 120.4 119.2
5363 (5314) 2649 (2621) 1773 (1756)
5303 2648 1766
109.7 122.5 117.4 119.8 119.6 120.5 120.8 120.1 119.3 120.0 119.0 119.4 119.3 120.5
Rotational constants (MHz)b A 5687 (5651) B 2627 (2619) C 1797 (1790) a b
excitation spectra is also worthwhile because studying the differences between the simulated and measured excitation spectrum can provide a more detailed description of the excited state. For example, the frequencies that show up in the excitation spectrum can be used to correct the calculated force field of the excited state, while the differences between observed and predicted intensities can elucidate the importance of decay processes. In Fig. 3, comparison is made between the experimental S0 Resonance Enhanced MultiPhoton Ionization S1 (REMPI) (1 + 1) excitation spectrum recorded by Dopfer [23] (Fig. 3a) and the spectrum predicted by the present calculation (Fig. 3b). As expected, the agreement between the two spectra is quite good for the major in-plane modes (6a, 12, 1 and 18a). Similar to what was observed for the S1
Labeling according to Fig. 2. Experimental values in brackets [21,22].
3.2. S1
S0 REMPI excitation spectrum
As described before, the simulation of excitation spectra is intrinsically more difficult than the simulation of emission spectra. Apart from the fact that it is not always possible to obtain well-resolved emission spectra, simulation of
Fig. 2. Schematic representation of the change in geometry of phenol upon excitation to the S1 state.
Fig. 3. (a) Experimental S1 00 S0 excitation spectrum from [20]. (b) Franck–Condon simulation of the spectrum (calculated frequencies have been scaled with a factor of 0.956). (c). Simulated spectrum for the S1 geometry derived from the experimental intensities in the emission spectrum. Simulated spectra have been constructed by convolution of stick spectra with a Gaussian line profile with a width of 6 cm1.
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
00 ! S0 emission spectrum, the intensity of the transition to the fundamental of mode 12 is underestimated in the simulation. This is corrected if the simulation is performed with the S1 equilibrium geometry obtained from the fit to the emission spectrum (Fig. 3c). Despite the overall good agreement, there remain a number of eye-catching differences between the simulated and the experimental spectra. Firstly, the frequency of mode 16a is underestimated by 44 cm1. Because the first overtone of this mode is visible in the spectrum, this leads to a rather large difference in peak position between simulated and measured spectra (88 cm1). Secondly, a number of small peaks (all with intensities below 10% of the intensity of the 0–0 transition) are predicted for overtones and combinations of out-of-plane modes but cannot be found in the experimental spectrum. Both differences can be explained by an inaccurate description of the force field of the excited state for the out-of-plane modes. Similar to what was observed for the CAS-SCF calculations [24], we see in the present study that the frequency of almost all out-of-plane modes is underestimated by the TD-DFT calculation and that the differences with the experimental frequencies are much larger than for the in-plane modes. Additionally, inspection of the Duschinsky matrix shows that most out-of-plane modes are subject to severe mixing upon excitation. Predicted intensities for combination modes could differ from measured intensities if calculated and true mixing parameters differ slightly. Finally, the intensity of the transition to the mode that is assigned to 9a is severely underestimated (0.9%) in our simulation. Because the calculation does not reproduce the rather large decrease in frequency found experimentally for this mode [17] (1168 ! 975 cm1), one has to conclude that this mode is not described well by the calculation. 3.3. D0
S1 00
S0 REMPI-PES
Additional proof for the high quality of the S1 calculation can be found if we can project the excited state wavefunction onto another vibronic manifold than that of the ground state of the neutral molecule. This is done in a REMPI-PhotoElectron Spectroscopy (REMPI-PES) experiment. In such an experiment, the state under study – in this case the vibrationless level of S1 – is excited and subsequently ionized by one or more photons. The ionization process is assumed to occur vertically, leading to an ionic vibrational state distribution that is governed by the Franck–Condon factors between the vibrationless level of S1 and the various vibrational levels in D0. The quality of the calculated geometry and force field for D0 is usually of the same quality as the geometry and force field of S0. This makes the simulation of a REMPI-PES spectrum a sensitive tool for the quality of the excited state calculation. The simulation of the D0 S1 0 0 S0 REMPI-PES spectrum of Anderson et al. [25] is shown in Fig. 4. The equilibrium geometry of D0 that was employed for this
463
Fig. 4. Franck–Condon simulation of the D0 S1 00 S0 REMPI-PES. The simulated spectrum has been constructed by convoluting the stick spectrum with a Gaussian line profile with a width of 16 meV. The background in the experimental spectrum is modeled by a linear intensity decreasing toward higher photoelectron energies.
Table 2 Assignment of the dominant peaks in the photoelectron spectrum of the S1 00 level of phenol (energies in meV) Experiment b
52 64 101 121 150 207
Calculation
Modea
52 65 102 123 143, 148, 150 205
18b(5) 6a 12(1) 18a(1) b(9b), 9b(15), 9a 8a
a Modes are labeled according to Bist et al. In brackets the major contributing benzene mode as calculated by Grafton and Wheeler [18]. b Not explicitly labeled by Anderson et al.
simulation is given in Table 1. Although the resolution of the experimental spectrum is not as high as in the emission-detected spectra, the impressive agreement between the experimental and simulated spectra provides further confirmation of the accuracy of the S1 calculation.1 The interpretation of the spectrum by Anderson et al. is by far and large supported by our calculations. Table 2 lists the most important modes that contribute to the photoelectron spectrum. The strongest band in the spectrum is indeed associated with mode 8a (the quinoidal mode). On the other hand, although we agree with the conclusion that the OH bending mode is important – when the phenol becomes more quinone-like, the O atom becomes more
1
The origin of the background in the experimental spectrum is unclear, but is well modeled by a linear intensity decreasing toward high photoelectron energies.
464
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
sp2 hybridized and the ring-OH angle approaches 120° – we find that this mode is mainly found in modes b and 9b instead of 8a as previously assumed. 4. Conclusions We have shown that Franck–Condon calculations using the geometries and force fields obtained with TD-DFT yield simulations of the excitation, dispersed emission spectra, and ionization spectra of the first excited state of phenol of much higher quality than have been obtained with high-quality ab initio calculations such as CAS-SCF calculations. In most aspects, the calculated and measured spectra exhibit almost quantitative agreement, demonstrating that the predictive power of the calculations is large and therefore able to provide a reliable assignment of the various transitions. At a somewhat more detailed level, we have found that fitting the calculated spectra to the experimental intensities leads to only minor adjustments in the excited state geometry. We are therefore presently applying the same methods to interpret the results of experimental studies on the excited states of chromophores of PYP [26], for which such extensive knowledge as available for phenol is not present. Acknowledgments We gratefully acknowledge Dr. Markus Gerhards (Heinrich-Heine Universita¨t Du¨sseldorf) for providing us with the S1 S0 REMPI excitation spectrum from the Ph.D. Thesis of Dr. Otto Dopfer. This work has been supported by Netherlands Organization for Scientific Research (NWO).
References [1] E.V. Doktorov, I.A. Malkin, V.I. Man’ko, J. Mol. Spectrosc. 64 (1977) 302. [2] S. Schumm, M. Gerhards, K. Kleinermanns, J. Phys. Chem. A 104 (2000) 10648. [3] T.D. Crawford, H.F. Schaefer III, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, VCH Publishers, New York, 2000, p. 33. [4] J.F. Stanton, R.J. Bartlett, J. Chem. Phys. 98 (1993) 7029. [5] E.V. Gromov, I. Burghardt, H. Ko¨ppel, L.S. Cederbaum, J. Phys. Chem. A 109 (2005) 4623. [6] M. Dierksen, S. Grimme, J. Chem. Phys. 120 (2004) 3544. [7] F. Furche, R. Ahlrichs, J. Chem. Phys. 117 (2002) 7433. [8] O. Treutler, R. Ahlrichs, J. Chem. Phys. 102 (1995) 346. [9] M.V. Arnim, R. Ahlrichs, J. Comp. Chem. 19 (1998) 1746. [10] P. Deglmann, F. Furche, R. Ahlrichs, Chem. Phys. Lett. 362 (2002) 511. [11] P. Deglmann, F. Furche, J. Chem. Phys. 117 (2002) 9535. [12] A. Scha¨fer, C. Huber, R. Ahlrichs, J. Chem. Phys. 100 (1994) 5829. [13] D. Becke, J. Chem. Phys. 98 (1993) 5648. [14] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [15] D. Gruner, P. Brumer, Chem. Phys. Lett. 138 (1987) 310. [16] P.R. Callis, J.T. Vivian, L.S. Slater, Chem. Phys. Lett. 244 (1995) 53. [17] H.D. Bist, J.C.D. Brand, D.R. Williams, J. Mol. Spectrosc. 24 (1967) 413. [18] K. Grafton, R.A. Wheeler, J. Comp. Chem. 19 (1998) 1663. [19] A.P. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502. [20] W. Roth, P. Imhof, M. Gerhards, S. Schumm, K. Kleinermanns, Chem. Phys. 252 (2000) 247. [21] N.W. Larsen, J. Mol. Struct. 51 (1979) 175. [22] G. Berden, W.L. Meerts, M. Schmitt, K. Kleinermanns, J. Chem. Phys. 104 (1996) 972. [23] O. Dopfer, Dissertation, Technische Universita¨t Mu¨nchen, Mu¨nchen, 1994. [24] S. Schumm, M. Gerhards, W. Roth, H. Gier, K. Kleinermanns, Chem. Phys. Lett. 263 (1996) 126. [25] S.L. Anderson, L. Goodman, K. Krogh-Jespersen, Ali G. Ozkabak, R.N. Zare, C. Zheng, J. Chem. Phys. 82 (1985) 5329. [26] M. De Groot, W.J. Buma, in preparation.
A time-dependent density functional study of vibrationally resolved excitation, emission, and ionization spectra of the S1 state of phenol Mattijs de Groot, Wybren Jan Buma
*
Van’t Hoff Institute for Molecular Sciences, Faculty of Science, University of Amsterdam, Nieuwe Achtergracht 166, 127-129, 1018 WV Amsterdam, The Netherlands Received 2 November 2005; in final form 15 December 2005 Available online 2 February 2006
Abstract Franck–Condon simulations of excitation, dispersed emission and ionization spectra, using geometries and force fields obtained with Time Dependent Density Functional Theory, are reported for the lowest excited singlet state of phenol. The quality of the simulations is much higher than that previously obtained with ab initio complete active space self-consistent field (CAS-SCF) calculations, demonstrating the large predictive power of these calculations and their usefulness for the interpretation of experimental, vibrationally resolved spectra. Specifically, the shortening of the C–O bond length in the excited state and the activity of the 6a mode, both missing in the CAS-SCF results, are modeled accurately. Ó 2006 Elsevier B.V. All rights reserved.
1. Introduction One of the great challenges in the current design of new photoresponsive materials is to understand and thereby predict the changes in the electronic properties after electronic excitation that are ultimately responsible for their photophysical and photochemical behavior at a molecular level. These changes can be probed effectively by vibrationally resolved electronic spectroscopy. The changes in the vibrational frequencies are a direct probe for the force constants and thus the bond strengths, while transition intensities are a direct measure of the geometry changes that occur upon excitation and thus contain information about the geometry of the probed excited state. Franck–Condon (FC) simulations can be used to extract this information quantitatively from vibrationally resolved spectra. Ideally, a high-quality FC simulation makes the assignment of the peaks in a spectrum a trivial task and – in as far that minor differences still are present – allows *
Corresponding author. Fax: +31 20 525 6456/6422. E-mail address: [email protected] (W.J. Buma).
0009-2614/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.12.095
one to come even closer to the true geometry changes by using experimental peak intensities [1]. The necessary prerequisites for such simulations are reliable geometries and force fields for the ground and excited states. It has become clear in recent years that for the majority of ground state problems, calculations based on density functional theory are the method of choice. Calculations of the same quality for excited states have proven much more difficult. Up to recently, the most widely used methods for the calculation of optimized excited state geometries and force fields of medium to large size molecules were based on single excitation configuration interaction (CIS) and complete active space self-consistent field (CAS-SCF) methods. On the whole, it is found that the applicability of CIS is limited. CAS-SCF calculations, on the other hand, have been successfully used in various studies, but cannot be employed for larger molecular systems. Moreover, in some cases they fail to give even a good qualitative description of the excited state geometry. One of the notorious cases in this respect is phenol, where FC simulations based on the bare results CASSCF calculations hardly bear any resemblance to the
460
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
experimental spectra [2]. By clever use of ‘chemical intuition’ and additional experimental data available from rotationally resolved excitation spectra, the authors of this study nevertheless succeed to fit the data to an impressive degree, but the fact that the excited state force field is apparently described more accurately than the geometry is surprising and the success of this approach might be called somewhat fortuitous. More distressing, however, is the conclusion that in such cases the FC calculation loses its predictive power and that an easy assignment of peaks becomes impossible. This particularly holds for excitation spectra because here one needs excited state frequencies, which are always calculated less accurately than the ground state frequencies needed for emission spectra. The methods based on Coupled-Cluster (CC) theory [3], such as CC2 and CCSD for the ground state and their Equation-Of-Motion (EOM) analogs for excited states [4], are in principle capable of producing very accurate results. With the development of an analytical gradient technique, EOM–CCSD presently allows for geometry optimization of excited states [5], but because an analytical second derivative technique is not available yet, calculations of excitedstate force fields still represent a challenging computational task. One consequently has to conclude that these methods cannot yet be applied routinely to larger molecules. We have recently started an extensive program in which we study experimentally the photophysics and photochemistry of several chromophores of the Photoactive Yellow Protein (PYP), amongst else by performing high-resolution spectroscopy on chromophores seeded in molecular beams. During these studies, it rapidly became clear that reliable simulations of the vibronic structure are essential to extract the maximum amount of information from the measured spectra. Dierksen and Grimme [6] have shown that with the advent of analytical excited state gradients for TDDFT developed by Furche and Ahlrichs [7] the accurate calculation of geometries and force fields with TD-DFT has become possible and can in fact yield very high quality simulations of vibrational structure in spectra. Since a large amount of experimental and theoretical data is available for phenol, we have employed this molecule as a benchmark system to assess the performance of TD-DFT for interpreting the vibrationally resolved spectra of the molecular systems of interest to us. To this purpose, we have calculated equilibrium geometries and harmonic force fields of the molecule in the S0, S1, and D0 states and present here simulations of the vibrationally resolved S1 S0 fluorescence excitation, S1 ! S0 dispersed emission, and D0 S1 photoelectron spectra. The results demonstrate that an impressive agreement between theory and experiment can be obtained on the basis of the bare theoretical results, and validate the predictive capabilities of the method. 2. Computational details Density functional calculations were performed using the TURBOMOLE 5.7 suite of programs [8–11]. The calcula-
tions used a Gaussian AO basis set of triple-f quality augmented with a double set of polarization functions (defTZVPP [12]), the hybrid B3-LYP functional [13,14], a large grid for numerical quadrature (‘grid m3’ option), and convergence of the ground state energy and density matrix to at least 108 a.u.. The geometry of the molecule in ground and excited states was optimized using analytical DFT and TDDFT [7] gradients, respectively. At the equilibrium geometries, harmonic force fields were calculated using numerical (two-side) differentiation of the gradients for both ground and excited states. Vibronic transition intensities were calculated on the basis of the Condon approximation in which the electronic transition dipole moment is assumed to be independent of the nuclear geometry. Within this approximation, the calculation of the vibronic spectrum of the transition to a single electronic state is reduced to the calculation of the appropriate Franck–Condon factors. These factors were determined employing an in-house developed program that is based on the recursion relations derived by Doktorov et al. [1], which have been applied by numerous other groups as well [2,15,16]. Even though it is found that calculations based on geometries and force fields of good quality can reproduce experimental spectra quite reasonably, there generally remain still differences between the two. If we assume that the primary cause for these differences is that the calculation of the equilibrium geometry of the excited state is not accurate enough, i.e., we assume that the S0 equilibrium geometry and the S0 and S1 force fields are ‘correct’, we can reconstruct the excited state equilibrium geometry from the experimental intensities of the mi ða0 Þ01 transitions using the procedures described by Doktorov et al. [1]. It is obvious that this assumption is not completely valid, but it is justified because the ground state calculations are usually the most accurate ones and the influence of changes in the force fields on the spectrum is small compared to the influence of changes in the equilibrium geometries. In the following, we will employ the labeling of the normal modes as used by Kleinermanns and co-workers [2], which is based on the labeling of Bist et al. [17]. This labeling scheme is, however, somewhat misleading because it suggests that the modes of phenol can easily and uniquely be described in the basis of the benzene modes. Grafton and Wheeler [18] have shown that for most modes an expansion into benzene modes is indeed possible, but because the majority of the modes needs to be described by combinations of two or more benzene modes, the labeling becomes somewhat arbitrary, in particular when different phenol modes have the same benzene mode as the largest contributing mode. To obtain a better guide to the eye in comparing simulated and experimental spectra, we have not employed the frequency scaling factor that is usually employed for B3-LYP calculations (0.9614 [19]), but the average of the scaling factors derived for each of the totally symmetric modes present in the pertaining spectrum.
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
461
3. Results and discussion 3.1. S1 00 ! S0 emission spectrum Dispersed emission and excitation spectra both contain information on the geometry differences between ground and excited state. For several reasons dispersed emission spectra, and in particular that obtained after excitation of the electronic origin can, in general, be simulated with larger accuracy than the excitation spectrum. Firstly, peak intensities in the excitation spectrum are not only determined by transition probabilities between initial and final states, but also by the efficiency of (radiationless) decay processes in the excited state. Because this efficiency can vary between vibrational levels, this can influence the peak intensities in an excitation spectrum, but not the relative peak intensities in a dispersed emission spectrum. Secondly, in recording excitation spectra with acceptable signal-to-noise ratios for the weaker transitions, one is easily lured into increasing the intensity of the excitation light. This can cause saturation of the stronger transitions in the excitation spectrum, and thus an underestimation of their transition probability. Saturation of a particular S0 ! S1 transition, on the other hand, does not affect the relative intensities in the emission spectrum. Thirdly, the errors in the calculated ground state frequencies are usually smaller than the errors in the calculated excited state frequencies. Fig. 1 displays the experimentally obtained S1 00 ! S1 emission spectrum [20] (Fig. 1a) and the simulation based on the present calculations (Fig. 1b). While simulations based on ‘raw’ CAS-SCF calculations [2] do not even reproduce qualitatively the experimental spectrum, it is clear that the present calculations lead to an excellent agreement between experiment and theory, thereby giving evidence for the high quality of the equilibrium geometries and force fields calculated with (TD-)DFT. In particular, we notice that aspects where the CAS-SCF calculations failed – the shortening of the C–O bond length in the excited state and the activity of the 6a mode – are now predicted correctly. Also all other major features in the spectrum are simulated with great accuracy. The only exception is the transition to the (16b) [2] level, which is absent in the simulation. Roth et al. explain the intensity in this transition with the assumption of a Fermi resonance with the nearby fundamental of mode 1. The harmonic model used for the present simulation does not take this anharmonic effect into account, and is therefore not capable of reproducing it. The intensity of the transition to mode 12 is somewhat underestimated, but the correlation between the calculated and measured peaks (and therefore the assignment of the experimental peaks) is evident. An ‘experimental’ excited-state geometry can be obtained from the emission spectrum when the experimentally observed intensities of the transitions to the fundamentals of modes 6a, 12, 1, 18a and 7a are used as input for a Franck–Condon calculation and a back-transformation is
Fig. 1. (a) Experimental S1 00 ! S0 emission spectrum from [21] (figure taken from [2]). (b) Franck–Condon simulation of the spectrum (calculated frequencies have been scaled with a factor of 0.985). (c). Simulated spectrum for the S1 geometry derived from the experimental intensities in the emission spectrum. Simulated spectra have been constructed by convolution of stick spectra with a Gaussian line profile with a width of 6 cm1.
done [1]. Fig. 1c shows the S1 00 ! S0 emission spectrum that is obtained for the fitted geometry; the geometrical parameters for the calculated S0 and S1 states as well as the fitted S1 state are reported in Table 1. From this Table, it can be concluded that the adjustments in the S1 geometry for the fit are small compared with the differences between the calculated S0 and S1 geometrical parameters. The same Table reports a comparison between calculated and experimentally determined rotational constants. The differences between the two are similar for ground and excited states, the largest one being only 1.1%. These results thus show a quasi quantitative agreement between the calculated geometry changes upon excitation and those determined experimentally in rotationally resolved fluorescence excitation spectroscopy studies [21,22]. The ring expands and is elongated along the b-axis, while the C–O bond length is shortened as shown in Fig. 2.
462
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
Table 1 Geometry of phenol in S0, S1 and D0 state S0 a ˚ Bond lengths (A) O–H 0.962 O–C1 1.367 1.393 C1–C2 C1–C6 1.393 C2–C3 1.391 1.388 C6–C5 C3–C4 1.390 C5–C4 1.393 C2–H2 1.085 1.082 C6–H6 C3–H3 1.083 C5–H5 1.083 C4–H4 1.081
S1
S1 fit
D0
0.966 1.348 1.421 1.414 1.421 1.422 1.415 1.413 1.082 1.079 1.079 1.079 1.082
0.966 1.349 1.428 1.419 1.424 1.426 1.420 1.418 1.082 1.080 1.080 1.080 1.082
0.971 1.308 1.428 1.432 1.365 1.366 1.419 1.413 1.083 1.081 1.081 1.081 1.082
Angles (°) C1–O–H C2–C1–O C6–C1–O C1–C2–C3 C1–C6–C5 C2–C3–C4 C6–C5–C4 C2–C1–C6 C3–C4–C5 C1–C2–H2 C1–C6–H6 C2–C3–H3 C6–C5–H5 C3–C4–H4
109.6 120.0 116.3 117.6 117.8 119.1 119.0 123.8 122.7 120.2 119.3 120.7 120.7 118.5
109.5 119.8 116.2 117.7 117.8 118.5 118.4 124.0 123.7 120.2 119.0 121.4 121.3 118.0
113.8 122.6 116.0 118.6 118.9 120.1 119.8 121.4 121.2 119.7 118.3 120.4 120.4 119.2
5363 (5314) 2649 (2621) 1773 (1756)
5303 2648 1766
109.7 122.5 117.4 119.8 119.6 120.5 120.8 120.1 119.3 120.0 119.0 119.4 119.3 120.5
Rotational constants (MHz)b A 5687 (5651) B 2627 (2619) C 1797 (1790) a b
excitation spectra is also worthwhile because studying the differences between the simulated and measured excitation spectrum can provide a more detailed description of the excited state. For example, the frequencies that show up in the excitation spectrum can be used to correct the calculated force field of the excited state, while the differences between observed and predicted intensities can elucidate the importance of decay processes. In Fig. 3, comparison is made between the experimental S0 Resonance Enhanced MultiPhoton Ionization S1 (REMPI) (1 + 1) excitation spectrum recorded by Dopfer [23] (Fig. 3a) and the spectrum predicted by the present calculation (Fig. 3b). As expected, the agreement between the two spectra is quite good for the major in-plane modes (6a, 12, 1 and 18a). Similar to what was observed for the S1
Labeling according to Fig. 2. Experimental values in brackets [21,22].
3.2. S1
S0 REMPI excitation spectrum
As described before, the simulation of excitation spectra is intrinsically more difficult than the simulation of emission spectra. Apart from the fact that it is not always possible to obtain well-resolved emission spectra, simulation of
Fig. 2. Schematic representation of the change in geometry of phenol upon excitation to the S1 state.
Fig. 3. (a) Experimental S1 00 S0 excitation spectrum from [20]. (b) Franck–Condon simulation of the spectrum (calculated frequencies have been scaled with a factor of 0.956). (c). Simulated spectrum for the S1 geometry derived from the experimental intensities in the emission spectrum. Simulated spectra have been constructed by convolution of stick spectra with a Gaussian line profile with a width of 6 cm1.
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
00 ! S0 emission spectrum, the intensity of the transition to the fundamental of mode 12 is underestimated in the simulation. This is corrected if the simulation is performed with the S1 equilibrium geometry obtained from the fit to the emission spectrum (Fig. 3c). Despite the overall good agreement, there remain a number of eye-catching differences between the simulated and the experimental spectra. Firstly, the frequency of mode 16a is underestimated by 44 cm1. Because the first overtone of this mode is visible in the spectrum, this leads to a rather large difference in peak position between simulated and measured spectra (88 cm1). Secondly, a number of small peaks (all with intensities below 10% of the intensity of the 0–0 transition) are predicted for overtones and combinations of out-of-plane modes but cannot be found in the experimental spectrum. Both differences can be explained by an inaccurate description of the force field of the excited state for the out-of-plane modes. Similar to what was observed for the CAS-SCF calculations [24], we see in the present study that the frequency of almost all out-of-plane modes is underestimated by the TD-DFT calculation and that the differences with the experimental frequencies are much larger than for the in-plane modes. Additionally, inspection of the Duschinsky matrix shows that most out-of-plane modes are subject to severe mixing upon excitation. Predicted intensities for combination modes could differ from measured intensities if calculated and true mixing parameters differ slightly. Finally, the intensity of the transition to the mode that is assigned to 9a is severely underestimated (0.9%) in our simulation. Because the calculation does not reproduce the rather large decrease in frequency found experimentally for this mode [17] (1168 ! 975 cm1), one has to conclude that this mode is not described well by the calculation. 3.3. D0
S1 00
S0 REMPI-PES
Additional proof for the high quality of the S1 calculation can be found if we can project the excited state wavefunction onto another vibronic manifold than that of the ground state of the neutral molecule. This is done in a REMPI-PhotoElectron Spectroscopy (REMPI-PES) experiment. In such an experiment, the state under study – in this case the vibrationless level of S1 – is excited and subsequently ionized by one or more photons. The ionization process is assumed to occur vertically, leading to an ionic vibrational state distribution that is governed by the Franck–Condon factors between the vibrationless level of S1 and the various vibrational levels in D0. The quality of the calculated geometry and force field for D0 is usually of the same quality as the geometry and force field of S0. This makes the simulation of a REMPI-PES spectrum a sensitive tool for the quality of the excited state calculation. The simulation of the D0 S1 0 0 S0 REMPI-PES spectrum of Anderson et al. [25] is shown in Fig. 4. The equilibrium geometry of D0 that was employed for this
463
Fig. 4. Franck–Condon simulation of the D0 S1 00 S0 REMPI-PES. The simulated spectrum has been constructed by convoluting the stick spectrum with a Gaussian line profile with a width of 16 meV. The background in the experimental spectrum is modeled by a linear intensity decreasing toward higher photoelectron energies.
Table 2 Assignment of the dominant peaks in the photoelectron spectrum of the S1 00 level of phenol (energies in meV) Experiment b
52 64 101 121 150 207
Calculation
Modea
52 65 102 123 143, 148, 150 205
18b(5) 6a 12(1) 18a(1) b(9b), 9b(15), 9a 8a
a Modes are labeled according to Bist et al. In brackets the major contributing benzene mode as calculated by Grafton and Wheeler [18]. b Not explicitly labeled by Anderson et al.
simulation is given in Table 1. Although the resolution of the experimental spectrum is not as high as in the emission-detected spectra, the impressive agreement between the experimental and simulated spectra provides further confirmation of the accuracy of the S1 calculation.1 The interpretation of the spectrum by Anderson et al. is by far and large supported by our calculations. Table 2 lists the most important modes that contribute to the photoelectron spectrum. The strongest band in the spectrum is indeed associated with mode 8a (the quinoidal mode). On the other hand, although we agree with the conclusion that the OH bending mode is important – when the phenol becomes more quinone-like, the O atom becomes more
1
The origin of the background in the experimental spectrum is unclear, but is well modeled by a linear intensity decreasing toward high photoelectron energies.
464
M. de Groot, W.J. Buma / Chemical Physics Letters 420 (2006) 459–464
sp2 hybridized and the ring-OH angle approaches 120° – we find that this mode is mainly found in modes b and 9b instead of 8a as previously assumed. 4. Conclusions We have shown that Franck–Condon calculations using the geometries and force fields obtained with TD-DFT yield simulations of the excitation, dispersed emission spectra, and ionization spectra of the first excited state of phenol of much higher quality than have been obtained with high-quality ab initio calculations such as CAS-SCF calculations. In most aspects, the calculated and measured spectra exhibit almost quantitative agreement, demonstrating that the predictive power of the calculations is large and therefore able to provide a reliable assignment of the various transitions. At a somewhat more detailed level, we have found that fitting the calculated spectra to the experimental intensities leads to only minor adjustments in the excited state geometry. We are therefore presently applying the same methods to interpret the results of experimental studies on the excited states of chromophores of PYP [26], for which such extensive knowledge as available for phenol is not present. Acknowledgments We gratefully acknowledge Dr. Markus Gerhards (Heinrich-Heine Universita¨t Du¨sseldorf) for providing us with the S1 S0 REMPI excitation spectrum from the Ph.D. Thesis of Dr. Otto Dopfer. This work has been supported by Netherlands Organization for Scientific Research (NWO).
References [1] E.V. Doktorov, I.A. Malkin, V.I. Man’ko, J. Mol. Spectrosc. 64 (1977) 302. [2] S. Schumm, M. Gerhards, K. Kleinermanns, J. Phys. Chem. A 104 (2000) 10648. [3] T.D. Crawford, H.F. Schaefer III, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, VCH Publishers, New York, 2000, p. 33. [4] J.F. Stanton, R.J. Bartlett, J. Chem. Phys. 98 (1993) 7029. [5] E.V. Gromov, I. Burghardt, H. Ko¨ppel, L.S. Cederbaum, J. Phys. Chem. A 109 (2005) 4623. [6] M. Dierksen, S. Grimme, J. Chem. Phys. 120 (2004) 3544. [7] F. Furche, R. Ahlrichs, J. Chem. Phys. 117 (2002) 7433. [8] O. Treutler, R. Ahlrichs, J. Chem. Phys. 102 (1995) 346. [9] M.V. Arnim, R. Ahlrichs, J. Comp. Chem. 19 (1998) 1746. [10] P. Deglmann, F. Furche, R. Ahlrichs, Chem. Phys. Lett. 362 (2002) 511. [11] P. Deglmann, F. Furche, J. Chem. Phys. 117 (2002) 9535. [12] A. Scha¨fer, C. Huber, R. Ahlrichs, J. Chem. Phys. 100 (1994) 5829. [13] D. Becke, J. Chem. Phys. 98 (1993) 5648. [14] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [15] D. Gruner, P. Brumer, Chem. Phys. Lett. 138 (1987) 310. [16] P.R. Callis, J.T. Vivian, L.S. Slater, Chem. Phys. Lett. 244 (1995) 53. [17] H.D. Bist, J.C.D. Brand, D.R. Williams, J. Mol. Spectrosc. 24 (1967) 413. [18] K. Grafton, R.A. Wheeler, J. Comp. Chem. 19 (1998) 1663. [19] A.P. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502. [20] W. Roth, P. Imhof, M. Gerhards, S. Schumm, K. Kleinermanns, Chem. Phys. 252 (2000) 247. [21] N.W. Larsen, J. Mol. Struct. 51 (1979) 175. [22] G. Berden, W.L. Meerts, M. Schmitt, K. Kleinermanns, J. Chem. Phys. 104 (1996) 972. [23] O. Dopfer, Dissertation, Technische Universita¨t Mu¨nchen, Mu¨nchen, 1994. [24] S. Schumm, M. Gerhards, W. Roth, H. Gier, K. Kleinermanns, Chem. Phys. Lett. 263 (1996) 126. [25] S.L. Anderson, L. Goodman, K. Krogh-Jespersen, Ali G. Ozkabak, R.N. Zare, C. Zheng, J. Chem. Phys. 82 (1985) 5329. [26] M. De Groot, W.J. Buma, in preparation.