Given nowadays popularity of smartphones and many other potential applications, OLED materials and organic electronics are a very active field of research and industrial developments.

The following tutorial illustrates how advanced electroluminescence phenomena of OLED materials can be simulated and studied.
After a more general introduction to the topic, the approach presented below is largely based on the following recent publications:

OLED devices typically consist of a thin layer of a host material which is doped with an organic (or metalorganic) emitter material.
Applying an electric voltage across such semiconductor layers results in free electrons and electron holes, so-called charge carriers, travelling through the material.
When these charge carriers recombine, strongly bound electron-hole pairs, so-called excitons are formed, which are entirely localized on a single emitter molecule and thereby correspond to electronic excitations.

As free electrons (and holes) are equally likely to be created with spin up or spin down configurations, 75% of the formed excitons emerge in a triplet state (total spin \(S\) = 1) while 25% are found in a singlet state (\(S\) = 0).

This spin statistics negatively affects the quantum efficiency of an OLED material, hence on the proportion of electron-hole pairs that successfully undergo a radiative decay in to the ground state \(S_{\mathrm{0}}\) and emit a photon.
A direct de-excitation from \(T_{\mathrm{1}}\) to \(S_{\mathrm{0}}\) is spin-forbidden in a non-relativistic theory.
In the absence of any spin-orbit coupling effects in the dye, a radiative \(T_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) emission becomes very unlikely and light is only be emitted from a flourescent \(S_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) transition.

To address this problem a second generation of OLED emitters was developed involving OLED dyes containing heavy element atoms.
Electrons in the valence shells of such atoms (e.g. Iridium) are subject to significant spin-orbit coupling (SOC).
SOC effects in turn, permit \(T_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) transitions and facilitate so-called intersystem-crossings (ISC) e.g. from \(S_{\mathrm{1}}\) to \(T_{\mathrm{1}}\).
As a result the \(T_{\mathrm{1}}\) states become more populated.
From there the system undergoes a phosphorescent \(T_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) decay and emit light.

Fig. 3 Second generation OLED dyes rely on intersystem crossing to convert all exitons into a triplet state from which phosphorescent emissions occur.

While SOC-based emitters can reach quantum efficiencies of nearly 100%, the presence of long-living triplet states can cause chemical degradations.
This tendency increases with the energy difference to the ground state (and thus with the frequency of the emitted light).
As a result no stable, long-lasting OLED materials to deep blue light have been found as of yet.
Apart from that, SOC-based dyes involve relatively rare and expensive heavy elements which significantly increases their costs.

These disadvantages motivated the development of third generation OLED materials:
Dye molecules in which the states \(S_{\mathrm{1}}\) and \(T_{\mathrm{1}}\) are close in energy and (to some degree) vibrationally coupled, can exhibit a phenomenon known as thermally activated delayed flourescence (TADF).
Using emitter molecules specifically tailored to maximize TADF, this effect can be exploited to avoid the degradation problems of SOC-based emitters.
In addition, TADF occurs in purely organic dyes so that these third generation OLED emitters do not require expensive heavy metal elements.

TADF involves a reverse intersystem crossing (RISC) process, in which \(S_{\mathrm{1}}\) states are populated from \(T_{\mathrm{1}}\) levels.
This is followed again by a normal flourescent transition \(S_{\mathrm{1}} \rightarrow S_{\mathrm{0}}\) and allows for light emissions with nearly 100% quantum yields.

Fig. 4 The thermally assisted reverse intersystem crossing mechanism in third generation OLED materials converts exitons from triplet into singlet states from which a flourescent radiation process can occur.

RISC competes with other processes such as the \(S_{\mathrm{1}} \rightarrow T_{\mathrm{1}}\) intersystem crossing.
However, if non-radiative paths are negligible, \(k_{\mathrm{RISC}}\) represents the rate-determining factor for TADF.

TADF represents a rather complex process and \(k_{\mathrm{RISC}}\) has been shown to be strongly affected by various individual properties of a given dye material:

A small singlet-triplet gap \(\Delta E_{\mathrm{ST}} = E(S_{\mathrm{1}}) - E(T_{\mathrm{1}})\) is essential for the RISC process.
Indeed, the reverse intersystem crossing rate \(k_{\mathrm{RISC}}\) is very sensitive regarding \(\Delta E_{\mathrm{ST}}\), so that this quantity has to be predicted very precisely.
Because of the required accurate electronic structure descriptions of excited states and their tiny energy differences, computational studies of TADF processes can be quite challenging for practical applications.
As TADF dyes typically consist of hundreds of atoms, time-dependent density functional theory (TDDFT) often is the only computationally viable approach still accurate enough for modelling TADF processes.

To perform such TDDFT calculations, a suitable DFT approximation (i.e. exchange-correlation functional) needs to be chosen.
Indeed, this is nontrivial question in itself.
TDDFT results based on standard density functional approximations are inherently plagued by the self-interaction error.
For the charge-transfer (CT) excitations typically found in OLED emitters such self-interaction effects cause spurious artifacts in the energies and electronic structure of excited states.
While hybrid DFT methods ameliorate self-interaction errors, the excitation energies and electronic coupling parameters of TADF emitters are very sensitive regarding the amount of exact exchange incorporated in a given hybrid DFT approximation.

In range-separated hybrid (RSH) DFT methods the amount of exact exchange increases with the electron-electron distance.
In many RSH approximations this used to restore the correct long-range asymptotic behavior of the corresponding potential.
Their proper asymptotic behavior renders RSH methods potentially more accurate for the description of electronic excitations.
The range-separation of RSH methods is controlled by and additional parameter \(\omega\), which corresponds to the inverse distance of the midpoint of the switching function between short- and long-range terms.
The parameter \(\omega\) can be tuned for a given emitter material to improve the description further, say for example by minimizing the difference between the HOMO energy and the first vertical ionization energy.
This improves the accuracy compared to RSH descriptions with predetermined \(\omega\)-values.

Following Samanta et al., we use TDDFT calculations in the Tamm-Dancoff approximation (TDA) to describe the TADF process in 1,2,3,5-tetrakis(carbazol-9-yl)-4,6-dicyanobenzene (4CzIPN), a typical TADF emitter.

Fig. 5 1,2,3,5-tetrakis(carbazol-9-yl)-4,6-dicyanobenzene (4CzIPN)

Thereby we use the range-separated LCY-PBE functional, which employs a Yukawa potential to realize the switching function between long- and short-range terms.
The sharpness of the Yukawa potential is controlled by the parameter \(\gamma\), for which a value of 0.22 is used in the following.

Note

The value \(\gamma\) = 0.22 is the result of the RSH-tuning prescription of Karolewski et al. for the 4CzIPN molecule. Thereby the difference between the ionization potential and the HOMO energy is simultaneously minimized for a neutral and anionic system.
Please follow the advanced tutorial tuning the range separation for range separated hybrids.
Note, that \(\gamma\) is thereby optimized in the gas-phase. Prescriptions for tuning RSH functionals in solvents can be found in the publications of Joo et al. and Sun et al.

To model the TADF process of 4CzIPN properly, the geometry of the system needs to be relaxed for all three relevant electronic states, \(S_{\mathrm{0}}\), \(S_{\mathrm{1}}\), and \(T_{\mathrm{1}}\), respectively.
We can do this from GUI with ADF with the following settings:

A few more additional options are needed to define the excited states geometry optimizations.
In the case of the \(S_{\mathrm{1}}\)-optimization:

In the window with the previously created 4CzIPN_GSopt.adf file still open:

Properties → Excitations (UV/Vis), CD

Type of excitations: → SingletOnly

TDA: → Yes

Number of excitations: → 2

Calculate: → NTOs

Calculate: → Charge transfer descriptors

Properties → Excited State Geometry

Excitation: enter 1A

Spin type → Singlet

Save the input File → Save as... under e.g. 4CzIPN_S1opt.adf

Analogously for the \(T_{\mathrm{1}}\)-state:

Properties → Excitations (UV/Vis), CD

Type of excitations: → TripletOnly

TDA: → Yes

Number of excitations: → 2

Calculate: → NTOs

Calculate: → Charge transfer descriptors

Properties → Excited State Geometry

Excitation: enter 1A

Spin type → Triplet

Save the input File → Save as... under e.g. 4CzIPN_T1opt.adf

We are now ready to run all three optimization jobs.
Note that, for a system like 4CzIPN a excited states geometry optimization with RSH can take multiple hours.
For your convenience, the relaxed structures for the \(S_{\mathrm{0}}\), \(S_{\mathrm{1}}\), and \(T_{\mathrm{1}}\) states, respectively can be downloaded here:

After these optimizations are concluded, the respective minimum energies of the states can be extracted at the end of the logfile:

From there we can extract the following energy differences (neglecting vibrational contributions for the moment):

\(E_{\mathrm{S1}}-E_{\mathrm{S0}}\qquad\qquad\qquad\)

\(E_{\mathrm{T1}}-E_{\mathrm{S0}}\qquad\qquad\qquad\)

\(\Delta E_{\mathrm{ST}} = E_{\mathrm{S1}}-E_{\mathrm{T1}}\)

2.55 eV

2.39 eV

0.16 eV

The result is somewhat different from the results (\(\Delta E_{\mathrm{ST}}\) = 0.01 eV) of the original paper due to the different electronic structure description used therein.
Compared to the experimental value of \(\Delta E_{\mathrm{ST}}\) = 0.1 eV, the difference is however smaller and we will use our computed results of 0.16 eV in the subsequent calculations of this tutorial.

The output file provides the charge transfer descriptors, which are based on transition density matrices and were obtained from the prescription of Plasser and Lischka.

For the \(T_{\mathrm{1}}\) state we find a charge transfer (CT) character of 90% with an average electron-electron hole distance (R_HE) of 4.7 Å, while in the case of \(T_{\mathrm{1}}\) the CT character amounts to 94% with an average of 5 Å between the excited electron and electron hole.
Note that such high CT characters and electron-hole distances are a direct consequence of the design principles of typical OLED emitters featuring well-separated electron donor and acceptor groups.
In consequence, these excited states also exhibit a relatively small exchange energy, which in turn leads to a small \(\Delta E_{\mathrm{ST}}\) value.

Following Samanta et al. and to analyse the system further, we compute vertical excitations at the ground state geometry to obtain the absorption spectrum of 4CzIPN.
Regarding most aspects this is done analogously to the previous calculations.
The main differences are that we now treat the excitations in a spin-orbit description which allows for pure spin-states (i.e. no spin contamination).
Furthermore, the response of the solvation model is ignored as the solvent molecules cannot rearrange within the short timespan of a vertical excitation; a better approximation of this interaction requires the optical part of the solvents dielectric constant (see the corresponding section in the ADF Manual).

As discussed above, transitions between \(S\) and \(T\) states are formally forbidden within a non-relativistic framework.
Phosphorescence or intersystem crossings thus become only possible due to the presence of spin-orbit coupling (SOC) effects.
Compared to other interactions these couplings are typically small and can therefore be treated as a perturbation of the systems electronic structure mediated by \(\hat{H}_{\mathrm{SOC}}\), the SOC operator.
There exist various different approaches and approximation for \(\hat{H}_{\mathrm{SOC}}\).
See the recent review of Marian for an overview of spin-orbit methods in the context of excited states electronic structures and intersystem crossing rates.

In this case Fermi’s golden rule provides an expression for the reversed intersystem crossing rate

which will be used in the following.
\(|V_{\mathrm{SOC}}|^{2}\) is directly proportional to \(k_{\mathrm{RISC}}\) and is to be understood as coupling term of the matrix representation of \(\hat{H}_{\mathrm{SOC}}\).
We model the transition from \(T_{\mathrm{1}}\) to \(S_{\mathrm{1}}\), whereas \(T_{\mathrm{1}}\) actually consists of three individual states the system can assume; one for each possible value of the total angular momentum \(J=-1, 0, 1\).
\(|V_{\mathrm{SOC}}|^{2}\) is therefore an average of the three spin-orbit coupling matrix elements (SOCME) between \(S_{\mathrm{1}}\) and the \(T_{\mathrm{1}}\):

As \(|V_{\mathrm{SOC}}|^{2}\) is strongly affecting the reverse intersystem crossing rate, developing novel TADF emitters with larger SOCME values is an important part of many OLED related research activities.

The other important factor entering the \(k_{\mathrm{RISC}}\) rate is \(\rho_{\mathrm{FCWD}}\), the Frank-Condon-weighted density of states, which will be addressed in the corresponding section below.

The reverse intersystem crossing process itself occurs on a too short timescale for the systems geometry to readapt during the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition.
The factor \(|V_{\mathrm{SOC}}|^{2}\) in the above expression for \(k_{\mathrm{RISC}}\) is therefore computed at the optimized geometry of the \(T_{\mathrm{1}}\) state.

While \(|V_{\mathrm{SOC}}|^{2}\) can be interpreted as probability for the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition, \(\rho_{\mathrm{FCWD}}\) describes the thermokinetic barrier associated with this process.

Semiclassical Marcus theory (SCMT) provides a simple way to describe this thermokinetic barrier.
The potential energy surfaces of both states around their respective minimum structures are thereby described in terms of parabolic functions.

Fig. 6 Potential energy surfaces of the \(S_{\mathrm{1}}\) and \(T_{\mathrm{1}}\) states as quadratic functions centered at \(x_{\mathrm{S1}}\) and \(x_{\mathrm{T1}}\), respectively.
In Marcus theory the difference between \(E_{\mathrm{S1/S1}}\) and \(E_{\mathrm{S1/T1}}\) is used to estimate the thermokinetic barrier of the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition.

Assuming the same curvature of both states, i.e. \(E(T_{1}) = a\cdot x^{2}\) and \(E(S_{1}) = a\cdot (x_{} + x_{\mathrm{S1}} - x_{\mathrm{T1}})^{2} + \Delta E_{\mathrm{ST}}\), one finds the for the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) barrier height as the intersection of both parabolas:

The reorganization energy \(\lambda\) can thereby be considered as the energy required to bring the system in the \(S_{\mathrm{1}}\) to the minimum structure of the \(T_{\mathrm{1}}\) state.
Designating the energies of the \(S_{\mathrm{1}}\) state at \(x_{\mathrm{S1}}\) and \(x_{\mathrm{T1}}\) as \(E_{\mathrm{S1/S1}}\) and \(E_{\mathrm{S1/T1}}\), respectively, on can use

\[\lambda = E_{\mathrm{S1/T1}} - E_{\mathrm{S1/S1}}\]

to calculate \(k_{\mathrm{RISC}}\).
We can thereby set up the necessary calculations analogously to the vertical absorption calculations:

After running the calculations the required energies can be extracted again with ADFoutput.
Note, that the excited states energies are obtained from the ground state energy and the respective excitation energies.
These can be retrieved in ADFoutput as Total Bonding Energy under Properties → Bonding Energy Decomposition and as the first excitation with a significant oscillator strength under Response Properties → All Spin-Orbit Coupling Excitation Energies, respectively.

For 4CzIPN this results in a value of \(\lambda\) = \(E_{\mathrm{S1/T1}} - E_{\mathrm{S1/S1}}\) = 0.06 eV.
This result can now be used along with the previously obtained values \(\Delta E_{\mathrm{ST}}\) = 0.16 eV and \(|V_{\mathrm{SOC}}|^{2} =\) 8.2\(\cdot\)10-10 eV2 in the above expressions for \(\rho_{\mathrm{FCWD}}^{\mathrm{SCMT}}\) and \(k_{\mathrm{RISC}}\).

This yields a reverse intersystem crossing rate of \(k_{\mathrm{RISC}}\) = 2.3\(\cdot\)105 / s.
Given that the rather simplistic approximation for \(\rho_{\mathrm{FCWD}}\) and \(\lambda\), this result seems acceptable in comparison to the experimental reference of \(k_{\mathrm{RISC}}\) = 2.7\(\cdot\)106 / s.

Compared to that, Samanta et al. obtained higher rates for 4CzIPN (4.9\(\cdot\)106 – 1.8\(\cdot\)107/ s).
This difference can be rationalized by their usage of different DFT approximations, the significantly smaller value of \(\Delta E_{\mathrm{ST}}\) obtained thereby, and a different handling of the spin-orbit couplings.
We also note that Samanta et al. argued that for emitters like 4CzIPN \(\lambda\) can be replaced with an effective value between 0.1 and 0.2 eV.
When using \(\lambda\) = 0.1 eV instead of the previously calculated reorganization energy, the RISC rate rises to \(k_{\mathrm{RISC}}\) = 6.3\(\cdot\)105 / s.

While the semiclassical often provides acceptably accurate results, it relies on some strong assumptions and the real situation turns out more intricate:
First, the potential energy surfaces around both minima are described by the entirety of vibrational normal modes.
Each normal mode defines a separate harmonic oscillator system which, in a quantum mechanical description, leads to a separate set of vibrational eigenfunctions (quanta).
The Franck-Condon principle then states that the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition requires a significant overlap between nuclear wavefunctions of both electronic states.

Fig. 7 The Franck-Condon principle describes \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition rate in terms of the overlap between the nuclear vibronic wavefunctions of the initial and final state, respectively.

One can distinguish between two limiting cases: the quanta low frequency vibrations from shallow modes and from rearranging solvent molecules (see below) typically form an essentially continuous spectrum.
Pairs of overlapping quanta from low frequency modes will therefore occur at energies close to the intersection of the two parabolic potential curves, which is leads to results equivalent to those of semiclassical Marcus theory.
Within semiclassical Marcus theory the effect of these modes can be described in terms of a Marcus reorganization energy \(\lambda_{\mathrm{M}}\).
If the transition is dominated by high frequent modes, nuclear tunnelling effects and the overlap of individual pairs of quanta are prevalent.
The effect of such modes consists mainly in a renormalization of the coupling parameters.

An extension of SCMT, Marcus-Levich-Jortner theory includes includes both of these limiting cases.
The high frequency modes can thereby be projected onto the transition path which summarizes them into a single effective mode with frequency \(\omega_{eff}\) and gives rise to the following expression:

Thereby, \(S_{eff}\) denotes the Huang-Rhys factor of this effective mode, which is a measure for the strength of the electron-phonon coupling.
\(S_{eff}\) and the corresponding frequency \(\omega_{eff}\) are obtained from the Huang-Rhys factors and frequencies of the individual high frequency normal modes.

In some cases, the DFTB method can be used as a more approximate but efficient alternative to compute the normal modes.
The FCF program of the Amsterdam Modeling Suite can then be used to calculate the \(S_{i}\) factors which then lead to the \(S_{eff}\) and \(\omega_{eff}\) values.
The following procedure illustrates how this can be achieved with the Amsterdam Modeling Suite.

Warning

The calculation of Huang-Rhys factors crucially depends on the quality of the normal modes.
Their computation from TD-DFTB excited states gradient must be carefully tested to determine their validity for a given dye system.

After completing both calculations, the binary result files dftb.rkf containing the normal modes are found in the respective directories <jobname>.results.
Create a new directory and copy these rkf files naming them e.g. S1_optFreqDFTB.rkf and T1_optFreqDFTB.rkf.
Based on these results the program-tool fcf can be invoked to compute \(S_{eff}\) and \(\omega_{eff}\) values using the following bash script:

We are interested in the \(T_{\mathrm{1}}\rightarrow S_{\mathrm{1}}\) transition and are therefore using the results from the first state.
The mode-resolved Huang-Rhys factors \(S_{i}\) are obtained as the square of the above Electron-Phonon coupling parameters.
Using the above expressions for \(S_{eff}\) and \(\omega_{eff}\), we then sum over all modes with \(\omega_{i} >\) 1000 cm-1.
This results in \(S_{eff}\) = 0.837 and \(\omega_{eff}\) = 1498 cm-1.

These quantities are then used in the above expression for \(\rho_{\mathrm{FCWD}}^{\mathrm{SCMT}}\) (with \(\lambda_{\mathrm{M}}\) = \(\lambda\) = 0.06 eV) whereas the summation is stopped at \(n\) = 5.
With \(|V_{\mathrm{SOC}}|^{2} =\) 8.2\(\cdot\)10-10 eV2 this results in a reverse intersystem crossing rate that is essentially identical to the result from semiclassical Marcus theory of \(k_{\mathrm{RISC}}\) = 2.3\(\cdot\)105 / s.
This supports the finding of Samanta et al. that quantum tunneling effects are negligible in the reverse intersystem crossing process of 4CzIPN and that a classical treatment suffices for its description.

The approach presented above included solvent effects implicitly, as the COSMO model was used in the calculations.

Besides examining vibrational modes that simultaneously involve both, the OLED emitter and explicit solvent molecules, and pursuing analogously to the approach for intramolecular vibrations discussed above,
many more explicit treatments have been suggested in the literature.
Outersphere vibronic contributions can for example be obtained in terms of dielectric response functions (Rühle et al.), polarizable force fields (McMahon and Troisi), or from QM/MM approaches (Norton and Brédas).