Organic Emitters Showing Excited-States Energy Inversion: An Assessment of MC-PDFT and Correlation Energy Functionals Beyond TD-DFT

: The lowest-energy singlet ( S 1 ) and triplet ( T 1 ) excited states of organic conjugated chromophores are known to be accurately calculated by modern wavefunction and Time-Dependent Density Functional Theory (TD-DFT) methods, with the accuracy of the latter heavily relying on the exchange-correlation functional employed. However, there are challenging cases for which this cannot be the case, due to the fact that those excited states are not exclusively formed by single excitations and/or are affected by marked correlation effects, and thus TD-DFT might fall short. We will tackle here a set of molecules belonging to the azaphenalene family, for which research recently documented an inversion of the relative energy of S 1 and T 1 excited states giving rise to a negative energy difference ( ∆ E ST ) between them and, thereby, contrary to most of the systems thus far treated by TD-DFT. Since methods going beyond standard TD-DFT are not extensively applied to excited-state calculations and considering how challenging this case is for the molecules investigated, we will prospectively employ here a set of non-standard methods (Multi-Conﬁgurational Pair Density Functional Theory or MC-PDFT) and correlation functionals (i.e., Lie–Clementi and Colle–Salvetti) relying not only on the electronic density but also on some modiﬁcations considering the intricate electronic structure of these systems.


Introduction
The violation of Hund's rule in molecules [1], analogously to atoms, is commonly ascribed to an inversion of the excitation energies of the lowest states of spin-singlet (S 1 ) or spin-triplet (T 1 ) multiplicity. In common situations, the energy difference between S 1 and T 1 excited states, that is ∆E ST = E(S 1 ) − E(T 1 ) is positive, contrary to what happens if Hund's rule is altered (in that case, ∆E ST would be negative). Note that the negative sign contradicts the fact that the exchange energy (K) is normally thought to be positive, historically ∆E ST ≈ 2K, thus, implying that the lowest singlet excited state lies energetically above the lowest triplet excited state and not the opposite.
That exchange energy is known to be of the order of hundreds of meV for common organic chromophores, and this becomes a key parameter for photophysics and related applications [2][3][4]. However, it has been demonstrated that strong correlation effects can decrease the ∆E ST value [5,6] and even invert the energies of S 1 and T 1 excited states due to a more marked stabilization of the former vs. the latter state [7], although very few molecules are discovered up to now displaying such an excited-state energy inversion. Additionally, fast environmental effects (whenever they are reliably introduced) could also lead to negative ∆E ST values.
Among those environmental effects, we mention thermal fluctuations of molecular conformations or microscopic electronic polarization effects in amorphous films of carbazole derivatives [8], thus, opening a whole world of future studies and applications around this unexpected issue. However, for solvation effects, current implementations of continuum solvation models should be employed with caution since it could lead to spurious excitedstate energy inversion [9].
The physical origin of this inversion and its possible practical implications should not merely be considered as an academic questions. As an example of practical use, due to the spin statistics, triplet excitons (dark) are known to be formed upon a 3:1 ratio with respect to singlet excitons (bright) thus limiting the efficiency of electroluminiscent processes. Therefore, the mentioned energy inversion could be further exploited in the, e.g., recovery of triplet excitons created upon electroluminiscence to increase the internal quantum efficiency or quantum yields of Organic Light-Emitting Diodes (OLEDs).
Many other applications in related fields (photocatalysts, covalent organic frameworks, liquid crystals, etc.) have also been envisioned and recently reviewed [10]. However, from an experimental point of view, the range of disclosed molecules showing that excitedstate energy inversion is very limited, and this goes back to the discovery in the 1980s of some azaphenalene molecules candidates [11,12] for such a violation of Hund's rule. Those initial molecules were also extended to other (not-yet-synthesized) candidates after a massive screening of compounds recently performed [13], showing that the topic is still open and under active investigation.
Therefore, theoretical methods based on one-electron excitations (i.e., Time-Dependent Density Functional Theory or TD-DFT) are questioned in its current implementation to recognize that excited-state energy inversion, due to the lack of inclusion of higher-thansingle excitations into their formulation. In contrast, wavefunction methods have been shown in recent studies [14][15][16][17][18][19][20][21] to provide reasonably accurate values for that ∆E ST energy difference, although at a cost certainly higher than TD-DFT.
That limitation of TD-DFT is known to occur independently of the underlying exchangecorrelation functional and basis sets chosen. However, since excited-state wavefunction methods can capture double and higher-order excitations by definition, depending of the truncation done for the excitation operator, those methods are able to predict the excitedstate energy inversion while concomitantly providing accurate individual (i.e., S 1 and T 1 ) energies for the excited-states involved.
Based on these findings, our goal here is to investigate if methods going beyond standard (TD-)DFT could predict that excited-state energy inversion and thus compete in accuracy with wavefunction results. To assess that, we will employ methods merging wavefunction and correlation functionals, in the hope of including both kind of correlation effects (dynamical and non-dynamical) for any of the electronic states involved. These results will also serve to confirm the key role played by marked correlation effects, as well as to invigorate more research of DFT methods out of the most commonly found implementations.

Choice of the Target Systems
The set of systems selected is exclusively motivated by the previous experimental and theoretical findings mentioned above. The set of azaphenalene molecules shown in Figure 1 is known to display negative ∆E ST values at various wavefunction levels, from pioneering studies [14][15][16] later extended to related and/or larger systems [17][18][19][20][21], which clearly constitutes a challenge for any theoretical method.
Interestingly, chemical substitution of the heptazine core -C 6 H 7 H 3 -with chlorine -C 6 H 7 Cl 3 -, cyano -C 6 H 7 (CN) 3 -, or p-methoxyphenylene -C 6 N 7 (p-C 6 H 4 OCH 3 ) 3 -groups preserved the negative ∆E ST value, showing that the chemical structure of the core is indeed responsible as well as that more potential and synthetically viable molecules could soon be theoretically disclosed and/or experimentally achieved. From the experimental point of view, some studies that were performed in the 1980s also concluded with a real possibility of having an excited-state inversion for some of these azaphenalenes [11,12]. More recently, some heptazine derivatives have been successfully employed as emitters in real devices [22,23], with an exceptionally high quantum yield reported and explained by the (partial) conversion of triplet into singlet excitons possibly assisted by a negative ∆E ST value.
Note also that other related cores could be also potential candidates, after the conclusions reached by a massive screening of thousands of potential azaphenalene candidates by Pollice et al. [13]. However, we are more interested in assessing the reliability of theories going beyond (TD-)DFT and we will, thus, restrict this work to the compounds shown in Figure 1 for which reference results are available as well.

Physical Meaning of Reduced Density Matrices
While the electronic density, or first-order reduced density matrix, is given by: and thus integrates over the number of electrons N, ρ(r)dr = N, the corresponding spinless second-order reduced density matrix integrates to the total number of interacting electron pairs: and represents the probability density of finding a particle at point r 1 and simultaneously another particle at point r 2 . The explicit form is given by: or better its reduced form ξ 2 (r 1 , r 2 ; r 1 , r 2 ) = ∑ s 1 ,s 2 = γ 2 (x 1 , x 2 ; x 1 , x 2 )| s i =s i . Finally, the diagonal element or ρ 2 (r 1 , r 2 ) would account from any correlation effect arising from interparticle interaction, as it can be easily seen from the electron-electron mean value as a function of this new magnitude: Note also that ρ(r) and ρ 2 (r 1 , r 2 ) are related through: In the following, we will denote ρ 2 (r) = ρ 2 (r 1 , r 2 )| r 1 =r 2 as the function at the twoelectron coalescence point, whose modelling has been extensively pursued in the past [24][25][26], as well as its integration into excited-state formalisms [27,28], as the next step for the description of electronic structure beyond the use of merely the electronic density ρ(r).

Theories Going beyond (TD-)DFT
The methods included in this study can be categorized in popular language as methods going beyond (standard) DFT, in the sense that they are based on the on-top second-order reduced density matrix ρ 2 (r), and not only on the first-order density ρ(r), or post-MCSCF, in the sense that a Multi-Configurational Self-Consistent Field (MCSCF) calculations needs to be done first from which the magnitude ρ 2 (r) is obtained.
Note that the on-top second-order reduced density matrix represents the probability that two opposite-spin electrons are found at point r and integrates to the total number of interacting pairs. Multiconfiguration Pair-Density Functional Theory (MC-PDFT [29][30][31][32]) can be thus viewed as a post-MCSCF method that evaluates the energy of any state with on-top pair-density function theory. Basically, for a MCSCF wavefunction, |Ψ MCSCF = ∑ µ C µ |Ψ µ , one can obtain the total electronic energy as: (6) with the different terms being the kinetic, potential, Coulomb, and exchange-correlation energies, respectively, the latter relying on a modification of a common DFT functional to be employed together with Equation (1). The most extensively tested on-top density functional is called tPBE and will, thus, be used here consequently.
On the other hand, one can also directly employ a functional explicitly depending upon the design of the on-top second-order reduced density matrix, such as the Colle-Salvetti (CS [33][34][35]) correlation functional. We will denote, in the following, the CS expression as a two-body correlation functional, cast as E c [ρ(r), ρ 2 (r)], in contrast with conventional or one-body functionals commonly used for standard DFT calculations, or simply E c [ρ(r)]. Note that the famous Lee-Yang-Parr (LYP [36]) correlation functional is a reformulation of the Colle-Salvetti expression to avoid the explicit dependence on ρ 2 (r) at the price of neglecting its use with, e.g., MCSCF wavefunctions. For this case, the total electronic energy is calculated by a two-step procedure, after adding the correlation (mainly dynamic) energy to the energy calculated by the underlying MCSCF procedure, which already includes the non-dynamic (or static) correlation energy, withT,V Ne ,V ee the kinetic, nuclear-electron, and electron-electron operators.
Another not-so-common approximation is given by the Lie-Clementi (LC [37,38]) correlation functional, with an explicit dependence on the natural (fractional) occupation numbers, if a MCSCF calculation is done first. A modified density is built, such as depending on the density built from the natural orbitals,ρ(r), and the corresponding natural orbital occupation numbers (n i ). That density is, thus, inserted into the reparameterized correlation functional of Gombas et al., generally denoting this class of functionals as E c [ρ m (r)]. Interestingly, those orbitals with n i < 2 do not contribute to the correlation energy as much as those doubly occupied (n i = 2), thus, describing both ground-and excited-states independently. Therefore, MC-PDFT, E xc [ρ(r), ρ 2 (r)], and E c [ρ m (r)] exchange-correlation functionals will be all based here on a Complete Active Space Self-Consistent Field (CASSCF) wavefunction with an active space of N electrons housed in M orbitals, or simply (N, M), to incorporate non-dynamic (or static) correlation effects in a consistent way. The active space chosen, (6,6) or (12,12), is indeed based on the occupancy (and degeneracy) or molecular orbitals found at the uncorrelated level.
Note that: (i) MC-PDFT could be instead used with many MCSCF wavefunction types, such as GVB, CASSCF, RASSCF, CAS-CI, and RAS-CI, as is the case for the other two-body functionals too. We will, however, limit this work to the same wavefunction type for both schemes for the sake of coherence. (ii) Since these theories can be applied to any state of interest, independently of its spin, there is no need to invoke a linear-response regime as it happens for TD-DFT. (iii) Any of these methods will incorporate all correlation effects, be they static or dynamic, thus, possibly disentangling the importance of any of these contributions into the final results.

Computational Details
The ground-state (S 0 ) geometry of all the compounds was optimized by the B97-3c method [39], without any imaginary frequency obtained. The energy difference between the lowest-energy spin-singlet (S 1 ) and spin-triplet (T 1 ) excited-states is denoted as ∆E ST , which is normally positive unless for an energy inversion of the S 1 and T 1 energies, thus, giving rise to ∆E ST < 0. The def2-SVP and def2-TZVP basis sets [40] are used for all the calculations, with the auxiliary def2/JK and def2-TZVP/C basis sets [41] to reduce the computational cost.
For some control TD-DFT calculations, we will employ the ωB97, ωB97X [42], and ωB97X-2 [43], which form a set of range-separated exchange-correlation functionals belonging, respectively, to non-hybrid (ωB97), hybrid (ωB97X), and double-hybrid (ωB97X-2) rungs to infer if the addition of exact-exchange (for a hybrid) or pertubation-like (for a doublehybrid) brings any difference to the results.
We used the following quantum-chemical packages for the calculations performed here: ORCA 5.0 [44] for the (standard) TD-DFT with hybrid and double-hybrid functionals, GAMESS [45] for the MC-PDFT method with the tPBE functional, and an in-house program [46,47] (interfaced with GAMESS) for the E c [ρ m (r)] and E c [ρ(r), ρ 2 (r)] calculations employing the Lie-Clementi and the Colle-Salvetti correlation functionals.

Reference Results Available
We compare, in the following, the results obtained here not only with previously applied wavefunction methods but also with respect to the experimental information available: Leupin et al. obtained [11] for MAP a S 1 ← S 0 and T 1 ← S 0 excitation energies of 0.972 and between 0.972 and 0.984, respectively, and thus with a ∆E ST energy difference possibly negative. For 4AP, Leupin et al. obtained [12] excitation energies below 2.39 for S 1 ← S 0 (although the final value could not be completely determined experimentally) and 2.29 eV for T 1 ← S 0 , respectively, again not excluding a ∆E ST < 0 value depending on how low S 1 ← S 0 was in reality.
In all cases, all these correlated wavefunction methods were able to predict a ∆E ST < 0 value for all the MAP, 4AP, 5AP, and 7AP systems. Choosing one of these methods as reference, necessarily motivated by the completeness of the values found in the literature [17] together with the high accuracy demanded, and the SCS-CC2/def2-TZVP results are listed next for cross-comparison along the study.

TD-DFT Calculations
First, we also illustrate here if a standard TD-DFT calculation is able or not to provide a ∆E ST < 0 value and (if so) with what accuracy. For this purpose, we will choose not a selection of different functionals but a set of models of increasing complexity (i.e., semilocal, hybrid, and double-hybrid functionals). Additionally, since range-separation was before invoked as an accurate tool for the modelling of organic emitters of this type [49], Table 1 gathers the results obtained by the ωB97, ωB97X, and ωB97X-2 range-separated exchange-correlation functionals.
Strikingly, ωB97 and ωB97X are unable to predict the excited-state energy inversion, providing values for individual excitation energies not differing greatly between both methods, which is not often the case since those values are known to depend on the amount of (short-range) exact-exchange introduced (0 % in ωB97 vs. ≈16 % in ωB97X) and/or the range-separation parameter (ω = 0.4 in ωB97 vs. ω = 0.3 in ωB97X).
Furthermore, we note how the S 1 ← S 0 excitation energies are overestimated (severely in the case of 7AP) by both methods with respect to experimental or SCS-CC2 reference results. Note that many previous publications [14][15][16][17][18][19][20][21] already demonstrated how this was the case for all exchange-correlation functionals assessed up to now, which is also confirmed here.
We additionally assessed if the use of the ωB97X-D3 or ωB97X-D4 models, which are reparameterized to include the dispersion correction, would lead to any difference: the ∆E ST value only changed only 0.01 eV with respect to the value calculated by the original ωB97X.
We next analysed the performance of the ωB97X-2 double-hybrid functional, which, in its extension to excited states [50,51], includes a contribution from double excitations, thus, going beyond the single excitations introduced by TD-DFT routinely. We can observe how this method actually predicts ∆E ST < 0 values between −0.4 and −0.7 eV, roughly speaking that are, thus, too large with respect to the reference results. An additional concern arises from the inspection of individual S 1 ← S 0 and T 1 ← S 0 excitation energies, since the method seems to severely underestimate (by up to 1 eV) the former values progressively as a function of the N atoms introduced into the chemical structure.
The overestimation of the T 1 ← S 0 excitation energies is slightly attenuated with respect to the S 1 ← S 0 ones but again with values deviating too much with respect to the reference results. Shortly speaking, although this method is able to provide negative values for ∆E ST , the results appear to be affected by a systematic error. The use of double-hybrid functionals has been recently and more systematically examined [52], with some of the models assessed being promising enough to display accurate individual excitation energies: we thus refer the reader to that study for further information and confirmation about the key role played by double excitations into the final values.

MC-PDFT Calculations
We will inspect the CASSCF results shown in Table 2 to first observe the effect of using both basis sets, def2-SVP and def2-TZVP, for this set of calculations. For the CASSCF(6,6) results, going from def2-SVP to def2-TZVP implies a slight increase of the S 1 ← S 0 and T 1 ← S 0 excitation energies, with the exception of the latter for the 7AP molecule, but asymmetrically, with the corresponding ∆E ST values altered significantly. For the CASSCF (12,12) case, the variations for 7AP are also significant and deviate from the SCS-CC2/def2-TZVP reference values.
For MAP, the CASSCF(6,6)/def2-TZVP calculation already provided close results to reference SCS-CC2/def2-TZVP results, with S 1 ← S 0 and T 1 ← S 0 excitation energies differing by 0.15 and 0.10 eV, respectively, and thus leading to a negative ∆E ST value of −0.18 compared to −0.22 eV as reference. However, a larger active space is not definitively giving any advantage here, as it was also found before [17], stabilizing too much the S 1 (T 1 ) state and leading consequently to an overly negative (positive) ∆E ST value with the def2-SVP (def2-TZVP) basis set. For 4AP, the agreement is not so close to CASSCF (6,6), with both excitation energies largely overestimated (by 0.3-0.7 eV) independently of the basis set chosen. A larger active space, CASSCF (12,12), seems beneficial only with the def2-SVP basis set, which appears to indicate a not so balanced treatment of correlation effects in the absence of a dynamical correlation correction.
This overestimation was also found for 5AP and 7AP, particularly striking for the latter and again independently of the active space fixed. Overall, it seems that the CASSCF results do not suffice to lead to accurate and robust results by themselves, although negative ∆E ST values are mostly obtained. Previous publications show how the need of dynamic correlation effects (i.e., NEVPT2) as a further step to obtain more accurate and robust results [17,18]. The use of the tPBE correlation functional together with Equation (7) is presented next in Table 3, again for both basis sets (def2-SVP and def2-TZVP) and both active spaces of the underlying CASSCF calculation. For MAP, the CASSCF(6,6) + tPBE results are considerably accurate with both basis sets, not only for the target ∆E ST energy difference but also for the individual excitation energies. The use of the larger CASSCF (12,12) active space instead increases both excitation energies, especially the former, and reverse the sign of ∆E ST .
For 4AP and 5AP, the CASSCF(12,12) + tPBE results are very accurate with the def2-TZVP basis set and with respect to the SCS-CC2 reference values. For 7AP, the CASSCF(12,12) + tPBE results are also relatively accurate, leading to a negative ∆E ST value with both basis sets and correcting the overestimation of values found at the CASSCF (12,12) level with the def2-TZVP basis set. It thus appears that the addition of a (modified) correlation functional is qualitatively beneficial; however, more research is still needed to confirm the application of MC-PDFT to other chromophores and related systems. Table 3. Excited-state energies and associated ∆E ST energy difference (all in eV) calculated with MC-PDFT.

Basis Set
Molecule

Lie-Clementi (LC) and Colle-Salvetti (CS) Calculations
These two functionals are here applied with the CASSCF(6,6) active space and the def2-SVP basis set to avoid a known (and long-standing) problem with these methods, i.e., the double counting of the dynamical correlation energy, which might be minimized using the smallest admissible active space in the underlying CASSCF calculations [53,54]. Table 4 presents these pioneering results, from which interesting features can be observed: (i) a negative ∆E ST is provided in all cases, contrarily to some of the former methods and the isolated CASSCF(6,6) calculations of Table 2, which clearly underlines the major role played by dynamical correlation effects for describing at least quantitatively these states; (ii) the cost-effective Lie-Clementi functional overestimates the individual excitation energies, systematically increasing them with respect to the CASSCF(6,6) values; and (iii) the Colle-Salvetti functional, relying on the introduction of the ρ 2 (r) variable, provides closer values to reference results, again with the exception of 7AP for which these are still overestimated.
Given the variety of methods tested, Figure 2 reports the calculated ∆E ST energy difference for a better comparison. Contrary to previous TD-DFT applications, and with very few exceptions, the figure clearly shows how it is possible to obtain negative values (i.e., inverted S 1 and T 1 excitation energies) with these methods thanks to the combination of non-dynamical and dynamical correlation effects. Taking into account the use here of the def2-SVP basis set, we can next compare these results with the high-quality values available in literature: (i) For MAP, NEVPT2(6,6)/def2-SVP results [18] gave S 1 ← S 0 and T 1 ← S 0 excitation energies of 1.102 and 1.319 eV with the CASSCF(6,6) + CS results differing by only 0.13 and 0.05 eV, respectively. (ii) For 5AP, EOM-CCSD/cc-pVDZ results [13] are 2.251 and 2.329 eV for S 1 ← S 0 and T 1 ← S 0 excitation energies, respectively, with the CASSCF(6,6) + tPBE results of 2.153 and 2.437 eV differing by only 0.10 and 0.11 eV, respectively. (iii) For 7AP, DLPNO-NEVPT2(6,6)/def2-SVP results [13] led to S 1 ← S 0 and T 1 ← S 0 excitation energies of 2.552 and 2.906 eV with the CASSCF(12,12) + tPBE results being the closest ones but still differing by 0.30 and 0.47 eV, respectively.

Conclusions
The field of (TD-)DFT has so impressively advanced over the recent decades thanks to the two-fold and concurring efforts of continuously merging developments and applications. As a corollary, the latter would have not been possible without major advances from developments, often considered part of basic but completely needed Science. In this regard, new methods have historically been fostered by providing cost-effective yet accurate expressions and implementations for wide applications or, on the other hand, by tackling cutting-edge applications at the frontier of knowledge to move the field forward. In other words, inaccurate results are often needed to question why (TD-)DFT behaves as it does and how it can be rigorously and systematically improved.
Therefore, we attempted to continue building that interface between both worlds (developments and applications) by selecting a long-standing chemical problem of revisited interest: the energy inversion of the lowest spin-singlet and spin-triplet excited states of azaphenalene compounds intended to be used as organic emitters or photocatalysts. For that purpose, knowing from previous works than TD-DFT was not a reliable path to adequately address this issue, we applied methods employing not only the electronic density as the sole ingredient of an exchange-correlation functional but also other more involved magnitudes, e.g., the on-top pair density.
Overall, the use of the latter into the MC-PDFT scheme or as part of the explicit formulation of the Colle-Salvetti correlation functional offers an attractive way to overcome the limitations found for (TD-)DFT, although at a higher computational cost. However, we are also aware that further research is needed to benchmark these non-conventional methods as well as to reduce their computational cost and scaling with the system size, for which more challenging applications will be also welcome in the near future.

Concluding Remarks: A Personal Note
We would like to contribute with this article, as part of the Special Issue in honour of Professor Karlheinz Schwarz on the occasion of his 80th birthday, to celebrate the outstanding role played by Professor Schwarz in the field of Density Functional Theory (DFT) through his scientific career [55]. The authors met him as part of the International Scientific Committee of the "International Conference on Density-Functional Theory and its Applications", which is likely one of the longest-lived events to exist in the fields of theoretical and computational Chemistry and Physics after Paris Thanks to the strong activity of Heinz promoting DFT worldwide, we had the opportunity to enjoy, as local organizers of the last edition (see Figure 3), his compromise and illusion with this series of conferences. This work is, thus, our small recognition to his talented, vibrant, active, and kind figure, promoting DFT worldwide alongside an outstanding scientific career.