Assessment of the Potential Energy Hypersurfaces in Thymine within Multiconfigurational Theory: CASSCF vs. CASPT2

The present study provides new insights into the topography of the potential energy hypersurfaces (PEHs) of the thymine nucleobase in order to rationalize its main ultrafast photochemical decay paths by employing two methodologies based on the complete active space self-consistent field (CASSCF) and the complete active space second-order perturbation theory (CASPT2) methods: (i) CASSCF optimized structures and energies corrected with the CASPT2 method at the CASSCF geometries and (ii) CASPT2 optimized geometries and energies. A direct comparison between these strategies is drawn, yielding qualitatively similar results within a static framework. A number of analyses are performed to assess the accuracy of these different computational strategies under study based on a variety of numerical thresholds and optimization methods. Several basis sets and active spaces have also been calibrated to understand to what extent they can influence the resulting geometries and subsequent interpretation of the photochemical decay channels. The study shows small discrepancies between CASSCF and CASPT2 PEHs, displaying a shallow planar or twisted 1(ππ*) minimum, respectively, and thus featuring a qualitatively similar scenario for supporting the ultrafast bi-exponential deactivation registered in thymine upon UV-light exposure. A deeper knowledge of the PEHs at different levels of theory provides useful insight into its correct characterization and subsequent interpretation of the experimental observations. The discrepancies displayed by the different methods studied here are then discussed and framed within their potential consequences in on-the-fly non-adiabatic molecular dynamics simulations, where qualitatively diverse outcomes are expected.


Introduction
Over the last decade, a great deal of work has been done and yet plenty of controversy still surrounds the characterization of the main intrinsic photochemical mechanisms governing the photoinduced processes occurring in the DNA nucleobases [1][2][3][4][5][6][7][8][9][10][11]. Insight into their photochemical pathways provides unique information to ascertain the effect of light irradiation on the cellular system and is tightly related to the ability possessed by the genetic material to dissipate the excess of energy gained upon absorption in a harmless manner. This feature, known as photostability [12,13], is ascribed to the ultrafast decay of the excited state population to the ground state in a non-radiative manner mediated by conical intersections (CIs) between the initially accessed bright 1 (ππ*) and the ground state, justifying the sub-ps deactivation registered experimentally [3] and explaining the reduced damage registered in DNA despite our constant UV-light exposure [14,15].
Several components in different time ranges have been detected over the years in time resolved experiments on DNA nucleobases: in the sub-ps regime, a bi-exponential decay was characterized for nucleobases, nucleosides or nucleotides [3,24,60,62], being mainly ascribed to ultrafast reactions along the initially accessed 1 (ππ*) PEH, whereas a single component has also been observed in the few-ps timescale and being ascribed to dark nπ* states [55], and a longer component closer to the ns regime has been attributed to the eventual involvement of triplet states [55,56]. Theoretically, most studies agree in the initial involvement of the 1 (ππ*) state for the sub-ps channels, having characterized a series of key structures in its PEH that rationalize its concrete role. Figure 1 displays these structures for thymine, the system under study in the present work. They correspond to two different excited state minima, one planar [11] and one partially twisted or "boat-like" found in more correlated computations or solvated environments [24,34,63,64], and an "ethene-like" CI [65,66] connecting the 1 (ππ*) and ground states. The CI is usually employed to rationalize the ultrafast decay of thymine [8] and the extremely low fluorescence quantum yields recorded experimentally [62], and was initially reported in the pioneering study of Perun et al. [22]. The involvement of either dark nπ* or triplet states are usually ascribed to non-adiabatic interactions along the deactivation routes in these systems, which populate nπ* states that can in turn funnel the population towards the triplet manifold via intersystem crossing [67] given the large spin-orbit couplings featured [30,68,69], in agreement with El-Sayed rules [70].
In this paper we focus on the assessment of different multiconfigurational electronic structure theory approaches and geometry optimization techniques to benchmark the PEHs of the lowest-lying excited states in thymine and calibrate the extent to which the theoretical method employed can influence the model put forth to explain the photoinduced phenomena in DNA. To this end, widely used CASSCF approaches are compared to more expensive but currently feasible dynamically correlated CASPT2 methods, elucidating whether more correlated approaches are essential to achieve a proper description of the photoinduced phenomena. This comparison on static computations draws out the importance of the electron dynamic correlation to describe the PEHs and subsequent spectroscopic magnitudes for understanding the photoinduced processes occurring in thymine, and by extension other DNA or RNA nucleobases.

Results
The results are presented in three sections. First, calculations making use of the CASPT2//CASSCF protocol are discussed, focusing on the most important and initially accessed 1 (ππ*) state. In this protocol, the geometries are optimized with the CASSCF method and the energies are computed with the CASPT2 method at the optimized structures. Next, results from CASPT2//CASPT2 computations are given. In this case, the geometries are also determined with the highest level of theory. A comparison between both methodologies is then drawn. The results obtained are summarized and discussed together with previous theoretical and experimental studies present in the literature, isolating the concrete role played by dynamic correlation in the PEH topology of the different photo-activated states and establishing its influence on our subsequent interpretation of the experimental evidence.

CASPT2//CASSCF
Minimum energy paths (MEPs) and unconstrained geometry optimizations are employed to map the PEH of the 1 (ππ*) state in thymine at the CASSCF level and for different one-electron basis sets (see Computational Details), being summarized in Table 1. As can be seen, different behaviors might be inferred depending on the electronic correlation added through the basis set and active space employed. The smallest basis set, 6-31G, yields a planar minimum in the PEH of the 1 (ππ*) excited state in all cases, whereas the larger and more flexible 6-31G* gives rise to the 1 (ππ*/S0)CI crossing when either the optimization thresholds are enhanced or the MEP algorithm is employed.

Results
The results are presented in three sections. First, calculations making use of the CASPT2//CASSCF protocol are discussed, focusing on the most important and initially accessed 1 (ππ*) state. In this protocol, the geometries are optimized with the CASSCF method and the energies are computed with the CASPT2 method at the optimized structures. Next, results from CASPT2//CASPT2 computations are given. In this case, the geometries are also determined with the highest level of theory. A comparison between both methodologies is then drawn. The results obtained are summarized and discussed together with previous theoretical and experimental studies present in the literature, isolating the concrete role played by dynamic correlation in the PEH topology of the different photo-activated states and establishing its influence on our subsequent interpretation of the experimental evidence.

CASPT2//CASSCF
Minimum energy paths (MEPs) and unconstrained geometry optimizations are employed to map the PEH of the 1 (ππ*) state in thymine at the CASSCF level and for different one-electron basis sets (see Computational Details), being summarized in Table 1. As can be seen, different behaviors might be inferred depending on the electronic correlation added through the basis set and active space employed. The smallest basis set, 6-31G, yields a planar minimum in the PEH of the 1 (ππ*) excited state in all cases, whereas the larger and more flexible 6-31G* gives rise to the 1 (ππ*/S 0 ) CI crossing when either the optimization thresholds are enhanced or the MEP algorithm is employed.  Opt 1 (ππ*) min 1 (ππ*/S 0 ) CI Opt Tight 1 1 (ππ*/S 0 ) CI 1 (ππ*/S 0 ) CI MEP 1 (ππ*/S 0 ) CI 1 (ππ*/S 0 ) CI 1 Convergence thresholds of 10 −8 and 10 −6 for the energy change and the norm of the gradient, respectively. Figure 2 shows the computed MEPs by using different basis sets, increasing the amount of electron correlation considered. As can be seen, the MEP computed with the 6-31G basis set leads to the planar minimum 1 (ππ*) min in the 1 (ππ*) hypersurface (see Figure 1), as previously reported by Perun et al. [22]. Using the same active space but increasing the flexibility of the one-electron basis set employed yields a very different result though (cf. Figure 2). In this case, the system decays in a barrierless manner towards the ring-puckering ethene-like crossing 1 (ππ*/S 0 ) CI . The same result holds true after adding extra correlating orbitals to the active space [CASSCF (10,11)] (see Figure S1), thus validating the previous result obtained with the 6-31G* basis set and the CASSCF(10,8) methodology. Active space selection and validation becomes essential as the size of the system under study increases and computations become impractical [64], and has also been documented to affect the outcome in photoinduced dynamics [71]. This illustrates the importance of the electronic correlation in order to give a proper description of the photochemical pathways actually accessible in a given molecular system [72]. Moreover, the technique employed to optimize the excited state seems to also play a role [8]. Geometry optimization of the 1 (ππ*) state with the 6-31G* basis set yields the planar minimum 1 (ππ*) min (cf. Table 1), identical to the one found with the CASSCF/6-31G MEP, while the MEP at the CASSCF/6-31G* decays to 1 (ππ*/S 0 ) CI as shown in Figure 2. In other words, for thymine, determination of the different stationary points independently (geometry optimization for the ground state and the excited state using standard convergence thresholds, and CI search) at the CASSCF/6-31G* level, and connecting those regions of the hypersurface subsequently by linear interpolation of internal coordinates (LIIC), offers a view of the photochemistry through a well-defined minimum on the excited state hypersurface. A similar behavior is observed when the geometry optimizations are carried out at the CASPT2 level, with the sole difference that the minimum found along the hypersurface can be slightly twisted or "boat-like" 1 (ππ*) min-twisted (cf. Figure 1) [17,34,63]. The conclusions derived from those studies seem to be dependent on the strategy used, highlighting the additional insights provided by more sophisticated MEP techniques to map PEHs, while stressing the need for computing experimental observables derived from accurate PEHs to be able to unambiguously compare face-to-face with the experiment [17,73].
To exert a more accurate treatment of the PEH, several ANO-type basis sets have been employed yielding qualitatively analogous results (see Figure 2 and Supplementary Materials), pointing to the direct decay towards the ground state as the main pathway to be followed upon light irradiation. Upon inspection of Table 1, it can be seen that, once extra polarization functions are added on the C,N,O atoms, as in 6-31G*, ANO-S 321/21, ANO-L 321/21, and ANO-L 432/21 basis sets, the resulting structures after geometry optimization/MEP with CASSCF(10,11) converge to the same 1 (ππ*/S 0 ) CI structure with the exception of the CASSCF(10,11)/ANO-S 321/21 optimization, in which a planar minimum is obtained. A steady increase of electron correlation seems to converge towards the CI. The effect of the convergence thresholds on the PEHs has also been analyzed. While normal thresholds lead to minima in the rather flat 1 (ππ*) PEH, the decrease of those values lead towards the CI, its resulting structures are in agreement with all other 1 (ππ*/S 0 ) CI structures found and are described below. Careful analysis shows a rather planar 1 (ππ*) PEH that might explain the different results yielded by different thresholds and optimization procedures. Due to the planarity of the surface, and depending on the threshold employed during the optimization procedure, the presence of a given minimum may be warranted or not. This particularly flat topology of the PEH is not intrinsic to CASSCF/CASPT2, having been reported with other methods [24,73], and is in agreement with the broad emission band recorded and with its small fluorescence quantum yield [62].
The geometrical parameters defining the structures optimized at the different levels of theory are provided in the SI (Table S46). All results obtained employing basis sets with polarization functions are in agreement with those obtained by using the largest ANO-L 432/21. Deviations are small, as confirmed by the root mean square (RMS) averaged over all bond distances. In contrast, the smaller 6-31G basis set shows larger discrepancies. The C5-C6 and C5-C7 bond lengths and the dihedral describing the out of plane motion are the main geometrical distortions along the pyramidalization process towards the ring-puckering 1 (ππ*/S 0 ) CI , which strongly resembles those found for ethene [65]. The differences obtained for the geometry-optimized and MEP-resulting structures at the CASSCF (10,8) and CASSCF (10,11) levels are also shown in the SI (Table S47). Similar trends are found, supporting the fact that a steady increase of the correlation in the model leads directly towards the conical intersection.
In general, no qualitative changes on the topography of the PEH along the MEPs are discernible by going from the 6-31G* to ANO-L 432/21 (cf. Figure 2, Figures S1,S2), even though the vertical transition energies at the Franck-Condon region are overestimated on the former as compared to the results obtained with the ANO-type basis sets [49]. The computed values are 5.18, 4.83, 4.83, and 4.72 eV, with the 6-31G*, ANO-S 321/21, ANO-L 321/21, and ANO-L 432/21 basis sets, respectively, in agreement with the experimental band maximum placed at around 4.8 eV [3] and more correlated coupled cluster approaches [74]. Other properties highly dependent on the basis set, like the dipole moment, do change more prominently in the Pople-type basis set even though the mapped PEH appears to be qualitatively similar. This might have some implications in external perturbations to the wave function such as those given for QM/MM schemes [75][76][77], where solvation is described as dipole-driven electrostatic interactions between the solvent and the solute, and thus a more accurate characterization of the dipole moment may be important.
Additional computations have been carried out to explore the correlation effects in the less involved and dark 1 (n O π*) state. CASSCF (14,10) has been employed in this case, adding both n O lone-pair orbitals as described in the Computational Details to exert an even treatment to that given to the 1 (ππ*) state, and only geometry optimizations have been carried out because its direct population upon absorption is discarded given its negligible oscillator strength. This active space has been shown to provide analogous results for the 1 (ππ*) state as those previously depicted. A minimum on the 1 (n O π*) PEH has been obtained, characterized by an elongated C4=O bond of~1.37 Å for all basis set considered, and a C6C5C4O dihedral of −165.5 • for the ANO-L 432/21 basis set, differing in 1-2 degrees depending on the basis set employed thus yielding qualitative akin results (see SI). Their overall RMS compared to ANO-L 432/21 is~0.004, showing how covalent states such as 1 (n O π*) show a lesser dependence on the dynamic correlation included in the model than that found for the 1 (ππ*) state (cf . Table S48) [63]. It can therefore be concluded that 6-31G* provides a reasonable description of the PEH and will be the basis set used to test the adequacy of the CASPT2//CASSCF protocol against CASPT2//CASPT2. Table 2 compiles the different results obtained by means of geometry optimizations/MEP calculations at the CASPT2//CASPT2 level of theory. As can be seen, the threshold in the optimization procedure plays an akin role to what has been shown at the CASSCF level, prompting towards 1 (ππ*/S 0 ) CI when a tighter threshold is employed. The computed MEPs show a somewhat different trend compared to their CASSCF counterparts, in which the size of the hypersphere (r; see definition in Section 4) [78] employed along the path seems to influence the fate of the calculation. This appears to be inherent to the planarity of the 1 (ππ*) state PEH as discussed above, and the hypersphere employed seems to be the deciding factor in this particular case, small r values being retained in the shallow minimum encountered along the flat PEH whereas larger r values skip the minimum and go straight to the CI. The results here strongly emphasize the planarity of the PEH from the immediate surroundings of the FC to the CI, being even flatter than at CASSCF, and agree with the low fluorescence quantum yield as well as the broad emission band registered [62], and with its ultrafast sub-ps deactivation reported experimentally [3,54,79]. Table 2. CASPT2(10,8)/6-31G* geometry-optimization and MEP geometries obtained for thymine. r denotes the hypersphere radius employed in the calculation (in a.u.) [78].

Type of Calculation Method Structure
Opt The CASPT2 structures obtained in this section and their main parameters are analyzed in more detail in Table 3. As can be seen, the dihedral angle provides a good estimate to assess the resulting structure, as 1 (ππ*) min 1 (ππ*) twisted-min , and 1 (ππ*/S 0 ) CI feature values of~0, 30, and 105 • -110 • , respectively, and the latter is similar to the~120 • angle obtained at CASSCF. Bond distances differ from CASSCF (cf . Table S47): C6-C5 bond distances are shorter for the planar minimum, while they appear to be slightly elongated for the CI geometry at the CASPT2 level. CASPT2 C5-C7 bond lengths, on the other hand, appear to be slightly shorter than the~1.55 Å computed at CASSCF for the 1 (ππ/S 0 ) CI . Interestingly, MS-CASPT2 results seem to be closer to CASSCF values, and this supports the CASSCF description of the CI given the necessity to employ multistate treatments to recover the right dimensionality of the CI [80]. Regarding the new minimum identified in this section, 1 (ππ*) twisted-min , C6-C5 and C5-C7 bond distances are~1.45 and~1.49 Å, respectively.  The CASPT2 method has additionally been benchmarked against MS-CASPT2 (see Table 2), a technique used previously to tackle the photochemical behavior of thymine [32,34,41]. The results here obtained are rather homogeneous, yielding similar estimates at both levels of theory. Regarding the MEP calculations, the MS-CASPT2 method appears to require slightly larger r values to reach the 1 (ππ*/S 0 ) CI , a result that does not appear in the geometry optimizations independently of the threshold employed. The 1 (ππ*/S 0 ) CI structure obtained at the MS-CASPT2 level of theory resembles the one resulting from CASSCF calculations, whereas the CASPT2 CI is characterized by shorter bond length distances and a slightly lower dihedral angle. The differences obtained between CASPT2 and MS-CASPT2 relate to those previously reported for the penta-2-4-dieniminium cation (PSB-3) [81], where slight deviations are shown between these two different zeroth-order Hamiltonians. These could be partially caused by overestimated off-diagonal terms in the MS-CASPT2 effective Hamiltonian, as their definition can be lacking in some cases where the states included in the multistate procedure are of strongly different character, giving rise to an unphysical overestimation. Larger active spaces are demonstrated to palliate this effect but cannot be presently considered in this work to properly validate the MS-CASPT2 result, so we expect future results employing larger active spaces to be placed in between the estimates here given for the CASPT2 and MS-CASPT2 methods [81]. More accurate estimates are expected from the extended multistate (XMS) approach of Granovsky [82] and its related implementations [83][84][85][86][87][88].
Representative MEPs yielding the different key structures on the 1 (ππ*) state hypersurface are depicted in Figure 3. No planar structure is obtained by employing the MEP strategy at the CASPT2 level. Instead, the presence of a twisted minimum along the decay path, previously reported in the literature [34,63], appears to be the main structure hampering the barrierless route towards the 1 (ππ*/S 0 ) CI . This relates to what has also been found for PSB-3, where the dynamically correlated CASPT2 displays a more stable minimum under the planar minimum described by CASSCF [63]. Lower values of r favor this situation in which thymine finds itself to be stranded in a shallow minimum. Two different behaviors can be distinguished by inspection of Figure 3. In CASPT2/r = 0.09 (see Figure S3) and MS-CASPT2/r = 0.18 the MEPs lead towards a minimum, which is different from the planar minimum reported above at the CASSCF level. CASPT2/r = 0.30, on the other hand, exhibits a barrierless pathway towards the CI with the ground state. This finding is also obtained at the MS-CASPT2/r = 0.30 level with a slight difference in the dihedral angle and the C6-C5 and C5-C7 bond distances at the CI geometry.
As in CASPT2//CASSCF, additional computations have been carried out to characterize the 1 (n O π*) minimum, employing CASSCF (14,10). The minimum has been computed at the CASPT2 and MS-CASPT2 levels employing the 6-31G* basis set. Shorter C4=O distances (~1.365 Å) are obtained as compared to those previously reported with the CASSCF method, and a nearly planar structure displaying a C6C5C4O dihedral of −178 • (see SI), contrary to CASSCF findings. It is worth noting that both CASPT2 and MS-CASPT2 converge to the same structure, the choice of zeroth-order Hamiltonian exerting a lesser influence in this case as opposed to what has been previously reported for 1 (ππ*), which can be related to the different nature of the states: zwitterionic (ππ*) and covalent (nπ*). As in CASPT2//CASSCF, additional computations have been carried out to characterize the 1 (nOπ*) minimum, employing CASSCF (14,10). The minimum has been computed at the CASPT2 and MS-CASPT2 levels employing the 6-31G* basis set. Shorter C4=O distances (~1.365 Å) are obtained as compared to those previously reported with the CASSCF method, and a nearly planar structure displaying a C6C5C4O dihedral of −178° (see SI), contrary to CASSCF findings. It is worth noting that both CASPT2 and MS-CASPT2 converge to the same structure, the choice of zeroth-order Hamiltonian exerting a lesser influence in this case as opposed to what has been previously reported for 1 (ππ*), which can be related to the different nature of the states: zwitterionic (ππ*) and covalent (nπ*).

Discussion
CASSCF and CASPT2/MS-CASPT2 results obtained in the previous sections are here compared to assess the accuracy of the usually employed CASSCF method for geometry optimizations versus the CASPT2 method.
As can be seen in Figure 4, which summarizes the results displayed in Figures 2 and 3 and Tables 1 and 2, the results of the MEP [89] calculations at the CASPT2//CASSCF/6-31G* and CASPT2//CASPT2/6-31G* levels of theory exhibit in this particular case similar potential energy profiles that validate the CASPT2//CASSCF approach as an adequate protocol to predict the ultrafast nature of the decay process. However, it has been shown that the extremely flat 1 (ππ*) PEH seems to have a certain dependence on the dynamic correlation in its description, featuring a different minimum and displaying distinct gradients along the PEH that are expected to yield more noticeable differences in non-adiabatic molecular dynamics simulations. In such simulations, either the CASSCF or the CASPT2 methods are used, but not a combination of both as in the static CASPT2//CASSCF protocol [89][90][91][92][93]. Whereas the mentioned small discrepancies are not expected to alter our vision of the photoinduced decay (ultrafast process) and the interpretation of the experimental lifetimes from a static standpoint, it is important to note that they may influence spectroscopic observables arising from the minima characterized, as magnitudes like excited state absorptions show a strong dependence on the molecular structure [94].
CASSCF and CASPT2/MS-CASPT2 results obtained in the previous sections are here compared to assess the accuracy of the usually employed CASSCF method for geometry optimizations versus the CASPT2 method.
As can be seen in Figure 4, which summarizes the results displayed in Figures 2 and 3 and Tables 1 and 2, the results of the MEP [89] calculations at the CASPT2//CASSCF/6-31G* and CASPT2//CASPT2/6-31G* levels of theory exhibit in this particular case similar potential energy profiles that validate the CASPT2//CASSCF approach as an adequate protocol to predict the ultrafast nature of the decay process. However, it has been shown that the extremely flat 1 (ππ*) PEH seems to have a certain dependence on the dynamic correlation in its description, featuring a different minimum and displaying distinct gradients along the PEH that are expected to yield more noticeable differences in non-adiabatic molecular dynamics simulations. In such simulations, either the CASSCF or the CASPT2 methods are used, but not a combination of both as in the static CASPT2//CASSCF protocol [89][90][91][92][93]. Whereas the mentioned small discrepancies are not expected to alter our vision of the photoinduced decay (ultrafast process) and the interpretation of the experimental lifetimes from a static standpoint, it is important to note that they may influence spectroscopic observables arising from the minima characterized, as magnitudes like excited state absorptions show a strong dependence on the molecular structure [94]. Some differences are encountered for the geometry of the 1 (nOπ*) minimum. Whereas CASSCF predicts a partially twisted structure, both CASPT2 and MS-CASPT2 yield completely planar geometries more similar to the 1 (nOπ*/ππ*)CI crossing mediating the 1 (nOπ*) non-adiabatic population (see SI). The structural similarities featured between the 1 (nOπ*) minimum and 1 (nOπ*/ππ*)CI, as well as their energetic proximity, according to the CASPT2//CASSCF and CASPT2//CASPT2 protocols, may promote a re-crossing back to the 1 (ππ*) manifold [25] that may explain the reduced 1 (nOπ*) yield registered experimentally [55]. This would drastically reduce the role of the 1 (nOπ*) state in the photo-process, in contrast to what has been reported in MRCIS and CASSCF non-diabatic molecular dynamics studies [26,93,95], where a vast fraction of the excited Some differences are encountered for the geometry of the 1 (n O π*) minimum. Whereas CASSCF predicts a partially twisted structure, both CASPT2 and MS-CASPT2 yield completely planar geometries more similar to the 1 (n O π*/ππ*) CI crossing mediating the 1 (n O π*) non-adiabatic population (see SI). The structural similarities featured between the 1 (n O π*) minimum and 1 (n O π*/ππ*) CI , as well as their energetic proximity, according to the CASPT2//CASSCF and CASPT2//CASPT2 protocols, may promote a re-crossing back to the 1 (ππ*) manifold [25] that may explain the reduced 1 (n O π*) yield registered experimentally [55]. This would drastically reduce the role of the 1 (n O π*) state in the photo-process, in contrast to what has been reported in MRCIS and CASSCF non-diabatic molecular dynamics studies [26,93,95], where a vast fraction of the excited state population ends up trapped in the 1 (n O π*) state in the sub-ps timescale. This incidence is related to the likely overestimated role of the 1 (n O π*) PEH in MRCIS/CASSCF, where this state is proposed to be the main responsible of the ultrafast deactivation, and points towards an uneven treatment of the 1 (n O π*) and 1 (ππ*) states, which are described quite differently with methods lacking dynamic correlation like CASSCF [63,81], and that has motivated studies in the literature to assess the importance of the dynamic correlation in the overall photochemical picture [19,72]. This is further assessed in Figure 5, where an energy diagram of the involvement of the 1 (n O π*) state at both CASSCF//CASSCF and CASPT2//CASPT2 levels is shown, together with the energy profile arisen from the CASPT2//CASSCF protocol. Figure 5 depicts pronounced differences. The main ones are: i) in the Franck-Condon (FC) region, a large energy gap between the bright 1 (ππ*) and dark 1 (n O π*) states appears with the CASSCF//CASSCF approach, while at the CASPT2//CASPT2 level, the states appear to be near degenerate, and ii) the twisted CASSCF 1 (n O π*) minimum appears to be much lower in energy with respect to 1 (n O π*/ππ*) CI than its planar CASPT2 counterpart. This artificial stabilization is corrected at the CASPT2//CASSCF level (cf. Figure 4). The former difference (i) will have important incidences in the early time events of molecular dynamics simulations employing the CASSCF method, where kinetic energy is built along the pathway from the initially accessed 1 (ππ*) state towards the 1 (n O π*/ππ*) CI that may artificially increase the population on the 1 (n O π*) state, whereas the latter discrepancy (ii) is expected to ensure the population trapping in the 1 (n O π*) state. On the other hand, at the CASPT2//CASPT2 level, both 1 (ππ*) and 1 (n O π*) states are degenerate at the FC region and both 1 (n O π*/ππ*) CI and the 1 (n O π*) minimum appear within half an eV, which is expected to increase the re-crossing back to the 1 (ππ*) manifold [25] and thus reduce the population ending in the 1 (n O π*) state [55].
state population ends up trapped in the 1 (nOπ*) state in the sub-ps timescale. This incidence is related to the likely overestimated role of the 1 (nOπ*) PEH in MRCIS/CASSCF, where this state is proposed to be the main responsible of the ultrafast deactivation, and points towards an uneven treatment of the 1 (nOπ*) and 1 (ππ*) states, which are described quite differently with methods lacking dynamic correlation like CASSCF [63,81], and that has motivated studies in the literature to assess the importance of the dynamic correlation in the overall photochemical picture [19,72]. This is further assessed in Figure 5, where an energy diagram of the involvement of the 1 (nOπ*) state at both CASSCF//CASSCF and CASPT2//CASPT2 levels is shown, together with the energy profile arisen from the CASPT2//CASSCF protocol. Figure 5 depicts pronounced differences. The main ones are: i) in the Franck-Condon (FC) region, a large energy gap between the bright 1 (ππ*) and dark 1 (nOπ*) states appears with the CASSCF//CASSCF approach, while at the CASPT2//CASPT2 level, the states appear to be near degenerate, and ii) the twisted CASSCF 1 (nOπ*) minimum appears to be much lower in energy with respect to 1 (nOπ*/ππ*)CI than its planar CASPT2 counterpart. This artificial stabilization is corrected at the CASPT2//CASSCF level (cf. Figure 4). The former difference (i) will have important incidences in the early time events of molecular dynamics simulations employing the CASSCF method, where kinetic energy is built along the pathway from the initially accessed 1 (ππ*) state towards the 1 (nOπ*/ππ*)CI that may artificially increase the population on the 1 (nOπ*) state, whereas the latter discrepancy (ii) is expected to ensure the population trapping in the 1 (nOπ*) state. On the other hand, at the CASPT2//CASPT2 level, both 1 (ππ*) and 1 (nOπ*) states are degenerate at the FC region and both 1 (nOπ*/ππ*)CI and the 1 (nOπ*) minimum appear within half an eV, which is expected to increase the re-crossing back to the 1 (ππ*) manifold [25] and thus reduce the population ending in the 1 (nOπ*) state [55].

Computational Details
CASSCF and CASPT2 methods have been used in this study as implemented in MOLCAS 7.6 [96,97]. Two different types of basis sets were used: i) the segmented 6-31G and 6-31G* basis sets, and ii) two generally contracted basis sets of atomic natural orbital (ANO) type with C,N,O(10s,6p,3d)/H(7s3p) [98]  Multiconfigurational wave functions have been built using in the CAS space the whole valence π occupied and virtual orbitals of thymine, giving rise to the CASSCF (10,8). Further calculations with three extra π virtual orbitals have been performed to assess the convergence of the active space for such computations, resulting in the CASSCF(10,11) level [11]. No symmetry restrictions were imposed to the molecule as required to avoid any symmetry breaking problems and to allow complete freedom during the optimization procedures. Additional computations including the n O lone pair orbitals in the active space were also carried out and labeled as CASSCF (14,10). Two roots were averaged for CASSCF (10,8) and CASSCF (10,11), whereas five states were included for CASSCF (14,10), accounting for both low-lying n O π* states. Second-order perturbative corrections have been computed on top of the CASSCF wave functions, maintaining all core electrons frozen during the perturbation step and making use of the zeroth-order Hamiltonian as originally implemented [43,48], and employing a 0.2 a.u. imaginary level shift to avoid weakly-interacting intruder states [100]. Additional computations at the multistate CASPT2 (MS-CASPT2) level have also been performed [101] for the sake of comparison with previous results present in the literature [34].
Geometry optimizations and minimum energy paths (MEPs) have been obtained at both CASSCF and CASPT2 levels. CASSCF MEPs have been computed with analytical gradients [102] and subsequent CASPT2 calculations have been performed on the converged geometries along the path in order to include dynamic correlation, as usually performed nowadays in photochemical studies following the CASPT2//CASSCF protocol [46]. CASPT2 optimizations and MEPs, on the other hand, have been computed, making use of numerical gradients. MEPs have been extensively used in this work, making use of mass weighted coordinates. This technique follows (if present) a steepest descent path, in which each step is built by the minimization of the energy on a hyperspherical cross section of the PEH of a predefined radius centered on the geometry optimized in the previous step, providing a powerful tool to study the photophysics and photochemistry of molecular systems [78]. Different optimization thresholds and hypersphere radii (r) for the MEPs have been computed to assess their influence in the PEH topography. For geometry optimizations, a standard convergence threshold of 10 −5 and 10 −3 (in a.u.) for the energy change and the norm of the gradient have been used, respectively. Additionally, a tighter threshold characterized by an energy change of 10 −8 and a norm of the gradient of 10 −6 has also been employed. On the other hand, several values for the hypersphere radius (r) defining the MEP have been computed, ranging from 0.09 to 0.3 (in a.u.). CIs are characterized at the CASSCF level with the GAUSSIAN 09 package [103] and employing the projection technique of Bearpark et al. [104], while CASPT2/MS-CASPT2 crossings are directly taken from the minimum energy crossing points found along the MEP.

Conclusions
We can conclude that the flat PEH of the 1 (ππ*) bright low-lying excited state of thymine, widely discussed in the present study, presents a couple of shallow minima on the deactivation path towards the 1 (ππ*/S 0 ) CI , which might account for the sub-ps to few-ps signals registered experimentally. The different minima exhibit small barriers, within the intrinsic error of the methods employed, and thus the decay is predicted to be ultrafast, the fluorescence quantum yield is expected to be small, and the emissive feature broad, in accordance with the experimental evidence. The 1 (n O π*) state, on the other hand, displays a slightly twisted minimum with the CASSCF method and a planar minimum with the CASPT2. In the latter case, the minimum of the 1 (n O π*) state is placed at the vicinities of the 1 (n O π*/ππ*) CI crossing close to the Franck-Condon region which might facilitate the re-crossing towards the 1 (ππ*) state and thus account for the reduced 1 (n O π*) yield registered experimentally. Overall, CASPT2//CASSCF and CASPT2//CASPT2 (or MS-CASPT2//MS-CASPT2) appear to yield on static grounds a similar qualitative picture of the photoinduced phenomena in the sub-ps timescale. However, this is expected to change in molecular dynamics simulations where the overly steep CASSCF surfaces/energy gradients have been reported to predict too short decay times compared to those recorded. The study supports the qualitative picture obtained on static grounds of the CASPT2//CASSCF protocol, given the energetic discrepancies between CASSCF and CASPT2 are carefully assessed while calling for dynamically correlated CASPT2 gradients for an even-footing description of the photoinduced events on the different states involved in on-the-fly non-adiabatic molecular dynamics simulations. Novel spectroscopic techniques with enhanced temporal and spectral resolution should help overcome these ambiguities over the coming years [105][106][107].
Supplementary Materials: The following are available online at www.mdpi.com/1420-3049/21/12/1666/s1: Cartesian coordinates of the minima and minimum energy crossing points obtained within the study, absolute energies of relevant structures, and additional MEP profiles are given.