Aiming at an accurate prediction of vibrational and electronic spectra for mediumtolarge molecules: An overview

Consiglio Nazionale delle Ricerche, Istituto di Chimica dei Composti OrganoMetallici (ICCOM-CNR), UOS di Pisa, Area della Ricerca CNR, Via G. Moruzzi 1, Pisa 56124, Italy Scuola Normale Superiore of Pisa, Classe di Scienze, Piazza dei Cavalieri 7, Pisa 56126, Italy International Centre for Quantum and Molecular Structures, College of Sciences, Shanghai University, 99 Shangda Road, Shanghai 200444, China Correspondence Julien Bloino, Consiglio Nazionale delle Ricerche, Istituto di Chimica dei Composti OrganoMetallici, Area della Ricerca CNR, Via G. Moruzzi, 1, Pisa, Italy. Email: julien.bloino@pi.iccom.cnr.it; Alberto Baiardi: alberto.baiardi@sns.it; Malgorzata Biczysko: biczysko@shu.edu.cn Funding Information This work was supported by Italian MIUR (FIRB 2012: “Progettazione di materiali nano-eterogenei per la conversione di energia solare,” protocollo: RBFR122HFZ and PRIN 2012,”STAR: Spettroscopia e Tecniche computazionali per la ricerca Astrofisica, atmosferica e Radioastronomica,” protocollo: 20129ZFHFE). J.B. and M.B. also acknowledge the COST CMTS-Action CM1405 (MOLIM: MOLecules In Motion). Support from Gaussian, Inc. is greatly appreciated. Abstract In this tutorial review, we present some effective methodologies available for the simulation of vibrational and vibrationally resolved electronic spectra of medium-to-large molecules. They have been integrated into a unified platform and extended to support a wide range of spectroscopies. The resulting tool is particularly useful in assisting the extensive characterization of molecules, often achieved by combining multiple types of measurements. A correct assessment of the reliability of theoretical calculations is a necessary prelude to the interpretation of their results. For this reason, the key concepts of the underlying theories will be first presented and then illustrated through the study of thiophene and its smallest oligomer, bithiophene. While doing so, a complete computational protocol will be detailed, with emphasis on the strengths and potential shortcomings of the models employed here. Guidelines are also provided for performing similar studies on different molecular systems, with comments on the more common pitfalls and ways to overcome them. Finally, extensions to other cases, like chiral spectroscopies or mixtures, are also discussed.

provide a simple interface, lowering the barrier to the broad adoption of state-of-the-art methods able to produce accurate data for the problems at hand. [11,[15][16][17] This can be achieved by automatizing most aspects of the calculations. Due to the inherent complexity of some theoretical models, such a goal may not be attainable. In this case, the necessary parametrization to the problem of interest should be made as simple as possible. A second requirement is to be able to choose the most appropriate methodology based on the size and topology of the molecular system under study, the target accuracy, and the complexity of the environment and its possible effects on the system.
In this article, we will focus on cost-effective methods, suitable for the simulation of electronic or vibrational spectra of medium-to-large molecular systems, ranging from half-a-dozen to more than a 100 atoms, and how they can be integrated within a unified platform, that we will call the virtual spectrometer. [11,[15][16][17] While those methods are primarily intended for semirigid molecules, their extension to support also more flexible systems will be discussed.
To understand the choice of the methods described here, a brief review of the theoretical context is necessary. First of all, it should be stressed out that the theoretical data used to support experimental spectroscopic studies are still mainly related to the direct outcome of standard computational approaches, such as geometry optimizations and harmonic frequency computations. They provide a model for the potential energy surface (PES) and can be complemented by additional equilibrium properties, such as the vertical excitation (VE) energies and transition moments, which are now available for a large number of electronic structure quantum mechanical (ESQM) approaches and in several computational packages. [18][19][20][21][22][23][24][25] However, such data are based on an approximated description of spectroscopic phenomena, for instance by neglecting anharmonic effects in vibrational transition energies and intensities, or the vibrational effects associated to electronic transitions. Indeed, the actual PESs are never harmonic and the vibrational effects are present even at the absolute zero. Thus, (i) molecular structures observed experimentally do not correspond to the computed equilibrium geometries (r eq ) but include vibrational effects (r 0 ), (ii) harmonic vibrational energies do not match experimental fundamental transitions (they usually overestimate), and intensities of overtones and combination bands are null, (iii) all electronic spectra band-shapes are the result of a potentially large number of vibronic transitions, which may not be clearly visible. Here, we will discuss theoretical models which take into account all these important effects related to the vibrations of the molecules and allow a direct comparison between simulated and experimental spectra. [26][27][28][29] To fulfill the accuracy (for structures and frequencies) and interpretability (for intensities) requirements of vibrational spectra, it is mandatory to go beyond the so-called double-harmonic approximation accounting for both mechanical and electric (and more generally property-related) anharmonic effects. For very small molecules, it is possible to compute ro-vibrational energy levels using full variational approaches [14,30,31] while simplified, less expensive methodologies are necessary for larger molecular systems (see also reference 32 and references therein), either based on variational [33][34][35][36] or perturbative [37][38][39][40][41][42][43][44][45][46][47][48] schemes. Here, we will employ the framework of the generalized second-order vibrational perturbation theory (GVPT2), which permits the computation of thermodynamic properties, vibrational energies, and transition intensities (for infrared (IR), Raman, vibrational circular dichroism and Raman optical activity) from the vibrational ground state to fundamentals, overtones and combination bands up to the three quanta. [17,29,[49][50][51][52] Computations of vibrationally resolved electronic spectra are even more challenging since at least two PESs corresponding to different (initial and final) electronic states need to be taken into account. For this case as well, very accurate variational anharmonic calculations based on a ro-vibronic Hamiltonian in internal coordinates [53] or by means of discrete variable representation (DVR)-based approaches, [54] and considering several interacting electronic states beyond the Born-Oppenheimer (BO) approximation, are only feasible for the smallest molecules. [54][55][56][57][58]  For larger molecular systems, most theoretical models are based on the harmonic description of PESs, using either the time-independent (TI) [59] or time-dependent (TD) [60,61] formalisms, with the possible inclusion of anharmonic corrections to the band positions, [62] or the anharmonic description of some specific vibrational degrees of freedom. [63][64][65] The TI approach is based on a sum-over-states formalism, where the spectrum is obtained as an ensemble of distinct transitions between the vibrational initial and final states. This is the method of choice when high-resolution spectra are available or band identification is required, but can suffer from convergence issues related to the band-shape, especially when temperature effects are of interest. In such circumstances, the TD route is better suited since it exploits the properties of the Fourier transform to achieve fully converged spectra including also temperature effects, without additional computational cost. In this work, we will employ both strategies, through an approach supporting both TI and TD formulations [61,66,67] and based on algorithms allowing to take into account Franck-Condon (FC), Herzberg-Teller (HT), Duschinsky and anharmonic effects, with the consequent availability of both highresolved and fully converged spectra at finite temperatures.
While the approaches described above to simulate vibronic spectra are effective for semirigid systems, their reliability decreases as the structural changes involved in the electronic transition become important, for instance with large amplitude motions (LAMs) present. In fact, those deformations are poorly described in terms of the Cartesian-based normal modes computed at the harmonic level. A full anharmonic treatment of vibronic transitions is usually feasible only for small-size systems, [68,69] so more affordable methods must be developed to tackle large-size systems as well. A less expensive, but still effective solution is offered by the use of internal coordinates, which are better-suited to represent LAMs. In fact, even if in vibrational spectroscopies the difference between Cartesian and internal coordinates emerges only at the anharmonic level, [70,71] this is not true for vibronic spectroscopy. In this case, two different structures (the equilibrium geometries of the electronic states involved in the transition) must be taken into account and the definition of the Duschinsky transformation changes already at the harmonic level. [72] As shown by some of us, [73,74] the definition of the Duschinsky transformation in internal coordinates leads to more reliable results for flexible systems. An additional advantage of using internal coordinates is the reduction of the couplings between the large-amplitude mode(s) and the other vibrational degrees of freedom. For simple deformations, the large-amplitude motion can be associated to a well-defined single mode, which is uncoupled from the other ones. It is then possible to develop hybrid schemes, where the monodimensional LAM is treated using high-accuracy, anharmonic variational approaches, and the other modes at the harmonic level. [63] Based on those considerations, we present here a new method to simulate vibronic spectra of flexible systems undergoing monodimensional, large-amplitude deformations. Within this approach, the reaction path Hamiltonian (RPH) [75][76][77] is used to describe the motion along the LAM, and the vibrational energies are determined through a variational method, which relies on a DVR-based approach. [78,79] Within the so-called adiabatic approximation, the couplings between the LAM and the other degrees of freedom can be neglected, [80,81] with the latter treated using the TI approach described above. Such an hybrid model can be straightforwardly extended to support multiple LAMs if the coupling between them is negligible. A more complete treatment would require the use of the so-called reaction surface hamiltonian (RSH) model, [82][83][84] and will not be discussed here.
The simulation of the vibrational and vibrationally resolved electronic spectra of thiophene (T) and 2,2 0 -bithiophene (2T) will cover cases ranging from semirigid systems and small PES changes on electronic transition up to floppy systems with a well-defined large-amplitude degree of freedom. Those molecules have been extensively studied experimentally [85][86][87][88][89][90] and used as benchmarks to test theoretical models. [91][92][93][94][95][96] The importance of a direct comparison between theoretical and experimental spectra have been recently demonstrated from a joint study on the phosphorescence of 2,2 0 -bithiophene. [95] Indeed, previous works showed discrepancies between theory and experiment, [97] which was found out to be caused by inaccurate measurements. Once redone, a nearly quantitative agreement was reached for the phosphorescence spectrum of 2T. In this respect, there are still experimental spectra of T and 2T, which require further theoretical investigations. For instance, the recently recorded high-resolution photo-absorption Fourier transform spectrum of thiophene [89] calls for a more detailed assignment of its vibrational structure. Thiophene and its oligomers and derivatives are receiving also significant attention as important "building blocks" molecules, being the simplest heteroaromatic cycles of interest in astrochemistry [98,99] or prototypes of p-conjugated systems for technological applications. [90,[100][101][102] From a computational point of view, the smallest oligomer, which is 2T, is an ideal test-case to check the reliability of the DVR-based approach described above. In fact, the dihedral angle between the two thiophene moieties can change significantly on electronic excitation, [103] and therefore an internal coordinates-based treatment is required to reproduce the experimental vibronic spectra. [86] The manuscript is organized as follows. The first section presents key aspects of the theoretical background used for the simulation. This will pave the way to a more practical discussion in the next section on the actual simulation of vibrational and vibronic spectra, including important features or possible shortcomings one should keep in mind when analyzing the results. This will be illustrated with the examples of thiophene and bithiophene, where application of the models described here leads to a direct vis-a-vis comparison with experimental results. Some further discussion concerning chiroptical spectra and rotorlike molecules as well as final remarks on the current capabilities of such a virtual spectrometer and possible strategies for future improvements will conclude this presentation.

| T HE O RE TI CA L BA CK GR O U N D
An extensive presentation of the theoretical backgrounds underlying both anharmonic vibrational and vibrationally resolved electronic spectroscopies bear the risk of obscuring important aspects, which need to be known when assessing the reliability of simulations. For this reason, emphasis will be put on specific parts of the methodologies implemented in the virtual spectrometer. Interested readers can refer to the BLOINO ET AL. | 3 papers cited in the following for more in-depth and technical discussions on each feature presented here.
Before starting, we should emphasize that, to support easily multiple spectroscopies and facilitate future extensions, the tool described and used here relies on general formulations and methodologies. However, since we will only need some spectroscopies, significant simplifications can be made. Their impact on the overall discussion will be indicated where necessary.
Let us consider the most general integral of an arbitrary operator H between two molecular states, j W I i and j W F i, Within the BO approximation, the molecular wavefunction can be written as the product of a nuclear (w) and an electronic (/) components, so the previous integral can be rewritten, where the orthonormality of the electronic wave function has been used to simplify the denominator (h / I j/ F i5d IF ). Note that a shortcut has been taken in the notation regarding the subscripts, which referred to the molecular state. A more accurate representation would require the use of different indexes for each state, at the cost of the overall readability. For this reason, and until explicitly needed, this simplified notation will be kept.
The nuclear wavefunction can be further factorized if the Eckart-Sayvetz conditions are assumed met, which can be considered a good approximation for semirigid molecules. [104,105] In this case, it is possible to factor out the translational and rotational components, and only keep the vibrational wavefunction, which will be at the center of our discussion here. Within this framework, Equation 1 becomes, Within the harmonic approximation, the vibrational wavefunctions can be written as products of one-dimensional harmonic oscillators corresponding to the normal modes of vibrations. While such a calculation is rather straightforward for vibrational spectroscopies, where only one electronic state is involved (j / F i5j / I i), this is not the case for electronic spectroscopies since the initial and final vibrational states are generally not expressed in the same basis of normal coordinates.
In the following, except if stated otherwise, we will assume that the property or quantity of interest and the associated electronic transition dipole moment (h / I jĤj/ F i) are known, with the latter noted H e for simplicity.

| Vibrational spectroscopy
While calculation of vibrational spectra at the harmonic level has become routine and is commonly used to interpret experimental spectra, it suffers from important shortcomings, among which a systematic overestimation of vibrational transitions, worsening at higher quanta, and the fact that all nonfundamental bands have null intensities. A cost-effective way to improve this is offered by the second-order vibrational perturbation level of theory (VPT2). [37] As a reminder, the starting point is the harmonic Hamiltonian, where x i is the harmonic wavenumber, q i and p i are, respectively, the dimensionless normal coordinate and conjugate momentum associated to mode i, and N the total number of normal modes. The potential energy operatorV is expanded as a Taylor series up to the fourth order, @ 3V @q i @q j @q k q i q j q k 1 1 24 @ 4V @q i @q j @q k @q l q i q j q k q l (4) By comparison with the definition of the harmonic oscillator, the perturbative orders of the Hamiltonian operator are, @ 4V @q i @q j @q k @q l q i q j q k q l 5 1 24 The coupling between the vibrational and rotational wavefunc- q i p j q k p l B eq s is the rotational constant at equilibrium geometry and f ij;s the Coriolis constant coupling modes i and j along the rotation axis I s . In order for this formula to be valid, the molecule must be in Eckart orientation, in which the vibro-rotational couplings are minimized and the inertia tensor is diagonal.
The most common strategies to carry out the development and obtain the vibrational energies are based on the Rayleigh-Schr€ odinger perturbation and Van-Vleck contact transformation [107] theories. The result is an analytic formula, sufficient to compute the energy of any vibrational state j m i (in cm 21 ), where v m i is the number of quanta associated to mode i in state m, and e 0 is the zero-point vibrational energy (ZPVE), v is the matrix containing the anharmonic contributions, whose elements are defined as, Here, the problem of normal mode degeneracies has been ignored. A recent, more general formulation, supporting doubly degenerate modes (symmetric-and linear-top molecules) can be found in reference 51.
An alternative form of Equation 7, sometimes used in the literature, relies on the v 0 matrix, also noted G 0 . [45,108] Its major drawback, however, is the presence of resonances in v 0 , which disappear when recast into e 0 . [108] The use of v 0 is nevertheless interesting when dealing with transition states. [49] The implementation of Equations 7-10 is straightforward and the VPT2 model has been available in multiple computational chemistry programs and standalone codes. [18,22,43,109] Practical discussions on the calculations themselves will be deferred to the next section.
An important aspect of VPT2 calculations is the problem of resonances, which can be encountered in a large number of cases. The most critical situation comes from the so-called Fermi resonances (FRs, type I: x i % 2x j and type II: x i % x j 1x k , used in the following as the general form), in which case some denominators in Equations 9 and 10 can become very small or almost negligible, leading to an inaccurate definition of the anharmonic correction. This represents a failure of the standard VPT2 model, which must be treated accordingly. The most common strategy is to identify the resonant states and remove the corresponding terms from the calculation of v, thus e m . This procedure is far from being straightforward since the definition of the resonance is not unequivocal. For this reason, several approaches have been proposed in the literature [45,46,49,110,111] and only the ones used in this work will be described. The so-called deperturbed VPT2 (DVPT2) relies on a two-step procedure in which the states close in energy are first identified (low frequency difference) and the resonance magnitude is then estimated from the deviation of the VPT2 term from a model variational treatment. The latter step is commonly known as Martin's test. [110] The overall procedure can be summarized as, Step 1 : jx i 2ðx j 1x k Þj1x 122 x Step 2 : jh v11 i jĤ 0 jv11 j 11 k ij2 k ijk whereĤ 0 is the contact transformed hamiltonian, which is used here for consistency during the variational correction. [47] Here and in the following, harmonic vibrational states will be represented as vectors of To avoid the definition of thresholds, an alternative strategy is to replace all potentially resonant terms with nonresonant equivalents.
Such an approach was for instance proposed by Kuhler, Truhlar, and Isaacson as the degeneracy-corrected PT2 (DCPT2). [112] One important drawback of this approach is its inaccuracy with respect to the original VPT2 formulation far from resonance. However, this can be significantly reduced by introducing an ad hoc transition function to recover the VPT2 term in this case. This hybrid approach is referred to as HDCPT2 (hybrid DCPT2-VPT2). [49] An important difference between DCPT2 and DVPT2 in presence of resonance is that the anharmonic correction will become truncated in the latter case, since some terms The final vibrational energies are obtained after diagonalization of the full variational matrix. The overall procedure (DVPT2 1 variational correction) will be called generalized VPT2 (GVPT2  [113] (DDRs), need to be treated variationally, which is done concurrently with the correction to the FR. Similarly to the latter, a two-step procedure is employed to identify the resonances. First, the states close in energy are selected.
Then, within this ensemble, those for which the associated variational terms have a significant magnitude are defined as resonant. This can be summarized as follows, Step 1 : jDxj1D A-B x Step 2 : jh v 0 jĤ 0 jv 00 ij ! K A-B where "A-B", v 0 , and v 00 depend on the type of DDR, generally identified by the number of annihilated and created quanta (ex: 1-1, 2-2).
For the sake of completeness, another type of DDR, 1-3, needs to be taken into account when 3-quanta transitions are included in the variational correction. In the fundamentals region (below 4000 cm 21 ), those transitions generally contribute marginally to the band-shapes, [114] especially for spectroscopies like IR and Raman, and thus will be ignored in this work. Extension to 3-quanta transitions is nevertheless straightforward. [47,52] The variational terms have been recently reviewed by Rosnik and Polik, who fixed some errors present in previous works. [47,115] They are reported in the Supporting Information (Equations S1-S10), with the notation adapted to the present work. A potential issue of those variational terms is that they are not immune to FRs themselves. This problem has been already noted by several authors. [47,115,116] The most common strategy is to discard the resonant terms, either systematically by assuming all of them to be resonant, or after a preliminary identification step. Here, an alternative approach, derived from the HDCPT2 scheme is employed, to replace a priori all potentially resonant terms, assuming the following transformation to be accurate, Considering for instance one of the potentially resonant term in the correction associated to the 1-1 DDR (Equation S2 Supporting Information), we apply the following equivalency, leading to the corrected term, To get accurate band-shapes, the transition moments for the properties of interest has to be computed at the anharmonic level as well.
Here, only IR and Raman will be used, so that only the electric dipole (l) and the polarizable tensor (a) are needed. It is important to note that the latter is generally computed in the far-from-resonance regime in electronic structure calculations. For the resonant (or pre-resonant) regime, a different formalism needs to be followed, which will not be presented here. Discussions on the vibronic calculations of resonance Raman spectra can be found in references 73, 117 and 118, and references therein for instance. Both quantities considered here (l and a) are functions of the normal coordinates, so there is no need for a gen-eral formulation such as the one adopted in references 50 and 52. For this reason, the more compact notation, used in references 32 and 119, will be followed. At variance with vibrational energies, two states are involved in the definition of the transition moments (instead of one) so that different analytical formulas are required depending on the initial and final levels. Equations for transitions from the ground state to states with up to 2 quanta are reported in the Supporting Information (Equations S11-S12).
In addition to FRs, intensity calculations can also be directly impacted by DDRs as well. Since they are connected to the mechanical anharmonicity, the method adopted here is to use the identification procedure from the energy calculations to generate the set of resonant Finally, the Raman and IR band-shapes will be computed as, where ðx 0 Þ is the molar absorption coefficient and "dr=dX " the first differential cross section, with x 0 the wavenumber of the incident light. for small molecules [68,69,120] or for some specific degrees of freedom. [64] Considering the size of the systems of interest here, the harmonic approximation will be used to compute the transition energies and intensities. We will see later how some of the anharmonicity can be partially recovered in a cost-effective way.
The first issue to address is the definition of a common basis set to compute the integral h H i I;F . Here, the affine transformation proposed by Duschinsky to express one set of (mass-weighted) normal coordinates with respect to the other will be adopted, [121] Q I 5JQ F 1K where J is the Duschinsky matrix and K the shift vector, whose definitions depend on the type of reference coordinates. For the sake of simplicity, we will assume the latter to be the Cartesian coordinates so, where L x is the transformation matrix from mass-weighted Cartesian to normal coordinates, M is the diagonal matrix of atomic masses, and X eq I and X eq F are the equilibrium Cartesian coordinates of the initial and final states, respectively. Extension to internal coordinates is rather straightforward [74] and only requires modified definitions of J and K, as well as the derivatives of the electronic transition moments for the property of interest. For the sake of readability, the Cartesian coordinates will be used as reference in following. To ensure the generality of the discussion, J will be assumed nonorthogonal, which is the case in internal coordinates.
The definition of K requires the knowledge of the initial and final states' equilibrium geometries. This leads to another issue related to the harmonic approximation of the PESs. Indeed, the electronic transition can be characterized by a structural change, so that the two PESs will differ and their minima will be shifted. The latter can put a strain on the reliability of the harmonic representation, which worsens with the "distance" from the minimum. For this reason, two strategies can be adopted to compensate this problem. The first one, known as the adiabatic model, maximizes the overall description of each PES, which is computed about its actual minimum. This represents the best definition of the initial-and final-state equilibrium structures and frequencies (within the harmonic description), so that the highest accuracy for the band positions in the vibronic spectrum should be achieved. An alternative way is to emphasize the correctness of the preeminent bands in the overall band-shapes, thus the most intense transitions. This requires a better description of the PES about the vertical transition (highest overlap), so that this approach is often referred to as the vertical model. Within each strategy, two reasonings can be adopted. One is to emphasize accuracy, so that the force constants (hessian matrix) are computed in each electronic state. They will be referred to as Adiabatic Hessian (AH) and Vertical Hessian (VH) depending on the model.
Another one is to prioritize computational cost and feasibility. In this case, the two PESs are assumed to be equal, so that only one set of force constants needs to be calculated, generally the less expensive one, commonly in the ground state. They will be labeled based on the main quantities, which need to be known from the other state, namely Adiabatic Shift (AS), since the geometry optimization is still necessary, and Vertical Gradient (VG), since only the forces at the reference geometry are required (the latter is also known as the linear coupling model in the literature [122] ). The definitions of J and K differ based on the model. [59,74] The most important aspects are that J is the identity matrix in the approximated models (AS and VG), K is the same for the adiabatic models (AS and AH), and K is obtained from an extrapolation of an harmonic PES for vertical models (VH and VG). By comparing the definitions of K in the adiabatic and vertical models, it is possible to obtain an estimate of the unknown equilibrium geometry for the vertical models, which can provide a visual aid in assessing the reliability of those models.
The definition of a transformation between the two sets of initialand final-state normal coordinates is necessary but not sufficient to compute the transition moment h H i I;F . Indeed, the electronic transition moment (H e 5h / I jĤj/ F i) has a dependence on the nuclear coordinates, which is not known analytically. As a consequence, an approximate form is adopted. The most common way was proposed by Franck and formalized by Condon. [123][124][125][126] In practice, it is assumed that the electronic transition moment does not depend on the nuclear configuration and thus is constant. This approximation was later refined by Herzberg and Teller, who proposed to include a linear variation with respect to the normal coordinates. [127] Those models can be generalized as a Taylor expansion of the electronic transition moment about one equilibrium geometry, where the zeroth-order term corresponds to the FC approximation and the first order to the HT one. The choice of the reference state for the expansion depends on the model used to describe the transition. It is generally the initial one for vertical models and the upper (excited) one for adiabatic models, mainly for practical reasons. Note that in the following, the development up to the first order will be labeled FCHT, which is sometimes simply defined as the HT approximation in the literature.
Thus, the definition of the model for the description of the electronic transition (AH, AS, VH, and VG) and the approximation of the transition property (FC, FCHT) will pave the way to the calculation of the transition moments. The band-shape is then obtained by computing the intensity of each transition separately. This corresponds to the sum-over-states or TI approach. An alternative way to get the bandshape is through an implicit inclusion of all possible transitions, within the so-called path integral or TD formalism. To understand better the differences between those two methods, let us write the formula of the total sum-over-states intensity of the one-photon absorption (OPA) or emission (OPE) spectrum for a given absorption or emission wavenumber (x 0 ), respectively, where a is a constant factor, which depends on the type of transition, where the single and double overbars were used to represent the initial and final states, respectively, to avoid confusion with other subscripts. Hence, the original integral h v I jl e jv F i can be written as a combination of integrals between vibrational states of the initial and final electronic states, also called FC integrals. Several approaches have been proposed in the literature, based either on analytic [128][129][130][131][132][133] or recursion [134][135][136][137] formulas. The latter have proven to be more efficient, especially when dealing with medium-to-large systems since they only require the implementation of a few equations, generally the integral between the vibrational ground states (h 0 I j0 F i) and one or more recursion relations. The approach used here is the one proposed by Ruhoff based on the original work by Sharp and Rosenstock. [128,136] The formulas needed to compute recursively any FC integrals are listed in the Supporting Information (Equations S13-S15).
Let us remark that for the special case where the PESs are equal (approximation used in AS and VG), significant simplifications can be made to the recursion formulas. In particular, if temperature effects are negligible, so that only the vibrational ground state of the initial state is At the FC level, the total intensity is then, factors [138] gives an estimate of the intensity of any overtone when the coupling between them is negligible.
Note that, in theory, the summations over the initial and final states are infinite. At low temperature, the ground state can be assumed to be the only one populated, so that the summation over I can be discarded.
However, some screening is still necessary to limit the number of transitions to a practical level. To this end, various strategies have been pro- However, the performance of the sum-over-states approach is generally low when temperature effects become important and the number of initial states to take into account increases strongly.
This problem can be overcome with the TD approach. The principle is to recast Equation 18 within the time domain by replacing the Dirac function by its Fourier transform, and the Boltzmann population with where Z is the canonical vibrational partition function. The total intensity can then be rewritten as follows, dt Tr l e e 2sF HF l e e 2sIHI À Á e i DE2x0 where Tr l e e 2sF HI l e e 2sIHI À Á is referred to as the electronic transition An important assumption done in the TI and TD formulations is the fact that the PESs are harmonic. As discussed before, improvement of the band-shape by mean of a full anharmonic treatment is not conceivable due to the steep increase of the computational cost. However, an alternative way is to focus only on the band positions, and more specifically on the fundamental energies. Those energies can be computed at the VPT2 level as shown before. While the harmonic approximation is still assumed valid, this provides a systematic improvement to the transition energies, even for higher quanta, with the computational cost of a standard VPT2 calculation. One prerequisite remains the availability of analytic harmonic force constants, which may not be the case for all methods, for instance most implementations of TD-density functional theory (DFT). Another hurdle may be that the cost is too high to build the anharmonic force constants for one of the PESs. In such cases, where the VPT2 fundamental wavenumbers can only be known for one state, it is possible to use a simple scheme based on the Duschinsky matrix to derive scaling factors to be applied to each harmonic wavenumber of the other state. [62] More formally and assuming that the final state's anharmonic wavenumbers are unknown, then the extrapolation formula is, where m i is the anharmonic fundamental wavenumber (m i 5e 1i 2e 0 ), and J is assumed orthogonal. Extension to the nonorthogonal case can be found in reference 74.

| One-dimensional large amplitude motion
The vibronic models described above assume the molecule to be sufficiently rigid and the vibrations not associated to any LAM. However, this is not always the case and the flexibility of the system needs to be properly taken into account. In some cases, the main structural change can be associated to a single coordinate, and therefore a block-based approach, where the LAM is treated at the anharmonic level, whereas the harmonic approximation is used for the other modes, can be devised. However, in order for this block-based approach to be effective, the coupling between the LAM and the other vibrational degrees of freedom must be negligible, and this condition is met only rarely using Cartesian-based normal modes as the reference coordinate set.
At variance, the coupling can be strongly reduced using curvilinear, internal coordinates, which make them better suited to treat LAMs.
The extension of the theoretical framework described above to support internal coordinates in vibronic spectroscopy has been presented recently [74] for both TI and TD-based algorithms. Without going into details, this generalization requires the expression of the quantities needed to simulate the vibronic band-shape in terms of a general set of curvilinear, internal coordinates. The harmonic frequencies do not depend on the coordinate system, and therefore, at the FC level, only the Duschinsky transformation needs to be generalized. This can be done following, for example, the approach proposed by Reimers. [72] Inclusion of HT effects requires also the expression of the derivatives of the transition dipole moment in terms of internal coordinates, as discussed, for example, in reference 74. To apply this framework to largesize systems, a fully automatized procedure to select the nonredundant set of internal coordinates is needed. In our discussion, delocalized internal coordinates (DICs) will be used, [142] which are well-suited for the study of large-size systems since they can be built using only the molecular topology. Moreover, they are more reliable than the widely used Z-matrix coordinates when dealing with nonsimple deformations associated to the electronic transition. [73] As already remarked above, the use of internal coordinates is crucial to reduce the coupling between the LAM and the other vibrations.
A good way to check if the coupling is sufficiently low is to inspect the Duschinsky matrix computed for the full system. This procedure can be straightforwardly automatized and has been described in references 74 and 143. The principle is to build a consistent set of normal modes (same number in the initial and final states) based on their overlap given by the squared elements of the Duschinsky matrix (note that the rows and columns need to be normalized in internal coordinates). The vibration corresponding to the deformation (largest element of the shift vector) is first selected and its projection on each normal coordinate of the other state is calculated. If the highest overlap is above a given threshold, the mode is considered uncoupled from the rest of the system (to be valid, the limit should be set high, for instance 0.9). Once this condition is met, it is possible to factorize the FC integrals as follows: The factorization given in Equation 21 can be used to develop the hybrid scheme outlined above. Indeed, the one-dimensional FC integrals associated to the LAM coordinate (h v LAM j v LAM i) can be computed at the anharmonic level, while the other integrals are treated at Let us present now the anharmonic model used to describe the LAM. As a first step, the LAM must be expressed as a linear combination of internal coordinates. The appropriate linear combination can be determined, for example, using a potential energy distribution analysis. [144] Then, the PES along the LAM coordinate is computed Miller [145] in terms of the gradient of the PES along the MEP. Within this framework, each Cartesian geometry is expressed in function of s and the displacements of the other "N 2 1" modes. If the deformation along those modes is negligible, the distance between two geometries (x 1 and x 2 , corresponding to two values of the IRC parameter s 1 and s 2 ) can be approximated as the corresponding distance in mass-weighted Cartesian coordinates, as expressed in the following: The previous equations can be used to build the Hamiltonian matrix Further extension to higher number of flexible degrees of freedom has been proposed only for small systems. [147] In any case, the extension of the model presented in this work to multidimensional problems is not straightforward, and therefore will not be discussed here.
The c matrix obtained by diagonalization of H s can be used to express the vibrational eigenfunctions as linear combinations of the DVR basis set as follows, where j v LAM i is the wavefunction for the LAM associated to the vibrational quantum number v. In vibronic spectroscopy, the secular equation where the orthogonality of the DVR basis sets has been implicitly employed.

| Combining One-dimensional large amplitude motion and vibronic calculations
The procedure described above can be used to compute the FC factors involving the large-amplitude motion. Indeed, within the harmonic approximation and the sum-over-states formalism presented in Equations 18 and 21 can be used to combine the LAM and the FC factors involving the other "N 2 1" modes and thus obtain the overall vibronic spectrum.
where the I and F refer to combinations of the "LAM" and "N 2 1" initial and final states, respectively. The FC approximation has been assumed, so that the transition dipole moment is a constant factor. It should be noted that the inclusion of HT effects within this framework is not straightforward and will not be studied here. Indeed, the first-order terms of the Tay where: The relation given in Equation 30 shows that the final vibronic spectrum is not simply the sum of the two vibronic spectra, but is given by the convolution of the spectra of the LAM motion (I LAM ) with the vibronic spectrum of the "N 2 1"-dimensional system (I N21 ).
The development version of the suite of quantum chemical programs Gaussian was used for the vibrational and vibronic simulations. [152] To choose the most suitable combination of exchange-correlation functional and basis set, a preliminary group of the currently most popular functionals were selected (B3LYP, [153] CAM-B3LYP, [154,155] PBE1PBE (PBE0), [156] LC-xPBE, [157] M06-2X, [158,159] x-B97X, [160] and B2PLYP [161][162][163] ) and used to compute the vibrational and vibronic spectra, with thiophene as the reference molecule. The double-f basis set SNSD, [17,164] which has been extensively validated for several spectroscopic properties in the ground [17,165] and excited electronic states [67,166] has been used, except for the double-hybrid B2PLYP functional, since the latter has been shown to require a larger basis set. [26,163,167,168] In this case, a modified version of maug-cc-pVTZ, [169] where d functions on hydrogens have been removed, was employed.
Since the dispersion correction is expected to play an important role in the conformational stability of bithiophene [170] Grimme's type dispersion corrections, [171] mainly its D3 formulation [172] in conjunction with the Becke-Johnson (BJ) damping [173] were used in all DFT calculations (except for M06-2X-D3 and xB97XD). This preliminary assumption has been confirmed by comparing the rotational profiles of bithiophene computed employing dispersion-corrected DFT-D approaches with standard DFT models.
Whenever suitable, bulk solvent effects have been taken into account by means of the polarizable continuum model [174] within its integral equation formalism (IEF-PCM), [175] with the default set of

| VPT2 calculations
The anharmonic data needed for the VPT2 calculations were obtained through numerical differentiation along the normal coordinates The overlap was then checked based on the squared matrix elements J ik 2 . A threshold of 90% for each coordinate was required to consider the two sets consistent.
Raman intensities were computed from frequency-dependent polarizabilities, [177] using an incident frequency matching experiment (514.5 nm for the Raman spectrum [87] ). The validity of the far-fromresonance regime was confirmed by comparison of x 0 with the energy of the first excited electronic state. Temperature effects were only included within the TD framework, at 298 K.

| One-dimensional scan
The one-dimensional potential energy profiles along the SACAC 0 AS 0 dihedral angle h of 2,2 0 -bithiophene (see Figure 1 describing the rotation between the two thiophene rings of 2T have been selected to assess the reliability of the different DFT approaches. The equilibrium structures have been compared with the so-called semiexperimental (SE) equilibrium geometry (r eq SE ), [179][180][181][182] accurate within 0.001 Å and 0.1 degrees for bond lengths and angles, respectively, taken from a database providing geometries for a large number of small and medium-size molecules. [181,182] Moreover, the accuracy of the theoretical vibrationally averaged structures have been evaluated by comparing the theoretical and experimentally measured rotational constants (B 0 ). [183] The computed values have been obtained by correcting the equilibrium rotational constants (B eq ) with the vibrational contributions DB vib obtained from VPT2 computations. [37,39] The mean absolute errors (MAEs) and maximum absolute deviations (jMAXj) between each DFT functional and the reference data are reported in The good performance of double-hybrid functionals have been already shown for other semirigid systems. [114,181,182,184,185] We note that for thiophene, the dispersion correction does not introduce structural changes; however, its inclusion shall be more important for flexible systems, [170,186,187] as will be shown for bithiophene. Very accurate structures have been obtained with LC-xPBE/SNSD too, and all other DFT functionals deliver CAC and CAH bond lengths accurate within 0.01 Å, while larger discrepancies are observed for the SAC bonds, which are overestimated by 0.013 up to 0.028 Å. The trends observed for the SAC bond lengths are consistent with the results for B 0 , with a better accuracy for CAM-B3LYP(-D3BJ) compared to B3LYP(-D3BJ).
Comparison between theoretical anharmonic and experimental [87] gas phase fundamental wavenumbers, reported in Figure 3, shows clearly that the best results, with MAE and jMAXj as low as 5.  with D3BJ. Based on those results and taking into consideration our own experience, we can conclude that the D3 and D3BJ dispersion corrections can be recommended for the anharmonic vibrational computations, leading to systematic improvements even for semirigid molecules, which becomes more significant for instance in case of hydrogen-bonded [188,189] or stacked systems. [190] Interestingly, the worst performance, with mean and maximum errors of 28 cm 21 and 80 cm 21 , is observed with LC-xPBE-D3BJ/SNSD, which, however, gave very good geometry structures. This demonstrates that the good accuracy achieved at the bottom of the well does not guarantee a reliable prediction of the PES curvature, so that both aspects need to be taken into account in spectroscopic studies. Another important aspect to consider is the relative performance of harmonic frequencies and anharmonic corrections taken separately, in particular with reference to hybrid models where the harmonic part can be computed at a higher level of theory, for instance coupled cluster (CC) or B2PLYP in conjunction with basis sets of at least triple-f quality. [163,184] In this regard, we will use the harmonic wavenumbers (x) and anharmonic corrections (Danh) computed at the B2PLYP-D3BJ/maug-cc-pVTZ (B2D3) level as the reference (Figure 3). For CAM-B3LYP-(D3BJ) and PBE0-D3BJ, the lower accuracy for the fundamental transitions can be fully related to the harmonic part, which is also the case for LC-xPBE-D3BJ. At variance, for M06-2X-D3 and xB97XD, the large discrepancies with experiment are due to large errors in both harmonic frequencies and anharmonic corrections, further highlighting that these functionals are less suitable for vibrational computations. [17,163] The good accuracy of hybrid schemes combining B2PLYP-D3BJ/maug-cc-pVTZ harmonic contributions with B3LYP-D3BJ/SNSD or CAM-B3LYP-D3BJ/SNSD anharmonic corrections is confirmed by comparison with experimental data in Figure 3, with both models performing on par with the full B2PLYP-D3BJ/maug-cc-pVTZ anharmonic results, but at a significantly lower computational cost.
For the computation of vibrationally resolved electronic spectra we can distinguish two important factors whose combination governs the overall accuracy of the results: (i) the excitation energy and electronic transition moments and (ii) the shape of the initial and final PESs.
The first group determines the position of the absolute energy range and the intensity of the electronic transition, while the second group influences the magnitude of the overlap integrals and spectral bandshape. It can be noted that the two sets of data can be computed with different methods if needed, leading to hybrid ESQM/ESQM' models where the vertical properties, computed at one, usually higher level of theory [26] are combined with PESs obtained from less expensive computations. In the framework of DFT/TD-DFT, it can be also convenient to combine two different functionals which perform best for each of the two aspects. The accuracy of DFT models for the prediction of excitation energies is often performed by comparing the VE energies with experimentally observed k max , even though the latter depends on the band-shape, made of vibronic transitions. For this reason, a more reliable validation can only be achieved by simulating the band-shape from the convolution of vibronic transitions, including, where necessary, solvent and temperature effects. [191] Figure 4 compares the TD AHjFC spectra computed using different DFT functionals, at 298K, with the experimental photoionization spectra of thiophene. [89] The good match between the theoretical envelope and experimental band at about 9 eV suggests that indeed the latter is related to the 2 A 2 X 1 A 1 transition as postulated in earlier works. [88,89] [89] and theoretical photoelectron spectra of thiophene. TD AHjFC calculations were done for the 2 A 2 X 1 A 1 transition at T 5 298 K with different DFT functionals, using Lorentzian distribution functions with HWHMs of 1000 cm 21 for the broadening. The vertical bar corresponds to the vertical transition energy at the CAM-B3LYP-D3BJ level.  be concluded that all DFT functionals lead to the reliable prediction of the 2 A 2 X 1 A 1 transition for thiophene. However, this is not always the case and, for instance, the correct description of Rydberg and charge transfer states with TD-DFT is more challenging. [192] In this respect, CAM-B3LYP [193,194] presents a significant improvement over the standard B3LYP. intensity band at about 2600 cm 21 , which is due to "hot" transitions.
It can be conluded that all DFT models predict reliably the band-shape of the medium-resolution spectra, while for the high-resolution spectra, where the accuracy of the band position becomes more critical, the same observations as for the vibrational energies would apply.
The spectral properties of 2,2 0 -bithiophene are strongly connected of the potential energy profile of the rotation between the two thiophene rings. The potential energy scan along the SACAC 0 AS 0 dihedral angle is shown in Figure 6, with particular emphasis on the region close to the two global minima of trans-2T. Following the results obtained for thiophene and previous studies, [114,170] B2PLYP-D3BJ will be used as reference. All methods predict trans-2T to be more stable than To conclude this validation procedure, and following this preliminary study on the potential energy along the SACAC 0 AS 0 dihedral angle, the S 1 S 0 OPA spectrum of bithiophene was computed and compared to experiment. [195] Due to the system size, full vibronic calculations at the B2PLYP level are unfeasible and only the most successful functionals up to now, that is B3LYP, CAM-B3LYP, and their dispersion corrected D3BJ versions were considered. The TD AHjFC model has been applied to all vibronic simulations to obtain a fully converged band-shape, and DICs [74] have been used to provide a correct description of the large-amplitude distortion along the SACAC 0 AS 0 dihedral angle. Furthermore, on the base of the previous analysis, all computations have been performed for the minima associated to the trans conformer, which is the most stable. To reproduce the experimental conditions, [195] solvent effects (n-hexane) have been included by mean of PCM and the broadening was simulated with Gaussian distribution functions with HWHMs of 100 cm 21 . As shown in Figure 7, In conclusion, considering every aspect influencing the spectra calculations, it is rather clear that, as could be expected, the doublehybrid B2PLYP-D3BJ with maug-cc-pVTZ outperforms the other considered models. However, that is also the most expensive approach and, due to the MP2 contributions, it has larger basis set requirements too, so its applicability might be limited to the equilibrium structure and harmonic calculation of the PES for larger molecular systems.
More importantly, several properties, which are widely available for

| Vibrational spectra
The IR and Raman spectra computed with different types of approximations are compared to their experimental counterparts recorded in gas phase [85,87,196] in Figures 8 and 9, respectively. The most common approach, which is also the starting point for more accurate techniques, is the harmonic approximation ("Harm" in the figures). At this level, all bands are shifted, with MAE of 43 cm 21 , the differences being more pronounced at higher energies resulting in jMAXj over 140 cm 21 . A common refinement to improve the overall agreement is to apply a FIG URE 8 Experimental [196] and theoretical IR spectra of thiophene. Spectra from top to bottom correspond to the raw harmonic spectrum (Harm), the harmonic spectrum with scaled wavenumbers (Scaled), the anharmonic spectrum with harmonic intensity (GVPT2 freq.), the full anharmonic band-shape (GVPT2), and the experimental spectrum (Exp). The theoretical spectral lineshapes have been convoluted by means of Lorentzian distribution functions with HWHMs of 2 cm 21 . All spectra have been shifted vertically for better visual comparison.
FIG URE 9 Experimental [87] and theoretical Raman spectra of thiophene. Spectra from top to bottom correspond to the raw harmonic spectrum (Harm), the harmonic spectrum with scaled wavenumbers (Scaled), the anharmonic spectrum with harmonic intensity (GVPT2 freq.), the full anharmonic band-shape (GVPT2), and the experimental spectrum (Exp). The theoretical spectral lineshapes have been convoluted by means of Lorentzian distribution functions with HWHM of 2 cm 21 . All spectra have been shifted vertically for better visual comparison. is that, depending on the complexity of the system, the correction can be inconsistent, and thus, the reliability of a band assignment based on the scaled energies may be questionable. A final hurdle to their application is that scaling factors need to be derived for each computational model (DFT functional 1 basis set), [197] which means an "ad hoc" parametrization whenever a new ESQM approach is to be chosen. [198][199][200] A more consistent approach requires a proper inclusion of anharmonic effects, here by mean of GVPT2 computations ("GVPT2 freq"). As SNSD. [26,114,184,201] Up to now, focus has been given on improving the vibrational states' energies, hence the band positions. To get band-shapes comparable to experiment, it is necessary to treat correctly the transition intensities as well. As a matter of fact, all nonfundamental transitions have null intensity at the harmonic level, so several bands are missing FIG URE 10 Experimental [196] and theoretical GVPT2 IR spectra of thiophene in the 1000-2000 cm 21 range, along with the assignment of the most intense fundamental (green) and nonfundamental (blue) transitions. "n m " represents the final vibrational state with n quanta associated to mode m. The theoretical line-shape has been convoluted by means of Lorentzian distribution functions with HWHM of 2 cm 21 . favored for the assignment. In fact, it has been already shown that some corrections to previous assignments might be necessary based on the fully anharmonic computations. [202] The detailed list of all experimentally observed IR and Raman transition of thiophene compared to theoretical results is reported in Table 1, in most cases confirming previous experimental assignments. Some corrections can be proposed in the spectral ranges where several overtones or combination bands match closely the observed experimental energy. Indeed, in such situations it is very difficult to perform the assignment based solely on the TABLE 1 GVPT2 fully anharmonic wavenumbers (m, in cm 21 ), IR intensities (IR, in kmÁmol 21 ), and Raman activity (RA, in Å 6 ) compared to gas phase experimental data [87]    match between the estimated and observed band positions (even with selection rules based on symmetry when possible). For instance (see Table 2 21 . In all those cases, once the intensity of the nonfundamental transitions is known, the assignment is straightforward. As an example, the band at 2190 cm 21 observed in the IR spectra can be assigned to 1 5 1 19 with significant IR intensity, and the band at 2127 cm 21 observed in Raman spectra to 1 15 1 17 , which shows a strong Raman intensity. Such a procedure has been already applied for the detailed assignment of spectra, also in the near IR region, where only nonfundamental transitions are present. [52,114,203] On the whole, we can conclude that for thiophene, GVPT2 com-

| Electronic spectra
To illustrate some of the difficulties, which can be faced when simulating electronic spectra, three types of transitions have been considered, each one accompanied by structural changes of increasing magnitude, depicted in Figure 12, namely the photoionization 2 A 2 X 1 A 1 and valence electronic 1 A 1 X 1 A 1 transitions of thiophene and the S 1 S 0 transition of 2,2 0 -bithiophene.

| Small structural changes: 2
The 2 A 2 X 1 A 1 photoelectron spectrum of thiophene corresponds to the smallest structural changes on electronic excitation (ionization), mainly related to the 60.55 Å elongation of the C2AC3/C4AC5 bonds and shortening of the C3AC4 one (see Figure 1), resulting in an inversion of the bond order, in agreement with previous theoretical investigations. [88] Nevertheless, both ground-and ionic-state equilibrium structures are planar and overlap strongly as shown in Figure 12, representing a suitable example to discuss the performance of the different theoretical approaches available within the FC principle. The spectral  2 GVPT2 fully anharmonic wavenumbers (m, in cm 21 ), IR intensities (IR, in kmÁmol 21 ), and Raman activity (RA, in Å 6 ) for all transitions up to the 2 quanta in the 2100-2200 cm 21 range, compared to gas phase experimental data [87]  band-shapes obtained with different models for the description of the initial and final PESs (computed at the CAM-B3LYP-D3BJ/SNSD level) are presented in Figure 13. The vibronic transitions were all computed within the FC approximation using the TI formulation, with a total convergence of the spectrum of nearly 100%. Lorentzian broadening functions were used for the convolution with larger (500 cm 21 ) and smaller (75 cm 21 ) HWHMs, to simulate the low and medium-resolution spectra. The low-resolution spectra are presented in an absolute energy scale, highlighting the shifts between the bands. The differences are related to the position of the 0-0 transition, which also include a contribution from the ZPVE of the initial and final states. AH represents the most accurate description in this case. Two factors can influence the value of the 0-0 transition energy. In vertical models (VH and VG), the PES minimum is extrapolated by assuming a parabolic form. Moreover, for VH, the frequencies are computed outside the minimum, which may induce some further approximations. A second aspect is that the vibrational contribution is null when the PESs are assumed equal (AS and VG). In the present case, these differences are rather small, and all models predict the bands to be blue-shifted (by 0.03-0.06 eV; 230-500 cm 21 ) with respect to the reference AH model. The most shifted band-shape is obtained for VH, which can be related to the extrapolation of the PES minimum. VG performs better, most likely because of some error compensation between the extrapolation and the approximation on the excited-state PES, and has a shift similar to AS, which differ from AH only by the lack of contribution from the ZPVEs to the 0-0 transition energy. It is interesting to note that it is not possible to assess a priori if the neglect of the difference between the ZPVEs, as well as the PES extrapolation, will lead to a red-or blueshift of the 0-0 transition as it depends on the system and changes in the structure or the vibrational frequencies on excitation. However, in the present case, all models lead to very similar low-resolution spectra.
In particular, all predict the same band position and envelope, in good agreement with the experimental band at about 9 eV, as already pointed out.
The medium-resolution spectra are reported with energies relative to the 0-0 transition, and shifted on the vertical scale, to facilitate the qualitative comparison between band-shapes obtained with the different models used before. To provide further insight to the PES changes between the two electronic states, the HR factors, [138] which are closely related to the shift vector as shown in the theoretical section, and the Duschinsky matrices (J) are reported in Figures 14 and 15, respectively. The largest HR factors are associated to the same normal modes, indicating that all approaches reproduce the same trends for the geometrical changes on ionization. The main difference between the vertical and adiabatic models is related to mode 11 (t 6 in Table 3), which corresponds to an in-plane ring deformation. However, more important effects, related to the changes in the normal modes, should 14 Huang-Rhys factors for the 2 A 2 X 1 A 1 transition of thiophene, simulated within the TI framework based on Cartesian coordinates with the FC approximation and the VG, VH, AS, and AH models at the CAM-B3LYP-D3BJ/SNSD level. Internal coordinates (AH Int) and the special case "AH J 5 I" are also considered. other models a single, more intense band. These differences indicate clearly that the mode mixing can play an important role even for systems not undergoing significant structural changes on excitation, [204] and should not be readily discarded if accurate band-shape or high resolution spectra are required. Conversely, for low-resolution spectra, the VG model represents a cost-effective way to account realistically for the band envelopes. It should be noted that, for spectra encompassing several electronic states, the band envelopes can vary strongly between transitions, with the overall band-shape made of their sum, which can lead to destructive effects in chiroptical spectroscopies, so the relative band-widths and intensities need to be correctly taken into account for a comparison with experiment. [28,[205][206][207] Table 3. Before analyzing the theoretical and experimental spectra presented in Figure 16, it should be noted that, TABLE 3 Harmonic and anharmonic fundamental wavenumbers (in cm 21 ) in the ground (X 1 A 1 ) and excited ( 1 A 1 ) electronic states of thiophene It is interesting to note that spectra computed with Cartesian and DICs are essentially identical at lower energies, which also correspond to smaller displacements from the equilibrium structures, and diverge at higher energies, so associated to regions more distant from the PESs minima.
The good agreement with experimental findings, which, to the best of our knowledge, also represents the first direct comparison between simulated vibrationally resolved spectra and high-resolution measurements on this molecule, [89] allows us to suggest some alternative interpretations. First of all, the absorption band at 5.156 eV was attributed until now to the origin of the 1 A 1 X 1 A 1 transition, [209] with the weak bands below this energy tentatively associated to the transition to the 1 B 2 state. Our simulation predicts that the intensity of the 0-0 transition is very low, and corresponds to a band at the limits of the experimental measurements (at about 4.9 eV). Thus, no other, lower-lying excited electronic state is needed to describe the weak bands below 5.1 eV, which are in fact related to the 1 A 1 X 1 A 1 transition. These results are in agreement with the most accurate theoretical predictions showing that 1 A 1 is the lowest excited electronic state, [89] as well as with suggestions that the vibronic structure of the transition to the 1 B 2 state might not be observable due to its unbound character. [210] In a similar way, a previous study based on a fulldimensional AHjFC model allowed us to correct an assignment of the low-intensity 0-0 transition for the phenyl radical. [204] Finally, we note  Table 3). These results confirm further that the direct simulation of vibrationally resolved electronic spectra is a valuable help for the proper analysis of experimental results.

| Large amplitude changes:
The preliminary study on 2T for the definition of the most suitable ESQM models can now be used to analyze more in detail the vibronic properties of this system. To this end, we will focus on the OPA and OPE spectra for transitions involving the S 0 and S 1 states. As shown in Figure 12, trans-2T is planar in the S 1 state, and nonplanar in S 0 , which corresponds to values of the dihedral angle S1AC2AC 0 2 AS 0 1 of 180 and $155 , respectively. Therefore, the electronic excitation is accompanied by a large-amplitude distortion and, similarly to the 1 A 1 X 1 A 1 spectrum of thiophene, the choice of the coordinate system used in  Table 3). the simulation is expected to be relevant to reproduce correctly the vibronic spectrum. [74,211] The theoretical OPA spectra, computed at the TD AHjFC level using Cartesian and DICs, are reported in Figure 17, together with the experimental spectrum, taken from reference 195. The TD approach has been used in all cases since, for flexible systems such as bithio-phene, the convergence of the prescreening scheme used in TI computations can be slow, especially when mode-mixing effects are relevant, such as when cartesian coordinates are employed. Conversely, the TD approach ensures a full convergence of the band-shape also in those cases. As expected, the spectra obtained with the two kinds of coordinate systems are very different, the shape being significantly broader  with Cartesian coordinates than with DICs. Another important aspect is that, while the spectra have been normalized to facilitate the comparison of the band-shape, the intensity is also halved with Cartesian coordinates. The spectrum is computed at the FC level, therefore the difference between Cartesian and internal coordinates is directly con-  [74] for 2T, the main source of inaccuracy is the definition of J. The impact of mode mixing is highlighted in Figure 17, where the spectra computed at the ASjFC level, therefore neglecting the mode-mixing effects, and at the AHjFC level are compared. The AS spectra are significantly less broadened than the AH ones, even though the same Gaussian broadening functions with HWHMs of 100 cm 21 have been used. Therefore, mode-mixing effects are relevant, since they enhance the intensity of a large number of vibronic transitions, at the expense of the intensity maximum, which is lowered. Anyway, even at the AS level, the agreement between the theoretical spectrum based on Cartesian coordinates and the experimental one [195] is still unsatisfactory since, for example, the intensities of the lower energy transitions are overestimated. Conversely, the agreement is better both at the AH and AS levels if DICs are used. Following our previous works, the couplings effects between the final-states modes, which contribute to the intensity of the combination bands, has been evaluated for each kind of coordinate system from the magnitude of the elements of the Sharp and Rosenstock C matrix. [74] The graphical representation of C for 2T in Cartesian and DICs, reported in Figure 18, confirms that, if the latter are used, the torsional mode (corresponding to mode 1 in the plot of To improve further the accuracy of the theoretical results, the DVR-based anharmonic approach presented in the theoretical section has been used to describe the torsional motion along the dihedral angle h between the two rings. Among the various high-resolution spectra available in the literature, [86,195] the laser-induced fluorescence (LIF) spectrum of 2T [86] will be used as a reference. In addition to our reference ESQM model, CAM-B3LYP-D3BJ, which is expected to provide the most accurate results, the low-barrier PES obtained from the B3LYP computations (see Figure 6) has been selected to show some effects of the barrier height on the vibrational energies and the overall spectroscopic properties. The PES of the S 0 electronic state obtained through a relaxed scan along the angle h, is reported at both levels in . [86] To compute the FC factors for the LAM, the PES along the torsional angle h was also computed at the TD-CAM-B3LYP-D3BJ/SNSD level for the S 1 state, with the results reported in Figure 20. The shape of the S 1 PES is significantly different from S 0 due to the presence of a single planar equilibrium structure. Therefore, the vibrational levels are nearly equally spaced (the computed spacing between the first two vibrational levels is 98 cm 21 ), with a spacing larger than for S 0 , since the PES of S 1 is significantly steeper.
The OPE spectra along the torsional motion are compared to the experimental LIF measurements [86] in Figure 21 tions, but the computed ones tend to be too low at higher quanta.
Except for the lowest energy transition, the relative intensities are also qualitatively comparable to experiment, reproducing the trend in the FIG URE 21 Theoretical S 1 ! S 0 OPE spectra trans-2,2 0 -bithiophene computed at the CAM-B3LYP-D3BJ (green line) level from the ground (upper panel) and the first excited (lower panel) vibrational levels of the S 1 electronic state. The experimental LIF spectrum, taken from reference 86, is also reported.

FIG URE 22
Theoretical S 1 ! S 0 OPE spectra trans-2,2 0 -bithiophene computed at the TI AHjFC level of theory using delocalized internal (green lines) and Cartesian coordinates (red lines) for the reduced-dimensionality system, where the torsional mode is neglected (dashed lines), and for the hybrid DVR-harmonic scheme (solid line). The experimental spectrum has been taken from reference 90.  Figure 19). In fact, as shown in the lower panel of Figure 21, the reproduction of the band position is in this case poor, and does not allow to assign univocally all the bands. Those results show that the DVR-based approach, when coupled with an accurate description of the PES, provides a significantly better characterization of large-amplitude motions. However, this model may need to be refined on several aspects, such as the support of multidimensional problems or the inclusion of coupling terms between the LAM and the other modes, [81] to further increase its accuracy.
The results shown above have been used to simulate the full OPE spectrum of 2T using the hybrid DVR-harmonic approach presented in the theoretical section, and the results are reported in Figure 22. The vibronic spectra have been computed at the TI AHjFC level for the "N 2 1" harmonic modes, explicitly excluding the mode associated to the LAM, which has been treated at the DVR level. The vibronic spectrum for the "N 2 1" modes has been convoluted with the progressions due to the LAM to obtain the full spectrum. As remarked in the previous section, the TD approach would be better suited to treat flexible systems such as 2T. However, to convolute the "N 2 1" spectrum with the LAM-specific one, the intensity of the individual vibronic transitions needs to be known, and therefore the TI approach must be employed.
In this case, the convergence is high both in cartesian coordinates (85%) and in DICs (99%), and this ensures the reliability of the TI results. As shown in Figure 22, the agreement is satisfactory for the spectrum computed using DICs, whereas it is poor with Cartesian coordinates. Those results show that, if the LAM is strongly coupled to the other "N 2 1" modes, such as when Cartesian coordinates are used, the accuracy is unsatisfactory independently of the method used to describe the LAM. A more accurate reproduction of the experimental spectrum is obtained using DICs, since in this case the coupling between the LAM and the other modes is negligible, increasing the accuracy of the hybrid DVR-harmonic approach. The spectrum obtained using DICs has been compared also to the one obtained including only the "N 2 1" harmonic block, and neglecting the vibronic transitions involving the LAM. As shown in Figure 22, the two spectra are nearly superimposable. As a matter of fact, to reproduce the experimental broadening, Gaussian functions with HWHMs of 200 cm 21 have been used. This value is much larger than the separation between the vibrational levels of the LAMs, and therefore the convolution of the vibronic transitions involving the LAM is equivalent to an overall broadening of the vibronic spectrum obtained for the "N 2 1" harmonic block. Such a conclusion is however only true for low-resolution spectra, since, when high-resolution spectra are simulated (such as for the LIF spectrum described above), an explicit inclusion of the vibronic transitions of the LAM becomes mandatory.

| Extension to new cases studies: Some illustrative examples
The protocol presented above can be applied to a wide range of molecular systems and spectroscopies. Here, we will illustrate some additional kinds of studies, regarding chiral spectroscopies and rotor-like molecules, which have one or more low-energy torsional modes.
The computational protocol to follow for electronic or vibrational chiral spectroscopies is very similar to what has been described above.
An important difference, however, is that they depend on multiple properties. In addition to requiring some specific theoretical and algorithmic developments, mentioned before (see references 16 29  n ECD spectrum of 1A conformer of (Z)-8-Methoxy-4-Cyclooctenone in the 340-240 nm range, compared to experiment in hexane solution taken from reference 212. The theoretical spectra have been simulated with the TI AHjFC and AHjFCHT models at the TD-CAM-B3LYP/aug-N07D/C-PCM level. [67] The line-shapes have been convoluted by means of Lorentzian distribution functions with HWHMs of 250 cm 21 .
anharmonic calculations or select the most intense vibronic transitions for instance, become critical as a validation done a posteriori on the results is generally not obvious. To illustrate another interesting peculiarity, we will focus on vibronic spectroscopy, and in particular electronic circular dichroism (ECD). When defining the computational model, a common way to choose between FC and FCHT is to look at the oscillator strength from the electronic calculations. For fully allowed transitions, FC can be sufficient, while for weakly allowed or dipole-forbidden ones, FCHT is necessary. However, such consideration does not hold for chiral spectra. The ECD intensity is related to the rotatory strength R I;F 5=fh v I jl e jv F ih v F jm e jv I ig, with m the magnetic dipole. At the FC level, the transition moments are constant and the ECD spectrum is a scaled OPA spectrum. This is not anymore true at the FCHT level since the products between the properties and their derivatives, which can be positive or negative, will scale differently the contributions of the HT (dl e =dQ i Á dm e =dQ j ) terms and mixed FC-HT terms (l e Á dm e =dQ i and dl e =dQ i Á m e ). An important aspect is that, even for fully allowed transitions, the FC contributions can be negligible if the transition moments of the electric and magnetic dipoles moments are near orthogonal (l e ðQ eq Þ Á =fm e ðQ eq Þg % 0). As an example, let us consider the case of (Z)28-Methoxy-4-Cyclooctenone (MCO) [67] where one broad band is observed in the 340-250 nm range of the experimental OPA spectrum in hexane solution, but a more complex pattern is present in the ECD one. [212] The OPA spectrum is satisfactorily described as the p Ã n transition of the 1A conformer computed at the TI AHjFC level (full convergence), but this is not the case for ECD, as shown in Figure 23. At this level, a single, broad negative band is present, compared to the two, positive and negative, narrower bands from experiment. By including HT contribution (AHjFCHT), a correct pattern is obtained, with the right sign change in the band-shape. Those results indicate that, for MCO, inclusion of the HT contributions was necessary to reach at least a qualitative agreement with experimental data. Because of the sensitivity of chiral spectroscopies, vibrational contributions are often more critical that for the nonchiral ones. While HT effects can be less important than in the case presented here, with the experimental band having a single sign, they can still modulate the shape and improve the overall agreement. [16,67] The situation can become even more complex when multiple electronic states are involved since vibronic band-shapes are very likely to overlap, canceling or enhancing each other, so that a proper representation of their envelope asymmetry is still necessary. [16,206,207] Other studies based on chiral spectroscopies, including vibrational ones, can be found in references 29, 213, and 214.
Switching to problems related to flexible molecules, an interesting case is the presence of torsional modes with quite low vibrational energies, which behave as an hindered rotor. The problem is very similar to the LAM seen for bithiophene, but the conformers connected by the torsional coordinate are separated by a low energy barrier. Because of the nature of the vibration, it must be separated from the rest of the vibrational treatment. This is done by mean of a hindered rotor identification algorithm. [215] The modes identified as coupled to the torsion are then ignored from the rest of the system and considered inactive in the VPT2 computations. Combination of these two approaches leads to the final hindered-rotor anharmonic-oscillator model (HRAO), necessary to compute correctly the thermodynamic quantities of flexible systems. [49,216] Another important aspect is that, due to the low-energy barrier, multiple conformers generally need to be included in the simulation of the spectra. A good illustration is the case of glycine, for which the three lowest energy conformers (Ip/ttt, IIn/ccc, and IIIp/tct, see Figure 24) have been concomitantly identified through their IR vibrational features in the low-temperature matrix environment. [201,217,218] In this condition, the identification of each conformer is based on the analysis of the intense vibrational transitions that remarkably differ in frequency from one conformer to another, for instance the C@O stretch, which are highlighted in Figure 24. The simulation of the overall spectra, which can be directly compared with experiment, is then performed by applying weighted contributions of each conformer according to their abundance obtained from the theoretical Boltzmann population at experimental conditions, such as the temperature. In this respect glycine represents an interesting case, since the harmonic approximation provides semiquantitative results for enthalpies and free energies at low temperatures, but at higher temperatures, relevant for the comparison with experiment, the entropy of the IIIp/tct conformer is strongly overestimated at the harmonic level, predicting this conformer to be the most stable one in contrast to experimental evidence.
However, the full HRAO approach is able to provide relative free energies and theoretical Boltzmann population at 410 K for all conformers in excellent agreement with experimental estimates (i.e., Ip 5 77%, IIn 5 18%, IIIp 5 6%, see reference 201 for the details). The overall simulated spectrum obtained by weighting the contribution of each conformer according to its abundance, is also depicted in Figure 24 and Theoretical IR spectra of the most stable conformers of glycine computed at the CC/DFT level [201] : single contributions from Ip/ttt, IIn/ccc, and IIIp/tct, and their sum weighted with relative abundances at 410 K (Ip-IIn-IIIp) compared to the experimental IR spectrum recorded in low-temperature Ar matrix. [217] Red, blue, and green arrows mark the bands assigned to the Ip/ttt, IIn/ccc, and IIIp/tct conformers, respectively. IR spectra line-shapes have been convoluted with Lorentzian functions with HWHMs of 1 cm 21 .
compared to its experimental counterpart. [217] It is demonstrated that the most intense transitions are indeed related to the most stable conformer, Ip/ttt, but finger-marks for IIn and IIIp are also clearly present.
As pointed out before, at variance with simple methodologies based on the double-harmonic approximation, fully anharmonic spectra have non-null intensities for overtones and combination bands. It is thus possible to recognize low-intensity features related to nonfundamental transitions of the most abundant conformer from the fundamental transitions of the less abundant ones. [219] The same computational procedure can be easily applied also to sets of isotopologues to facilitate an unequivocal identification of several conformers present in a given sample. [219] This strategy can be even extended to combine complementary vibrational spectroscopies, for example, IR and Raman, and can be useful for the identification and recognition of molecular species, like free radicals, present in experimental mixtures. [17]

| C O NC LU S I O N S
This study, focusing on thiophene and bithiophene, has shown that computational methods still commonly used for medium-to-large molecules, such as purely electronic vertical models for UV-vis spectra or the harmonic approximation in vibrational spectra, can fall short in providing a reliable description of experimental results. Thanks to the availability of cost-effective models targeted for systems of dozen of atoms or larger, a systematic and significant improvement of the overall simulations can be achieved. While implementations had been scarce in the past, limiting their adoption, this is not anymore the case as general-purpose programs have started to propose more advanced methodologies for the quantum treatment of nuclear motion. As a result, it is now possible to define more reliable computational protocols for the simulation of spectra. While their use can be facilitated with a careful automatization of the whole process, most models, due to their intrinsic limitations, are not without their own caveats and one should be cautions of the assumptions on which they are built, thus, their applicability. As shown in the case of bithiophene, those limitations, once correctly identified, can be overcome, by transforming the system in well-defined subsystems through the use of hybrid schemes.
An important aspect of such strategies is to maintain the computational cost affordable for medium-to-large molecules. In practice, this means that the largest part of the N-dimensional problem to solve must remain describable with the cost-effective methods presented here, the rest being represented by more refined, albeit expensive, models.
Such an approach can be very interesting as it provides a way to include progressively new, more accurate methodologies to tackle more complex molecular systems in a cost effective way. An important hurdle remains the subdivision of the system and the subsequent combination of the results, which is not necessarily straightforward for all cases, and the extensions of such hybrid computational models to other cases, like FCHT vibronic calculations or anharmonic vibrational spectroscopies for instance, will be the focus of future works.

ACKNOWLEDGMENT
The high performance computer facilities of the DREAMS center