On the performance of DFT/MRCI for singlet–triplet gaps and emission energies of thermally activated delayed fluorescence molecules

Abstract

This work investigates the performance of the density functional theory multireference configuration interaction (DFT/MRCI) method for the donor–acceptor and multi-resonance thermally activated delayed fluorescence (TADF) emitters of the recent STGABS27 benchmark set [L. Kunze, A. Hansen, S. Grimme and J.-M. Mewes, J. Phys. Chem. Lett., 2021, 12, 8470–8480]. Comparing the accurate experimental singlet–triplet energy gaps (ΔE_ST) and fluorescence energies (E_EM) to values computed with DFT/MRCI reveals a robust performance without large or systematic errors. Specifically in the vertical approximation without a solvation model, DFT/MRCI achieves mean absolute deviations (MADs) for singlet–triplet gaps and emission energies of 0.06 eV and 0.21 eV, respectively. Surprisingly, these values do not improve systematically when geometric relaxation and state-specific solvation effects are included. Apparently, part of these effects are absorbed in the parameterization of DFT/MRCI and attempting to include them explicitly via a ROKS+PCM reaction field leads to an imbalanced treatment. As a result, the simplest approach of running calculations in the vertical approximation in gas phase turns out to be the most accurate. Albeit less accurate and more computationally demanding than state-specific orbital-optimized DFT, DFT/MRCI has the advantage that all low-lying excited states are obtained in a single calculation, including transition properties between them. At the same time, the aforementioned performance for the ΔE_ST and E_Em values is achieved without molecule-specific or state-specific adjustments like optimal tuning that is often necessary for time-dependent DFT. Hence, we conclude that DFT/MRCI is particularly useful during the initial stage of computational investigations of TADF emitters to screen for the ΔE_ST and identify the relevant states, whose energies can then be refined with accurate state-specific DFT methods like ROKS or (Δ)UKS with MADs for ΔE_ST below 0.03 eV.

---

1 Introduction

Thermally activated delayed fluorescence (TADF) emitters have received great attention in the context of organic light-emitting diodes (OLEDs) and other optoelectronic devices. The delayed fluorescence is enabled by a small energy gap ΔE_ST between the lowest excited singlet (S₁) and triplet states (T₁), which facilitates efficient reverse intersystem crossing (rISC). Since singlet and triplet excitons can be harvested this way for emission from the S₁ state, the theoretical maximum fluorescence quantum yield can, in principle, increase from 25% in pure fluorescence emitters to 100% in TADF emitters.^1–4 Hence, the design goal underlying the numerous and structurally diverse TADF emitters proposed in recent years is a vanishing ΔE_ST, which ensures a minimal barrier in the thermally-driven rISC process. One way to achieve a small singlet–triplet gap on the order of the thermal energy (k_BT ≈ 0.025 eV) is by reducing the magnitude of the respective exchange-type integral that leads to the energetic differences of singlet and triplet states with otherwise identical spatial electron configurations. In the literature, minimizing this exchange-type integral is sometimes misleadingly set equal to reducing the overlap between the highest occupied (HOMO) and lowest unoccupied molecular orbital (LUMO), which is, however, always zero for canonical orbitals. This design principle was prominently exploited by Adachi and coworkers.^5,6

Spatially separated donor (D)–acceptor (A) pairs offer a convenient way to design molecules by exploiting low-lying charge transfer (CT) states, which minimize the singlet–triplet gap.^7–9 In the resulting DA-TADF emitters, the lowest excited singlet CT state is nearly degenerate with the respective triplet CT state, often with locally excited (LE) triplets in close energetic proximity. The latter are typically of ππ* or nπ* character with more spatial overlap of the electron–hole pair than in the CT states. Nearby LE triplets lead to mixing with the triplet CT state, and thus, to small but finite spin–orbit coupling (SOC) with the pure singlet CT state (SOC would be zero between pure CT states), which is sufficient to explain the observed rISC rates.^10,11 A disadvantage of the DA emitter design is that the emission originates from a singlet CT state, whose strong coupling to the dielectric environment renders the emission undesirably broad and, because of the small spatial overlap of the electron–hole pair, also slow (low oscillator strength f_osc). Both effects negatively impact the overall luminescence quantum yield.^10–14

So-called multi-resonance (MR-)TADF emitters based on alternating electron-donating and electron-accepting atoms (e.g., nitrogen and boron, respectively) in a conjugated π-system offer an alternative.^15–19 Compared to DA-type emitters, MR-TADF emitters have reduced structural flexibility and lower excited-state polarity leading to sharp emission spectra with high quantum yields. However, this comes at the cost of a generally larger singlet–triplet gap compared to donor–acceptor architectures, which, so far, prevented their commercial application as TADF-enabled emitters. Nevertheless, due to their very favorable sharp emission, DABNA derivates are used commercially as fluorescent blue emitters,²⁰ and the development of advanced MR emitters with a small ST gap is a hot research topic.²¹ Because of their relevance, the last two molecules in the STGABS27 set are MR-TADF emitters. Further development of the MR-TADF-based emitter concept aims at more efficient singlet-to-triplet conversion processes by enabling inverted singlet–triplet gaps (INVEST).^22–25

With the growing interest in TADF emitter materials, computational methods have become more important for the virtual screening of candidates.^26–33 However, an accurate description of both the strong CT states in earlier DA-type emitters and the electronic structure of MR- or INVEST-type emitters, which involves double excitations,³⁴ pose challenges for established theoretical methods.^25,35,36 In the case of CT states, commonly used linear response methods based on (hybrid) density functional theory (TD-DFT³⁷) or its Tamm–Dancoff-approximated variant (TDA-DFT³⁸), lead to unreasonably underestimated CT energies due to failures of the local density functional approximation,^16,31,39 and suffer from a lack of accurate excited-state solvation models.⁴⁰ As a result, most protocols to obtain accurate ΔE_ST values with TD-DFT often rely on some kind of error cancellation, which comes at the cost of strongly underestimated emission energies E_Em and wrong state characters.^40,41 Recently, some progress was reported concerning this issue by Khatri and coworkers who applied dielectrically-screened range-separated hybrid functionals to the DA-TADF emitters in the STGABS27 benchmark. In doing so, they achieved a mean absolute deviation of 0.06 eV for ΔE_ST and reasonable emission energies (the latter were not evaluated).⁴² However, for MR-TADF and INVEST emitters, the issues of TD-DFT are more fundamental. For these systems, moving to double-hybrid density functionals (DHDFs) with perturbative second-order corrections is another option that has been explored and found to repair TD-DFT for MR-TADF⁴³ and inverted singlet–triplet gap emitters from the INVEST set,⁴⁴ but this comes at an increased computational cost and steeper scaling with system size.

Realizing that DHDFs improve the agreement, it is no surprise that correlated wave-function-theory-based methods like linear response second-order approximate coupled cluster singles and doubles (LR-CC2⁴⁵), second order algebraic diagrammatic construction (ADC(2)^46,47), and especially their spin-scaled variants^48,49 (SCS-CC2, SCS-ADC(2)) are also accurate for all kinds of TADF systems, including MR-TADF and INVEST.^11,43,50 However, their steep scaling with system size renders routine calculations on emitters with often more than 100 atoms impractical, especially in combination with the desirable triplet-zeta basis sets.³⁶ Another technical issue of wave-function-theory-based methods, especially for DA-TADF emitters, arises from the lack of complete solvation models. Although state-specific equilibrium, and non-equilibrium solvation models are available for ADC(2) and even ADC(3),^51–53 excited-state gradients have only recently been introduced and are again prohibitively expensive for routine application.⁵⁴

It should be emphasized that in the aforementioned linear response and configuration interaction (CI)-type approaches, it is the linear combination coefficients (amplitudes) of the individual excited determinants that are variationally optimized to yield the many-particle wave function of the excited-state. Unless combined with a state-specific solvation treatment, the underlying orbitals of the reference determinant remain unaltered. A conceptually different approach is the one followed in state-specific orbital-optimized DFT (OO-DFT) methods. Here, the many-particle wave function comprises only one (e.g., UKS)⁵⁵ or (a set of) a few determinants (e.g., OSS-ROKS)⁵⁶ with fixed values for the amplitudes. In the case of a few determinants, these are generally chosen to be the coupling coefficients to form a suitable spin-pure open-shell configuration state function (CSF). Instead of optimizing the amplitudes, it is the molecular orbital coefficients that are variationally optimized in OO-DFT until a self-consistent field (SCF) is achieved. Consequently, these methods are similar and, to some degree, equivalent to ground-state DFT, the only difference being that they impose an approach for the many-particle wave function that uses an excited-state instead of the ground-state aufbau occupation. Different flavors of OO-DFT are in use, and we will review the most relevant ones here, due to their increasing relevance for excited-state calculations. The most widely used variant is the one that combines a so-called ΔSCF treatment with an unrestricted Kohn–Sham determinant (in short (Δ)UKS). ΔSCF methods^55,57,58 use schemes to enforce the excited-state occupation inside an, otherwise, unaltered ground-state SCF procedure. The outcome of a ΔUKS treatment to describe a singly excited singlet state is an open-shell determinant which is used as a proxy to compute the singlet-state energy. Yet, the converged open-shell determinant is never a spin eigenfunction, irrespective of additional spin contamination due to symmetry breaking in (Δ)UKS. Still, in practice, this is of less concern for the computation of energies and excited-state geometries. ΔSCF approaches manage to converge quite robustly for well-separated excited-states, such as the lowest excited-state or core excited-states, but convergence issues can arise, e.g., when other possible solutions are nearby.

To ameliorate the aforementioned shortcomings of ΔUKS, a restricted open-shell Kohn–Sham (ROKS^56,59) approach for open-shell singlet (OSS) states was first developed by Frank et al. and later picked up by van Voorhis and coworkers. This variant of OO-DFT starts from Ziegler's spin purification formula⁶⁰ and minimizes the resulting OSS energy expression using a restricted set of orthonormalized orbitals, hence the term OSS-ROKS. This not only leads to a pure spin eigenfunction, but to a variational minimum solution corresponding to the energetically lowest OSS that can be represented by exactly two determinants.

OO-DFT methods were shown to provide an excellent description of states that are inherently difficult for TD-DFT, such as X-ray and CT states.^61–63 Likewise, such OO-DFT methods were shown to be very successful in computing the singlet–triplet gaps and emission energies of TADF emitter molecules.^25,36,41,64 For various types of purely organic TADF emitters, ΔUKS (and partially also ROKS) with a range-separated hybrid functional and coupled to a solvation model recovers experimental or high-level theoretical singlet–triplet gaps with sub-kcal mol⁻¹ precision (mean absolute deviation well below 0.05 eV), with most errors below the thermal energy at ambient conditions of k_BT ≈ 0.025 eV. The good performance of these methods stems from the state-specific orbital optimization and the resulting natural inclusion of state-specific solvation effects, which are more difficult to incorporate in linear response and CI-type approaches.^11,64 Because of this natural inclusion of solvation effects, OO-DFT methods are among the few methods that afford excited-state geometry optimizations in solvent, which is especially important for CT states due to the failure of TD-DFT+LR-PCM.^40,64

However, their performance is naturally limited by the chosen many-particle wave function and its ability to describe only one excited-state of interest at a time. While OSS-ROKS is restricted to singly excited OSS configurations, ΔUKS yields metastable solutions that are no spin eigenfunctions. Though the issue of ΔSCF becoming unstable and collapsing to a lower-lying (non-targeted) state is relevant for higher-lying and nearly degenerate states, this is not a practical issue for calculation of the well-separated HOMO–LUMO-dominated excited-states in many TADF emitters in STGABS27. From a pragmatic point of view, the main downsides of OO-DFT include the difficulty to compute absorption and emission intensities (transition properties in general), the issue of state targeting, and the lack of oversight over other states as each calculation yields only a single (excited) state. These disadvantages become particularly relevant when multiple energetically close states that are dominated by orbitals further away from the HOMO and LUMO, e.g., in TADF emitters with more than one donor or acceptor moiety.

Particularly with respect to the latter point, another viable but less explored option are multireference CI methods with modified Hamiltonians to effectively capture dynamic correlation in a truncated CI expansion. This way, other important excited-states aside from just the lowest one can be identified to gain an overview and consider in a more thorough subsequent investigation with the aforementioned state-specific methods. This is of particular relevance if multiple adiabatic states of different character (local vs. CT excitations) are energetically close and can switch order during excited-state relaxation, which was found to be the case in bridged DA-type TADF emitter molecules.³² In the “real-life” scenario, no a priori information about the number and nature of involved states is available. Thus, an efficient and sufficiently accurate CI-type method is desirable to compute multiple low-lying excited-states, generally, starting from the ground-state geometry.

As such, we explore DFT in conjunction with multireference configuration interaction (DFT/MRCI) in the present work.^65,66 By explicitly including multiply excited configurations, DFT/MRCI promises to account for both the orbital relaxation effects required to properly describe CT states and the description of doubly excited states, which is challenging even for truncated coupled-cluster methods like CC2.^67,68 Admixture of doubly-excited determinants was found to be important in some MR-TADF emitters.^25,35 Unlike conventional (MR)CI, the reference configuration of DFT/MRCI is not constructed from (uncorrelated) Hartree–Fock (HF) but from correlated DFT orbitals. Thus, electron correlation is introduced through the underlying exchange–correlation (XC) functional, an empirical scaling of CI-matrix elements, and naturally through the MRCI expansion. Large fractions of Fock exchange, which are employed for computing the Kohn–Sham orbitals of the reference configuration and appropriately scaling the two-electron integrals in DFT/MRCI, effectively alleviate the CT failure observed for TD(A)-DFT.⁶⁶ DFT/MRCI was originally parametrized for the BHLYP⁶⁹ functional (50% Fock exchange) and only recent work explored the use of other functionals for this method.⁷⁰ The problem of “double counting” of electron correlation in DFT/MRCI is reduced by scaling of CI matrix elements.⁶⁵ The damping of energetically distant CSFs leads to a sparse MRCI Hamiltonian. Additionally, CSFs with a mean energy above a certain threshold are discarded,⁶⁶ thus making the DFT/MRCI method computationally much cheaper than conventional HF-based counterparts and applicable to large systems with over 100 atoms. From past experience, DFT/MRCI represents a robust CI scheme to be used in a nearly “black-box” manner to describe multiple excited-states of various molecular systems with proper spin symmetry. It has already proven to be useful in explaining the properties of different TADF materials.^13,71–75 One of us has successfully used it in describing the “hot exciplex” mechanism in anthracene-bridged DA-type TADF emitters.³² Here, the treatment of higher-lying and near degenerate excited-states is crucial, something that could not be achieved easily with ΔUKS methods. However, a benchmark study involving the by now established STGABS27 benchmark set has not been carried out so far. With this work, we wish to provide a thorough benchmarking of DFT/MRCI for computing the singlet–triplet gaps and emission energies of the STGABS27 TADF emitters. In outlining crucial computational parameters, we present an efficient protocol for describing the excited-states of these emitters.

2 Benchmark set

For the assessment in this work, we consider the STGABS27 benchmark set (see Fig. 1 for the molecular structures) proposed by Kunze et al.⁶⁴ It contains accurate ΔE_ST values based on the temperature-dependent measurement of the TADF rate.⁷⁶ Due to the exponential dependence of the TADF rate on ΔE_ST, we consider here ±0.05 eV as an upper limit for a viable prediction of this quantity, even though previous theoretical calculations indicate an even smaller experimental error. Recently, three of us complemented the data with experimental emission energies to ensure a physically consistent description of the whole TADF process and recognize purely error-cancellation-based computational protocols.⁴¹ While the experimental uncertainty of E_Em is limited to about 0.20 eV by the assumption of vertical transitions, previous theoretical calculations indicate this to be a rather conservative upper bound.⁴¹ The STGABS27 benchmark set comprises 27 TADF emitters, grouped according to their structural motifs in Fig. 1.

Fig. 1 The STGABS27 benchmark set categorized by structural motives: linear-, branched-, spiro donor–acceptor (DA), and multi-resonance (MR) emitters. Highlighted in orange is a compact representative subset, named STGABS15, which involves smaller systems of the benchmark set. System (12) is shown in light grey and was excluded from our assessment.

15 emitters show a linear alignment of the DA moieties, eight are branched, two feature a spiro motif, and two exhibit multi-resonance TADF, respectively. We excluded the largest emitter (12), which consists of 139 atoms, from our assessment due to computational costs and challenges in the convergence of the configuration selection⁶⁶ in the DFT/MRCI calculation. Based on the remaining molecules, a structurally diverse subset of smaller emitters with a minimum of one molecule per structural group, named STGABS15, was assembled (see Fig. 1). This subset was used to explore the influence of computationally more demanding settings for DFT/MRCI.

3 Theory

Theoretical aspects of the DFT/MRCI method are briefly reviewed in this section for a restricted closed-shell reference. The working equations that describe the (off-)diagonal matrix elements of the Hamiltonian are given below. We will mostly adapt the notation used in ref. 66 and refer the interested reader to this article or the original work ref. 65 for further details on the technical aspects of this method. The DFT/MRCI Hamiltonian matrix elements can be divided into diagonal and off-diagonal elements between CSFs, which are defined by their spatial occupation number vector (ONV) w and the spin-coupling pattern ν.

3.1 Diagonal Hamiltonian matrix elements

In analogy to conventional (MR)CI based on a restricted HF reference, the diagonal Hamiltonian matrix elements between CSFs sharing the same w and ν read:whereand we use Mulliken notation to denote the two-electron integrals:with basis functions ϕ and electron coordinates r_i. Eqn (1) contains a HF-like matrix element 〈νw|Ĥ − E^HF|νw〉 and is modified with effective one- and two-electron corrections. n_exc is the number of excitations for a given ONV w with respect to the restricted Kohn–Sham (KS) anchor configuration. The occupation numbers of MO i in the anchor configuration are given by w̄_i. The indices i and k denote occupied orbitals in the reference configuration. We highlight here, that quantities denoted with “HF” are not explicitly formed from HF, but from the KS orbitals according to eqn (1) and (2). Hence, F^HF_ii is an effective one-electron matrix element constructed from the KS orbitals, but using the ab initio Hamiltonian leading to Fock matrix elements in the KS canonical orbital basis.⁶⁶ Similarly, the E^HF is not a true Hartree–Fock energy but computed using the summed diagonal Fock matrix elements of the occupied canonical KS orbitals and subtracted from this, the double counting of the electron–electron interaction. The effective difference in electron–electron interactions relative to the reference configuration due to excitations in DFT/MRCI appears through the Coulomb ΔE_J and exchange-type ΔE_X correction terms. In the original, multiplicity-specific Hamiltonian,⁶⁵ this term reads:

Here, the indices a and i denote created or annihilated electrons relative to the anchor configuration, respectively. p_J and p[N_o] are empirical parameters in DFT/MRCI determined by a global fit to minimize the root mean square error (RMSE) on vertical excitation energies for singlet and triplet states. Initially, this set comprised experimental data (e.g., from photoelectron loss spectroscopy)^65,77 and was extended for redesigns of the Hamiltonian, thereby also including theoretical reference data.^78–81 For clarity, we note here already that the multiplicity-specific parametrization of the original Hamiltonian has been removed in latter flavors accompanied with adjusted terms in the (off-)diagonal matrix elements. Redesigned expressions of eqn (5) were developed elsewhere,^66,81 but are not shown here.

3.2 Off-diagonal Hamiltonian matrix elements

The off-diagonal matrix elements in the DFT/MRCI Hamiltonian can further be divided according to equal or different (primed variables) ONVs and spin-parts of the CSFs, which leads to the following cases (eqn (6) and (7)):

〈νw|Ĥ|ν′w′〉 = 〈νw|Ĥ^CI|ν′w′〉·p₁ξ(ΔE_ww′;p₂) (7)

In the original publication, elements 〈νw|H^DFT|ν′w〉, i.e., with different spin-part and equal ONVs (eqn (6)), remained unmodified to conventional CI. In newer Hamiltonians,⁸¹ this was adjusted to account for an additional parameter p_X. In DFT/MRCI, “double counting” of electron correlation is counteracted by damping energetically distant CSFs with different ONVs as displayed in eqn (7). p₁ξ(ΔE_ww′;p₂) is an energetic damping function that operates with the effective CSF energy obtained from the sum of KS orbital energies ΔE_ww′ and the parameters p₁ and p₂. Again, with different flavors of the Hamiltonian, the employed energetic damping function was varied (see ref. 66 and 78).

The damping and, related to this, discarding of energetically distant CSFs critically determines the computational efficiency of DFT/MRCI. So-called “short” and “standard” parameter sets were optimized and, associated with them, are fixed energy cutoff values of 0.8 E_h and 1.0 E_h, respectively. Depending on the energy, an effective selection criterion to discard CSFs is employed,⁶⁶ ultimately determining the size of the truncated CI expansion and the computational cost of a DFT/MRCI calculation.

4 Computational details

The ground and excited state geometries used in this study were taken from ref. 64 without modifications. Since no analytical implementation of nuclear gradients is available for the DFT/MRCI method, different levels of theory were used for the optimized geometries. In principle, nuclear gradients can be obtained using numerical differentiation techniques. However, due to the extreme computational costs associated with this task,⁶⁶ this was not attempted here. The S₀ geometries from ref. 64 were calculated in the gas phase using KS-DFT with the composite method PBEh-3c.⁸² This level of theory combination has been used in the past for modeling DA-type TADF emitters.³² For the S₁ and T₁, the implicitly solvated ROKS/UKS-optimized geometries were selected due to their excellent performance on the STGABS27 set. Specifically, these geometries were optimized with the optimally-tuned (OT^83,84) LC-ωPBE^85,86 functional with the Q-Chem program (version 5.4)⁸⁷ in combination with D3 dispersion correction with Becke–Johnson damping (D3(BJ)^88,89) and a def2-SVP⁹⁰ basis set. During the geometry optimization of the excited states, implicit solvation in the integral equation formalism polarizable continuum model (IEF-PCM^91,92) with a dielectric constant of ε = 3 had been used. This dielectric constant was chosen to emulate the non-polar organic solvents (mostly toluene) and thin-film environments used in the experimental measurements for the STGABS27 (see ref. 40 and 64). In the following, we will use the short-hand notation PCM or state-specific (SS-)PCM, implying that the IEF-PCM formalism has been used.

The anchor configuration for DFT/MRCI calculations was computed with Turbomole (version 7.6).^93–96 The iterative MRCI procedure in DFT/MRCI works as an iterative MRCISD procedure (i.e., including single and double excitations from the reference space), which is considered converged when the targeted states in the Davidson procedure only contain contributions from CSFs that are also part of the reference space. Our initial CI reference space considered all single and double excitations in a (4,4) active space. Molecular orbitals (MOs) with orbital energies ε outside the interval −3.0 E_h < ε < +2.0 E_h of the BHLYP/def2-SV(P) anchor configuration are kept frozen for the MRCI, if not stated otherwise. We employed the resolution-of-the-identity approximation for Coulomb integrals (RI-J^97,98 and RI-C for ERIs in the CI Hamiltonian) with the corresponding auxiliary basis set during the BHLYP/def2-SV(P) reference calculation.⁹⁹ The command line define (cefine)¹⁰⁰ tool was used to generate the required input files. If not stated otherwise, we converged the CI space towards the three lowest excited-states for all calculations using “short settings”. This implies reducing the default selection threshold used to damp or remove energetically distant CSFs from the CI space from 1.0 E_h (“standard settings”) to 0.8 E_h. Accordingly, we used the corresponding “short” parametrization for the DFT/MRCI Hamiltonian.

Since there exists no state-specific solvation model for DFT/MRCI, we employed gas-phase as well as ground-state-solvated orbitals with the conductor-like screening (COSMO)¹⁰¹ implicit solvation model as implemented in Turbomole. State-specific solvation was emulated by reading in the converged excited state reaction field (RF) from a ROKS+PCM calculation performed in Q-Chem (version 5.4)^41,87,102 at the OT-ωB97M-V¹⁰³/def2-SVP level of theory. In the latter, the excited-state equilibrium RF (ε = 3) was included as point charges in the calculation of the reference BHLYP/def2-SV(P). It should be emphasized that DFT/MRCI is not a fully variational method and that solvation effects included in the reference calculation translate to the final DFT/MRCI energies through the direct use of canonical KS-DFT orbital energies in the CI Hamiltonian (see Section 3).

We also tested including non-equilibrium solvation effects for the emission energies via the respective perturbative state-specific correction (ptSS-PCM⁵¹) from the aforementioned ROKS+PCM calculation (see also SI, Sections S2 and S3).⁴¹ Complementary information on further geometry refinements, computational protocols, and raw data are provided in the SI.

5 Results and discussion

5.1 Hamiltonian and DFT/MRCI settings

As a first task, we assess the performance of the available DFT/MRCI parametrizations, including the original Hamiltonian by Grimme and Waletzke,⁶⁵ as well as the four Hamiltonians proposed by Marian and coworkers (R2016,⁸¹ R2017,⁸⁰ R2018,⁷⁹ and R2022⁷⁸). We will not review particular details of the Hamiltonians but refer, instead, to the original papers and this review article that covers all but the most recent version.⁶⁶ For this initial assessment, we will simply compare singlet–triplet gaps computed in gas-phase by each DFT/MRCI scheme in the vertical approximation at the S₀ PBEh-3c geometry. While we recognize that this is an approximation, given that the experiment refers to adiabatic gaps in solution, this approach serves as a simple means to compare the different DFT/MRCI schemes and select one that we will consider in more detail. Furthermore, this vertical treatment is actually quite representative of the “real-life” use case of DFT/MRCI. That is, if nothing is known about the excited states a priori, there are no relaxed geometries of particular excited states available, and one will start the investigation by computing multiple low-lying excited at the ground state geometry. Fig. 2 displays the deviations for ΔE_ST computed in the vertical gas phase approximation, while statistical measures are collected in Table 1.

Fig. 2 Violin plot for the difference between computed and experimental S₁–T₁ gaps for the available DFT/MRCI Hamiltonians. For all computations, short settings (0.8 E_h selection threshold and adjusted parameters, “short parametrization”) and a BHLYP/def2-SV(P) anchor configuration were used. The computed energy differences corresponds to vertically computed S₁–T₁ gaps at PBEh-3c S₀ gas phase geometries.

Table 1 Statistical measures (in eV) of the computed vertical and adiabatic ST gaps for the STGABS27 benchmark set using DFT/MRCI with different Hamiltonians, short settings, and a BHLYP/def2-SV(P) anchor configuration. MD: mean deviation, MAD: mean absolute deviation, AMAX: absolute maximum deviation, RMSE: root mean square error, SD: corrected standard deviation (see Table S1 for definitions). Structures for adiabatic gaps were optimized with ROKS/UKS-OT-LC-ωPBE-D3(BJ)/def2-SVP+SS-PCM (ROKS-OT-wB97M-V/def2-SVP+SS-PCM for ^a), while vertical calculations were conducted at the ground state PBEh-3c geometries

Hamiltonian	Vertical singlet–triplet gaps
	MD	MAD	AMAX (system)	RMSE	SD
a The values refer to the OT-ωB97M-V/def2-SVP+PCM level of theory.b The values refer to the OT-LC-ωPBE/def2-SVP+SS-PCM level of theory.
R2022^78	0.06	0.07	0.22 (10)	0.09	0.09
R2018^79	0.05	0.06	0.20 (10)	0.08	0.09
R2017^80	0.09	0.10	0.24 (10)	0.11	0.08
R2016^81	0.10	0.10	0.25 (24)	0.12	0.09
Original^65	0.06	0.08	0.26 (24)	0.10	0.11

---

	Adiabatic singlet–triplet gaps
	MD	MAD	AMAX (system)	RMSE	SD
R2022	0.12	0.12	0.35 (21)	0.15	0.13
R2018	0.10	0.11	0.34 (21)	0.14	0.13
R2017	0.14	0.15	0.35 (21)	0.17	0.12
R2016	0.15	0.15	0.33 (21)	0.17	0.11
Original	0.10	0.12	0.36 (8)	0.16	0.15
R2018 (ES eq.)	−0.07	0.14	0.57 (2)	0.19	0.21

---

Other methods	Adiabatic singlet–triplet gaps
	MD	MAD	AMAX (system)	RMSE	SD
ROKS^a^64	0.00	0.02	0.06 (25)	0.03	0.03
TDA-DFT^b^40	0.10	0.14	0.50 (21)	0.21	0.19

---

All four Hamiltonians by Marian et al. provide similar results concerning their standard deviation (SD), with only small differences in the systematic positive shift (see also Fig. S1). Overall, the R2018 Hamiltonian shows the smallest absolute maximum deviation (AMAX = 0.20 eV) and a mean deviation (MD = 0.05 eV) closest to zero for the computed ΔE_ST. The original Hamiltonian has the largest statistical errors (SD = 0.11 eV), while the MD remains almost on par (MD = 0.06 eV) with R2018. The larger spread in ΔE_ST may be related to the fact that the original Hamiltonian utilizes different parameter sets for singlet and triplet states, respectively. Later Hamiltonians unify the treatment of all spin-states initially motivated by the goal of a consistent treatment for non-covalent photoexcited dimers with four-fold open-shell orbital configurations.⁶⁶ Moreover, it is worthwhile to mention that beyond conceptional redesigns in later Hamiltonians, the fit set has been modified and expanded, which should generally enable a more robust parametrization.⁶⁶ Given the small differences in errors for ΔE_ST, the treatment of CT states in DFT/MRCI appears rather insensitive to the specific Hamiltonian.

We next used the S₁ geometries (cf. Section 4) to compare the performance of the different DFT/MRCI Hamiltonians for the calculation of emission energies. Inspection of emission energies (see Fig. S2 and Table S3) shows that for E_Em, the original Hamiltonian performs best (MAD = 0.14 eV), see Table S3 for more details. Most newer Hamiltonians (R2016, R2017, R2018) show overall similar performance to the original Hamiltonian with slightly larger MADs. For R2022, the MAD increases to 0.36 eV. Furthermore, the AMAX in E_Em for all Hamiltonians is due to the branched emitter p-AC-DBNA (16) and ranges between 0.43 eV (original) and 0.76 eV (R2022). The similar SD between all Hamiltonians shows that the multiplicity-specific parametrization of the original Hamiltonian does not lead to a consistent improvement of E_Em over the entire set. We will inspect trends regarding E_Em in more detail later and investigate how solvation effects influence the predictions obtained using DFT/MRCI. At this point, it should be made clear that the excited state minimum geometries are optimized with ROKS/UKS using OT-LC-ωPBE-D3(BJ)/def2-SVP in implicit solvation (IEF-PCM). However, no implicit solvation correction was employed in any of the DFT/MRCI single-point calculations discussed so far.

We tested, furthermore, the influence of the chosen basis set as well as the effect of selecting the more computationally demanding “standard” instead of the “short” parameter settings. We present this data and discussion in the SI. Increasing the basis set can decrease the ΔE_ST by up to 0.05 eV for DA-type emitters, but we find that, overall, increasing the basis set and moving away from the “short” settings does not lead to a significant improvement. Given the much lower computational demand, we continue with the “short” settings in combination with the def2-SV(P) basis set.

5.2 Comparison of adiabatic and vertical energies

After considering initially only singlet–triplet gaps in the vertical approximation, we now want to address structural effects in more detail and move on to the physically more sound adiabatic gaps. Since no analytical gradients are available for the DFT/MRCI method,⁶⁶ we employ the aforementioned ROKS/UKS+PCM-optimized excited state geometries for the S₁ and T₁ states. Fig. 3 compares the computed adiabatic gaps with the vertical ones for the R2018 Hamiltonian, while Table 1 collects again the statistical measures.

Fig. 3 Difference between experimental and computed ΔE_ST as obtained with DFT/MRCI using the R2018 Hamiltonian, short settings and a BHLYP/def2-SV(P) anchor configuration. Vertical gaps were obtained using PBEh-3c S₀ gas phase geometries (blue crosses) whereas the ROKS/UKS OT-LC-ωPBE-D3(BJ)/def2-SVP+PCM excited S₁ and T₁ geometries (dark red circles) were used to determine adiabatic gaps.

Notably, the deviations are generally larger for adiabatic gaps, both due to a more systematic overestimation (MD = 0.10 eV vs. 0.06 eV) and an increased statistical error (SD = 0.13 eV vs. 0.09 eV). While this contrasts the general improvements of adiabatic gaps observed with state-specific theories like ROKS/UKS,⁶⁴ it is reminiscent of the results reported for TDA-DFT.⁴⁰ As for TDA-DFT, we suspect that part of the error arises from using two distinct geometries for S₁ and T₁, which hampers stable error compensation observed with the single S₀ structure. In the case of adiabatic TDA-DFT calculations, this loss of error compensation necessitates the explicit treatment of state-specific solvation effects for even qualitatively correct ΔE_ST values.⁴⁰ The systematically more positive MD of the adiabatic gaps (see above and SI, Fig. S1) suggests that the lack of solvation for the more polar singlet state may also be an issue for DFT/MRCI, which we address below. A further explanation might be a bias in the DFT/MRCI parametrization for equilibrium structures, i.e., from using computed vertical excitations, while the experimental reference includes vibronic effects and may even be tilted towards adiabatic energies since the electron loss spectroscopy is a much slower process than photoexcitation.^65,70,77,81

While TDA-DFT can also improve on this set in performance based on fortuitous error cancellation, that requires a functional with only a negligible admixture of Fock exchange (≈10%).⁴⁰ Consequently, vertical emission and absorption energies are red-shifted by up to 1.0 eV. Out of the box, DFT/MRCI shows a rather good agreement with the experimental values (ΔE_ST and E_Em), while slightly overestimating the E_Em compared to experiment due to the missing solvation (see below and see Fig. S2 and Table S3 for emission energies).

Some systems in Fig. 3 show large discrepancies between the vertical and adiabatic ΔE_ST in the gas phase. One particular example is the linear emitter TPA-Ph2CN (8), where ΔE_ST amounts to 0.14 eV in the vertical approximation and 0.41 eV calculated using the excited state minimum geometries. We exemplarily checked how different ground state conformations affect the computed vertical ΔE_ST (Fig. S9), and found that this may be a minor source of error. ΔE_ST of system (8) is somewhat sensitive to changes in the dihedral angle between the D–A moieties θ_DA. ΔE_ST computed with DFT/MRCI ranges from 0.04 eV to 0.14 eV, i.e., a similar range that was found previously in DA emitters with OO-DFT by one of us.¹¹ This angular dependence also became evident from a comparison of different DFT composite methods used to compute the ground state geometries from which vertical ΔE_ST (Fig. S7) were computed. For this system, as well as systems 19–24, the variance of the ΔE_ST for different geometries was the most noticeable. A comparison of vertically computed ΔE_ST at the PBEh-3c ground-state and the adiabatically computed ΔE_ST at the OT-LC-ωPBE-D3(BJ)/def2-SVP+PCM excited-state geometries shows that the adiabatic gap becomes notably larger as a result of a decreased θ_DA (cf. Fig. S10). In most cases, changing the theory level among the various efficient “3c” methods to describe the ground-state geometry has a subordinate impact on ΔE_ST. Such methods are obvious and common choices to be used with DFT/MRCI for (vertical) calculation of excitation energies.^32,104,105

Another source of error is the absence of a proper solvation treatment, which generally leads to overestimated adiabatic ΔE_ST. The errors for (8) can be reduced if an approximate state-specific excited-state solvation treatment is introduced (cf. Fig. 5 below). Previously, Marian et al.⁷⁵ analyzed how the excited S₁ and T₁ states are altered with θ_DA for (8) and related compounds, and how this impacts different excited state properties. They could show that the use of explicit solvent molecules (toluene) reduces the adiabatic ΔE_ST as determined using DFT/MRCI as well. Given that there are other cases, where the adiabatic ΔE_ST computed in gas-phase with DFT/MRCI perform worse than the vertical ΔE_ST, e.g. (21) or (23), we decided to look further at the effect of including state-specific implicit solvation effects.

5.3 Effect of continuum solvation

Having discussed the choice of Hamiltonian and the effect of the chosen geometry, we now turn to the role of solvation. The high polarity of CT excited states leads to large electrostatic solvation effects, which lend themselves ideally to an implicit treatment.^40,41,64,102 In this work, we consider two approaches for combining DFT/MRCI with implicit continuum solvation: (i) BHLYP ground state solvation with the COSMO model as a zeroth order approximation for excited state solvation and (ii) approximate state-specific solvation with the reaction field from a ROKS-OT-ωB97M-V¹⁰³/def2-SVP+PCM calculation (see Section S3) of the targeted excited states.⁶⁴ It should be noted that, in DFT/MRCI, the implicit solvation effects enter in two ways: the potential from the reaction field surface charges directly enters the Kohn–Sham orbital energies in eqn (1). Indirectly, also the canonical DFT orbital shapes change, which affects the Fock matrix elements and electron repulsion integrals.

Starting again with the singlet–triplet gaps, we note that the approach (i) based on the ground-state solvation field leaves the vertical ΔE_ST values for most systems unchanged (Fig. S8). The less polar charge density distribution of the ground state appears to induce only a weak and uniform reaction field, which leads to only small but similar effects for S₁ and T₁. This observation aligns well with previous studies suggesting that polar CT states require a fully state-specific account of solvation.^{41,52,106,107}

To analyze this further, we compare in Fig. 4 the deviations of ΔE_ST from the experiment in gas phase computed vertically (a) and adiabatically (b) as well as with the approximate state-specific reaction field approach (c) from the aforementioned ROKS+PCM calculation. Furthermore, we characterize both, the singlet and triplet excited states, by the change of the static dipole moment relative to the electronic ground state χ_n (see SI, Section S5 for details) as either LE or CT states. Consequently, the ST-gaps fall into four categories depending on the state character of the involved singlet and triplet states, which allows us to disentangle the effects of the geometry, solvation, and state character a bit further.

Fig. 4 Difference between computed and experimental ΔE_ST for (a) vertical gaps (PBEh-3c S₀ gas phase geometries), (b) adiabatic gaps (ROKS/UKS-OT-LC-ωPBE-D3(BJ)/def2-SVP+PCM geometries) and (c) adiabatic gaps that include the excited-state reaction fields (from ROKS-OT-ωB97M-V/def2-SVP+PCM) are presented. Different state characters of the low-lying singlet and triplet excited states (CT or LE) are visualized. The states are characterized based on the dipole moment norm relative to the ground state (see SI, Section S5 for details). The nominal S₁–T₁ energy gap is marked for all cases (black box). Singlet–triplet ΔE_ST resulting from CT–CT (blue circles), CT–LE (purple diamonds), LE–CT (green triangles) and LE–LE (yellow crosses) are differentiated. Vertical grey lines are drawn for different structural groups of the emitters for visual guidance.

A comparison of the excited-state characters in Fig. 4a and b shows that several ΔE_ST categorizations change already due to the different geometries used. For oTE-DRZ (9), the CT–LE vertical gap changes to both singlet and triplet being characterized as CT states. Other LE–LE gaps are changing similarly to yield CT–CT gaps between the lowest excited states (e.g. 16, 18, 19), in alignment with the state character present in the theory level used for the geometry optimization. While these changes agree well with the state characters assigned in the corresponding ROKS+PCM calculations (see Table S4), the agreement of the adiabatic ΔE_ST with the experimental references is worsened. To make sense of this, it should be emphasized that the geometries were optimized with PCM solvation, while the (adiabatic) DFT/MRCI calculations performed at those are still gas-phase calculations. The worse ΔE_ST from these gas-phase DFT/MRCI calculations indicate that, though a geometry-induced change in state character is observed, a state-specific electronic relaxation effect is still missing in these calculations.

When including state-specific reaction fields from the selected ROKS+PCM to model the excited-state solvation (Fig. 4c), comparison to the adiabatic ΔE_ST obtained in the gas phase demonstrates that this approximate state-specific solvation treatment leads to improvements in previously bad cases (e.g., 8, 16, 21, 22, and 24) but, even more so, leads to several unphysically inverted ΔE_ST values in other cases (notably, 2, 3, 5, 7, 11, 13–16, 18, and 19). The more polar S₁ CT state appears to become overstabilized in this setting, which overcorrects the generally slightly too positive gaps observed in the adiabatic gas phase calculation (cf. Fig. 3 and 4b). This effect is most pronounced for TMCz-BO (2), where the ROKS+PCM reaction field of the CT singlet state leads to a strong stabilization of the corresponding state in DFT/MRCI. For the LE triplet state reaction field from ROKS+PCM, this does not nearly have the same effect for the DFT/MRCI-computed triplet state, which causes a completely unreasonable gap of −0.56 eV. To a lesser degree, similar behavior is observed for other systems, which can be rationalized from the fact that the CT character is typically less pronounced in the triplet state minimum compared to the singlet minimum structure. Ultimately, DFT/MRCI is able to capture the excited-state character in agreement with ROKS+PCM with and without solvation at the adiabatic excited state minimum geometries. While including state-specific solvation from ROKS+PCM leads to improvements in some cases, it overall leads to more severe outliers for the ΔE_ST values. Due to the empirical and not fully variational structure of the DFT/MRCI Hamiltonian (eqn (1)), it is difficult to pinpoint this behavior to one cause. Of course, one is the mixing of two theory levels here: the state-specific reaction field could not be generated with the existing DFT/MRCI method directly. Hence, we used a different, well-performing theory level (ROKS-OT-ωB97M-V/def2-SVP+PCM) for generation of the reaction field and used it to converge the BHLYP reference wavefunction in its presence. This is expected to lead to some imbalance, as the potential energy surfaces for any two theory levels may not be exactly parallel. This may be amplified by the empirical nature of the DFT/MRCI Hamiltonian: its design and parametrization have been partially adjusted to match experimental data in solution based on gas-phase calculations. When explicitly taking into account state-specific solvation effects via the reaction field, double counting of solvation effects might be taking place. Surprisingly, the vertical approximation in the gas-phase still performs best, likely due to a well-working error cancelation. As seen in Fig. 4, DFT/MRCI then produces different state characters compared to the adiabatic setting. It should be noted that this does not imply a change in the order of the adiabatic states within each spin manifold. Instead, the LE and CT states appear to be the same adiabatic state with dominant HOMO–LUMO contributions. Geometry relaxation or the presence of the reaction field then leads to the change in state character of that state.

Since the solvation energy of a CT excited and a less-polar ground state differs much more than that of two rather similar CT states, the effect of state-specific solvation for E_Em is expected to be more pronounced than for the previously discussed ΔE_ST. Fig. 5 shows emission energies including the dielectric continuum again via the state-specific reaction field from the ROKS-OT-ωB97M-V/def2-SVP+PCM calculation, next to the gas-phase and ground state (COSMO) solvation results. Additionally, we include non-equilibrium solvation effects (relaxation of only the fast polarization component during emission) via a ptSS-PCM correction only for E_Em, but not for ΔE_ST.

Fig. 5 Emission energies computed using DFT/MRCI with the R2018 Hamiltonian and short settings at the ROKS S₁ geometries. The DFT/MRCI anchor configuration (BHLYP/def2-SV(P)) was computed in gas phase (blue), using ground state equilibrium solvation conditions with the COSMO implicit solvation model (denoted “GS eq.”, orange) or the S₁ excited state reaction field without (denoted “ES eq.”, light purple) and with ptSS-PCM correction (“ES eq. + ptSS”, purple) obtained from ROKS-OT-LC-ωPBE/def2-SVP+PCM. An uncertainty of 0.20 eV for the experimental references is assumed and represented as grey shading.

We find that emission energies computed in gas phase mostly overestimate the experimental reference (Table 2, MD = 0.17 eV). As discussed before, this is expected due to missing stabilization of the excited CT singlet state from a matching reaction field. The introduction of ground-state COSMO solvation further increases this trend on average (MD = 0.20 eV), as now, the ground state becomes stabilized instead of the excited state. Including the approximate state-specific reaction field from the chosen ROKS+PCM theory level (purple) reduces the emission energies far below the experimental reference (MD = −0.42 eV). Hence, the aforementioned imbalance due to using the empirical DFT/MRCI Hamiltonian in combination with that excited-state reaction field, which led to larger but less systematic errors observed for ΔE_ST (Fig. 4), is now reflected clearly in systematically underestimated E_Em values. Only after also including the fast relaxation to the ground-state reaction field by the approximate ptSS-correction from the same ROKS calculation, most emission energies are within the experimental uncertainty, but still underestimated in almost all cases (MD = −0.15 eV). Part of the remaining underestimation with the ptSS treatment is likely due to similar incompatibilities of the solvation treatment and theory levels that we observed for ΔE_ST, so far. Still, the results align with previous studies of the solvent effects in the emission energies and the accuracy is comparable to the best TDA-DFT based approaches.⁴¹ Overall, this leads to similar errors in the E_Em values as the gas-phase calculations with both having opposite systematic shifts. It is noteworthy that for DA systems with a dimethylacridine donor group, we always observe an overestimated emission energy when using DFT/MRCI gas-phase calculations (molecules (3), (5), (11), and (14)–(18)) with molecules (16) and (18) being the strongest outliers. Particularly, the inclusion of solvation effects through the ptSS-PCM correction improves the emission energies for these systems with little or no underestimation tendency.

Table 2 Statistical measures (in eV) of the errors in the computed emission energy for the STGABS27 benchmark set using DFT/MRCI with different environment descriptions, short settings, and a BHLYP/def2-SV(P) anchor configuration compared to experimental references. Gas phase values, zeroth order ground state solvation conditions (“GS eq.”) described using COSMO, approximate excited-state solvation conditions (“ES eq.”) employing a ROKS-OT-LC-ωPBE/def2-SVP+PCM reaction field are compared. For the latter, values obtained by adding a ptSS-PCM correction (“ES eq. + ptSS”) are given as well. MD: mean deviation, MAD: mean absolute deviation, AMAX: absolute maximum deviation, RMSE: root mean square error, SD: corrected standard deviation (see Table S1 for definitions). Excited state S₁ structures were optimized at the ROKS/UKS-OT-LC-ωPBE-D3(BJ)/def2-SVP+SS-PCM (ROKS-OT-wB97M-V/def2-SVP+SS-PCM for ^a) level of theory

	MD	MAD	AMAX (system)	RMSE	SD
a The values refer to the OT-ωB97M-V/def2-SVP+PCM level of theory.b The values refer to the OT-LC-ωPBE/def2-SVP+SS-PCM level of theory.
Gasphase	0.17	0.21	0.60 (16)	0.26	0.37
GS eq.	0.19	0.24	0.70 (16)	0.30	0.37
ES eq.	−0.42	0.44	0.92 (19)	0.52	0.46
ES eq. + ptSS	−0.15	0.21	0.59 (19)	0.27	0.40
ROKS^a^41	0.04	0.13	0.32 (2)	0.16	0.15
TDA-DFT^b^41	−0.07	0.20	0.59 (12)	0.27	0.27

---

The observed impact of the state-specific solvation treatment for E_Em and ΔE_ST suggests that future redesigns of the DFT/MRCI may benefit from including a well-founded state-specific solvation treatment already during the parametrization. This way, the appreciated advantages of DFT/MRCI (multiple-state treatment, spin-pure states) could be leveraged when studying solvent-specific polarity differences for various excited states. As of now, the gas-phase treatment in combination with the R2018 Hamiltonian appears to perform best for the ΔE_ST gaps in the STGABS27, while leading to somewhat overestimated but fairly acceptable E_Em values. We could demonstrate DFT/MRCI's intriguingly good performance for the ΔE_ST gaps in the vertical approximation. In a “real-life” scenario, the lowest excited state and its minimum geometry are not known. In that case, DFT/MRCI can be used as a screening method to identify the lowest excited states and potentially small ΔE_ST gaps, without introducing molecule- or state-specific tuning as in (TD-)DFT methods. Given its ability to simultaneously describe multiple excited states, including doubly excited and charge transfer states, it can be used to identify the relevant adiabatic states involved in the TADF process, which are not limited to the lowest excited states,³² before moving to a more accurate state-specific OO-DFT scheme.

5.4 Comparison to other methods

Last, we want to compare DFT/MRCI to other well-performing computational methods for the STGABS27 set. Here, we consider the results from previous works using ROKS (OT-ωB97M-V/def2-SVP+PCM)⁶⁴ and TDA-DFT (OT-LC-ωPBE/def2-SVP+SS-PCM).^40,41 Both methods take into account state-specific solvation effects, an optimally-tuned density functional, and, in case of emission energies,⁴¹ a non-equilibrium solvation correction. We compare them with the best-performing and much simpler DFT/MRCI approach using the vertical approximation in gas phase. In Fig. 6, the DFT/MRCI data acquired in this work is compared to the theoretical estimates using ROKS and TDA-DFT published recently.^40,64

Fig. 6 Comparison of vertical ΔE_ST and adiabatic E_Em computed using DFT/MRCI(R2018, short) with a BHLYP/def2-SV(P) anchor configuration (blue) in the gas phase, adiabatic TDA-DFT (OT-LC-ωPBE+SS-PCM, red line, upper triangles), and ROKS (OT-ωB97M-V+PCM, black line, crosses). Reference data for ΔE_ST and E_Em from ref. 40, 41 and 64. Uncertainties of 0.05 eV and 0.20 eV are given as grey shading.

The first thing to note from Fig. 6 is that state-specific ROKS+PCM clearly performs best, as mentioned already in the introduction. The computed ΔE_ST values are excellent with only two errors above 0.05 eV and it is also the best performing method for the emission energies. The accuracy of this OO-DFT scheme is due to its native account for orbital relaxation and the consistent state-specific treatment of solvation effects.

Albeit less accurate than OO-DFT, DFT/MRCI shows an overall robust performance for both ΔE_ST and E_Em. Particularly in comparison to the best-performing TDA-DFT-based approach, which also solves a multi-state eigenvalue problem, the performance of DFT/MRCI for ΔE_ST in the vertical gas-phase approximation is surprisingly good.

Overall, using DFT/MRCI without any solvation treatment offers a comparably state-independent and robust description of the low-lying excited states in the investigated TADF emitters.

6 Conclusions

We studied the performance of DFT/MRCI for the energies of the lowest singlet and triplet states in TADF emitters. The results show that the DFT/MRCI (R2018 Hamiltonian) provides accurate singlet–triplet gaps (MAD of 0.06 eV), while other DFT/MRCI Hamiltonians presented by the Marian group perform similarly or slightly worse. Specifically, the R2018 Hamiltonian provides slightly better results for ΔE_ST using the computationally more efficient short settings of DFT/MRCI and a small def2-SV(P) basis set.

Notably, the most accurate description of ΔE_ST was obtained using DFT/MRCI with ground state geometries in the gas phase, i.e., the vertical approximation. Moving to adiabatic ΔE_ST gaps computed for ROKS+PCM-optimized structures and including dielectric relaxation through a ROKS+PCM reaction field, did not systematically improve the agreement with experiment compared to the vertical gas-phase approximation. In contrast, for some systems of the STGABS27 benchmark set, we observed a severe overstabilization of charge-transfer singlet states compared to the more local-in-character triplet states. Overall, this leads to underestimated E_Em values (MD of −0.15 eV, RMSE of 0.27 eV) and, in some DA systems, to erroneously inverted ΔE_ST gaps. Hence, despite being physically motivated, our findings suggest that the current DFT/MRCI (R2018) Hamiltonian design does not benefit systematically from including a state-specific solvation compared to a gas-phase treatment, indicating that some solvation effects are included in the parametrization.

Generally, our analysis of the low-lying excited singlet and triplet states showed that DFT/MRCI identifies the S₁–T₁ adiabatic energy gap to originate from CT–CT excited states in most cases, in line with results obtained using state-specific approaches. Emission energies are slightly overestimated using DFT/MRCI with an anchor configuration computed in gas phase (MD of 0.17 eV, RMSE of 0.26 eV). Still, the method performs reasonably accurate in this setting. Aside from that, we also briefly commented on the fact that structural changes, e.g., from considering different conformers may be an additional source of error when computing ΔE_ST for the studied TADF emitters.

Summarizing, we do find the best performance of DFT/MRCI for ΔE_ST gaps and the E_Em energies when computed in the vertical gas-phase approximation. In this configuration, it outperforms Tamm–Dancoff-approximated TD-DFT for ΔE_ST without the need for optimally-tuning and state-specific solvation, and approaches the accuracy of state-specific orbital-optimized DFT variants like ROKS/UKS+IEF-PCM (though never reaches it, cf. Tables 1 and 2).

It should be pointed out that DFT/MRCI has a significantly higher computational cost than TD-DFT and ΔUKS. Combined with the def2-SV(P) basis set a molecule with 90 atoms (system 16) requires already about 12 h of computation time on 16 CPU cores compared to minutes for ΔUKS (see SI). However, it offers great potential for the initial assessment of the excited states with no a priori knowledge about the molecules, because of its robust performance for local excitations, charge transfer excitations, and also double excitations⁶⁶ and the ability to simultaneously calculate multiple low-lying excited states and their transition properties. Therefore, we see its primary role as a robust and fairly reliable excited state exploration tool that can be applied at the ground state geometry, before moving to a more accurate molecule and state-specific treatment of the excited states preferably via orbital-optimized DFT.

With respect to future investigations, assessing its performance for more MR-TADF and the novel generation TADF emitters with inverted singlet–triplet gaps (INVEST) should be carried out, as the involved states were reported to have smaller structural relaxation and partial double excitation character.^25,36

Author contributions

M. Pauls: conceptualization (supporting); data curation (lead); formal analysis (lead); investigation (lead); methodology (equal); validation (lead); visualization (lead); writing – original draft (equal); writing – review & editing (equal); T. Froitzheim: conceptualization (supporting); formal analysis (supporting); investigation (supporting); methodology (equal); writing – original draft (equal); writing – review & editing (equal); funding acquisition (supporting); A. Torgashov: formal analysis (supporting); investigation (supporting); methodology (supporting); Jan-Michael Mewes: conceptualization (lead); formal analysis (supporting); investigation (supporting); writing – original draft (supporting); writing – review & editing (equal); Stefan Grimme: writing – original draft (supporting); writing – review & editing (equal) funding acquisition (supporting); Christoph Bannwarth: conceptionalization (lead); formal analysis (supporting); investigation (supporting); writing – original draft (supporting); writing – review & editing (equal); funding acquisition (lead).

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information is available, including the used statistical measures, further computational details, explanations of the state-specific reaction fields, basis set convergence, and timings, as well as details on the assignment of state characters. Output files from the calculations are also provided. See DOI: https://doi.org/10.1039/d5cp01869b

Acknowledgements

The authors thank Lukas Kunze for providing geometries from the STGABS27 benchmark set. Furthermore, we thank Gereon Feldmann for reading the manuscript and helpful discussion. We appreciate Dr Martin Kleinschmidt's support and technical advice for questions regarding the DFT/MRCI program. Christoph Bannwarth acknowledges funding by the Ministry of Culture and Science of the German State North Rhine-Westphalia (MKW) via the NRW Rückkehrprogramm. Thomas Froitzheim is most grateful to the “Fonds der Chemischen Industrie (FCI)” for financial support under a Kekulé scholarship.

Notes and references

1. C. W. Tang and S. A. VanSlyke, Appl. Phys. Lett., 1987, 51, 913–915.
2. Z. Yang, Z. Mao, Z. Xie, Y. Zhang, S. Liu, J. Zhao, J. Xu, Z. Chi and M. P. Aldred, Chem. Soc. Rev., 2017, 46, 915–1016.
3. H. Nakanotani, Y. Tsuchiya and C. Adachi, Chem. Lett., 2021, 50, 938–948.
4. T. Huang, W. Jiang and L. Duan, J. Mater. Chem. C, 2018, 6, 5577–5596.
5. A. Endo, K. Sato, K. Yoshimura, T. Kai, A. Kawada, H. Miyazaki and C. Adachi, Appl. Phys. Lett., 2011, 98, 083302.
6. H. Uoyama, K. Goushi, K. Shizu, H. Nomura and C. Adachi, Nature, 2012, 492, 234–238.
7. Q. Zhang, J. Li, K. Shizu, S. Huang, S. Hirata, H. Miyazaki and C. Adachi, J. Am. Chem. Soc., 2012, 134, 14706–14709.
8. Y. Tao, K. Yuan, T. Chen, P. Xu, H. Li, R. Chen, C. Zheng, L. Zhang and W. Huang, Adv. Mater., 2014, 26, 7931–7958.
9. S. Hirata, Y. Sakai, K. Masui, H. Tanaka, S. Y. Lee, H. Nomura, N. Nakamura, M. Yasumatsu, H. Nakanotani and Q. Zhang, et al., Nat. Mater., 2015, 14, 330–336.
10. F. B. Dias, T. J. Penfold and A. P. Monkman, Methods Appl. Fluoresc., 2017, 5, 012001.
11. J.-M. Mewes, Phys. Chem. Chem. Phys., 2018, 20, 12454–12469.
12. C. M. Marian, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 187–203.
13. C. M. Marian, J. Phys. Chem. C, 2016, 120, 3715–3721.
14. I. Lyskov and C. M. Marian, J. Phys. Chem. C, 2017, 121, 21145–21153.
15. T. Hatakeyama, K. Shiren, K. Nakajima, S. Nomura, S. Nakatsuka, K. Kinoshita, J. Ni, Y. Ono and T. Ikuta, Adv. Mater., 2016, 28, 2777–2781.
16. S. Madayanad Suresh, D. Hall, D. Beljonne, Y. Olivier and E. Zysman-Colman, Adv. Funct. Mater., 2020, 30, 1908677.
17. J.-M. Teng, Y.-F. Wang and C.-F. Chen, J. Matter. Chem. C, 2020, 8, 11340–11353.
18. K. Shizu and H. Kaji, Commun. Chem., 2022, 5, 53.
19. K. R. Naveen, P. Palanisamy, M. Y. Chae and J. H. Kwon, Chem. Commun., 2023, 59, 3685–3702.
20. S. Seifermann, Organic molecules for optoelectronic devices, US Pat., US11545632B2, 2023.
21. X.-F. Luo, X. Xiao and Y.-X. Zheng, Chem. Commun., 2024, 60, 1089–1099.
22. G. Ricci, J.-C. Sancho-García and Y. Olivier, J. Mater. Chem. C, 2022, 10, 12680–12698.
23. P. de Silva, J. Phys. Chem. Lett., 2019, 10, 5674–5679.
24. L. Tučková, M. Straka, R. R. Valiev and D. Sundholm, Phys. Chem. Chem. Phys., 2022, 24, 18713–18721.
25. L. Kunze, T. Froitzheim, A. Hansen, S. Grimme and J.-M. Mewes, J. Phys. Chem. Lett., 2024, 15, 8065–8077.
26. Y. Bu and Q. Peng, J. Phys. Chem. C, 2023, 127, 23845–23851.
27. K. Zhao, Ö. H. Omar, T. Nematiaram, D. Padula and A. Troisi, J. Mater. Chem. C, 2021, 9, 3324–3333.
28. J. M. Jacob, P. K. Samanta and M. K. Ravva, New J. Chem., 2023, 47, 10552–10563.
29. P. Ulukan, E. E. Bas, R. B. Ozek, C. Dal Kaynak, A. Monari, V. Aviyente and S. Catak, Phys. Chem. Chem. Phys., 2022, 24, 16167–16182.
30. K.-L. Woon, P. A. Nikishau and G. Sini, Adv. Theory Simul., 2022, 5, 2200056.
31. J. Sanz-Rodrigo, Y. Olivier and J.-C. Sancho-Garcia, Molecules, 2020, 25, 1006.
32. A. L. Schleper, K. Goushi, C. Bannwarth, B. Haehnle, P. J. Welscher, C. Adachi and A. J. Kuehne, Nat. Commun., 2021, 12, 6179.
33. F. Majer, L. Roß, A. L. Respondek, C. Bannwarth and A. J. C. Kuehne, ChemPhotoChem, 2024, 8, e202400141.
34. R. Pollice, B. Ding and A. Aspuru-Guzik, Matter, 2024, 7, 1161–1186.
35. Sanyam, R. Khatua and A. Mondal, J. Phys. Chem. A, 2023, 127, 10393–10405.
36. L. Kunze, A. Hansen, S. Grimme and J.-M. Mewes, J. Phys. Chem. Lett., 2025, 16, 1114–1125.
37. E. Runge and E. K. Gross, Phys. Rev. Lett., 1984, 52, 997.
38. S. Hirata and M. Head-Gordon, Chem. Phys. Lett., 1999, 314, 291–299.
39. A. Dreuw and M. Head-Gordon, Chem. Rev., 2005, 105, 4009–4037.
40. T. Froitzheim, S. Grimme and J.-M. Mewes, J. Chem. Theory Comput., 2022, 18, 7702–7713.
41. T. Froitzheim, L. Kunze, S. Grimme, J. M. Herbert and J.-M. Mewes, J. Phys. Chem. A, 2024, 128, 6324–6335.
42. R. Khatri and B. D. Dunietz, J. Phys. Chem. C, 2024, 129, 436–446.
43. G. Ricci, E. San-Fabián, Y. Olivier and J.-C. Sancho-Garca, ChemPhysChem, 2021, 22, 553–560.
44. A. Derradji, D. Valverde, É. Brémond, Á. J. Pérez-Jiménez, Y. Olivier and J. C. Sancho-Garca, J. Phys. Chem. C, 2024, 128, 18313–18327.
45. O. Christiansen, H. Koch and P. Jørgensen, Chem. Phys. Lett., 1995, 243, 409–418.
46. J. Schirmer, Phys. Rev. A: At., Mol., Opt. Phys., 1982, 26, 2395.
47. A. Trofimov and J. Schirmer, J. Phys. B: At., Mol. Opt. Phys., 1995, 28, 2299.
48. S. Grimme, J. Chem. Phys., 2003, 118, 9095–9102.
49. A. Hellweg, S. A. Grün and C. Hättig, Phys. Chem. Chem. Phys., 2008, 10, 4119–4127.
50. D. Hall, J. C. Sancho-Garca, A. Pershin, G. Ricci, D. Beljonne, E. Zysman-Colman and Y. Olivier, J. Chem. Theory Comput., 2022, 18, 4903–4918.
51. J.-M. Mewes, Z.-Q. You, M. Wormit, T. Kriesche, J. M. Herbert and A. Dreuw, J. Phys. Chem. A, 2015, 119, 5446–5464.
52. J.-M. Mewes, J. M. Herbert and A. Dreuw, Phys. Chem. Chem. Phys., 2017, 19, 1644–1654.
53. B. Lunkenheimer and A. Köhn, J. Chem. Theory Comput., 2013, 9, 977–994.
54. C. Hättig and A. Pausch, J. Phys. Chem. A, 2025, 129, 6155–6169.
55. G. M. Barca, A. T. Gilbert and P. M. Gill, J. Chem. Theory Comput., 2018, 14, 1501–1509.
56. I. Frank, J. Hutter, D. Marx and M. Parrinello, J. Chem. Phys., 1998, 108, 4060–4069.
57. A. T. B. Gilbert, N. A. Besley and P. M. W. Gill, J. Phys. Chem. A, 2008, 112, 13164–13171.
58. D. Hait and M. Head-Gordon, J. Chem. Theory Comput., 2020, 16, 1699–1710.
59. T. Kowalczyk, T. Tsuchimochi, P.-T. Chen, L. Top and T. Van Voorhis, J. Chem. Phys., 2013, 138, 164101.
60. T. Ziegler, A. Rauk and E. J. Baerends, Theor. Chim. Acta, 1977, 43, 261–271.
61. D. Hait and M. Head-Gordon, J. Phys. Chem. Lett., 2021, 12, 4517–4529.
62. D. Hait, K. J. Oosterbaan, K. Carter-Fenk and M. Head-Gordon, J. Chem. Phys., 2022, 156, 201104.
63. N. Bogo and C. J. Stein, Phys. Chem. Chem. Phys., 2024, 26, 21575–21588.
64. L. Kunze, A. Hansen, S. Grimme and J.-M. Mewes, J. Phys. Chem. Lett., 2021, 12, 8470–8480.
65. S. Grimme and M. Waletzke, J. Chem. Phys., 1999, 111, 5645.
66. C. M. Marian, A. Heil and M. Kleinschmidt, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1394.
67. M. Schreiber, M. R. Silva-Junior, S. P. A. Sauer and W. Thiel, J. Chem. Phys., 2008, 128, 134110.
68. C. Bannwarth, J. K. Yu, E. G. Hohenstein and T. J. Martínez, J. Chem. Phys., 2020, 153, 024110.
69. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
70. T. S. Costain, V. Ogden, S. P. Neville and M. S. Schuurman, J. Chem. Phys., 2024, 160, 224106.
71. I. Lyskov and C. M. Marian, J. Phys. Chem. C, 2017, 121, 21145–21153.
72. H. Miranda-Salinas, A. Rodriguez-Serrano, J. M. Kaminski, F. Dinkelbach, N. Hiromichi, Y. Kusakabe, H. Kaji, C. M. Marian and A. P. Monkman, J. Phys. Chem. C, 2023, 127, 8607–8617.
73. N. Lüdtke, J. Föller and C. M. Marian, Phys. Chem. Chem. Phys., 2020, 22, 23530–23544.
74. X.-F. Song, Z.-W. Li, W.-K. Chen, Y.-J. Gao and G. Cui, Inorg. Chem., 2022, 61, 7673–7681.
75. J. M. Kaminski, T. Böhmer and C. M. Marian, J. Phys. Chem. C, 2024, 128, 13711–13721.
76. M. N. Berberan-Santos and J. M. Garcia, J. Am. Chem. Soc., 1996, 118, 9391–9394.
77. S. Grimme, Chem. Phys. Lett., 1996, 259, 128–137.
78. D. R. Dombrowski, T. Schulz, M. Kleinschmidt and C. M. Marian, J. Phys. Chem. A, 2023, 127, 2011–2025.
79. A. Heil, M. Kleinschmidt and C. M. Marian, J. Chem. Phys., 2018, 149, 164106.
80. A. Heil and C. M. Marian, J. Chem. Phys., 2017, 147, 194104.
81. I. Lyskov, M. Kleinschmidt and C. M. Marian, J. Chem. Phys., 2016, 144, 034104.
82. S. Grimme, J. G. Brandenburg, C. Bannwarth and A. Hansen, J. Chem. Phys., 2015, 143, 054107.
83. T. Stein, L. Kronik and R. Baer, J. Am. Chem. Soc., 2009, 131, 2818–2820.
84. R. Baer, E. Livshits and U. Salzner, Annu. Rev. Phys. Chem., 2010, 61, 85–109.
85. M. A. Rohrdanz and J. M. Herbert, J. Chem. Phys., 2008, 129, 034107.
86. M. A. Rohrdanz, K. M. Martins and J. M. Herbert, J. Chem. Phys., 2009, 130, 054112.
87. E. Epifanovsky, A. T. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al., J. Chem. Phys., 2021, 155, 084801.
88. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
89. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
90. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
91. S. Miertuš, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117–129.
92. J. Tomasi, B. Mennucci and E. Cancès, J. Mol. Struct., 1999, 464, 211–226.
93. TURBOMOLE V7.6 2020, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007, available from https://www.turbomole.com.
94. F. Furche, R. Ahlrichs, C. Hättig, W. Klopper, M. Sierka and F. Weigend, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 91–100.
95. R. Ahlrichs, M. Bär, M. Häser, H. Horn and C. Kölmel, Chem. Phys. Lett., 1989, 162, 165–169.
96. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński and J. M. Yu, J. Chem. Phys., 2020, 152, 184107.
97. F. Neese, J. Comput. Chem., 2003, 24, 1740–1747.
98. O. Vahtras, J. Almlöf and M. Feyereisen, Chem. Phys. Lett., 1993, 213, 514–518.
99. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
100. https://github.com/grimme-lab/cefine, (last accessed 01.01.2025).
101. A. Klamt and G. Schüürmann, J. Chem. Soc., Perkin Trans. 2, 1993, 799–805.
102. J. M. Herbert, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2021, 11, e1519.
103. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2016, 144, 214110.
104. P. Freund, M. Pauls, D. Babushkina, T. Pickl, C. Bannwarth and T. Bach, J. Am. Chem. Soc., 2025, 147, 1434–1439.
105. N. Pflaum, M. Pauls, A. Kumar, R. J. Kutta, P. Nuernberger, J. Hauer, C. Bannwarth and T. Bach, J. Am. Chem. Soc., 2025, 147, 13893–13904.
106. R. Cammi, S. Corni, B. Mennucci and J. Tomasi, J. Chem. Phys., 2005, 122, 104513.
107. C. A. Guido, D. Jacquemin, C. Adamo and B. Mennucci, J. Chem. Theory Comput., 2015, 11, 5782–5790.

---