Thermally Activated Delayed Fluorescence (TADF)
***********************************************

ADF Spectroscopy UV/Vis Organic electronics OLEDs Vertical absorption
Phosphorescence Intersystem crossing Spin-orbit coupling Marcus theory

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.

   [image]

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 publications:

P.K. Samanta, D. Kim, V. Coropceanu, J.-L. Brédas *Up-Conversion
Intersystem Crossing Rates in Organic Emitters for Thermally Activated
Delayed Fluorescence: Impact of the Nature of Singlet vs Triplet
Excited States*, J. Am. Chem. Soc. 139, 4042-4051 (2017).

A broader overview on the topic can be found for example in

Y. Olivier, B. Yurash, L. Muccioli, G. D’Avino, O. Mikhnenko, J. C.
Sancho-García, C. Adachi, T.-Q. Nguyen, D. Beljonne *Nature of the
singlet and triplet excitations mediating thermally activated delayed
fluorescence*, Phys. Rev. Materials 1, 075602 (2017).

and

Y. Olivier, J.-C. Sancho-García, L. Muccioli, G. D’Avino, D. Beljonne
*Computational Design of Thermally Activated Delayed Fluorescence
Materials: The Challenges Ahead*, J. Phys. Chem. Lett. 9, 6149-6163
(2018).

Note:

  This tutorial illustrates the modeling of TADF on the example of a
  realistic dye molecule. Expect several hours of computation time for
  each single-point evaluation.


General Remarks on Modelling OLED Emitters
==========================================


Electronic Structure of OLED Materials
--------------------------------------

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,
traveling 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).

   [image]Mechanism of first generation OLED emitters involving a
   fluorescent decay path with a 25% maximum theoretical quantum
   yield.

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 into 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 emitted from a
fluorescent 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.

   [image]Second generation OLED dyes rely on intersystem crossing to
   convert all excitons 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
fluorescence (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 fluorescent transition S_{\mathrm{1}}
\rightarrow S_{\mathrm{0}} and allows for light emissions with nearly
100% quantum yields.

   [image]The thermally assisted reverse intersystem crossing
   mechanism in third generation OLED materials converts excitons from
   triplet into singlet states from which a fluorescent 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:

* singlet-triplet gap \Delta E_{\mathrm{ST}} = E(S_{\mathrm{1}}) -
  E(T_{\mathrm{1}})

* spin-orbit interactions

* vibrational couplings between different electronic states

* solvent effects

Each of these aspects will be addressed in the following.


Computational Description of TADF 1: Electronic Structure
=========================================================

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 modeling TADF processes.

To perform such TDDFT calculations, a suitable DFT approximation (i.e.
exchange-correlation functional) needs to be chosen. Indeed, this is a
nontrivial question in itself as 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 is 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 an 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.


Excited States Geometry Optimizations
-------------------------------------

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.

   [image]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
  see the Python example 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:

   Download the 4CzIPN structure "here"
   Start AMSinput
   Load the structure with **File → Import Coordinates...**

   Select the following in the **Main** panel:
   **Task → Geometry Optimization**
   **XC Functional → Range Separated → LCY PBE**
   Click on [image: MoreBtn][image] right of the **XC Functional** field
   Enter **Gamma** "0.22"
   Return to the **Main** panel
   Make sure that **Relativity → Scalar**
   **Basis set → TZP**
   **Frozen core → None**
   **Numerical quality → Good**
   Click on the **Model** panel and select **Solvation**
   **Solvation method → COSMO**
   **COSMO solvent → Toluene**
   Save the input **File → Save as...** under e.g. *4CzIPN_GSopt.ams*
   Run the calculation **File → Run**

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.ams* file still open:
   Click in the panel bar on **Properties** and select **Excitations (UV/Vis), CD**
   **Type of excitations → SingletOnly**
   **TDA → Yes**
   Enter the **Number of excitations** "2"
   Tick the **Calculate → Charge transfer descriptors** box
   Click on **Properties** panel and select **Excited State Geometry**
   Enter at the **Excitation** "1A"
   **Spin type → Singlet**
   Save the input **File → Save as...** under e.g. *4CzIPN_S1opt.ams*
   Run the calculation **File → Run**

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

   Click on the **Properties** panel and select **Excitations (UV/Vis), CD**
   **Type of excitations → TripletOnly**
   **TDA → Yes**
   Enter the **Number of excitations** "2"
   Tick the **Calculate → Charge transfer descriptors** box
   Click on the **Properties** panel and select **Excited State Geometry**
   **Excitation** enter "1A"
   **Spin type → Triplet**
   Save the input **File → Save as...** under e.g. *4CzIPN_T1opt.ams*
   Run the calculation **File → Run**

We are now ready to run all three optimization jobs. Note that, for a
system like 4CzIPN an 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:

See also:

  * "TADF_4CzIPN_GS.xyz" 4CzIPN ground state structure

  * "TADF_4CzIPN_S1.xyz" 4CzIPN S_{\mathrm{1}} structure

  * "TADF_4CzIPN_T1.xyz" 4CzIPN T_{\mathrm{1}} structure

After these optimizations are concluded, the respective minimum
energies of the states can be extracted at the end of the logfile. For
the ground state,

   <Nov12-2020> <03:56:22>   Bond Energy         -33.76473746 a.u.

For excited states the end of the logfile looks like this (the T1
state in this case):

   <Nov13-2020> <04:18:25>  Total excited state energy:                 -33.67697839 Hartree

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

+-----------------------------------+-----------------------------------+-----------------------------------+
|                                   | Hartree                           | eV                                |
+-----------------------------------+-----------------------------------+-----------------------------------+
| E_{\rm{GS}}                       | -33.764737                        | -918.79                           |
+-----------------------------------+-----------------------------------+-----------------------------------+
| E_{\rm{S1}}                       | -33.671259                        | -916.24                           |
+-----------------------------------+-----------------------------------+-----------------------------------+
| E_{\rm{T1}}                       | -33.676978                        | -916.40                           |
+-----------------------------------+-----------------------------------+-----------------------------------+
| E_{\rm{S1}} - E_{\rm{GS}}         | 0.093478                          | 2.54                              |
+-----------------------------------+-----------------------------------+-----------------------------------+
| E_{\rm{T1}} - E_{\rm{GS}}         | 0.087759                          | 2.39                              |
+-----------------------------------+-----------------------------------+-----------------------------------+
| \Delta E_{\rm{ST}} = E_{\rm{S1}}  | 0.005719                          | 0.16                              |
| - E_{\rm{T1}}                     |                                   |                                   |
+-----------------------------------+-----------------------------------+-----------------------------------+

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 according to
the prescription of Plasser and Lischka, which are obtained from
transition density matrices. After opening the results with AMSoutput,
these descriptors can be found by navigating to **Section → Exciton
Energies**:

   [image]The charge transfer descriptors for the T_{\mathrm{1}} state

For the T_{\mathrm{1}} state we find a charge transfer (CT) character
of 82% with an average electron-electron hole distance (R_HE) of 4.6
Å, while in the case of S_{\mathrm{1}} the CT character amounts to 91%
with an average of 5.0 Å 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 contributes most to the small \Delta
E_{\mathrm{ST}} value computed above from TDDFT excitation energies.


Vertical Absorption
-------------------

Following Samanta et al. and to analyze the system further, we compute
vertical excitations at the ground state geometry to obtain the
absorption spectrum of 4CzIPN. Note, that these results do not enter
the calculation of the RISC rate but rather serve as an additional
assessment of the accuracy of the description employed here. 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).

   Open a new instance of AMSinput
   Load the previously optimized GS structure with **File → Import Coordinates...**
   In the **Main** panel:
   **Task → Single Point**
   **XC Functional → Range Separated → LCY PBE**
   Click on [image: MoreBtn][image] right of the **XC Functional** field
   Enter **Gamma** "0.22"
   Return to the **Main** panel
   Make sure that **Relativity → Scalar**
   **Basis set → TZP**
   **Frozen core → None**
   **Numerical quality → Good**
   Click on the **Model** panel and select **Solvation**
   **Solvation method → COSMO**
   **COSMO solvent → Toluene**
   Click on the **Details** panel and select **COSMO**
   Untick **Include response in TDDFT**
   Click on the **Properties** panel and select **Excitations (UV/Vis), CD**
   **Type of excitations → Spin-Orbit (Perturbative)**
   **TDA → Yes**
   Enter the **Number of excitations** "5"
   Tick the **Calculate → Charge transfer descriptors** box
   Save the input **File → Save as...** under e.g. *4CzIPN_GS_VA.ams*
   Run the calculation **File → Run**

After completing this calculation, we open the results with AMSoutput
and search for the results with **Response Properties → All Spin-Orbit
Coupling Excitation Energies**:

   [image]

Recalling the spin-forbidden transitions appear with near vanishing
oscillator strengths in an SOC treatment, we note the following from
these results:

* Excitation no. 8 is the lowest with a relevant oscillator strength
  f, it corresponds to S_{\mathrm{0}} \rightarrow S_{\mathrm{1}}.

* The S_{\mathrm{0}} \rightarrow S_{\mathrm{1}} excitation energy of
  2.77 eV compares well with the experiment value, 2.85 eV.

* At the GS geometry the T_{\mathrm{1}} and the T_{\mathrm{2}} state
  lie below S_{\mathrm{1}}.


Computational Description of TADF 2: Spin-Orbit Coupling
========================================================

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 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

   k_{\mathrm{RISC}} = \frac{2\pi}{\hbar}\left|V_{\mathrm{SOC}}\right
   |^{2}\times\rho_{\mathrm{FCWD}}

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}}, where
T_{\mathrm{1}} actually consists of three individual states the system
can assume, the three spin sublevels, which are, for example,
eigenfunctions of S_z, with eigen values m_S=-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}}:

   \left|V_{\mathrm{SOC}}\right|^{2} = \frac{1}{3}\left|\left\langle 
   S_{\mathrm{1}}\left|\hat{H}_{\mathrm{SOC}}\right|T_{\mathrm{1}}\ri
   ght\rangle\right|^{2} = \frac{1}{3}\!\!\!\sum\limits_{\
   m_S=-1,0,1}\!\!\!\left\langle S_{\mathrm{1}}\left|\hat{H}_{\mathrm
   {SOC}}\right|T_{\mathrm{1}}^{m_S}\right\rangle\left\langle S_{\mat
   hrm{1}}\left|\hat{H}_{\mathrm{SOC}}\right|T_{\mathrm{1}}^{m_S}\rig
   ht\rangle^{*}

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.


Calculating Spin-Orbit Couplings
--------------------------------

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.

   Open a new instance of AMSinput
   Load the previously optimized T1 structure with **File → Import Coordinates...**
   In the **Main** panel:
   **Task → Single Point**
   **XC-Functional → Range Separated → LCY PBE**
   Click on [image: MoreBtn][image] right of the **XC-Functional** field
   Enter **Gamma** "0.22"
   Return to the **Main** panel
   Make sure that **Relativity → Scalar**
   **Basis set → TZP**
   **Frozen core → None**
   **Numerical quality → Good**
   Click on the **Model** panel and select **Solvation**
   **Solvation method → COSMO**
   **COSMO solvent → Toluene**
   Click on the **Details** panel and select **COSMO**
   Make sure that the **Include response in TDDFT** box is ticked
   Click on the **Properties** panel and select **Excitations (UV/Vis), CD**
   **Type of excitations → Spin-Orbit (Perturbative)**
   **TDA → Yes**
   Enter the **Number of excitations** "5"
   Tick the **Calculate → Charge transfer descriptors** box
   Save the input **File → Save as...** under e.g. *4CzIPN_SOCME.ams*

After the calculation is completed, open the output file with
AMSoutput. In the AMSoutput the SOCME results can be found with
**Response Properties → Spin-orbit couplings**:

   [image]

The above root mean square coupling element |\langle
S_{\mathrm{1}}|\hat{H}_{\mathrm{SOC}}|T_{\mathrm{1}}\rangle| is then
divided by three to represent the average over the three triplet
states (see above). From the |\langle
S_{\mathrm{1}}|\hat{H}_{\mathrm{SOC}}|T_{\mathrm{1}}\rangle| value of
0.4 cm^-1 we then obtain |V_{\mathrm{SOC}}|^{2} = (0.4 cm^-1)^2 / 3 =
8.2\cdot10^-10 eV^2.


Computational Description of TADF 3: Vibrations
===============================================

Recall the aforementioned expression for the reverse intersystem
crossing rate:

   k_{\mathrm{RISC}} = \frac{2\pi}{\hbar}\left|V_{\mathrm{SOC}}\right
   |^{2}\times\rho_{\mathrm{FCWD}}

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.


Marcus Theory
-------------

Classical Marcus theory (CMT) 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.

   [image]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
T_{\mathrm{1}}\rightarrow S_{\mathrm{1}} barrier height as the
intersection of both parabolas:

   \lambda^{\ddagger} = \frac{\left(\Delta E_{\mathrm{ST}} +
   \lambda\right)^{2}}{4\lambda}

Using this barrier height in an Arrhenius-type expression yields then
the RISC rate:

   \rho_{\mathrm{FCWD}}^{\mathrm{CMT}} = \frac{1}{\sqrt{4 \pi \lambda
   k_{\mathrm{B}} T}} \exp\!\left(\!-\frac{\left(\Delta
   E_{\mathrm{ST}} + \lambda\right)^2}{4 \lambda k_{\mathrm{B}} T}
   \right)

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, one 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:

   Open a new instance of AMSinput
   Load the previously optimized S1 structures with **File → Import Coordinates...**
   In the **Main** panel:
   **Task → Single Point**
   **XC-Functional → Range Separated → LCY PBE**
   Click on the [image: MoreBtn][image] right of the **XC-Functional** field
   Enter **Gamma:** "0.22"
   Return to the **Main** panel
   Make sure that **Relativity → Scalar**
   **Basis set → TZP**
   **Frozen core → None**
   **Numerical quality → Good**
   Click on the **Model** panel and select **Solvation**
   **Solvation method → COSMO**
   **COSMO solvent → Toluene**
   Click on the **Details** panel and select **COSMO**
   Untick **Include response in TDDFT**
   Click on the **Properties** panel and select **Excitations (UV/Vis), CD**
   **Type of excitations → Spin-Orbit (Perturbative)**
   **TDA → Yes**
   Enter the **Number of excitations** "5"
   Tick the **Calculate → Charge transfer descriptors** box
   Save the input **File → Save as...** under e.g. *4CzIPN_S1_S1.ams*

Analogously for the S1/T1 calculation we set:

   Load the previously optimized T1 structures with **File → Import Coordinates...**
   Repeat the above steps
   Save the input **File → Save as...** under e.g. *4CzIPN_S1_T1.ams*

After running the calculations the required energies can be extracted
again with AMSoutput. Note, that the excited states energies of
E_{\mathrm{S1/S1}} and E_{\mathrm{T1/S1}} are obtained as sum of the
ground state energy and the respective excitation energy. These can be
retrieved in AMSoutput 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. So,
E_{\mathrm{S1/S1}} and E_{\mathrm{S1/T1}} can be calculated by
addition of its first singlet excitation energy to its bond energy.
E_{\mathrm{S1/T1}} = -918.591 eV (E_{\mathrm{GS/T1}}) + 2.437 eV
(excitation with significant oscillator strength) = -916.154 eV.
E_{\mathrm{S1/S1}} = -918.501 eV (E_{\mathrm{GS/S1}}) + 2.271 eV
(excitation with significant oscillator strength) = -916.230 eV.

For 4CzIPN, this results in a value of \lambda = E_{\mathrm{S1/T1}} -
E_{\mathrm{S1/S1}} = -916.154 eV -  -916.230 eV = 0.076 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\cdot10^-10 eV^2 in the above expressions for
\rho_{\mathrm{FCWD}}^{\mathrm{CMT}} and k_{\mathrm{RISC}} to yield at
a temperature of T = 298 K:

   \frac{2\pi}{\hbar}\left|V_{\mathrm{SOC}}\right|^{2} =
   \mbox{7.826}\cdot\mbox{10}^{6} \mbox{eV/s}\\
   \rho_{\mathrm{FCWD}}^{\mathrm{CMT}} =
   \mbox{6.696}\cdot\mbox{10}^{-3} /\mbox{eV}

and thus a reverse intersystem crossing rate of k_{\mathrm{RISC}} =
5.24\cdot10^4 / s. This result underestimates the crossing rate
compared to the experimental reference of k_{\mathrm{RISC}} =
2.7\cdot10^6 / s. The reason for this might be due to small
discrepancies in either \Delta E_{\mathrm{ST}} or \lambda which get
amplified by the exponential term in
\rho_{\mathrm{FCWD}}^{\mathrm{CMT}}.

Compared to that, Samanta et al. obtained higher rates for 4CzIPN
(4.9\cdot10^6 – 1.8\cdot10^7/ 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}} = 8.8\cdot10^5 / s.


Franck-Condon Principle and Marcus-Levich-Jortner Theory
--------------------------------------------------------

While the classical Marcus theory 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.

   [image]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 leads to results equivalent to those
of classical Marcus theory. Within classical 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
frequency modes, nuclear tunneling 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 CMT, Marcus-Levich-Jortner theory, 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:

   \rho_{\mathrm{FCWD}}^{\mathrm{MLJT}} = \frac{1}{\sqrt{4 \pi
   \lambda_{\mathrm{M}} k_{\mathrm{B}} T}} \sum\limits_{n =
   0}^{\infty}\exp\left(-S_{eff}\right)\frac{S_{eff}^{n}}{n!}
   \exp\!\left(\!-\frac{\left(\Delta E_{\mathrm{ST}} +
   n\hbar\omega_{eff} + \lambda_{\mathrm{M}}\right)^2}{4
   \lambda_{\mathrm{M}} k_{\mathrm{B}} T} \right)

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.

   S_{eff} = \sum\limits_{i}S_{i} \qquad \omega_{eff} =
   \sum_{i}\frac{\omega_{i}S_{i}}{S_{eff}}

Whereas modes with \omega_{i} > 1000 cm^-1 are typically included in
the summation.

Samanta et al. computed S_{i} factors from mode resolved
reorganization energies \lambda_{i} = k_{i}\Delta Q_{i}^{2}/2, which
in turn results from the frequencies k_{i} = \omega_{i}^{2} and \Delta
Q_{i}, the normal mode projection onto the transition path.

For a large emitter molecule like 4CzIPN however, the calculation of
the full Hessian matrix turns out to be computationally too expensive.


Effective Modes and Huang-Rhys Factors from DFTB and FCF
--------------------------------------------------------

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.

   Open a new instance of AMSinput
   Load the previously optimized S_{\mathrm{1}} structure via **File → Import Coordinates**
   In the panel bar click [image: ADFPanel][image] **→** [image: DFTBPanel][image]
   **Task → Geometry Optimization**
   Select the **Frequencies** box
   **Model → SCC-DFTB**
   **Parameter directory → DFTB.org/mio-1-1** (or similar)
   **Properties → Excitations (UV/Vis)**
   **Type of excitations → Singlet**
   Enter the **Number of excitations** "5"
   Calculate the excited state gradient for Excitation number "1"
   Save the input **File → Save as...**

Analogously, for the T_{\mathrm{1}} state:

   Open a new instance of AMSinput
   Load the previously optimized T_{\mathrm{1}} structure via **File → Import Coordinates**
   ...
   **Properties → Excitations (UV/Vis)**
   **Type of excitations → Triplet**
   Enter the **Number of excitations** "5"
   Calculate the excited state gradient for Excitation number "1"
   Save the input **File → Save as...**

See also:

  * "4CzIPN_optFreqDFTB_S1.run" DFTB optimization and normal modes on
    S_{\mathrm{1}}

  * "4CzIPN_optFreqDFTB_T1.run" DFTB optimization and normal modes on
    T_{\mathrm{1}}

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:

   $AMSBIN/fcf << eor
   STATE1 T1_optFreqDFTB.rkf
   STATE2 S1_optFreqDFTB.rkf
   NumericalQuality 0-0
   TRANSLATE
   ROTATE
   eor

Note, that we only couple the quanta 0 of each states as higher
vibronic quanta are expected to only have insignificant contributions
and because the number of pairs of vibronic wavefunctions to couple
quickly becomes overwhelming. Running this script yields the
frequencies and Electron-Phonon couplings for both states as shown in
the following example:

   =====================
   F I R S T   S T A T E
   =====================

                Frequency         Displacement  Electron-Phonon coupling  Huang-Rhys factor
                  cm^-1            a0 amu^1/2        (dimensionless)      (dimensionless)
               ----------          ----------          ----------          ----------
                12.923361          -11.749717            2.721999            7.409281
                14.307062            0.247000            0.060207            0.003625
                19.553068           -0.055138            0.015712            0.000247
    ...
              3074.945464           -0.044556            0.159220            0.025351
              3085.682184            0.014429            0.051651            0.002668
              3086.484483           -0.000025            0.000091            0.000000

    ...

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. We then sum these
squares of all modes with \omega_{i} > 1000 cm^-1 to get the Huang-
Rhys factor S_{eff} = \sum\limits_{i}S_{i}. The frequencies
corresponding to the Huang-Rhys factor, \omega_{eff}, can be
calculated by using these values: \omega_{eff} =
\sum_{i}\frac{\omega_{i}S_{i}}{S_{eff}}. This results in S_{eff} =
0.756 and \omega_{eff} = 1500 cm^-1.

These quantities are then used in the above expression for
\rho_{\mathrm{FCWD}}^{\mathrm{CMT}} (with \lambda_{\mathrm{M}} =
\lambda = 0.074 eV) whereas the summation is stopped at n = 5. With
|V_{\mathrm{SOC}}|^{2} = 8.2\cdot10^-10 eV^2 this results in a reverse
intersystem crossing rate that is essentially identical to the result
from classical Marcus theory of k_{\mathrm{RISC}} = 5.24\cdot10^4 / 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.


Computational Description of TADF 4: Solvent Effects
====================================================

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).
