1. Introduction

Metal-free, thermally activated delayed fluorescence (TADF) materials are characterized by luminescent molecules or complexes with small energy gap (ΔE_ST) between the lowest singlet (S₁) and triplet (T₁) excited states ^[1-4]. They can upconvert from the T₁ to S₁ by absorbing environmental thermal energy and then radiatively decay from the S₁. TADF emitters can be employed in OLEDs to convert both singlet and triplet excitons into light with a theoretical yield up to 100% ^[5-17]. Because OLEDs containing TADF emitters can be up to four times more efficient than OLEDs containing fluorescent molecules, and cheaper than OLEDs containing noble metal-based phosphor, they have garnered considerable research attention in recent years ^[18-28].

Small ΔE_ST value is the key to realizing an efficient TADF. We now know that two conditions must be met to ensure a small ΔE_ST: The first is a small exchange integral for the S₁ transition, which is roughly proportional to the degree of overlap between the initial and final state orbitals and is generally realized by an intermolecular or intramolecular charge-transfer (CT) ^[5-13]; the second is that the locally excited triplet state (³LE) be close to (or even exceed) the ³CT (Fig. 1) ^{[15, 18, 21]}. If theoretical calculations can accurately predict the transition energies of CT molecules, the time and cost of developing new TADF materials can be substantially reduced.

Owing to its ability to balance accuracy and efficiency, the linear response time-dependent density functional theory (LR-TDDFT), which was developed based on density functional theory, has become the most widely used method for theoretical interpretation of excited-state electron structure and dynamics for medium to large molecular systems ^[29-31]. The resolution of the LR-TDDFT equations performed at optimized ground- and excited-state geometries can reproduce the vertical absorption and emission energies, as well as the adiabatic energy (approximate to the zero-zero transition energy) under the Franck-Condon principle (Fig. 2) ^[32-34]. However, the TDDFT results for CT states are highly dependent on the selected exchange-correlation (XC) - especially the Hartree-Fock (HF) percentage contained therein - when hybrid-functionals are used ^{[31, 35-37]}. Consequently, the transition energies of different molecules must be calculated by using different XC-functionals. Moreover, when the singlet-triplet pairs of low-lying excited states are crossed in a TADF molecule, as shown in Fig. 1a, the transitions of S₁ and T₁ involving different orbitals must also be reproduced by different XC-functionals ^[37]. The challenge for theoretical chemists, to this effect, is to secure a systemic and universal strategy for the selection or correction of XC-functionals to facilitate accurate prediction of the transition energies of various TADF (and even CT) systems.

This brief review begins by introducing the theoretical and experimental error sources in the studies of transition energies of CT-type TADF molecules, then goes on to summarize current efforts to overcome these problems.

2. Error sources in theoretical prediction of TADF molecule transition energy 2.1. Drawbacks to existing theoretical models

The TDDFT scheme derived from linear response formalism is necessary for solving any non-Hermitian eigenvalue equation ^{[31, 35, 36]}:

$\left[ {\matrix{ A & B \cr {{B^ * }} & {{A^ * }} \cr } } \right]\left[ {\matrix{ X \cr Y \cr } } \right] = \omega \left[ {\matrix{ 1 & 0 \cr 0 & { - 1} \cr } } \right]\left[ {\matrix{ X \cr Y \cr } } \right]$

(1)

where A and B are the Hamiltonian matrix elements of the single (de)excitation determinants, X and Y are the single (de)excitation amplitudes, and ω is the excitation energy. When a hybrid- functional is used, matrices A and B can be expressed in Mulliken notation as follows:

$\begin{align} & {{A}_{ia,jb}}={{\delta }_{ij}}{{\delta }_{ab}}\left( {{\varepsilon }_{a}}-{{\varepsilon }_{i}} \right)+\left( ia\left| jb \right. \right)- \\ & {{C}_{HF}}\left( ij\left| ab \right. \right)+\left( 1-{{C}_{HF}} \right)\left( ia\left| {{f}_{XC}} \right|jb \right) \\ \end{align}$

(2)

${{B}_{ia,jb}}=\left( ia\left| bj \right. \right)-{{C}_{HF}}\left( ib\left| aj \right. \right)+\left( 1-{{C}_{HF}} \right)\left( ia\left| {{f}_{XC}} \right|bj \right)$

(3)

where ij are used for ground state occupied orbitals, a, b for vitual orbitals, δ is the Kronecker delta, ε is the orbital energy, f_xc is the exchange-correlation kernel, and C_HF is the percentage of Hatree- Fock exchange in the hybrid functional.

For a long-range CT transition, the overlap between the original and final orbitals decreases to zero. Thus, A and B can be simplified to:

${{A}_{ia,jb}}={{\delta }_{ij}}{{\delta }_{ab}}\left( {{\varepsilon }_{a}}-{{\varepsilon }_{i}} \right)-{{C}_{HF}}\left( ij\left| ab \right. \right)$

(4)

${{B}_{ia,jb}}=0$

(5)

The neglect of matrix B for a CT transition is well known as the Tamm-Dancoff Approximation (TDA) ^{[38, 39]}. The first term in the expression of matrix A is the energy difference between the virtual and occupied orbitals, while the second is a Coulomb-like term that is essential in reproducing the 1/R fall-off condition that results from electrostatic attraction between the separated hole (orbitals that lose electrons) and electron (orbitals that obtain electrons). The absence of a selection rule on C_HF generally prevents exact prediction of the transition energies. Additionally, the optimized geometries in excited states are generally more sensitive than those in ground state to C_HF due to the relatively flatter potential energy surfaces in the excited states ^[33], which makes predicting the transition energies (e.g., E_ve and E_AD) even more difficult on the basis of the excited state geometries.

To theoretically interpret the excitation process in a given environment, the solvation model must be introduced to the TDDFT calculation. Recent studies have indicated that the solvation of some TADF molecules cannot be satisfactorily reproduced by employing a polarizable continuum model (PCM) ^{[16, 40-42]}, while the effects of other solvation models such as CPCM ^{[43, 44]} SS ^{[45, 46]} COSMO ^[47], SCIPCM ^[48], SMD ^[49], or SM8 ^[50] remain unknown. In addition to the solvation model used, the size of the basis set ^{[37, 40, 51]} and weak interaction ^[40] also affect the TDDFT results (though the influence of either is relatively small).

2.2. Issues in spectrum analysis

Experimental spectroscopic data is the primary evidence utilized to evaluate the reliability of calculated transition energies. Although the E_va(S₁) , E_ve(S₁) , E_0-0(S₁) , E_0-0(T₁) , and ΔE_ST of TADF materials can be determined based on their adsorption, fluorescence, and low-temperature phosphorescence spectra in theory, the complicated experimental spectra do not readily provide reliable data - as discussed above, this makes developing a theoretical approach for the prediction of transition energies quite difficult.

The adsorption and emission spectra of organic molecules in gaseous state are narrow and sharp. Consequently, they allow level-to-level transitions to be easily identified (Fig. 3). In medium polar solvents (e.g. toluene) or organic semiconductor films, the spectrum of an LE transition may retain well-defined vibronic structures (e.g., anthracene and carbazole) while that of a CT transition is always broad and featureless (Fig. 3) ^[52]. According to the Franck-Condon principle, the maxima of CT adsorption and emission correspond to E_VA and E_ve, respectively, while the average of EVA and EVE or the crossing point between absorption and emission (normalized to the same height) corresponds approximately to the E_0-0 of the CT transition ^{[53, 54]}. A CT transition with slight overlap of the involved orbitals exhibits a weak absorption band which is easily covered by other higher transitions (e.g., S₂), making it difficult to determine E_va(S₁). An assumed mirror-image relationship between the first absorption band and the fluorescence spectrum (mirror-image rule) was recently used to determine E_va(S₁) from easily identified E_ve(S₁) ^[37], however, the premise of this rule is a small structure variation between the ground and excited state geometries, which cannot be guaranteed in CT transitions ^[55].

For the emission spectrum of a CT transition in medium polar medium, the intersection point (onset) between the base line and the tangent line of the high energy side (passing through the halfheight point) corresponds closely to the 0-0 energy (Fig. 3) ^[56]. In high polar solvents, however, there may be strong red-shift of the emission spectrum weakens the correspondence between the emission onset and the 0-0 energy.

The determination of E_0-0(T₁) from low-temperature phosphorescence spectra is a bit more complicated. Theoretically, similar to fluorescence spectra, the onset of a broad emission band corresponds to the E_0-0 of a ³CT transition while the highest energy peak of a vibrational emission band corresponds to the E_0-0 of a ³LE transition. However, different to the room temperature fluorescence which is mainly attributed to the radiative S₁→S₀ transition, T₂→S₀ transition as well as the TADF can contribute to the delayed spectra measured at low temperatures (Fig. 4a). Rigid and low-temperature (usually 77 K) environments slow the T₂→T₁ internal conversion but cannot inhibit the T₁→T₂ or T₁→S₁ upconversion ^{[12, 15, 56]}. To liminate the influence of delayed fluorescence, the testing temperature for any emitters with very small ΔE_ST (＜0.1 eV) should be extra low ^{[12, 15, 56]}. Comparison among the fluorescence and the phosphorescence spectra of molecular fragments enables us to judge the nature of the low temperature delayed spectra (Fig. 4b) ^{[8, 15]}, but even so, it is difficult to reliably estimate the E_0-0(T₁) of certain molecules from the spectra due to their low phosphorescence yield at low temperatures ^{[16, 57]}.

Estimating ΔE_ST by a different approach than emission spectra is necessary. In 2012, Gousih et al. ^[6] deduced the following equation for the reverse intersystem crossing (RISC) rate constant (kRISC) of TADF molecules:

${{k}_{PISC}}=\left( \frac{{{k}_{p}}{{k}_{d}}}{{{k}_{ISC}}} \right)\left( \frac{{{\phi }_{d}}}{{{\phi }_{p}}} \right)$

(6)

where k_p and k_d are the rate constants of the prompt and delayed fluorescence components, respectively, k_ISC is the intersystem crossing (ISC) rate constant from S₁ to T₁, and φ_p and φ_d are the photoluminescence quantum efficiencies of the prompt and delayed components. Because ln(k_RISC) is inversely proportional to ΔE_ST/T (where T is temperature, ) ΔE_ST can be obtained through linear fitting of ln(k_RISC) with a series of T.

Meanwhile, Nakagawa et al. ^[7] suggested estimating ΔE_ST using a Berberan-Santos plot from the following equation:

$\ln =\left[ \frac{\phi }{{{\phi }_{d}}}-\left( \frac{1}{{{\phi }_{T}}}-1 \right) \right]=\ln \left( \frac{{{k}_{p}}+{{k}_{nr}}}{{{k}_{RISC}}} \right)+\frac{\Delta {{E}_{ST}}}{RT}$

(7)

where φ_T is the triplet formation efficiency, k_p is the phosphorescence rate constant, k_nr is the non-radiative rate constant from T₁, and R is the ideal gas constant.

To estimate ΔE_ST using Eq. (6) or (7) , it is necessary to examine the quantum yields and emission decay characteristics of the samples at various temperatures between 140 K and 300 K, which is high cost and time-consuming. Very recently, Zhang et al. ^{[15, 16]} established a simple method for estimating the ΔE_ST of a TADF emitter at room temperature via the following approximate relationship:

${{k}_{d}}=\frac{1}{3}\left[ {{k}_{p}}\cdot \exp \left( -\frac{\Delta {{E}_{ST}}}{RT} \right) \right]$

(8)

This equation is simplified from a more general expression for the temperature dependence of the lifetimes of two excited states in thermal equilibrium when the gap between these two states ranges from 0.05 eV to 0.30 eV. According to this equation, if a molecule has a TADF lifetime as long as hundreds of microseconds, its ΔE_ST can be hardly lower than 0.1 eV.

3. Steps toward theoretical prediction of TADF 3.1. Traditional hybrid density functional methods

The properties of TADF materials are mainly determined by ΔE_ST as well as the energy and oscillator strength of the S₁→S₀ transition ^{[8, 12, 15]}. As mentioned above, choosing the right functional to calculate these parameters within TDDFT remains a problem. The most popular functionals used in TDDFT are generally hybrids (e.g., B3LYP ^[58], PBE0 ^[59], and M06-2X ^[60]) that perform well in terms of ground state calculations. Range- separated functionals (e.g., CAM-B3LYP ^[61], LC-GGA ^{[62, 63]} and vB97X-D ^[64]) in which the Coulomb operator is divided into short-range and long-range parts are also good candidates. In comparison with the theoretical energy of an LE transition, that of a CT transition is more sensitive to the Hartree-Fock percentage (HF%) in the XC functional of TDDFT. Previous studies have indicated that traditional hybrid functionals with low HF% such as B3LYP (20%HF) and PBE0 (25%HF) inherently underestimate the vertical transition energies of TADF molecules (Fig. 5) ^{[12, 65-70]} with few exceptions ^[51] (Table 1). By contrast, high HF% in the long-range parts of CAM-B3LYP (19%HF in short-range and 65%HF in long-range) and ωB97X-D (22%HF in short-range and 100%HF in long-range) generally leads to excess transition energy ^{[12, 41, 71, 72]} although quite accurate results can be obtained sometimes via CAM-B3LYP ^{[73, 74]}.

Very recently, Moral et al. ^[40] found that the double-hybrid functional B2GP-PLYP ^[75] performs well in calculating E_VA for several TADF molecules. For instance, the absolute difference between the theoretical and experimental values of E_VA(S₁) for ACRFLCN and 2CzPN are both smaller than 0.1 eV. The cost of a double-hybrid functional is typically one order of magnitude higher than that of a traditional hybrid functional, however, which limits the application of this method to further studies on TADF ^{[76, 77]}.

There are two main ways to calculate ΔE_ST: Taking the difference between E_VA(S₁) and E_VA(T₁) (ΔE_ST^vert), or taking the energy difference between the optimized S₁ and T₁ states ΔE_ST^{adia [32, 34, 41]}. As illustrated by the Franck-Condon principle energy diagram shown in Fig. 6, the first method is theoretically reliable only when S₁ and T₁ have nearly identical equilibrium geometries, i.e., the orbitals involved in S₁ and T₁ transitions are almost the same. Conversely, if the T₁ of a CT molecule is in a locally excited state, the first method tends to overestimate ΔE_ST because the relaxation energy of a CT state is always larger than that of a LE state (Fig. 6b). The ΔE_ST^vert calculated at the B3LYP level sometimes shows good agreement with the experimental ΔE_ST ^{[8, 12, 41, 65]} (Table 1) owing to an undervaluation of the S₁ (CT) transition energies and an error cancellation. Although the second method seems more reasonable in theory, especially when the zero-point energy is considered, the selection of correct XC-fuctionals for S₁ and T₁ geometry optimization is problematic. For instance, the ΔE_ST^adia of DTC-DPS, 2CzPN, 4CzIPN, and 4CzTPN calculated at the B3LYP level are all negative values representing a qualitative mistake ^[78].

3.2. Semiempirical descriptor-based optimal Hartree-Fock percentage (OHF) methods

There is most likely no exact and universal exchange- correlation kernel for predicting the energies of all transitions. In 2008, in effort to warn against possible TDDFT failures, Peach et al. adopted a descriptor ∧ to depict the degree of overlap between the initial and final orbitals in a transition and assess the suitability of B3LYP, PBE0, and CAM-B3LYP for transitions with different ∧ values ^[79]. A number of descriptors were then designed to depict the characteristics of various CT transitions [8084]. In 2013, Huang and Zhang et al. introduced the descriptor method to the TDDFT prediction of TADF and significantly improved the computational accuracy ^[37].

Huang and Zhang et al. calculated the E_VA(S₁) of all TADF molecules (Fig. 5) reported at that time via several different XC- functionals to find that the calculated E_VA(S₁) strongly depends on the HF%, almost regardless of the difference in the expression of said XC functionals. The continuous and monotonic dependence of E_va(S₁) on HF% allows for an optimal HF% (OHF) in the E_VA(S₁) calculation (Fig. 7a). There appears to be a highly linear correlation between OHF and CT amount q (OHF = 42q), where q is defined as the quantity of the charge on un-overlapped HOMO or LUMO (Fig. 7b). E_va(S₁) and E_ve(S₁), as well as the E_0-0(S₁) that approximates to their average, can be predicted according to this empirical correlation. Because it is not possible to judge the nature of T₁ based on configuration interaction description, the E_VA(T₁) values points corresponding to a CT or LE excited state were determined according to the change in E_VA(S₁)-E_VA(T₁) compared to HF% (Fig. 7a). Two semi-empirical formulas were then derived to calculate E_0-0(³CT) from E_0-0(¹CT), E_VA(¹CT), and E_VA(³CT), and to calculate E_0-0(³LE) from E_VA(³LE). After calculating E_0-0(S₁)and E_0-0(T₁) independently, the ΔE_ST of various TADF molecules were accurately reproduced (Table 1).

The corrected TDDFT method on the basis of OHF has been applied to the molecular design for TADF materials intended for OLED application. Highly efficient blue TADF emitters having extremely small ΔE_ST ^[15] and efficient red TADF emitters having both high transition rate (large f) and small ΔE_ST ^[16] have been successfully developed with the assistance of this method, in fact. Reducing the empirical content as well as number of processing steps would make the OHF method even more attractive.

In the last couple of years, more and more descriptors have been introduced into studies on TADF. For instance, Moral etal. observed a dependence of the experimental ΔE_ST on the calculated electron- hole distance Δr(NTO), an index defined as the distance between the centroids of the highest occupied and the lowest unoccupied natural transition orbitals (NTOs) ^[40]. Penfold pointed out that the traditional descriptor ∧ (marked O in his article) does not consider the (de)localization and spatial extent of related orbitals in a molecule and therefore fails to accurately indicate ΔE_ST; he suggested additionally using the HOMO-LUMO gap (ε_gap) to describe the orbital (de)localization degree and proposed a new descriptor ∧⋅ε_gap accordingly. The dependence of the experimental ΔE_ST on ∧⋅ε_gap has proven a better indicator than that on A ^[85]. Chen et al. established the following relationship regarding ΔE_ST：

$\Delta {{E}_{ST}}\propto \frac{\Lambda _{S}^{2}}{\Delta {{r}_{S}}}+\frac{\Lambda _{T}^{2}}{\Delta {{r}_{T}}}$

(9)

where ∧_S and ∧_T are the overlap integral between the dominating orbital pairs of S₁ and T₁ transitions, while Δr_S and Δr_T are the corresponding transition distances ^[86]. Accordingly, localized T₁ state may lead to large ΔE_ST. Like q, these recently proposed descriptors do have potential application terms of the proper selection of XC functionals for modeling ICT molecules.

3.3. Optimally tuned range-separated density functional

The optimal "tuning" approach based on the range-separated (RS) density functionals has attracted increasing research attention as it has proven to efficiently provide reliable prediction for the band gaps and charge-transfer excitation energies of organic molecular systems ^[87-94].The "non-empirical" nature of this tuning method lies in the way that the range-separated parameter v is determined not by fitting a large number of experimental data, but by satisfying the Koopmans theory to the maximum extent, i.e. ensuring the negative HOMO energy equals the vertical ionization potential.

The parameter ω in RS functionals exhibits significant system- dependence characteristics. By employing RS methodology, Penfold ^[85] and Sun et al. ^[78] obtained optimal v values based on LC-BLYP and LC-ωPBE functionals, respectively, for a series of TADF molecules. Their results indicate that the tuned v values (0.14-0.20 bohr^-1) decrease by half compared to their corresponding default values (0.33-0.40 bohr^-1). The calculated mean absolute errors (MAE) between the theoretical and experimental data (Fig. 8) further shows that the corresponding optimally-tuned RS functionals can significantly reduce error by up to one order of magnitude, with MAEs of 0.15 eV for the lowest singlet excitation energies and 0.07 eV for the singlet-triplet gaps, compared to their non-tuned versions or other traditional functionals.

Although the tuning method returns very accurate results, tuned RS functionals have also been shown to overestimate the vertical excitation energies of PhCz, NPh₃, and CBP (Fig. 5) possessing weak charge-transfer characteristics (Table 1). In addition, they may underestimate the significant charge-transfer molecules accompanying small CT transition distance, i.e., the excitation energies of 4CzIPN, 4CzPN, 4CzTPN, 4CzTPN-Me (Fig. 5), which were underestimated by roughly 0.3 eV (Table 1).

Although the universality of the tuning method needs further improvement, it is generally reliable and can simultaneously provide accurate predictions for lowest singlet and triplet excitation energies. To this effect, the prediction of singlet-triplet gaps is real instead of the result of error cancellation. Further, the success of optimally-tuned RS density functionals relies on a reasonable balance among the description of electronic localized and delocalized characteristics that are critically important for organic π-conjugated systems. The tuning approach can efficiently eliminate the so-called delocalization error ^{[95, 96]} and has correct asymptotic behavior within the long-range limit, which are not guaranteed for traditional density functionals. Finally, it is important to note that the traditional way to obtain optimal v value is based on a fitting approach and requires a large number of single-point calculations for cation, anion, and neutral systems. To further improve the efficiency of the tuning procedure, Sun et al. developed a "golden-ratio" method for searching the optimal ω parameters; their method can provide highly precise ω in relatively few steps ^{[78, 94]}.

4. Conclusions

In this review, we discussed several challenges and problemsolving strategies related to the TDDFT prediction of TADF. After excluding inappropriate ways of determining experimental and theoretical transition energies, the accurate prediction for vertical and zero-zero transition energies of TADF molecules appears to be highly dependent on the asymptotic behavior of the chosen XC- functionals. To a certain extent, the newly developed descriptor- OHF method and the optimally-tuned RS functionals eliminate the parameter selection problem in TDDFT, though their efficiency merits further improvement.

In addition to TADF materials for OLEDs, ICT compounds have shown considerable potential for application in organic photonics and electronics such as organic solar cells, sensors, probes, and nonlinear optics. A computer method that can accurately and rapidly predict the photophysical parameters of various ICT molecules would significantly promote innovative development in these fields.

Acknowledgment We are grateful to the Exploratory Research for Advanced Technology (ERATO) of Japan and the Key Special Program of the Ministry of Science and Technology of China (No. 2016YFB0401000) for financial support.