On-the-fly ab initio semiclassical evaluation of vibronic spectra at finite temperature

To compute and analyze vibrationally resolved electronic spectra at zero temperature, we have recently implemented the on-the-fly ab initio extended thawed Gaussian approximation [A. Patoz et al., J. Phys. Chem. Lett. 9, 2367 (2018)], which accounts for anharmonicity, mode-mode coupling, and Herzberg-Teller effects. Here, we generalize this method in order to evaluate spectra at non-zero temperature. In line with thermo-field dynamics, we transform the von Neumann evolution of the coherence component of the density matrix to the Schr\"{o}dinger evolution of a wavefunction in an augmented space with twice as many degrees of freedom. Due to the efficiency of the extended thawed Gaussian approximation, this increase in the number of coordinates results in nearly no additional computational cost. More specifically, compared to the original, zero-temperature approach, the finite-temperature method requires no additional ab initio electronic structure calculations. At the same time, the new approach allows for a clear distinction among finite-temperature, anharmonicity, and Herzberg-Teller effects on spectra. We show, on a model Morse system, the advantages of the finite-temperature thawed Gaussian approximation over the commonly used global harmonic methods and apply it to evaluate the symmetry-forbidden absorption spectrum of benzene, where all of the aforementioned effects contribute.


I. INTRODUCTION
2][3][4] The computational methods for simulating such spectra are, therefore, an essential tool in physical chemistry.
The most widespread is the global harmonic method, [5][6][7] which employs the harmonic approximation for both ground-and excited-state potential energy surfaces.][10][11][12] This approximation, however, neglects the effects of anharmonicity, which can significantly alter molecular spectra.Other quantum [13][14][15][16] and semiclassical 4,15,[17][18][19][20] methods do include anharmonicity effects on spectra, but at a substantial computational cost.Recently, we have been investigating the thawed Gaussian approximation, [21][22][23][24] an efficient semiclassical method that accounts partially for anharmonicity and requires no initial knowledge of the potential energy surface.6][27] Unfortunately, as a wavepacket propagation method, it has been limited to computing spectra in the zero-temperature limit, where only the ground vibrational state is populated initially.
To account for non-zero temperature, one typically employs the density matrix formalism, where a number of numerically exact [28][29][30] and approximate [31][32][33][34][35][36][37] approaches exist.9][40][41][42][43][44] Thermo-field dynamics 45,46 offers an alternative way to make wavefunction-based methods applicable at finite temperature: the problem, which seemingly requires the von Neumann equation for the density matrix, is mapped to a time-dependent Schrödinger equation with twice as many degrees of freedom.Recently, the thermo-field dynamics was employed in chemistry for solving the coupled electronicvibrational dynamics, [47][48][49][50] electronic structure, 51 and vibronic spectroscopy 11 problems.The application to vibronic spectroscopy, which is of central interest to this work, was, however, restricted to the global harmonic approximation.
Here, we combine the extended thawed Gaussian wavepacket propagation with the thermo-field dynamics in order to include both anharmonicity and finite-temperature ef-fects.Due to the favorable scaling of the thawed Gaussian approximation with the system's size, the new method adds nearly no additional cost to the original, zero-temperature approach.To illustrate the accuracy achieved by going beyond both global harmonic and zero-temperature approximations, we test the method on a set of Morse potentials with different degrees of anharmonicity and at different temperatures.Finally, we apply it to evaluate the spectrum corresponding to the symmetry-forbidden electronic transition ) of benzene and demonstrate that the simultaneous inclusion of Herzberg-Teller, anharmonicity, and finite-temperature effects is needed to reproduce the experimental spectrum.

II. THEORY
A. Extended thawed Gaussian approximation for zero-temperature spectra Before turning to vibrationally resolved electronic spectra at finite temperature, let us briefly describe the original, zero-temperature approach based on the extended thawed Gaussian approximation.
The absorption spectrum at zero temperature can be computed as the Fourier transform, 33,52,53 of the correlation function, where |1, g is the ground ("g") vibrational state of the ground ("1") electronic state, ω 1,g is its energy, Ĥ2 is the nuclear Hamiltonian corresponding to the excited (subscript "2") electronic state, and μ is the transition dipole moment μ21 = ˆ µ 21 • projected on the polarization of the external electric field.In other words, to compute the spectrum, one has to evolve the nuclear wavefunction |φ 0 = μ|1, g on the excited-state surface, which is, in general, a challenging task that scales exponentially with the number of atoms.
Different exact quantum 42,[54][55][56][57] and semiclassical 19,[58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74] methods were developed for solving the problem of wavepacket propagation.Sometimes, the region of the excited-state potential energy surface explored by the evolved wavepacket is fairly harmonic, meaning that it can be well approximated by a second-order Taylor expansion in nuclear coordinates about a reference geometry; we call this the global harmonic approximation.5,76 Moreover, in this case the excited-state surface is easily constructed from a single Hessian calculation.The prevalence of the global harmonic method in the vibronic spectroscopy literature 7,8,10,[75][76][77][78][79] testifies, on the one hand, to its applicability in a wide range of molecules and, on the other hand, to the absence of accessible alternatives that can account for anharmonicity effects.Often, due to the reduced resolution of electronic spectra, even such a crude approximation, which would nowadays be almost unacceptable for the simulation of vibrational (infrared) spectra, is considered appropriate.
To account for the anharmonicity effects on the spectrum at least approximately, we recommend using the simple and efficient semiclassical thawed Gaussian approximation. 21In contrast to many other exact or approximate quantum dynamics methods, this method is computationally feasible even for rather large molecules and can be employed in a "blackbox" fashion, i.e., it requires little human input.In particular, the thawed Gaussian approximation requires only local potential energy information along the classical trajectory (as described below) and, therefore, supports an on-the-fly implementation where the potential energy is provided by an ab initio electronic structure calculation.
Within the thawed Gaussian approximation, 21 the wavepacket is assumed to be a complex Gaussian function parametrized by the time-dependent D-dimensional real vectors q t and p t , D × D complex symmetric matrix A t , and complex number γ t ; D is the number of degrees of freedom.
The time dependence of the matrix A t implies that the width of the thawed Gaussian wavepacket changes with time, as opposed to the frozen Gaussian ansatz where the width remains constant.Wavepacket (3) solves exactly the Schrödinger equation where T (p) = 1 2 p T • m −1 • p is the kinetic energy and is the local harmonic approximation of the true potential energy V (q) about q t , if the timedependent parameters of ψ t satisfy the following equations of motion: In these equations, V (q t ) and V (q t ) denote, respectively, the gradient and Hessian of the potential energy evaluated at q t , m is the symmetric mass matrix, and L t = T (p t ) − V (q t ) is the Lagrangian.Note that due to the local harmonic approximation, Eq. ( 4) is a nonlinear Schrödinger equation because the potential V LHA (q) depends on the wavefunction through the parameter q t [Eq.( 5)], i.e., V LHA (q) ≡ V LHA (q; q t ) ≡ V LHA (q; ψ t ).
To construct the initial wavepacket, the ground-state potential energy surface V 1 (q) is assumed to be harmonic in the vicinity of its minimum q eq , i.e., the ground-state Hamiltonian is approximated as where K = V 1 (q eq ) is the symmetric force-constant matrix and ∂ q = ∂/∂q.In position representation, the lowest eigenstate ψ 0 (q) of the Hamiltonian (10) is a Gaussian (3) with parameters where . This initial wavefunction ψ 0 (q) is then evolved by solving differential equations ( 6)-( 9) with V = V 2 (the excited-state potential energy).
The thawed Gaussian wavepacket (3) is not suited to treat non-Condon effects, i.e., the effects due to the dependence of the transition dipole moment µ(q) on nuclear coordinates q.Within the Herzberg-Teller approximation-the simplest extension of the Condon approximation-the transition dipole moment is assumed to be a linear function 80 µ(q) = µ(q eq ) + µ (q eq ) T • (q − q eq ), (15)  where µ (q eq ) is the gradient of µ with respect to nuclear coordinates at the equilibrium geometry.Then, φ 0 (q) = µ(q)ψ 0 (q) is no longer a Gaussian wavepacket.Fortunately, the extended thawed Gaussian ansatz, 25,81 which is a special case of Hagedorn's "Gaussian times a polynomial" wavepacket, [82][83][84] solves the same Schrödinger equation [Eq.( 4)] as ψ t (q), provided that the Gaussian parameters evolve, as before, according to Eqs. ( 6)-( 9) and, in addition, Hence, with the extended thawed Gaussian approximation, one can include the Herzberg-Teller contribution at nearly no additional computational cost.

B. Vibrationally resolved electronic spectra at finite temperature
At non-zero temperature, the dipole-dipole correlation function needed in vibronic spectroscopy is where ρ = e −β Ĥ1 /Tr(e −β Ĥ1 ) is the vibrational density operator and β = 1/k B T .Note that in Eq. ( 19) we assumed that only the ground electronic state is populated in the thermal equilibrium, which is usually justified by the large energy gap between the ground and first excited electronic states.Because the time evolution in Eq. ( 19) involves two different Hamiltonians, an obvious classical analogue of Eq. ( 19) is missing, which, in turn, hinders the development of classical-like or semiclassical approximations for C(t).Here, we demonstrate that by transforming the problem to the one of wavepacket dynamics in an augmented space one can easily make use of the existing semiclassical methods for solving the time-dependent Schrödinger equation.
The correlation function can be re-written as = dq φ0 (q) * φt (q), ( where = e −iH 2 (q)t/ e iH 1 (q )t/ q|μρ 1/2 |q (25 q = (q, q ) T is a 2D-dimensional coordinate vector, and is a Hamiltonian in q coordinates.In Eq. ( 20), we used the relation [ρ, Ĥ1 ] = 0 and the cyclic property of the trace; in Eq. ( 21), we introduced the position representation in q and q coordinates; in going from (25) to (26), we used the fact that the two Hamiltonians H 1 (q ) and H 2 (q) commute because they act on different coordinates; finally, Eq. ( 22) follows from (21) because and (ρ 1/2 ) † = ρ1/2 .Equation ( 22) has a remarkable interpretation-the dipole-dipole correlation function C(t) for a D-dimensional system at finite temperature T can be thought of as a wavepacket autocorrelation function of φt (q) evolved with the Hamiltonian H(q) according to the Schrödinger equation, which describes an effective 2D-dimensional system at zero temperature.
The approach described here is, despite the explicit use of the position representation, equivalent to the thermo-field dynamics, as presented in Ref. 11.Indeed, the final result does not depend on the representation: where |q = |q |q is a general position state in the augmented direct-product Hilbert space and |q denotes a position state in the "fictitious" (or "tilde") Hilbert space.In Appendix A, we derive Eq. ( 31) using standard thermo-field dynamics notation and without invoking the position representation.
In principle, any known method for solving the time-dependent Schrödinger equation can be applied to obtain φt .However, the doubled number of coordinates adds a substantial, if not prohibitive, computational cost to the already large cost of zero-temperature calculations with exponentially-scaling exact quantum methods.In the following, we therefore employ the extended thawed Gaussian approximation, which scales favorably with the number of degrees of freedom.

C. Extended thawed Gaussian approximation for finite-temperature spectra
To solve Eq. ( 30) with the (extended) thawed Gaussian approximation, we must first identify φ0 and the local harmonic approximation to the potential energy If we assume, as in Sec.II A, that the ground-state surface V 1 is harmonic [Eq.( 10)], a gen- and a scalar See Appendix B for the derivation of Eqs.(32)- (35).The harmonic approximation for the ground-state potential energy surface is justified in the vicinity of its minimum and, therefore, for the construction of the equilibrium vibrational density matrix.6][87][88] Next, we assume μ to be diagonal in position representation, φ0 (q) = µ(q)ρ 1/2 (q), (36)   and employ the Herzberg-Teller approximation [Eq.( 15)] to obtain With these initial values, we propagate the time-dependent parameters qt , pt , Āt , and γt according to Eqs. ( 6)-( 9) and bt according to Eq. ( 18).The potential energy, its gradient, and its Hessian are given by while the where q t and q t are D-dimensional vectors composed of the first and second halves of coordinates of qt , i.e., qt = (q t , q t ) T .Interestingly, the classical equations of motion [Eqs.( 6) and (7)] for the parameters qt and pt are solved by propagating two independent trajectories in D spatial dimensions: the first trajectory evolves q t and p t with the excited-state Hamiltonian H 2 , while the second trajectory evolves q t and p t with the negative of the ground-state Hamiltonian, −H 1 , due to the negative signs of mass in Eq. ( 42) and gradient in Eq. (40).
Because the second trajectory is at a fixed point, i.e., at the minimum of the ground-state potential energy V 1 with zero momentum, it shows no dynamics.As a result, one requires only a single excited-state classical trajectory, to evolve the first D coordinates of qt , and Hessians of the excited-state potential energy surface along this trajectory, which is the same as in the original zero-temperature approach; no further potential energy evaluations are needed to account for the temperature effects.
Finally, let us note that an alternative approach to finite-temperature spectra with Gaussian wavepackets has been proposed in Ref. 12. There, the authors propose a similar scheme to directly evolve the coherence ρµ (t) = exp(−i Ĥ2 t/ )μρ exp(i Ĥ1 t/ ) in position representation with the doubled number of degrees of freedom.Then, the correlation function is evaluated simply as Their method, combined with the local harmonic approximation, is equivalent to ours and gives the same correlation function.In contrast to our approach, the method of Ref. 12 has so far been used to compute vibronic spectra only in systems described with globally harmonic potential energy surfaces, where it is equivalent to the global harmonic approximation for vibronic spectra, which is much simpler because analytical expressions for C(t) exist. 10,75Our approach based on thermo-field dynamics has the advantage of reducing the transition dipole autocorrelation function to the simpler and well-known expression (31) for the wavepacket autocorrelation, thus making it very easy to implement the finite-temperature treatment of vibronic spectra into the standard zerotemperature wavefunction-based codes, which typically contain procedures for computing the wavepacket autocorrelation.

III. COMPUTATIONAL DETAILS A. Morse potential
To test the accuracy of the proposed method, we construct a one-dimensional model system consisting of a ground-state harmonic potential and an excited-state Morse potential.
The ground-state surface is assumed harmonic to exclude the error (or error cancellation) due to using an approximate initial vibrational state-this is rarely an issue with zerotemperature methods because the harmonic approximation typically holds in the vicinity of the potential minimum but could affect the results at higher temperatures.In the current model, the error of the results obtained with thawed Gaussian approximation is only due to the anharmonicity of the excited-state potential energy surface.
A set of Morse potentials, was constructed by fixing the equilibrium position q 2 , minimum energy V 2 (q 2 ), and frequency )/m at q 2 and by varying the anharmonicity parameter χ.We set the minimum of the ground-state harmonic potential to zero (q 1 = 0) and its frequency to ω 1 = 1, while the excited-state Morse parameters were q 2 = 1.5, ω 2 = 0.9, and V (q 2 ) = 10.Mass was set to m = 1.The level of anharmonicity was tuned by changing the parameter χ in the range between 0.01 and 0.02, in steps of 0.001.
The exact spectrum was computed by evaluating Franck-Condon factors by numerical integration, which is feasible for this one-dimensional model system since both harmonic and Morse vibrational eigenfunctions are known analytically.The adiabatic harmonic model, which is constructed about the minimum of the potential energy surface, is the same for all constructed Morse potentials because it does not depend on χ.Since the (extended) thawed Gaussian approximation is exact for harmonic potentials, it was used to compute the adiabatic harmonic spectra.For both harmonic and thawed Gaussian dynamics calculations, time step was 0.1 and the total simulation time was 1000, i.e., 10000 steps in total were taken.Gaussian broadening with half-width at half-maximum of 0.1 was applied to all spectra.Spectra were evaluated at scaled temperatures T ω = 0, 0.5, and 1, where T ω = k B T / ω 1 = 1/β ω 1 (e.g., for an average molecular vibration of ω = 1000 cm −1 , T ω = 1 corresponds to the temperature T ≈ 1439 K).A constant transition dipole moment µ = 1 was used.
To compare reference (σ ref ) and approximate (σ) spectra, we used the spectral contrast angle θ, defined through its cosine as where σ 1 • σ 2 = dωσ 1 (ω)σ 2 (ω) is the inner product of two spectra and σ = √ σ • σ is the associated norm.In all calculations, the reference was the exact spectrum, while the approximate spectra were computed with the adiabatic global harmonic and thawed Gaussian approximations.

B. On-the-fly ab initio calculations
The S 1 ← S 0 absorption spectrum of benzene was computed with adiabatic harmonic, vertical harmonic, and thawed Gaussian approximations.In short, the adiabatic harmonic model is, as described above, obtained by the second-order Taylor expansion of the excitedstate potential energy surface about its minimum, while for the vertical harmonic model, the same expansion is performed about the ground-state minimum.
Density functional theory was used for the optimization and Hessian calculation of the ground electronic state, while its time-dependent version was employed for the excited-state optimization, energy, gradient, and Hessian calculations.We used the B3LYP functional with the ultrafine grid and 6-31+G(d,p) basis set, as implemented in the Gaussian09 89 package.For the thawed Gaussian propagation, we used a second-order symplectic integrator with a time step of 8 a.u.(≈ 0.2 fs) and 10000 steps in total.The Hessian of the potential energy was evaluated every four steps and interpolated in between, as done previously in Ref. 25.The ground-state surface was assumed to be harmonic.The gradient of the electronic transition dipole moment was computed numerically by the second-order finite difference method with a step of 10 −4 Å. 25 The computed correlation functions were multiplied by an exponential damping function e −t/τ with τ = 18000 a.u., resulting in a Lorentzian line shape with half-width at half-maximum of ≈ 12.2 cm −1 .To facilitate comparison with the experimental spectrum, computed spectra were shifted and scaled to match the experimental spectrum of Ref. 90 at its highest peak (data taken from the MPI-Mainz UV/VIS Spectral Atlas 91,92 ).
Finally, let us emphasize that the finite-temperature treatment of spectra requires no additional electronic structure evaluations, i.e., the same ab initio data could be reused to compute the benzene spectrum at any given temperature.We evaluated the spectra at zero temperature and at the temperature of the experiment (T = 298 K).

A. Morse potential
Thawed Gaussian and global harmonic spectra were compared with the exact result (see Fig. 1).Already for the system with weak anharmonicity (left panels, χ = 0.01), the thawed Gaussian approximation provides a more accurate spectrum than the harmonic method.The difference is seen mainly in the intensities of the high-frequency peaks.Since the adiabatic harmonic model describes well the region around the potential minimum, it can recover the positions and intensities of peaks corresponding to transitions between vibrational states with small quantum numbers.In contrast, the harmonic approximation breaks down for vibrational states with more quanta, resulting in incorrect intensities of highfrequency transitions.The effect of anharmonicity on the peak positions becomes significant for χ = 0.02 and even the thawed Gaussian approximation is inadequate.Nevertheless, it is still more accurate than the adiabatic harmonic model, which has no dependence on χ (harmonic spectra are, clearly, the same for different χ at a given temperature).In typical molecular systems, the peaks are often left unresolved due to the short excited-state lifetime or inhomogeneous broadening.Then, the intensities play an important role in recovering the overall shape of the spectrum, whereas even an error of tens of reciprocal centimeters in peak positions can be tolerated.
In contrast to the global harmonic method, the thawed Gaussian approximation can result in non-physical negative spectral features, which are due to the nonlinear character of the Schrödinger equation ( 4).This is a well-known disadvantage of the method and was discussed in more detail elsewhere. 23,24In the studied Morse system, a negative peak overlaps with the hot band around ω = 9, resulting in poor description of this spectral region at higher temperatures.In a way, the gain in accuracy in the high-frequency part of the spectrum is accompanied by a loss in accuracy in the frequency region below the 0-0 transition.
To compare the global harmonic and thawed Gaussian methods quantitatively, we measure the error of an approximate spectrum with the spectral contrast angle between the approximate and exact spectra (see Fig. 2).The thawed Gaussian approximation gives more accurate spectra than the harmonic approximation for all anharmonicities and at all temperatures studied.However, an interesting trend is observed: the harmonic approximation becomes more accurate as the temperature increases, whereas the thawed Gaussian approximation keeps the same degree of accuracy at all temperatures.The main reason for such behavior is closely related to the discussion above.As the temperature increases, the intensity of hot bands below the 0-0 transition grows and they become more relevant in measuring the error.Hence, the adiabatic harmonic method gains on accuracy, unlike the thawed Gaussian approximation, which always loses on accuracy in the low-frequency part of the spectrum.However, such behavior of the global harmonic method is not general; if the ground-state potential energy surface were anharmonic, high-temperature spectra would also reflect the effects neglected in the global harmonic models-those of ground-state anharmonicity on the initial density matrix.

B. Absorption spectrum of benzene
The symmetry-forbidden S 1 ← S 0 transition in benzene is a well-known example of the Herzberg-Teller effect, 1,2 where the spectrum arises only due to the coordinate dependence of the transition dipole moment, which is zero by symmetry at the equilibrium geometry.As such, it has been studied extensively both from the experimental 90,[93][94][95][96][97] and  45)], of the spectra computed with the thawed Gaussian approximation ("TGA", Sec.II B) or the adiabatic harmonic approach [Eq.( 44)], as a function of the anharmonicity parameter χ. Results are shown for three different temperatures theoretical 34,[98][99][100][101][102][103][104] points of view.The spectrum is a challenge for computational methods because it is highly resolved, exhibits Herzberg-Teller effects, and contains hot bands due to finite temperature.Although benzene is typically considered to be a rigid molecule, we have recently shown that the anharmonicity affects significantly the intensities of the peaks in the main progression of the spectrum. 25,26However, our previous work assumed zero temperature, therefore neglecting the weak hot bands present in the experimental spectrum.
Here, we complement our earlier result with the new finite-temperature extended thawed Gaussian method.First, we demonstrate [Fig.3(a)] the effect of non-zero temperature on the spectrum.Whereas the original, zero-temperature extended thawed Gaussian approximation neglects completely the weak, but non-negligible, hot bands, the finite-temperature approach reproduces all features of the spectrum.The inaccuracy in the frequencies of the peaks is most likely due to the electronic structure method used; we discuss this later.Nevertheless, Fig. 3(a) clearly shows the difference in the spectra computed without and with finitetemperature effects.
We argue that the benzene absorption spectrum is affected by the anharmonicity of the excited-state potential energy surface.This effect is best demonstrated by the difference in spectra based on two global harmonic models: if the potential energy surface were harmonic, the second-order expansion of the potential energy about any molecular geometry would result in the same spectrum.As shown in Fig. 3(b), in benzene, the adiabatic harmonic method is much more accurate than the vertical; in general, either of the two methods can be more appropriate. 22,23,105The extended thawed Gaussian approximation outperforms not only the vertical harmonic approach, whose spectrum is completely off, but also the adiabatic harmonic approximation, which fails to produce accurate peak intensities.
Figure 3(c) shows the importance of treating the Herzberg-Teller effect with the extended thawed Gaussian approximation.Since the transition is symmetry-forbidden, i.e., µ(q eq ) = 0, the spectrum computed within the Condon approximation [µ(q) ≈ µ(q eq )] vanishes, whereas the full, Herzberg-Teller treatment reproduces the experimental spectrum.
In computational chemistry, vibrational scaling factors, 106 which we denote by f , are often used to empirically correct for systematic errors in the vibrational frequencies computed with electronic structure methods.In vibronic spectroscopy, such scaling, applied to ground-and excited-state frequencies, can modify both peak positions and intensities. 10However, the effect on intensities is often weak; indeed, the adiabatic harmonic spectrum with scaled vibrational frequencies (red, dashed line in Fig. 4) exhibits almost perfect peak positions, but still the same errors in intensities as the adiabatic harmonic spectrum of Fig. 3(b).
For comparison-and for comparison only-we show an analogous, "corrected" spectrum computed with the extended thawed Gaussian approximation (blue, solid line in Fig. 4).
Since the simple procedure of scaling the vibrational frequencies is not applicable in this case, we scale directly the frequency axis by f , which corrects peak positions but leaves intensities unchanged.The results imply that the subtle anharmonicity effects on spectral intensities, described well with the on-the-fly semiclassical thawed Gaussian method, cannot be captured even with the corrected harmonic potential.i.e., that the position representation of state | φ0 introduced in this appendix is the function defined in Eq. ( 23) of the main text.Indeed, where we again used Eq.(A3) with Â being identity operator.

Appendix B: Derivation of the initial-state parameters
Here, we derive expressions (32)-( 35) for the Gaussian parameters of ρ 1/2 (q, q ).First, recall that the matrix element ρ(q, q ) of the thermal density operator in a one-dimensional harmonic oscillator Ĥ with mass m and frequency ω is 107 ρ(q, q ) = mω tanh(β ω/2) π Using Eq. (B2), we now derive the expression for the matrix element ρ(q, q ) in a D- where m and K are D × D symmetric mass and force-constant matrices, respectively, and q eq is the equilibrium position at which the potential energy has its minimum.

FIG. 3 .
FIG.3.Benzene S 1 ← S 0 absorption spectrum computed with the extended thawed Gaussian approximation ("Extended TGA") at 298 K (using the approach described in Sec.II C), compared with the experimental spectrum90,92 measured at 298 K and other approximate spectra simulations based on: (a) zero-temperature extended thawed Gaussian approximation ("Extended TGA 0 K") as described in Sec.II A, (b) adiabatic or vertical global harmonic models at 298 K (see Sec. III B), and (c) thawed Gaussian approximation, which assumes Condon approximation ["TGA (Condon)"].

FIG. 4 .
FIG.4.Benzene S 1 ← S 0 absorption spectra computed with the extended thawed Gaussian approximation ("Extended TGA") and adiabatic harmonic model, both at 298 K, compared with the experimental spectrum90,92 measured at 298 K.The adiabatic harmonic model was modified by scaling both ground-and excited-state frequencies by a constant f = 0.963, which was taken from Ref. 106 and is associated with the electronic structure method used (see Sec. III B).For the spectrum evaluated with the extended thawed Gaussian approximation, we applied the same scaling factor only to the values on the frequency axis.