FCclasses3: Vibrationally‐resolved spectra simulated at the edge of the harmonic approximation

Abstract We introduce FCclasses3, a code to carry out vibronic simulations of electronic spectra and nonradiative rates, based on the harmonic approximation. Key new features are: implementation of the full family of vertical and adiabatic harmonic models, vibrational analysis in curvilinear coordinates, extension to several electronic spectroscopies and implementation of time‐dependent approaches. The use of curvilinear valence internal coordinates allows the adoption of quadratic model potential energy surfaces (PES) of the initial and final states expanded at arbitrary configurations. Moreover, the implementation of suitable projectors provides a robust framework for defining reduced‐dimensionality models by sorting flexible coordinates out of the harmonic subset, so that they can then be treated at anharmonic level, or with mixed quantum classical approaches. A set of tools to facilitate input preparation and output analysis is also provided. We show the program at work in the simulation of different spectra (one and two‐photon absorption, emission and resonance Raman) and internal conversion rate of a typical rigid molecule, anthracene. Then, we focus on absorption and emission spectra of a series of flexible polyphenyl molecules, highlighting the relevance of some of the newly implemented features. The code is freely available at http://www.iccom.cnr.it/en/fcclasses/.


| INTRODUCTION
Electronic spectroscopy has ubiquitous applications in current research ranging from chemistry to materials science and biology. Different electronic spectroscopies investigate different properties of the system and in all cases the shape and position of the spectra arise from the interplay of intra-molecular and/or inter-molecular factors whose individual role cannot be easily extracted just from the inspection of the experimental trace. Thanks to recent progresses, computational spectroscopy stands nowadays as a mature and robust approach to complement the experimental measurements and help retrieving all possible information entailed in the spectroscopic signals. 1 It is well assessed that in order to describe the shape of the electronic spectra it is necessary to adopt theoretical models which explicitly address the coupling of the electronic transition with the nuclear motion (vibronic transitions). 2 The same is true for the simulation of the rates on nonradiative processes. The degree of complexity of the just formulated problem can be formidable, making it impossible finding a general solution.
Here we focus on the class of relatively large systems (dozens or hundreds of vibrational coordinates) which are rigid or with few flexible degrees of freedom (DoFs), so that harmonic approximation represents a useful framework to describe their PES, and whose electronic states involved in the radiative transition exhibit negligible or weak inter-state couplings. The field of the simulation of the vibronic spectroscopy of these systems has made remarkable progress in the last two decades thanks to the contributions of several groups, which is impossible to review here (see for instance the review papers in Refs. [2][3][4]. Our group has been active in this context and in this contribution, we present FCclasses3, the new and largely renovated and extended release of our code FCclasses distributed in 2009, which makes freely available all the computational methods we have developed in the last decade in the field of vibronic spectroscopy. As we mentioned above our computational approaches are based on the hypothesis that initial and final electronic states are not significantly coupled to any other state. Otherwise, the solution can be obtained through either time-independent approaches, like the Lanczos method [5][6][7] or numerical quantum dynamics propagation of the wavepacket over the coupled states. 5 Moreover, all PES are described within harmonic approximation, which is generally valid for stiff modes but inadequate to account for flexible and/or very anharmonic regions of the PES. The coupling triggering the transition between the initial and final states is generally considered small enough so that a Fermi Golden rule approach can be adopted. Finally, the nuclear dependence of electronic coupling terms, which can be expressed as a Taylor expansion, are usually described taking only the constant term (Franck-Condon, FC, approximation) or up to the linear term (usually called Herzberg-Teller, HT). The HT approximation introduces perturbatively the effect of inter-state couplings when they are weak and the states are well separated. Whereas most of the developments reported in literature focus on linear-spectroscopies, extension to some nonlinear spectroscopies like resonance Raman is straightforward, [8][9][10] while, vibronic computations of other signals like two-photon absorption, 11,12 circular-dichroism 12 or magnetic circular dichroism 13 within this framework are still doable but they require some additional approximation on the frequency-dependence of the transition amplitudes.
Adopting all the approximations sketched above, analytical expressions emerge for the intensities of the different spectroscopies and nonradiative processes, which can be formulated in terms of sum-over-states (what we will refer to as time independent [TI]) or Fourier-transform of correlation functions (what we will refer to as time-dependent [TD]). The basic general fundamental equations for both formulations have been known for decades. [14][15][16] However, the limited availability of second derivatives of the PES, which are required for harmonic models, made scarcely useful the development of efficient computational schemes for practical computations. The derivation of analytical formulas for first and second derivatives of the excited-state energy for different electronic structure methods, primarily the popular time-dependent density functional theory (TD-DFT) has completely changed the scenario in the last decades, stimulating, in the last 20 years, a fast flourishing of efficient methods for specific cases within both TI [17][18][19][20][21] and TD [22][23][24][25][26] formulations, and their widespread applications in literature. In the case of TI approaches, the main challenges to be overcome are related to the need to evaluate the intensity of the vibronic transitions between a vast number of possible initial and final states. This difficulty has called for the elaboration of effective filtering protocols to select only the most relevant transitions, and fast ways to compute them. Just to make an example, for a typical semi-rigid molecule with 100 normal modes millions or even billions of FC overlaps need to be evaluated for a spectrum convergence >90% (defined as the recovered fraction of the total intensity evaluated from analytical sums). Regarding TD approaches, the major developments have gone through deriving the exact analytical expression of the correlation function for different cases (under the harmonic approximation), including one photon spectroscopies with FC and HT terms, resonance Raman, 9,27,28 nonradiative process as internal conversion and intersystem crossing, 29,30 and two-photon absorption. 31,32 It is interesting to note that TI and TD formulations are somehow complementary. TI provides full information about the actual vibronic transitions responsible for the main spectral features, but can suffer from convergence problems. TD does not provide information on the individual vibronic transitions but it delivers fully converged spectra, that is, effectively accounting for all possible transitions.
Another essential point to consider are thermal effects, which tune the distribution of initially populated states. In TI formulations, this calls for an update of the filtering method to effectively account for the additional transitions. 33 On the contrary, in most cases TD formulations already account for temperature effects at no additional cost. It should be added that in some specific spectroscopies the electronic transition takes place from a nonthermal initial vibrational population, created by a infrared pre-excitation, leading to the so-called VIPER spectroscopy. This IR pre-excitation can be accounted for in both TI and TD formulations. 34 FCclasses3 is a program, written in Fortran90 with legacy parts in Fortran77, to compute electronic spectra and kinetic rates of nonradiative processes, including vibrational effects at the harmonic level.
The present code largely extends the capabilities of the program FCclasses, 35 which was based on an efficient algorithm for the calculation of one photon spectra adopting TI formulations. Concretely, among other features, this new version implements the full family of adiabatic and vertical harmonic models and it generalizes harmonic models allowing the quadratic expansion at arbitrary geometries using either Cartesian (valid only at stationary points) or valence internal coordinates. The number of supported properties has been notably increased from the ones included in the previous release, namely onephoton spectroscopies as one-photon absorption (OPA) and emission (EMI) and electronic circular dichroism (ECD). FCclasses3 can now compute additional one-photon chiral spectroscopic techniques: Circularly polarized luminescence (CPL) and magnetic circular dichroism (MCD) as well as two-photon spectroscopies, including absorption (TPA) and circular dichroism (TPCD), and it also supports resonance Raman (RR). Regarding nonradiative processes, it includes the calculation of internal conversion (IC) rates induced by nonadiabatic couplings and more generally of any transition between (diabatic) states driven by a constant scalar coupling (NRSC). Finally, the computation of all available properties can be performed with both TI implementations, exploiting and extending the TI algorithm present in the   previous version of the code, and new TD approaches, through purposely derived analytical correlation functions. As mentioned before, the simulation of vibronic spectra has attracted much interest in the past years. Beside the previous version of FCclasses (2.1), distributed since 2009, there are currently available several other codes devoted to such task such as FCBand 36 or ezSpectrum, 37 and it is now an add-on included in many electronic structure programs including Gaussian16, 38 Orca5.0, 39 Q-Chem, 40 or Molpro. 41 The important contributions in these fields by groups that implement their derivations in codes shared to different extents with the community should also be mentioned. The complete list is too vast to be listed here entirely, but, to name just a few, we can mention the contributions by Grimme, 17 Peluso, 42 Berger, 43 or Nooijen. 20 It is not the scope of the present article to make a detailed comparison of the different features of all the available codes. To the best of our knowledge, FCclasses3 stands among the most flexible and complete packages, and it includes a number of unique features, such as the possibility to define reduced dimensionality models through the availability of projectors in curvilinear internal coordinates, 44 the simulation of spectra starting from selected vibrationally excited states (which is the base of VIPER spectroscopy), 45 or the availability of some more exotic properties like the two-photon circular dichroism. 12 Moreover, we made a particular effort to give the user precise control of the different options in the different computational protocols and feed input data in several ways to ease its use.
All methods discussed up to now are formulated adopting the harmonic approximation, limiting the applicability to relatively rigid systems. The range of validity of the harmonic approximation is however somewhat dependent on the system of coordinates adopted to describe the normal modes. Namely, the use of linear coordinates, as Cartesian, can lead to severe issues when dealing with significant displacements along curvilinear DoFs, such as torsions. In these cases, Cartesian coordinates, artificially couple low-frequency modes (associated with such large-amplitude displacements) with high-frequency ones introducing large artifacts in the spectra. 46,47 When curvilinear coordinates (e.g., valence internal coordinates) are used instead, and this is possible in FCclasses3, such spurious couplings vanish. It has been shown that this extends the range of applicability of the harmonic approximation, preventing the occurrence of such artifacts. 48,49 It remains clear that, when significant anharmonicities of large amplitude motions come into play, they need to be dealt with specific methodologies beyond harmonic approximation.
FCclasses3 does not compute anharmonic vibronic spectra. The literature on aharmonic vibronic computations is too broad to be reviewed here. We just mention that focusing on TI approaches, anharmonic vibrational states in very few coordinates (and the corresponding FC overlaps) can be computed with straightforward variational approaches, and many groups have developed their own house-made codes for these tasks. When anharmonicity involves many dimensions, sophisticated approaches like the vibrational configuration interaction (VCI) or the vibrational coupled cluster (VCC), 50,51 can be applied but they require a quite different computational machinery; the interested user is referred to the implementations available in codes like Molpro 41 or MidasCpp. 52 Finally, also TD approaches can be profitably used for anharmonic systems, exploiting for instance the potentiality of the multiconfigurational timedependent Hartree method (MCTDH), 53 and its implementation in very versatile distributed codes like Quantics. 54 In this context, it is worthy to highlight that the use of curvilinear Hessian approach (Ad-MDjgVH) we proposed recently. 56 In this method, which is based on FCclasses3 and its interface with other codes, anharmonicity along soft coordinates is accounted for through molecular dynamics samplings, whereas stiff harmonic modes are treated at quantum vibronic level on effective PES that depend on the instantaneous position of the soft DoFs.
FCclasses3 program has been strongly re-structured, making it much more modular for helping its reading and modification. The main input style is completely renovated with a easy-to-follow keywordbased structure. Special care was put to make the code easily usable in combination with different codes, input sources, and computational protocols. Therefore, the code accepts inputs in both Cartesian and normal coordinates basis (or mixed forms), we establish a transparent and well documented format for providing the auxiliary input data, like geometries, and energy and transition amplitude derivatives, and we provide interface tools that automatically generate such data from many widespread quantum chemistry codes. Graphical tools for helping visualization and assignment of the results are also furnished. Section 4 is then devoted to some test applications, aiming at showing both the simple workflow to deal with the calculation of different properties and some key features to cope with scenarios where the harmonic approximations is challenged. Finally, Section 5 reports concluding remarks discussing possible perspectives for further development within the code.

| Harmonic PES models and normal mode analysis
In FCclasses3 both initial and final states are described by quadratic PES, that is, retaining up to the second term in the Taylor expansion around the reference geometry. The corresponding normal modes are adopted as a set of coordinates for each state to pose the vibronic problem. Normal coordinates are generally different for each electronic state, and they are related through the Duschinsky relation, 57 which in the form used in FCclasses3 reads: where Q 1 and Q 2 refer to the normal modes at the initial and final states, J is the so-called Duschinsky matrix and K are the normal mode displacements.
FCclasses3 can adopt both adiabatic (A) and vertical (V) approaches to build up the harmonic PES. In fact, in most common vibronic approaches the initial-state PES is always quadratically expanded around its own equilibrium geometry (it is however noteworthy that FCclasses3 provides some generalization, described below). In contrast, for the final state, the expansion is usually performed around either the equilibrium structure of that state (A) or the equilibrium structure of the initial state (V), that is, the so-called Franck-Condon point. In previous works, 58 we have designated such strategies as Adiabatic Hessian (AH) and Vertical Hessian (VH) models, respectively. Once normal modes are computed for the initial and final states, the Duschinsky matrix is evaluated as, The calculation of the normal mode displacements depends on the model adopted. For A models, it simply corresponds to the difference between the reference equilibrium geometries: where ξ is the set of coordinates adopted, either Cartesian or valence internal ones.
In the case of V models, the shift is evaluated finding the minimum of the quadratic final state PES, using both the gradient and the Hessian. It should be noted that, when V models are applied, normal modes for the final state are computed at a nonstationary point of the PES, so that in Cartesian coordinates it is not possible to rigorously separate rotations and vibrations. 59 Indeed, in this case, the use of curvilinear set of coordinates is required, such as valence internal coordinates. In FCclasses3, the minimum over the final state potential is first obtained in valence internal coordinates: where H s and g s are the Hessian and gradient for the final state in terms of valence internal coordinates. This displacement is then transformed to normal modes as indicated in Equation (3).
Simplified versions of adiabatic and vertical models can be generated. The most popular approximation imposes that all states share the same energy Hessian, leading to models that have many names in literature such as independent mode displaced harmonic oscillator (IMDHO) model 60 or linear coupling model (LCM) 11 and that we proposed to refer to as adiabatic shift (AS) and vertical gradient (VG) models so to make apparent the family of models they belong to (A or V). 58 Intermediate-quality models can be obtained keeping the same normal mode coordinates for both the initial and final states but updating the frequencies of the final state with the diagonal elements of the final state Hessian, leading to ASF and VGF models (where F stands for frequencies). A critical assessment of the relative quality of these models can be found in Ref. 58.
V and A models are schematized in Figure 1  The Hessian in curvilinear valence internal coordinates is obtained from that in Cartesian coordinates with the following expression: where g t s is the gradient in internal coordinates (thus vanishes at stationary points) and β is a rank-3 array whose elements correspond to the derivatives of the B matrix with respect to Cartesian coordinates.
Although the latter can be evaluated analytically, FCclasses3 implements numerical derivatives which revealed faster to compile and very stable and precise. 48 As mentioned above, internal coordinates are more convenient than Cartesian to apply the AH model when large displacements are found for curvilinear modes between the initial and final state minima. 47 Moreover, they are the only rigorous choice for VH model. 48 However, it should be noted that they also find some limitations.
Namely, the G matrix actually changes with the coordinates, which leads to additional nonharmonic terms in the Hamiltonian. 64 Such terms are ignored within FCclasses3, assuming that the G does not change drastically from the initial to final state minima. One way to analyze the significance of such terms is comparing the lineshape obtained with the VH model evaluating the G matrix at either the initial state geometry (default choice) or the extrapolated minimum of the final state. 48 This test can be performed within the code with a simple keyword in input.

| Reduced-dimensionality models in curvilinear internal coordinates
Valence internal coordinates further provide a suitable framework for projecting out some selected DoFs from the coordinate space used for vibronic calculations at a harmonic level. This strategy can be particularly useful for flexible coordinates related to anharmonic regions of the PES. It is indeed the basis of hybrid methods, which retain the harmonic description for the stiff modes while treating the anharmonic ones with more appropriate methodologies. 56,65,66 The projection to remove a curvilinear element has been already described in the literature, 67 and it can be formulated as follows 56 : where both the expressions both covariant, such as the gradient, g, and contravariant objects, such as an internal coordinate, s, are given.
Note also that the calculation of the modulus generally includes the metric tensor, that is, jsj 2 = s t G À1 s and jgj 2 =g t Gg.
Several coordinates can be removed by iteratively applying and updating the projector, as we recently proposed. 56 The actual protocol implemented in FCclasses3 is summarized in Figure S1.

| Supported spectroscopies and photophysical properties
FCclasses3 implements many different spectroscopies, namely OPA, EMI, ECD, CPL, MCD, TPA, TPCD and RR, and the computation of nonradiative rates driven by nonadiabatic couplings, IC, or any scalar coupling expanded up to the first-order in terms of the nuclear coordinates. In order to avoid a too long manuscript we give most of the details and derivations in the Supporting Information. Here in the manuscript, we only briefly discuss one-photon spectroscopies (OPA, EMI, ECD, and CPL) which can be described with common equations and allows us to introduce both TI and TD approaches.
The calculation of the spectrum, S(ω), which corresponds to the quantities experimentally recorded for each spectroscopy (molar absorptivity ε for OPA, its anisotropy Δε for ECD, number of emitted photons per unit of time I emi for EMI, …) can be formulated as, where L(ω) is usually referred as the lineshape and n is an integer whose value also depends on the type of spectroscopy. Concretely, n = 1 for absorption processes (OPA, ECD) and n = 3 for emission F I G U R E 1 Schematic representation of (A) adiabatic Hessian and vertical Hessian models and (B) generalized vertical Hessian models (gVH).
(EMI, CPL). C is a constant that depends on the type of spectroscopy.
As an example for OPA, when all variables are described in atomic units, C ≈ 703 leads to intensity as ε(M À1 cm À1 ). 68 The expression for the lineshape at temperature T within a TI formalism is given by the following sum-over-states, where ΔE is the adiabatic energy difference between the two PES in The equivalent expression within the TD formulation is given by the following Fourier transform: where g ' (t) is a damping function related to the broadening of the spectrum. Gaussian, Lorentzian and Voigt broadening profiles are implemented in FCclasses3 in combination with a TD calculation. χ is the transition dipole correlation function, which for OPA/EMI is given by, where τ f ¼ t=ℏ and τ i ¼ Àiβ À t=ℏ, and b H i and b H f are the Hamiltonians for the initial and final states.
Analytical expressions are available for both the vibronic transition dipole integrals, 14 ⟨v i jμ if jv f ⟩, and the correlation function 22,23,26,42,70 when the harmonic approximation is adopted and the electronic transition dipoles are described with FC and HT terms.
Finally, in the case of emission, the code also computes the radiative rate constant from the integration over the emission spectrum, 29 Whereas TD approaches are eigenstate-free approaches based on the fast evaluation of analytical expressions for time-correlation functions and their Fourier transform, TI approaches are based on the evaluation of the spectra and properties through an explicit sum over all the possible vibronic states. This is in principle a challenging problem, since their possible number increases so steeply with the dimension of the system that a brute-force calculation is out of the current possibilities for classical computers for the typical energy range covered by a spectrum and medium-size molecules. In fact such number is $10 20 for coumarin 153 (102 normal modes) for an energy window of $ 1 eV. 19 Such a challenge is so well recognized that vibronic computations were individuated as ideal test cases for showing quantum advantage in quantum computers. 71

| RUNNING FCCLASSES3
In this section, we review some details about the usage of the code, including the structure of the input file and the description of the tools provided with the package. More details about the implementation, for example, regarding the Fourier Transform, are given in the Supporting Information S1.

| Input and output files
The program options are controlled by an input text file structured in lines stating OPTION=VALUE assignments. The order of entries is flexible, and sensible defaults are given to many options to facilitate job preparation. An example input file is given for the first test case below (Listing 1).
Additionally, the user must also provide the information about the PESs (structure, energy, gradient and Hessian) of the initial (1) and final (2) states through the state files (STATE1_FILE, STATE2_FILE), and about the coupling term (in the case above, electronic transition dipole moments) through other additional files, specified in the input.
They all follow an intuitive format designed explicitly for FCclasses3 and described transparently in the manual to allow users to easily compose the needed data from any available source. In particular we also provide some tools (described below) to automatically generate these input files from the output files of different popular quantum-chemistry programs. It is worthy to specify that for HT computations

| An example of a rigid system: Anthracene
In this section, we focus on anthracene, a paradigmatic case where the harmonic approximation delivers an extremely good framework to simulate the electronic spectroscopy. 58 We use this example to illustrate the new easy-to-use workflow, from input preparation for the simulation of different properties, to the analysis of the stick transitions. As indicated above, the use of the pre-processing tools to extract the relevant data from the output of electronic structure program along with the generation of template inputs, allows the user to quickly set up all required files to run a calculation. All calculations in this section were carried out using DFT and TDDFT in the gas phase, adopting the B3LYP functional in combination with the 6-31G(d) basis set. Gaussian16 38 was used to perform all electronic structure calculations, except for the computation of the two-photon transition tensors, which were computed with Dalton 2016. 80 We start simulating the OPA spectrum. The required input to carry out the simulation of the spectrum at 300 K is shown in Listing 1.
Although almost all keywords adopt sensible defaults, we make The emission spectrum can be computed with a very similar input, setting PROPERTY = EMI and exchanging the data for the initial and final states. In addition to the emission spectrum, which for anthracene is nearly mirror symmetric with respect to absorption (see Figure S2), the calculation also provides the corresponding radiative rate (k r ) computed from Equation (11). For this transition, that is bright and correctly described at FC level, the rate is mainly dictated by the value of the transition dipole moment, and only slightly tuned by the vibronic shape. A value k r = 1.9 Â 10 7 s À1 is obtained, which is close to previous simulations 29 and experiments. 81 We now turn to explore the TPA spectrum of anthracene. Both S 1 (B 1u ) and S 2 (B 2u ) are TPA forbidden by symmetry. However, an experimental TPA signal in the spectral window associated to these states has been reported for anthracene in solution and assigned to the B 2u state. 82 Here, we simulate the TPA spectra for both S 1 and S 2 LISTING 1 Example input file structure for the calculation of a one photon absorption spectrum of anthracene at 300 K. Text after semicolons is ignored. were generated with a local script. 83 In order to ensure that exactly the same normal modes are used, we first run a FCclasses3 job that writes normal modes to files, which are subsequently read by the external script. Following this procedure, and taking into account that the tensor at the reference geometry is exactly zero, we lose the trace of the sign of the derivative. The sign, however, is not relevant within VG model and, for this reason, this is the model adopted to run the vibronic calculations. Energies, gradients, and Hessians are evaluated with the same functional and basis with Gaussian16. The input for FCclasses3 is similar to that of OPA. We simply need to change the keyword for the PROPERTY to TPA, and update the relevant settings for the state files, temperature (300 K) or harmonic model (VG). The two-photon transition tensor and its derivatives along normal modes are read from a file prepared from the finite differences calculation mentioned above.
The simulated TPA spectra including the symmetry forbidden transitions to S 1 (B 2u ) and S 2 (B 3u ) and the TPA active S 3 (B 1g ) are shown in Figure 3. HT effects are obviously required to reproduce the experimental observation of the symmetry forbidden bands, and the vibronic calculations included in FCclasses3 perform quite well in this sense. In this case, the limited resolution of the experimental measurements, and limited accuracy of the TD-DFT calculation, which does not account for potential multireference character, 84 does not allow a careful assessment of the simulated lineshape. Still, our results provide a correct semi-quantitative picture of the spectrum. Namely, the S2 (B 3u ) is much more intense than the S 1 band, which is consistent with the experimental evidences, although the ratio between both bands seems to be overestimated. 85 The TPA active transition leads to a more intense band, also in agreement with the experiment.
Note that, for the latter, only FC terms are included in the calculation as it is a TPA bright transition. Overall, taking into consideration the F I G U R E 2 (Top) Screenshot of the GUI tool to analyze TI results showing the analysis of the TI absorption spectra of anthracene. It includes the experimental spectrum (shifted by À0.30 eV). For the convolution we used a Lorentzian function with HWHM = 0.02 eV. (Bottom) The same plot exported to XMGRACE, where the TD spectra (at 0 and 300 K) have been added, computed with the same convolution scheme.
shift that was already required for the OPA spectrum (0.3 eV), the general reproduction of TPA spectrum is not far from the experiments, and actually supports the assignments of the main bands provided in Ref. 85.
We continue the spectroscopic analysis of anthracene by simulating its RR spectrum (at 0 K). The experimental spectrum in different solvents recorded using an excitation frequency in the UV region has been recently reported including also a detailed computational analysis. 86 Here, we will compare our gas phase results with those recorded in a nonpolar solvent, cyclohexane. The excitation frequency used in the experiments is nearly resonant with the peak of the brightest transition. According to our calculation in the gas phase, this state corresponds with S5 (B 2u ), which has a large oscillator strength (f = 1.9). We first simulate the OPA spectrum for the S5 state with the VH model, which results in vibronic profile with a maximum at the 0-0 band, consistent with the experiment, although appreciably narrower (see Figure S3 in the SI).
The simulation of the RR is carried out taking the energy of the maximum of the simulated spectrum ($5.1 eV) as incident frequency. We adopt the VH model and FC + HTi approximation we and account for fundamentals and overtones with two quanta for each mode whereas for combinations of two modes (1 + 1 quanta) we selected the modes with stronger fundamentals. For the computations, we use a damping function Γ = 150 cm À1 , following the simulations already reported for this system, 86 and convolute the resulting sticks with a Lorentzian function with HWHM = 5 cm À1 (along the Raman shift). Both TI and TD approaches are adopted, which deliver practically the same results (TI convergence is >99%; see Figure S4)  Table S1. As expected by symmetry, almost all fundamentals correspond to A g symmetry modes. The only exception is mode 55 (CC stretching), which belongs to B 3g and arises from the HT contribution. We also identify the over- In order to easily explore the 2D data (RR intensity at different incident and scattered frequencies), the tool fcc_RRinspector_PyQt5 can be used, which allows us to interactively inspect the RR spectra at F I G U R E 3 TPA spectrum of anthracene computed at 300 K with the VG model for S 1 , S 2 , and S 3 states. S 1 , S 2 (TPA forbidden) are computed with FC + HTi terms and applying a Gaussian broadening with HWHM = 0.02 eV. S 3 (TPA bright) is computed including FC terms only and a Gaussian broadening with HWHM = 0.05 eV. The experimental TPA spectrum adapted from Ref. 84 is also included. Simulated spectra are blue-shifted by 0.3 eV, consistently with the shift described for OPA in Figure 2.
For better visualization, the S 1 band is scaled by a factor of 40.

different incident frequencies (see movie in the Supporting
Information S1).
Finally, we turn to explore the nonradiative internal conversion rates of anthracene from S 1 , triggered by the nonadiabatic coupling with the ground state. The nonadiabatic coupling vector is computed at TDDFT level as implemented in Gaussian16. The input is very similar to that for emission, but setting PROPERTY = IC and specifying, instead of the ELDIP_FILE, the NAC_FILE, which contains the nonadiabatic coupling matrix elements in terms of Cartesian coordinates.
The profile for nonadiabatic rate vs the adiabatic energy obtained with the AH model is shown in Figure 5. The plot also includes the results with the TI method at 0 K, with a convergence of 90%. Since the rate constant is typically taken at adiabatic energy values corresponding to the decaying wing of the function in Figure 5, far from the 0-0 transition, a region which usually challenges the convergence of the TI calculation, the adoption of this implementation is encouraged only for interpretative purposes (e.g., assignments of the role of specific progressions).
In this section, we have illustrated the versatility of the program and its ease of use. All calculations require a minimal input from the F I G U R E 4 (Bottom) RR spectrum for anthracene with the incident frequency at the maximum of the bright UV transition (S5 in gas phase), computed with the VH model and FCjHTi dipole terms, and convoluting the stick transitions with a Lorentzian with HWHM = 5 cm cm À1 . All fundamental and overtone bands are accounted for, as the combination between the most intense fundamentals (8, 19, and 42). user, since many of the settings adopt sensible values, which can anyway be overridden with the corresponding keyword in input. The molecule adopted for this systematic analysis is rather rigid and represents a specially favorable case for the application of the harmonic approximation. In the next section, we will explore some systems where the harmonic approximation is not applicable in a straightforward way and discuss how can we approach such situations within FCclasses3.

| Vibronic spectroscopy of flexible systems
We turn now to show some of the new features of the code, which can be especially useful to treat flexible systems. Concretely, we approach two polyphenyl molecules: p-terphenyl and p-quaterphenyl (see inset in Figure 6), whose experimental spectra are obtained from Ref. 87. The simplest molecule from this series, that is, biphenyl is excluded from this study as TDDFT faces some issues to describe low lying excited states, so that multireference approaches such as CASPT2 88 or SAC-CI 89 would be required to deliver an accurate picture. The flexibility of the molecules is related to the inter-ring torsions and put in crisis the deployment of the harmonic approximation.
The interest of such systems in optoelectronics, for example, as building blocks for OLEDs, further motivates the accurate simulation of their optical properties. 90 The energy, gradient and Hessian of the ground and bright excited states are evaluated with DFT and TDDFT at the equilibrium geometries of each state, so to be able to adopt both adiabatic and vertical models. The PBE0 functional along with the 6-31G(d) basis set is adopted. The effect of the solvent (cyclohexane) has been accounted for with the IEFPCM model. 91 In Figure 6, tures have a remarkable effect on the spectral shape; indeed, the differences in the shape of the S 0 and S 1 PES is responsible for a marked symmetry-breaking between absorption and emission spectra. 92 Such a situation can be approached with different strategies within FCclasses3, adopting in all cases curvilinear internal coordinates. In the F I G U R E 6 Relaxed potential energy curves corresponding to one flexible dihedral of the different polyphenyls under study at PBE0/6-31G(d). Over the ground and excited (S 1 which is the bright H ! L transition). The profile for biphenyl is also included for comparison. Inset: Schematic representation of polyphenyl structure, where n = 1, 2, for p-terphenyl and p-quaterphenyl, respectively following, we discuss the appropriateness of the different models available within the code.
In Figure 7 we report absorption and emission spectra adopting can then be treated independently to account for the actual potential energy profile more appropriately. A similar strategy was already adopted by some of us to study the absorption and emission spectra of dithiophene, which shows a torsion with a similar change between S 0 and S 1. 55 In that work, the contribution of the torsion was evaluated both variationally, that is, computing the anharmonic energy profile, fitting it and then diagonalizing the vibrational Hamiltonian on a basis of sines and cosines, and classically, that is, obtaining the distribution of vertical energies over the torsion according to a Boltzmann population. The final spectrum was computed by convoluting the profile corresponding to the torsion and the lineshape for the remaining DoFs treated quantum mechanically under the harmonic approximation. We can follow the same approach for p-terphenyl and p-quaterphenyl. Since the contribution of the torsion at room temperature was shown to be correctly described classically in Ref. 55, we here adopt the classical approach only, which has also the benefit that the computation can be in practice completed running only FCclasses3.
The distribution that describes the torsion lineshape is computed by building a weighted histogram from the energy differences between the relaxed energy curves (Figure 6), and it is shown in Figure 8. We assume that the different torsional DoFs are independent within a given molecule, so that the total energy is computed by the sum of the individual (1D) relaxed scans. As shown in Figure S6, this is a good approximation for S 0 and still remains reasonable for S 1 . A cubic interpolation has been used to increase the sampling along the curves, and the weights for the histogram have been computed as where E i is the energy with respect to the bottom of the curve for the initial state. 93 The actual contribution of the torsion to F I G U R E 7 (Bottom) Absorption and emission spectra simulated at 300 K for p-terphenyl adopting different harmonic models: AH (both in Cartesian and internal coordinates) and VH. In the simulation of the emission spectra with VH model an imaginary frequency arose, which was arbitrarily turned real. (Top) Experimental absorption and emission lineshape spectral, adapted from Ref. 87.
the spectral shape is modeled with a Gaussian with the same width In summary, the total spectra are then computed by first building reduced vibrational spaces where the torsions are projected out with the iterative protocols implemented in the code. 56 In such a space the normal modes with imaginary frequencies present in the full space are removed (2 for p-terphenyl and 3 for p-quaterphenyl). Afterward, the contribution of the torsions is included with the convolution with a Gaussian with the computed HWHM. The results are shown in Figure 9, and confirm the good performance of the proposed simple approach. Namely, the spectral shapes generally well reproduced for both emission and absorption. For emission, some small changes when going from p-terphenyl to p-quaterphenyl are properly captured, as the relative intensity of the first (almost a shoulder) and second (most intense) peaks. For absorption, the overall width is slightly underestimated. This observation could be partly explained by the expected narrower distributions arising from classical sampling 44 or by the possible inaccuracies in the excited state PES around the FC point, which reveals critical to determine the spectral width. For F I G U R E 8 Distribution of transition energies at 300 K over the torsional degrees of freedom for p-terphenyl and p-quaterphenyl, both for S 1 ! S 0 (top) and S 0 ! S 1 (bottom).
F I G U R E 9 Absorption (top) and emission (bottom) lineshape spectra simulated with the VH model using PBE0. The torsions is removed from the calculation and the effect is recovered by convolution with a Gaussian having the same HWHM as the distributions showing in Figure 8.
instance, Figure S7 shows that using the profile computed for biphenyl to generate the distributions for all the dihedrals of p-terphenyl and p-quaterphenyl, the total width of the final spectra changes significantly.
In this section, we illustrated different strategies embedded in In order to help users to adopt the code for their own purposes, in this new release of the code we put special care in defining transparent standards for feeding information into the code, so to give the possibility to generate them from any source and exploiting different sets of coordinates (normal, Cartesian, internal). At the same time we remarkably increased the number of output data one can get, including for instance, energy and transition amplitude gradients (in normal coordinates), or Duschinsky matrix, for possible uses external to the code. For all the above reasons, we believe that FCclasses3 represents a useful tool that can be used also by nonexperts in the field. At the same time, it provides some advanced features, such as the generation of reduced-spaced models, that can be of interest in more Brownian oscillator models for simulating bi-dimensional electronic spectra. 96 The extended vibrational analysis tools embedded in FCclasses3 will also be exploited for working out novel advanced tools to investigate and interpret pure vibrational spectroscopies like infrared for complex multichromophoric systems. Moreover, we will also continue our efforts to improve the structure and modularity of the code. Indeed, in the long term, we plan to move towards a pluginbased paradigm, where the different parts of the code, like the ones dedicated to vibrational analysis, or projectors in curvilinear coordinates, and tools could be called independently from other programs.
Finally, regarding the tools, plans involve moving towards popular libraries, such as cclib, 97 and data formats, such as QCSCHEMA, 98 to extend the interfaces to more QM codes reducing maintenance efforts.

ACKNOWLEDGMENTS
The authors are indebted to FranciscoÁvila-Ferrer, Daniel Aranda, Na Lin, Yanli Liu, Daniele Padula, and Qiushuang Xu, who have been testing in the last decade the development version of the code. We also deeply thank all the colleagues with whom we have collaborated in these years and, in particular, those that stimulated and helped us to extend the capability of the codes to novel spectroscopies. They are too many to be mentioned here but their names are cited together with the relevant papers in the bibliography. Javier Cerezo thanks

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.