#  The role of spin polarization and dynamic correlation in singlet-triplet
gap inversion of heptazine derivatives

Daria Drwal  Institute of Physics, Lodz University of Technology, ul.
Wolczanska 219, 90-924 Lodz, Poland Mikulas Matousek  J. Heyrovský Institute
of Physical Chemistry, Academy of Sciences of the Czech Republic, v.v.i.,
Dolejškova 3, 18223 Prague 8, Czech Republic  Pavlo Golub  J. Heyrovský
Institute of Physical Chemistry, Academy of Sciences of the Czech Republic,
v.v.i., Dolejškova 3, 18223 Prague 8, Czech Republic  Aleksandra Tucholska
Institute of Physics, Lodz University of Technology, ul. Wolczanska 219,
90-924 Lodz, Poland Michał Hapka  Faculty of Chemistry, University of Warsaw,
ul. L. Pasteura 1, 02-093 Warsaw, Poland  Jiri Brabec  J. Heyrovský Institute
of Physical Chemistry, Academy of Sciences of the Czech Republic, v.v.i.,
Dolejškova 3, 18223 Prague 8, Czech Republic  Libor Veis  libor.veis@jh-
inst.cas.cz J. Heyrovský Institute of Physical Chemistry, Academy of Sciences
of the Czech Republic, v.v.i., Dolejškova 3, 18223 Prague 8, Czech Republic
Katarzyna Pernal<EMAIL_ADDRESS>Institute of Physics, Lodz University of
Technology, ul. Wolczanska 219, 90-924 Lodz, Poland

###### Abstract

The new generation of proposed light-emitting molecules for OLEDs has raised a
considerable research interest due to its exceptional feature—a negative
singlet-triplet (ST) gap violating the Hund’s multiplicity rule in the excited
S1 and T1 states. We investigate the role of spin polarization in the
mechanism of ST gap inversion. Spin polarization is associated with doubly
excited determinants of certain types, whose presence in the wavefunction
expansion favors the energy of the singlet state more than that of the
triplet. Using a perturbation theory-based model for spin polarization, we
propose a simple descriptor for prescreening of candidate molecules with
negative ST gaps and prove its usefulness for heptazine-type molecules.
Numerical results show that the quantitative effect of spin polarization is
approximately inverse-proportional to the HOMO-LUMO exchange integral.
Comparison of single- and multireference coupled-cluster predictions of ST
gaps shows that the former methods provide good accuracy by correctly
balancing the effects of doubly excited determinants and dynamic correlation.
We also show that accurate ST gaps may be obtained using a complete active
space model supplemented with dynamic correlation from multireference
adiabatic connection theory.

###### keywords:

gap inversion, Hund’s rule, adiabatic connection, spin polarization

Contributed equally. Faculty of Mathematics and Physics, Charles University,
Prague, Czech Republic Contributed equally.

##  1 Introduction

Detecting candidates for organic light-emitting diodes (OLEDs) through
screening of potential chromophores remains a challenge for quantum chemistry
methods. OLED molecules may rely on different light-emitting mechanisms,
ranging from pure fluorescence or phosphorescence (first and second generation
emitters) to more complex processes aimed at harvesting nonfluorescent triplet
excitons (third and fourth generation). 1, 2 Typically, organic molecules obey
Hund’s multiplicity rule3 and the triplet state T1 is lower than the first
excited singlet state S1. Hence, only a minor part of the excitons is
available for photon emission through radiative recombination of singlets.
This restricts internal quantum efficiency (IQE) to 25% when the distribution
of excitons between S1 and T1 states is most favorable 4. For molecules
characterized by small singlet-triplet gaps it is possible to thermally induce
reverse intersystem crossing (RISC), i.e. transition from lower-lying T1 (dark
state) to higher S1 (bright state). This process, known as thermally activated
delayed fluorescence (TADF), results in higher-rate fluorescent emission from
the singlet state. 5, 6, 7, 8, 9

Recent attention has focused on systems with inverted singlet-triplet gap
(INVEST) molecules. In INVEST emitters the S1 state lies below T1, so that the
relaxation from T1 takes place with no need for thermal induction and 100% IQE
of luminescence could, in principle, be achieved. 10, 11, 12, 13 It is worth
mentioning that aromatic chromophores with negative S${}_{1}-$T1 energy
difference (ST gap) are of interest in another area of active research, namely
in water-splitting photocatalysis. It has been shown that derivatives of
heptazine may act as efficient photocatalysts. Thanks to the negative ST
energy gap, the S1 state is exceptionally long-lived, as it does not suffer
from quenching through intersystem crossing (ISC) to T1.14

First calculations proving the possibility of ST gap inversion were carried
out independently in 2019 by de Silva for cycl[3.3.3]azine10 and Sobolewski,
Domcke and co-workers for heptazine. 15 Ref. 15 presented also the first
indirect experimental proof for the gap inversion and indicated that this
phenomenon could explain the high efficiency of OLEDs observed already in
201316 and 201417, 18 by Adachi and co-workers. A subsequent theoretical study
by Sobolewski and Domcke11 confirmed this hypothesis. Works of Sancho-Garcia
and co-workers19, 20 demonstrated ST gap inversion across N- and B-doped
triangulenes of different sizes. In 2022, Aizawa and co-workers21 proposed a
prescreening approach for INVEST candidates. Based on the elimination process,
two heptazine analogues were chosen for further, experimental evaluation,
which confirmed gap inversion.

The magnitude of the S1-T1 energy gap is directly related to the exchange
integral involving the highest occupied molecular orbital (HOMO) and lowest
unoccupied molecular orbital (LUMO). 22, 23, 24 A vanishing overlap between
frontier orbitals is a prerequisite for ST inversion, but it is not
sufficient. Several studies have demonstrated that accounting for electron
correlation via inclusion of double excitations is essential to stabilize the
singlet state with respect to triplet and obtain negative ST gaps of
chromophore molecules. 10, 24, 12 Results from single-reference response
methods suggest that double excitations contribute at a relatively low level
to the S1 state, estimated at ca. 10% by de Silva. 10 Time-dependent density
functional theory (TD-DFT) approaches are incapable of capturing double
excitations, and do not predict negative gaps, as corroborated by a range of
investigations. 10, 15, 25, 26, 27, 28

In pursuit of efficient strategies in designing novel INVEST materials,
Pollice et al.23 investigated all possible permutations of cyclazine and
heptazine derivatives with C–H replaced by a set of electron-donating and
electron-withdrawing substituents, and described them at the EOM-CCSD and TD-
DFT levels of theory. Aizawa and coworkers21 expanded on this work by
introducing 186 new substituents and finding almost 35 000 potential
candidates for INVEST molecules. Their further TD-DFT screening resulted in a
significant limitation of this set. Still, TD-DFT-based screening is
inherently limited, as the method cannot distinguish INVEST molecules from
cases in which the ST gap is small, but not negative. Having studied the
relationship between the structure and properties of selected INVEST
molecules, Olivier and co-workers13 formulated a set of design rules in which
the $C_{\rm 2v}$ point group of the triangulene core was identified as a
prerequisite for the gap inversion.

Advances in designing novel INVEST systems are hindered by two factors:
limited understanding of the mechanism behind gap reversal at the electronic
structure level and the lack of efficient descriptors that could account for
the main effects lowering S1 below T1. In one of the first theoretical works
devoted to the violation of Hund’s rule in closed-shell molecules, Kollmar and
Staemmler22 introduced a concept of dynamic spin polarization (sp), which
associates the ST gap inversion to the energetic effect exerted by a small
subset of doubly-excited configurations involving frontier orbitals. Numerical
investigations on conjugated hydrocarbons carried out by Koseki et al.29
showed that spin polarization may indeed lead to gap inversion in some
molecules. Spin polarization has been also indirectly identified as a key
effect in ST gap inversion of the heptazine molecule, without, however,
providing its quantitative measure.15

The main objective of our study is to elucidate the role of both spin
polarization and dynamic correlation energy in the mechanism of ST gap
inversion in heptazine-based molecules. Based on the spin polarization model,
we propose a simple descriptor for prescreening of molecules with a negative
ST gap. We also demonstrate the efficacy of selected multireference methods in
accurately predicting the magnitude of ST gaps in organic INVEST emitters. We
compare multireference adiabatic connection methods30, 31, 32, 33, 34, 35, 36
and the second-order $n$-electron valence state perturbation theory (NEVPT2)37
against Mukherjee's multireference coupled-cluster with non-iterative triple
excitations [Mk-MRCCSD(T)]. 38

##  2 Computational details

Our analysis is focused on six heptazine-based systems, see Fig. 1. Geometries
were optimized at the MP2 level using the Molpro39 program. All results were
obtained using the def2-TZVP40 basis set.

Complete active space self-consistent field (CASSCF) and NEVPT2 results were
obtained with the Molpro 39 program. We used the same software for restricted
Hartree-Fock (HF), CASSCF, and Kohn-Sham DFT calculations to obtain electron
integrals and one- and two-electron reduced density matrices needed for spin
polarization models and adiabatic connection, AC0 and ACn, calculations. The
latter were performed with the GammCor program. 41 Multireference coupled-
cluster Mk-MRCCSD(T) 38 energies for S1 and T1 states computed using the
CASSCF(2,2) natural orbitals were obtained with the NWChem program42 and set
as benchmark. CC2 results in def2-TZVP were taken from Ref. 24 while EOM-CCSD
ST gaps were calculated with the Orca code. 43

The ground state geometry of each molecule was used for both the S1 and T1
states. Except for system 4, calculations were performed employing the
$C_{2v}$ point group symmetry, which allowed for state-specific CASSCF
calculations for S1 and T1 states. System 4 is in a lower, $C_{s}$, point
group symmetry, so that two-state state-average CASSCF calculation is required
to access the S1 state.

Additional calculations were performed on extended test set including
heptazine-derivatives obtained by substituting C-H groups with N atoms in all
possible ways. Geometries of all molecules from the extended set were
optimized at B3LYP/cc-pVDZ level with no symmetry constraints using Orca code.
43

1

2

3

4

5

6

Figure 1: Structures of the studied molecules. The color codes are as follows:
N (blue), C (brown), and H (white).

##  3 Inverting S${}_{1}-$T1 energy gaps by including dynamic spin
polarization

Throughout the text $\Delta E_{\rm ST}$ will denote the singlet-triplet energy
gap defined as a difference of T1 energy subtracted from the energy of S1
state, namely

$\Delta E_{\rm ST}=E({\rm S_{1}})-E({\rm T_{1}})\ \ \ .$ (1)

Gap inversion would be equivalent to obtaining a negatively-valued $\Delta
E_{\rm ST}$ energy difference.

In the first approximation, the S1 and T1 states of the considered systems can
be described by singly-excited determinants, where one electron from the HOMO
($H$) orbital is excited to the orbital LUMO ($L$)

$\displaystyle\Psi_{S}^{0}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\left|\bar{H}L\right|-\left|H\bar{L}\right|\right)\
\ \ ,$ (2) $\displaystyle\Psi_{T}^{0}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\left|\bar{H}L\right|+\left|H\bar{L}\right|\right)\
\ \ .$ (3)

A notation $\left|\bar{H}L\right|$ stands for an $N$-electron Slater
determinant comprising $N/2-1$ doubly occupied orbitals and $H,L$ singly
occupied orbitals. $\bar{H}$, $\bar{L}$ and $H$, $L$ indicate orbitals with
$\alpha$ and $\beta$ spin components, respectively. Additionally, we use
$\\{pq\\}$ and $\\{pqrs\\}$ symbols to indicate which orbitals in open-shell
determinants are singly-occupied without explicitly specifying their spin
components. For example, both determinants $\left|\bar{H}L\right|$ and
$\left|H\bar{L}\right|$ would be denoted as $\\{HL\\}$.

As it is well known, the ST energy gap corresponding to a wavefunction
including only $\\{HL\\}$ singly excited determinants, Eqs. (2) and (3),

$\Delta E_{\rm
ST}^{0}=\left\langle\Psi_{S}^{0}|\hat{H}|\Psi_{S}^{0}\right\rangle-\left\langle\Psi_{T}^{0}|\hat{H}|\Psi_{T}^{0}\right\rangle=2\left\langle
HL|LH\right\rangle\geq 0\ \ \ $ (4)

is determined by the magnitude of the HL exchange integral (throughout the
text, two-electron integrals are written in the $\braket{12}{12}$ convention)
and is non-negative. Eq. (4) is a manifestation of the Hund’s rule satisfied
by wavefunctions $\Psi_{S}^{0}$ and $\Psi_{T}^{0}$.

In the Hund’s picture the ST gap is determined by the coupling of $H$ and $L$
electrons in the target spin state of the total wavefunction. This approach
ignores the interaction of the unpaired $\alpha$ and $\beta$ electrons with
the core electrons. However, Kollmar and Staemmler22 showed that including
these interactions is not merely a quantitative refinement, but may explain
the ST gap inversion in some systems. Their approach is based on extending the
first-order wavefunction with double excited determinants $\\{iaHL\\}$, where
core orbitals $i$ relax by excitations to an arbitrary virtual orbital $a$.
Including such determinants is referred to as ‘dynamic’ spin polarization to
distinguish if from the conventional ‘static’ spin polarization. While static
spin polarization leads to unequal spin components of electron density, the
dynamic spin polarization does not lead to nonzero spin density. Our goal is
to quantify the effect of dynamic spin polarization for heptazine-type
molecules and investigate if this mechanism is responsible for changing the
sign of $\Delta E_{\rm ST}^{0}$ gaps. We also investigate if sp-based gap
inversion model of Kollmar and Staemmler can be employed for prescreenig
molecules likely to act as INVEST systems.

We begin by following Ref. 22 and construct singlet functions by applying
$i\rightarrow a$ excitations to determinants in the $\Psi_{S}^{0}$ function,
which leads to two functions

$\displaystyle\Psi_{S}^{1}$
$\displaystyle=\frac{1}{\sqrt{12}}\left(2|\bar{i}\bar{a}HL|+2|ia\bar{H}\bar{L}|-|\bar{i}a\bar{H}L|-|\bar{i}aH\bar{L}|-|i\bar{a}\bar{H}L|-|i\bar{a}H\bar{L}|\right)\
\ \ ,$ (5) $\displaystyle\Psi_{S}^{1^{\prime}}$
$\displaystyle=\frac{1}{2}(|i\bar{a}H\bar{L}|-|i\bar{a}\bar{H}L|-|\bar{i}aH\bar{L}|+|\bar{i}a\bar{H}L|)\
\ \ .$ (6)

Analogously, three triplet $M_{S}=0$ functions can be generated from
$\Psi_{T}^{0}$ by considering $i\rightarrow a$ excitations and they read

$\displaystyle\Psi_{T}^{1}$
$\displaystyle=-\frac{1}{2}(|\bar{i}a\bar{H}L|-|\bar{i}aH\bar{L}|+|i\bar{a}\bar{H}L|-|i\bar{a}H\bar{L}|)\
\ \ ,$ (7) $\displaystyle\Psi_{T}^{2}$
$\displaystyle=-\frac{1}{\sqrt{2}}(|\bar{i}\bar{a}HL|-|ia\bar{H}\bar{L}|)\ \ \
,$ (8) $\displaystyle\Psi_{T}^{1^{\prime}}$
$\displaystyle=-\frac{1}{2}(|\bar{i}a\bar{H}L|+|\bar{i}aH\bar{L}|-|i\bar{a}\bar{H}L|-|i\bar{a}H\bar{L}|)\
\ \ .$ (9)

Recall that from the Epstein-Nesbet perturbation theory (PT), the drop in the
unperturbed energy, $E^{0}=\langle\Psi^{0}|\hat{H}|\Psi^{0}\rangle$, that
results from including in the wavefunction a configuration $\Psi^{\prime}$ and
corresponds to the second-order correction, is given by
$|\langle\Psi^{0}|\hat{H}|\Psi^{\prime}\rangle|^{2}/(E^{0}-E^{\prime})$, where
$E^{\prime}=\langle\Psi^{\prime}|\hat{H}|\Psi^{\prime}\rangle$. After
evaluating PT terms for functions given in Eqs. (5)–(9), see Supporting
Information for explicit expressions of Hamiltonian elements and energy
differences in terms of two-electron integrals, and subtracting triplet
contributions from the singlet ones, the following spin-free expression is
obtained

$\Delta E_{i,a}^{sp}=\frac{3}{2}\frac{\left(\left\langle
iH|Ha\right\rangle-\left\langle
iL|La\right\rangle\right)^{2}}{E_{S}^{0}-E_{S}^{1}}-\frac{1}{2}\frac{\left(\left\langle
iH|Ha\right\rangle-\left\langle
iL|La\right\rangle\right)^{2}}{E_{T}^{0}-E_{T}^{1}}-\frac{\left(\left\langle
iH|Ha\right\rangle+\left\langle
iL|La\right\rangle\right)^{2}}{E_{T}^{0}-E_{T}^{2}}\ \ \ $ (10)

in agreement with Ref. 22 (notice that contributions to $\Delta E_{i,a}^{sp}$
from states $\Psi_{S}^{1^{\prime}}$ and $\Psi_{T}^{1^{\prime}}$, not
considered in Ref.22, cancel each other). $\Delta E_{i,a}^{sp}$ represents a
contribution to the ST gap from spin polarization of two electrons occupying
the orbital $i$ by allowing their excitation to the orbital $a$. To ease its
interpretation, expression in Eq. (10) can be approximated by assuming a
common denominator given as Hartree-Fock orbital energy difference
$\varepsilon_{i}-\varepsilon_{a}$ for all three terms, leading to

$\Delta E_{i,a}^{sp}\approx 4\frac{\left\langle iH|Ha\right\rangle\left\langle
iL|La\right\rangle}{\varepsilon_{a}-\varepsilon_{i}}\ \ \ .$ (11)

It is now clear that spin polarization energy, $\Delta E_{i,a}^{sp}$,
pertaining to a pair of orbitals $(i,a)$, where $i<H$ and $a>L$, stabilizes
the singlet state with respect to triplet if orbitals $i$ and $a$ overlap
significantly with both $H$ and $L$ orbitals and, taking into account that
$\epsilon_{a}-\epsilon_{i}>0$, exchange integrals $\braket{iH}{Ha}$ and
$\braket{iL}{La}$ are of the opposite sign. The magnitude of the $\Delta
E_{i,a}^{sp}$ energy is dependent on the closeness of the $i$ and $a$ orbital
energy levels.

For the studied systems the $\left(H-1,H-2\right)$ and $\left(L+1,L+2\right)$
pairs of orbitals are degenerate, so the simplest model for the ST gap
accounting for major spin polarization would include two pairs of orbitals,
and the ST energy gap expression of such a “12” model would read

$\Delta E_{\rm ST}^{sp_{12}}=\Delta E_{\rm ST}^{0}+\Delta
E_{H-1,L+1}^{sp}+\Delta E_{H-2,L+2}^{sp}\ \ \ \ .$ (12)

For system 4, the model included also the contributions from the other two
combinations of orbitals $\Delta E_{H-1,L+2}^{sp}+\Delta E_{H-2,L+1}^{sp}$,
which are negligible in the other systems due to point group symmetry. To
fully account for spin polarization, considered pairs of orbitals $(i<H,a>L)$
must fulfill conditions leading to ST gap reversal, formulated below Eq. (11):
$i$) orbitals within each pair should strongly overlap with $H$ and $L$
orbitals, $ii$) the exchange integrals $\braket{iH}{Ha}$ and $\braket{iL}{La}$
are of the opposite signs. For heptazine derivatives the $\pi$ orbitals
satisfy these requirements and the all-$\pi$ sp-inclusive ST energy gap reads

$\Delta E_{\rm ST}^{sp_{\pi}}=\Delta E_{\rm ST}^{0}+\sum_{i<H}\sum_{a>L}\Delta
E_{i,a}^{sp}\ \ \ ,$ (13)

(orbitals $i$ and $a$ are of the $\pi$-type).

We consider two sets of orbitals in further analysis. The first set is given
by canonical HF orbitals. The second one corresponds to a ground state
CASSCF(14,14) wavefunction with the active space containing all the $\pi$
electrons and orbitals (2p${}_{\text{z}}$ orbitals on carbon and nitrogen
atoms plus one virtual orbital of the “double-shell” 3p${}_{\text{z}}$
character24). The CASSCF(14,14) natural occupation numbers for systems 2 and
4, used as illustrative cases, are given in Table 1. Notice that orbitals $H$
and $L$ pertain to $(N/2)$th and $(N/2+1)$th orbitals, respectively, of the
occupation numbers close to $0.5$. These orbitals are depicted in Figure 2
(see also Figures 1-12 in Supporting Information) together with four other
most strongly correlated, i.e. of the occupancies deviating most from $0$ and
$1$, orbitals.

Table 1: Natural occupation numbers for systems 2 and 4 from state-averaged CASSCF(14,14) calculations, for system 2 corresponding to the natural orbitals presented in Figure 2.  | Sys 2 | Sys 4
---|---|---
Orbital | S1 | T1 | S1 | T1
$H$ \- 6 | 0.995 | 0.993 | 0.990 | 0.989
$H$ \- 5 | 0.969 | 0.970 | 0.971 | 0.971
$H$ \- 4 | 0.967 | 0.968 | 0.970 | 0.971
$H$ \- 3 | 0.967 | 0.968 | 0.969 | 0.971
$H$ \- 2 | 0.928 | 0.933 | 0.936 | 0.939
$H$ \- 1 | 0.929 | 0.933 | 0.933 | 0.936
$H$ | 0.513 | 0.504 | 0.639 | 0.509
$L$ | 0.485 | 0.492 | 0.359 | 0.486
$L$ \+ 1 | 0.074 | 0.073 | 0.070 | 0.071
$L$ \+ 2 | 0.074 | 0.073 | 0.069 | 0.070
$L$ \+ 3 | 0.036 | 0.035 | 0.033 | 0.031
$L$ \+ 4 | 0.036 | 0.035 | 0.032 | 0.031
$L$ \+ 5 | 0.026 | 0.023 | 0.024 | 0.022
$L$ \+ 6 | 0.001 | 0.001 | 0.004 | 0.004
CI coefficients
$\\{HL\\}$ | 68.6% | 69.3% | 67.3% | 71.3%
$\\{(H-p)\ (L+p)\ H\ L\\}_{p=1,2}$ | 8.2% | 7.8% | 7.2% | 0.3%
$\\{i\ a\ H\ L\\}_{i,a\in\pi}$ | 15.0% | 14.6% | 13.4% | 5.6%

HOMO–2

HOMO–1

HOMO

LUMO

LUMO+1

LUMO+2

Figure 2: Natural orbitals of system 2 corresponding to state-averaged
CASSCF(14,14) calculations with one singlet and one triplet states.

The results for ST gaps obtained without spin polarization, $\Delta E^{0}_{\rm
ST}$, with partial and full spin polarization, $\Delta E_{\rm ST}^{sp_{12}}$
and $\Delta E_{\rm ST}^{sp_{\pi}}$ models, respectively, are presented in
Table 2. First, it can be seen that $\Delta E^{0}_{\rm ST}$ gaps are positive
and small, which indicates that for all considered systems the $H$ and $L$
orbitals overlap marginally both at the HF and CASSCF(14,14) levels of theory.
Accounting for dynamic spin polarization from only two pairs of orbitals,
$(H-1,L+1)$ and $(H-2,L+2)$, see Eq. (12), significantly reduces ST gaps.
Using canonical HF orbitals leads to negative ST gaps for systems 1 and 2,
while with CASSCF(14,14) natural orbitals all ST gaps turn negative.
Considering all $(i<H,a>L)$ pairs of $\pi$ orbitals, Eq. (13), removes another
$0.2-0.3$ eV from the ST gap, regardless of the employed orbitals. We conclude
that the PT-based model for dynamic spin polarization substantially lowers the
S1 state energy with respect to T1, possibly leading to sign reversal of the
ST gap. The combined sp effect of two pairs of orbitals, $(H-1,L+1)$ and
$(H-2,L+2)$, accounts for more than half of the ST gap reduction.

Table 2: S${}_{1}-$T1 energy gaps obtained without spin polarization, $\Delta E_{\rm ST}^{0}$, with spin polarization computed from two pairs of $\pi$ orbitals, $\Delta E_{\rm ST}^{{\rm sp}_{12}}$, from all occupied-virtual pairs of $\pi$ orbitals, $\Delta E_{\rm ST}^{{\rm sp}_{\pi}}$, and from CI with $\\{H\ L\\}$, $\\{H-1\ L+1\ H\ L\\}$, and $\\{H-2\ L+2\ H\ L\\}$ determinants, $\Delta E_{\rm ST}^{{\rm CI}_{12}}$. HF canonical and CASSCF(14,14) orbitals were used. The latter followed from state-averaged CASSCF calculations with one singlet and one triplet states. All values in eV. Orbitals | ST gap | $1$ | $2$ | $3$ | $4$ | $5$ | $6$
---|---|---|---|---|---|---|---
HF | $\Delta E_{\rm ST}^{0}$ | $0.25$ | $0.25$ | $0.49$ | $0.67$ | $0.68$ | $0.61$
$\Delta E_{\rm ST}^{sp_{12}}$ | $-0.16$ | $-0.18$ | $0.09$ | $0.29$ | $0.33$ | $0.22$
$\Delta E_{\rm ST}^{sp_{\pi}}$ | $-0.37$ | $-0.45$ | $-0.11$ | $0.09$ | $0.07$ | $0.00$
$\Delta E_{\rm ST}^{\rm CI_{12}}$ | $-0.21$ | $-0.25$ | $0.10$ | $0.27$ | $0.10$ | $0.08$
CASSCF | $\Delta E_{\rm ST}^{0}$ | $0.19$ | $0.19$ | $0.25$ | $0.32$ | $0.25$ | $0.28$
$\Delta E_{\rm ST}^{sp_{12}}$ | $-0.33$ | $-0.40$ | $-0.28$ | $-0.19$ | $-0.31$ | $-0.27$
$\Delta E_{\rm ST}^{sp_{\pi}}$ | $-0.60$ | $-0.73$ | $-0.54$ | $-0.46$ | $-0.61$ | $-0.55$
$\Delta E_{\rm ST}^{\rm CI_{12}}$ | $-0.29$ | $-0.36$ | $-0.25$ | $-0.14$ | $-0.28$ | $-0.24$

Spin polarization in the wavefunction picture is manifested by the presence of
doubly-excited determinants. To check if the simple PT-based model for sp
yields qualitatively correct predictions, we have performed CI calculations
constructing the S1 and T1 wavefunctions from $\\{H\ L\\}$, $\\{H-1\ L+1\ H\
L\\}$, and $\\{H-2\ L+2\ H\ L\\}$ singly- and doubly-excited determinants,
employing either HF or CASSCF(14,14) orbitals. The CI results, denoted as
$\Delta E_{\rm ST}^{\rm CI_{12}}$ in Table 2, are in close agreement with the
$\Delta E_{\rm ST}^{sp_{12}}$ model. This suggests that PT-based approach can
be used as a cost-saving alternative to CI calculations in prescreening for
molecules with inverted gaps.

ST gaps obtained from $\Delta E_{\rm ST}^{sp_{12}}$ and $\Delta E_{\rm
ST}^{sp_{\pi}}$ should be considered as basic approximations due to the fact
that there are only $\\{HL\\}$ and $\\{iaHL\\}$ determinants in the
wavefunctions. Thus, the model gaps cannot be quantitatively correct. A
comparison with Mk-MRCCSD(T) results, shown in the last column of Table 3,
reveals that ST gaps predicted by sp12 and spπ models employed with HF
orbitals, cf. Table 2, deviate on average by 0.26 and 0.15 eV from the Mk-
MRCCSD(T) benchmarks, respectively. With CASSCF(14,14) orbitals the respective
mean absolute errors amount to 0.14 and 0.43 eV. Such inaccuracies result from
the lack of dynamic electron correlation in the considered models. This
problem is addressed in the next section.

To further explore if the sp effect is decisive in inverting the ST gap and if
the model presented in Eq. (13) would be useful in screening for INVEST
candidates, we have applied it to an extended test set of molecules including
heptazine-derived systems obtained by substitution of C-H with N in all
possible ways, similar as in Ref. 23. The ST gaps have been computed using
Eqs. (10) and (13). Initial tests on systems 1-6 have shown, see Table 2, that
CASSCF(14,14) orbitals satisfy conditions required to invert the ST gaps via
the sp mechanism to a greater extent than HF orbitals and, unlike the latter,
they have led to obtaining negative gaps for all molecules. This suggests that
correlated orbitals are a better choice for INVEST screening. Guided by this
finding, we have used KS-DFT orbitals from ground-state BLYP44, 45, 46
calculations for tests on our extended set of molecules. In Figure 3, we
present ST energy gaps obtained without spin polarization, see Eq. (4), and
with spin polarization fully accounted for, cf. Eq. (13), compared with the
EOM-CCSD values. Evidently, employing only the HOMO-LUMO exchange integral for
INVEST prescreening does a poor job, as hardly any correlation between $\Delta
E_{\rm ST}^{0}$ and EOM-CCSD ST gap values can be seen. Shifting the $\Delta
E_{\rm ST}^{0}$ values by $-0.30$ eV and keeping only molecules with resulting
negative gaps would leave us with too large set: not only molecules with EOM-
CCSD-predicted negative gaps would be included, but also most of those with
positive gaps, see left panel in Figure 3. Thus, a criterion based on $\Delta
E_{\rm ST}^{0}$ values is not sufficiently selective in screening for INVEST
systems. A satisfactory correlation is obtained by employing the full spπ
model with KS-DFT orbitals and shifting the $\Delta E_{ST}^{sp_{\pi}}$
energies by $+0.2$ eV. This value is recommended in future INVEST prescreening
calculations, as it results in a good overlap between the subsets of molecules
with positive and negative gaps obtained from the sp model and EOM-CCSD, as
shown in Figure 3, right panel. Notice that the correlation is worse with HF
orbitals, see the right panel of Figure 13 in Supporting Information.

Figure 3: Left panel: ST energy gaps without spin polarization, Eq. (4). Right
panel: ST energy gaps with spin polarization from all occupied-virtual $\pi$
orbital pairs, Eqs. 13 and (10). Results obtained with BLYP ground state
orbitals vs. EOM-CCSD values for extended test set of molecules.

Another observation following from calculations on the extended set of
molecules is that the sp effect is approximately inversely proportional to the
HL exchange interaction, see Figure 4. Therefore, larger HL integrals
correspond to lower spin polarization, both effects being detrimental to
obtaining negative ST gaps. Apparently, sp for the considered isoelectronic
molecules is as sensitive as HL exchange integrals to the composition and
distribution of carbon and nitrogen atoms in heptazine derivatives.

In Ref.13 low HL exchange interaction has been associated with high symmetry
point groups $D/C_{3h}$, $C_{3v}$, $C_{2v}$. In particular, the existence of
the $\sigma_{v}$ plane was identified as critical in minimizing the HL
overlap. Since, as we have shown, minimization of the latter is accompanied by
maximization of the sp effect, in general it seems to be a good strategy to
account for symmetry while designing INVEST molecules. However, symmetry
cannot be used as the only criterion in high-throughput screening. Recall,
that system 4 does not possess $\sigma_{v}$ plane as symmetry element (Figure
1), but the corresponding ST gap is negative. In contrast, the molecule with
the largest HL exchange integral and most positive $\Delta E_{ST}^{sp_{\pi}}$
gap, cf. Figure 4, belongs to the $C_{2v}$ point group.

Results obtained on the extended set show that a simple spin polarization
model considered in this section used with KS-DFT orbitals can serve as a
computationally inexpensive predictor of INVEST molecules. Notice that for
systems 1-6 a common-denominator approximation adopted in Eq. (11) yields
energy gaps in a good agreement with those obtained if sp from all $\pi$
orbital-pairs is included as in Eq. (10), see Table 1 in Supporting
Information. Thus, further simplification of the model is achievable if Eq.
(11) is used instead of Eq. (10). Compared to predictors used in Ref. 23 based
on a double-hybrid DFT functional the proposed model is computationally more
efficient (it requires performing a ground state calculation with a semi-local
functional) and above all physically meaningful, as it explicitly includes the
sp effect responsible for the gap inversion.

Figure 4: Left panel: Contributions from spin polarization to ST energy gaps,
second term of the right-hand-side of Eq. (13), vs. ST energy gaps without
spin polarization, Eq. (4). Right panel: Molecules corresponding to the lowest
and highest H-L exchange interaction. Results obtained with BLYP ground state
orbitals.

##  4 The effect of dynamic correlation energy

CAS wavefunctions should effectively capture the sp effect. One expects that
expanding the active space from (2,2), two electrons on $H$ and $L$ orbitals,
via (6,6), six electrons on $H$, $L$, and two pairs of degenerate orbitals
$H-1$, $H-2$ and $L+1$, $L+2$, up to (14,14), all $\pi$ electrons on all $\pi$
orbitals, would lead to ST energy gaps that reflect the $\Delta
E^{0}_{ST}>\Delta E^{sp_{12}}_{\rm ST}>\Delta E^{sp_{\pi}}_{\rm ST}$ relation
observed for PT-based models, namely:

$\Delta E_{\rm ST}^{\rm CASSCF(2,2)}>\Delta E_{\rm ST}^{\rm
CASSCF(6,6)}>\Delta E_{\rm ST}^{\rm CASSCF(14,14)}\ \ \ .$ (14)

Indeed, CASSCF energy gaps presented in Table 3 agree with the predicted
inequality relations, Eq. (14), for each system (see the column denoted as
CAS). Going from CAS(2,2) to CAS(6,6) lowers the gap by as much as $0.6$-$1.9$
eV. Including all $\pi$ orbitals as active, in the CASSCF(14,14) model,
reduces the gaps by another $0.3$-$0.5$ eV. The $\Delta E_{\rm ST}^{\rm
CASSCF(14,14)}$ values are all negative and correspond well with the $\Delta
E_{\rm ST}^{sp_{\pi}}$ results obtained using CASSCF(14,14) orbitals, cf.
Table 2. We conclude that doubly excited determinants $\\{iaHL\\}$, where
$i<H$ and $a>L$ are $\pi$ orbitals, are the main contributors to the
wavefunction, after the leading singly excited $\\{HL\\}$ determinant. In
Table 1 we report sums of squares of pertinent CI coefficients for systems 2
and 4. The contribution of aforementioned double excitations is up to 15%.
This is far less than a contribution from $\\{HL\\}$ determinants amounting to
ca. 70%, but it determines ST gap inversion. It also stresses the importance
of having a small value of the HOMO-LUMO exchange integral to achieve gap
inversion.

Table 3: S${}_{1}-$T1 energy gaps in eV obtained from CASSCF, AC0, ACn, and
NEVPT2 methods for small, intermediate, and all-$\pi$ models for active
spaces, compared against CC2 and Mk-MRCCSD(T) values.

System | Active space | CAS | AC0 | ACn | NEVPT2 | CC224 | Mk-MRCCSD(T)
---|---|---|---|---|---|---|---
1 | (2,2) | $1.76$ | $0.70$ | $0.00$ | $-0.37$ | -0.13 | -0.18
(6,6) | $-0.18$ | $0.16$ | $0.01$ | $-0.08$
(14,14) | $-0.47$ | $-0.10$ | $-0.21$ | $0.07$
2 | (2,2) | $0.31$ | $-0.78$ | $-0.28$ | $-1.04$ | -0.24 | -0.28
(6,6) | $-0.31$ | $-0.02$ | $-0.14$ | $1.01$
(14,14) | $-0.62$ | $-0.21$ | $-0.34$ | $-0.06$
3 | (2,2) | $0.44$ | $-0.28$ | $0.02$ | $-0.59$ | -0.11 | -0.15
(6,6) | $-0.17$ | $0.08$ | $-0.01$ | $-0.04$
(14,14) | $-0.45$ | $-0.05$ | $-0.19$ | $0.09$
4 | (2,2) | $0.79$ | $-0.45$ | $-0.22$ | $-0.92$ | -0.08 | -0.04
(6,6) | $0.18$ | $-0.32$ | $-0.01$ | $-0.59$
(14,14) | $-0.31$ | $-0.04$ | $-0.12$ |  $-0.07$a
5 | (2,2) | $0.70$ | $-0.64$ | $-0.06$ | $-0.94$ | -0.14 | -0.15
(6,6) | $-0.12$ | $0.04$ | $-0.33$ | $-0.25$
(14,14) | $-0.50$ | $-0.06$ | $-0.20$ | $0.09$
6 | (2,2) | $0.72$ | $-0.67$ | $-0.06$ | $-0.98$ | -0.12 | -0.14
(6,6) | $-0.22$ | $0.08$ | $-0.08$ | $0.05$
(14,14) | $-0.50$ | $-0.01$ | $-0.15$ | $0.11$

a Due to problems with convergence of NEVPT2, the active space has been
reduced to (14,13).

A comparison between CASSCF(14,14) and Mk-MRCCSD(T) reference values, cf.
Table 3, reveals that CASSCF gaps are too negative with the mean unsigned
deviation exceeding 0.3 eV, see also Figure 5. To correct CASSCF for the
missing dynamic correlation energy, we use the multireference adiabatic
connection (AC) approach34, 35, 47, 48, 49, 50, 51, which recently has been
successfully applied to predicting singlet-triplet gaps of biradicals. 36 Two
AC variants are employed in this work: AC0 and ACn. The first one is based on
linearizing the adiabatic connection integrand.34, 51 ACn is free of such an
approximation and is expected to yield more accurate predictions.36 It is
worth mentioning that the computational cost of both AC0 and ACn scales with
the $5$th power of the system size, AC0 being more efficient than ACn due to a
smaller prefactor. AC methods rely only on 1- and 2-electron reduced density
matrices (1- and 2-RDMs, respectively), which makes them computationally more
efficient in treating large active spaces compared to canonical multireference
perturbation theory methods.50

AC0 (ACn) ST gaps are computed by adding dynamic correlation correction to the
CASSCF gap according to

$\Delta E_{\rm ST}^{\rm{AC}}=\Delta E_{\rm ST}^{\rm CASSCF}+E^{\rm AC}({\rm
S}_{1})-E^{\rm AC}({\rm T}_{1})\ \ \ ,$ (15)

where $E^{\rm AC}({\rm S}_{1})$ and $E^{\rm AC}({\rm T}_{1})$ are AC0 (ACn)
adiabatic connection correlation energies computed from CASSCF 1,2-RDM’s for
the S1 and T1 states, respectively. Inspection of the AC0 and ACn energy gaps
in Table 3 shows that accounting for dynamic correlation energy counteracts
the spin polarization effect. While spin polarization lowers the singlet state
energy more than the triplet, leading to a negatively-valued ST gap, dynamic
correlation shifts the S1 energy upward relative to T1. In the case of AC0
this effect is overestimated leading to too small gaps. AC0 gaps computed for
CASSCF(6,6) change the sign back to positive for systems 1, 3, 5, and 6.
Combining AC0 with CASSCF(14,14) retains the negative sign of the gaps, but
their magnitudes are underestimated compared with Mk-MRCCSD(T) values, see
also Figure 5. The best accuracy is obtained if ACn is applied together with a
CASSCF(14,14) model—the resulting gaps values agree up to 0.04 eV with the Mk-
MRCCSD(T) reference. AC correlation energy computed for CASSCF(2,2) models
yields poor results. This is expected, since the CASSCF(2,2) does not contain
doubly excited determinants, which AC approximations cannot amend. Ideally,
the CASSCF model should include all of the spin polarisation, leaving only
dynamic correlation energy for adiabatic connection. ACn combined with
CASSCF(14,14) proves that this strategy leads to accurate predictions.

In addition to AC, in Table 3 we present results of NEVPT2 calculations. 52
Although our previous studies show that NEVPT2 is typically as accurate as AC
methods, in particular AC034, 47, for heptazine-based systems the performance
of NEVPT2 is inferior. The NEVPT2 energy gap values do not show a systematic
improvement when enlarging the CAS(n,n) model. Quite contrary, the behavior is
erratic and the gaps are of wrong sign, except for systems 2 and 4, even if
CAS(14,14) model is employed, cf. Table 3. NEVPT2 gaps deviate on average by
0.2 eV from the reference values, see Figure 5.

Figure 5:  S${}_{1}-$T1 energy gaps for systems 1-6, panel (a). Mean unsigned
deviations of the gaps with respect to Mk-MRCCSD(T) values, panel (b). CASSCF,
AC0, ACn, and NEVPT2 results correspond to (14,14) active space.

The presented results show that two factors are equally important for accurate
prediction of ST gap inversion: accounting for a modest contribution of double
excitations and proper treatment of dynamic correlation. It is therefore not
surprising that single reference methods like CC2 or ADC(2)10, 24 are capable
of predicting ST gaps as accurately as multireference approaches. The CC2
energy gaps taken from Ref. 24 (see Table 3 and Figure 5) stay in a good
agreement with the multireference CCSD(T) results deviating from the latter by
only $0.03$ eV. Finally, we notice that the energy gaps predicted by both
single and multireference methods are quite insensitive to the basis set used,
compare the results from Table 3 and Table 2 in Supporting Information.

##  5 Summary and Conclusions

We have investigated sources of S${}_{1}-$T1 energy gap inversion in
heptazine-based molecules, focusing on the effect of dynamic spin
polarization. 22 We have found that spin polarization, which is equivalent to
including doubly excited determinants involving HOMO, LUMO and any two of the
other $\pi$ orbitals in the wavefunction expansion, drives the gap inversion.
Our findings are summarized in Figure 6. While ignoring spin polarization
leads to positive ST gaps (green bars in Figure 6), accounting for this effect
alone results in gaps which are too negative (red and blue bars in Figure 6).
For a quantitative ST gap prediction, it is essential to consider both spin
polarization and accurately treat dynamic correlation energy. Dynamic
correlation closes the gap, exerting an opposite effect to sp (see pink bars
in Figure 6). Comparing several multiconfigurational wavefunction approaches
that account for dynamic correlation, we recommend AC-based methods as having
the best accuracy/cost ratio. Among them, ACn method combined with
CASSCF(14,14) was able to reproduce accurate ST gaps for all studied heptazine
derivatives.

Figure 6:  S${}_{1}-$T1 energy gaps for systems 1-6. $0$: no spin
polarization, Eq.(4); spπ: PT-based model accounting for spin polarization,
Eq.(13); ${\rm CASSCF}$: spin polarization from CASSCF(14,14); ACn: spin
polarization at CAS(14,14) level and dynamic correlation via adiabatic
connection, Eq.(15).

In the investigated systems, the contribution of double excitations in
wavefunctions varies from 6% to 15%. This relatively weak doubly-excited
character of S1 and T1 states can be captured with single reference coupled
cluster methods. Indeed, we reported a good agreement between CC2 and
multireference approaches [ACn and Mk-MRCCSD(T)].

We have proposed a simple model for selecting INVEST systems which accounts
for two critical factors that favor negative gaps: vanishing HL exchange
integral and spin polarization. In this approach, ST gap is approximated via a
PT-based expression involving orbital energies and exchange integrals. When
KS-DFT orbitals are employed, the model is capable of efficient and accurate
prescreening for INVEST candidates, as we have demonstrated on a test set of
$\sim 100$ heptazine derivatives. We have found a previously unknown relation
between the sp contribution to the ST gap and the magnitude of the HL exchange
integral, which are inversely-proportional to each other. It implies that the
quantitative contribution to the gap from the effect of spin polarization is
negligible if the H-L exchange interaction is relatively large. Having the
advantage of being low-cost to compute, the proposed model for the
S${}_{1}-$T1 energy gap can be used in high-throughput calculations aimed at
searching for INVEST systems. It can be also used to build models for the
prediction of ST gap inversion by machine learning algorithms.

##  6 Supporting Information

Expressions: Hamiltonian matrix elements of doubly excited states, energy
differences of doubly excited states. Tables: ST energy gaps from sp models,
ST energy gaps in cc-pVDZ and cc-pVQZ basis sets, singlet- and triplet-state
energies. Figures: CASSCF(14,14) orbitals for S1 and T1 states, ST energy gaps
from sp models with HF orbitals vs. EOM-CCSD for extended set of molecules.

##  7 Acknowledgement

This work was supported by the Czech Science Foundation (Grant No. 23-04302L);
the National Science Center of Poland (grant no. 2019/35/B/ST4/01310); Lodz
University of Technology (Internal Grant FU2N - Fundusz Udoskonalania
Umiejętności Młodych Naukowców from Excellence Initiative – Research
University); the Charles University Grant Agency (Grant No. 218222); the
Center for Scalable and Predictive methods for Excitation and Correlated
phenomena (SPEC), which is funded by the U.S. Department of Energy (DOE),
Office of Science, Office of Basic Energy Sciences, the Division of Chemical
Sciences, Geosciences, and Biosciences. This work was also supported by the
Ministry of Education, Youth and Sports of the Czech Republic through the
e-INFRA CZ (ID:90254).

## References

  * Tyan 2011 Tyan, Y.-S. Organic light-emitting-diode lighting overview. _J. Photonics Energy_ 2011, _1_ , 011009–011009
  * Hong et al. 2021 Hong, G.; Gan, X.; Leonhardt, C.; Zhang, Z.; Seibert, J.; Busch, J. M.; Bräse, S. A Brief History of OLEDs—Emitter Development and Industry Milestones. _Adv. Mater._ 2021, _33_ , 2005630
  * Hund 1925 Hund, F. Zur Deutung verwickelter Spektren, insbesondere der Elemente Scandium bis Nickel. _Z. Phys._ 1925, _33_ , 345–371
  * Li et al. 2006 Li, Z.; Li, Z. R.; Meng, H. _Organic light-emitting materials and devices_ ; CRC press, 2006
  * Endo et al. 2011 Endo, A.; Sato, K.; Yoshimura, K.; Kai, T.; Kawada, A.; Miyazaki, H.; Adachi, C. Efficient up-conversion of triplet excitons into a singlet state and its application for organic light emitting diodes. _Appl. Phys. Lett._ 2011, _98_ , 42
  * Uoyama et al. 2012 Uoyama, H.; Goushi, K.; Shizu, K.; Nomura, H.; Adachi, C. Highly efficient organic light-emitting diodes from delayed fluorescence. _Nature_ 2012, _492_ , 234–238
  * Goushi et al. 2012 Goushi, K.; Yoshida, K.; Sato, K.; Adachi, C. Organic light-emitting diodes employing efficient reverse intersystem crossing for triplet-to-singlet state conversion. _Nat. Photonics_ 2012, _6_ , 253–258
  * Nakanotani et al. 2014 Nakanotani, H.; Higuchi, T.; Furukawa, T.; Masui, K.; Morimoto, K.; Numata, M.; Tanaka, H.; Sagara, Y.; Yasuda, T.; Adachi, C. High-efficiency organic light-emitting diodes with fluorescent emitters. _Nat. Commun._ 2014, _5_ , 4016
  * Hosokai et al. 2017 Hosokai, T.; Matsuzaki, H.; Nakanotani, H.; Tokumaru, K.; Tsutsui, T.; Furube, A.; Nasu, K.; Nomura, H.; Yahiro, M.; Adachi, C. Evidence and mechanism of efficient thermally activated delayed fluorescence promoted by delocalized excited states. _Sci. Adv._ 2017, _3_ , e1603282
  * de Silva 2019 de Silva, P. Inverted singlet–triplet gaps and their relevance to thermally activated delayed fluorescence. _J. Phys. Chem. Lett._ 2019, _10_ , 5674–5679
  * Sobolewski and Domcke 2021 Sobolewski, A. L.; Domcke, W. Are heptazine-based organic light-emitting diode chromophores thermally activated delayed fluorescence or inverted singlet–triplet systems? _J. Phys. Chem. Lett._ 2021, _12_ , 6852–6860
  * Pios et al. 2021 Pios, S.; Huang, X.; Sobolewski, A. L.; Domcke, W. Triangular boron carbon nitrides: An unexplored family of chromophores with unique properties for photocatalysis and optoelectronics. _Phys. Chem. Chem. Phys._ 2021, _23_ , 12968–12975
  * Ricci et al. 2022 Ricci, G.; Sancho-García, J.-C.; Olivier, Y. Establishing design strategies for emissive materials with an inverted singlet–triplet energy gap (INVEST): a computational perspective on how symmetry rules the interplay between triplet harvesting and light emission. _J. Mater. Chem. C_ 2022, _10_ , 12680–12698
  * Ehrmaier et al. 2020 Ehrmaier, J.; Huang, X.; Rabe, E. J.; Corp, K. L.; Schlenker, C. W.; Sobolewski, A. L.; Domcke, W. Molecular design of heptazine-based photocatalysts: effect of substituents on photocatalytic efficiency and photostability. _J. Phys. Chem. A_ 2020, _124_ , 3698–3710
  * Ehrmaier et al. 2019 Ehrmaier, J.; Rabe, E. J.; Pristash, S. R.; Corp, K. L.; Schlenker, C. W.; Sobolewski, A. L.; Domcke, W. Singlet–triplet inversion in heptazine and in polymeric carbon nitrides. _J. Phys. Chem. A_ 2019, _123_ , 8099–8108
  * Li et al. 2013 Li, J.; Nakagawa, T.; MacDonald, J.; Zhang, Q.; Nomura, H.; Miyazaki, H.; Adachi, C. Highly Efficient Organic Light-Emitting Diode Based on a Hidden Thermally Activated Delayed Fluorescence Channel in a Heptazine Derivative. _Adv. Mater._ 2013, _25_ , 3319–3323
  * Li et al. 2014 Li, J.; Zhang, Q.; Nomura, H.; Miyazaki, H.; Adachi, C. Thermally activated delayed fluorescence from 3n $\pi$* to 1n $\pi$* up-conversion and its application to organic light-emitting diodes. _Appl. Phys. Lett._ 2014, _105_ , 013301
  * Li et al. 2014 Li, J.; Nomura, H.; Miyazaki, H.; Adachi, C. Highly efficient exciplex organic light-emitting diodes incorporating a heptazine derivative as an electron acceptor. _Chemical Communications_ 2014, _50_ , 6174–6176
  * Ricci et al. 2021 Ricci, G.; San-Fabián, E.; Olivier, Y.; Sancho-García, J.-C. Singlet-triplet excited-state inversion in heptazine and related molecules: assessment of TD-DFT and ab initio methods. _ChemPhysChem_ 2021, _22_ , 553–560
  * Sanz-Rodrigo et al. 2021 Sanz-Rodrigo, J.; Ricci, G.; Olivier, Y.; Sancho-Garcia, J.-C. Negative singlet–triplet excitation energy gap in triangle-shaped molecular emitters for efficient triplet harvesting. _J. Phys. Chem. A_ 2021, _125_ , 513–522
  * Aizawa et al. 2022 Aizawa, N.; Pu, Y.-J.; Harabuchi, Y.; Nihonyanagi, A.; Ibuka, R.; Inuzuka, H.; Dhara, B.; Koyama, Y.; Nakayama, K.-i.; Maeda, S., et al. Delayed fluorescence from inverted singlet and triplet excited states. _Nature_ 2022, _609_ , 502–506
  * Kollmar and Staemmler 1978 Kollmar, H.; Staemmler, V. Violation of Hund’s rule by spin polarization in molecules. _Theor. Chim. Acta_ 1978, _48_ , 223–239
  * Pollice et al. 2021 Pollice, R.; Friederich, P.; Lavigne, C.; dos Passos Gomes, G.; Aspuru-Guzik, A. Organic molecules with inverted gaps between first excited singlet and triplet states and appreciable fluorescence rates. _Matter_ 2021, _4_ , 1654–1682
  * Tučková et al. 2022 Tučková, L.; Straka, M.; Valiev, R. R.; Sundholm, D. On the origin of the inverted singlet–triplet gap of the 5th generation light-emitting molecules. _Phys. Chem. Chem. Phys._ 2022, _24_ , 18713–18721
  * Ricci et al. 2021 Ricci, G.; San-Fabián, E.; Olivier, Y.; Sancho-García, J.-C. Singlet-triplet excited-state inversion in heptazine and related molecules: assessment of TD-DFT and ab initio methods. _ChemPhysChem_ 2021, _22_ , 553–560
  * Sanz-Rodrigo et al. 2021 Sanz-Rodrigo, J.; Ricci, G.; Olivier, Y.; Sancho-Garcia, J.-C. Negative singlet–triplet excitation energy gap in triangle-shaped molecular emitters for efficient triplet harvesting. _J. Phys. Chem. A_ 2021, _125_ , 513–522
  * Bhattacharyya 2021 Bhattacharyya, K. Can TDDFT render the electronic excited states ordering of Azine derivative? A closer investigation with DLPNO-STEOM-CCSD. _Chem. Phys. Lett._ 2021, _779_ , 138827
  * Ghosh and Bhattacharyya 2022 Ghosh, S.; Bhattacharyya, K. Origin of the Failure of Density Functional Theories in Predicting Inverted Singlet–Triplet Gaps. _J. Phys. Chem. A_ 2022, _126_ , 1378–1385
  * Koseki et al. 1985 Koseki, S.; Nakajima, T.; Toyota, A. Violation of Hund’s multiplicity rule in the electronically excited states of conjugated hydrocarbons. _Can. J. Chem._ 1985, _63_ , 1572–1579
  * Harris and Jones 1974 Harris, J.; Jones, R. O. The surface energy of a bounded electron gas. _J. Phys. F: Met. Phys._ 1974, _4_ , 1170
  * Langreth and Perdew 1977 Langreth, D.; Perdew, J. Exchange-correlation energy of a metallic surface: Wave-vector analysis. _Phys. Rev. B_ 1977, _15_ , 2884
  * Gunnarsson and Lundqvist 1976 Gunnarsson, O.; Lundqvist, B. Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism. _Phys. Rev. B_ 1976, _13_ , 4274
  * Teale et al. 2010 Teale, A. M.; Coriani, S.; Helgaker, T. Accurate calculation and modeling of the adiabatic connection in density functional theory. _J. Chem. Phys._ 2010, _132_ , 164115
  * Pernal 2018 Pernal, K. Electron Correlation from the Adiabatic Connection for Multireference Wave Functions. _Phys. Rev. Lett._ 2018, _120_ , 013001
  * Pernal 2018 Pernal, K. Exact and approximate adiabatic connection formulae for the correlation energy in multireference ground and excited states. _J. Chem. Phys._ 2018, _149_ , 204101
  * Drwal et al. 2022 Drwal, D.; Beran, P.; Hapka, M.; Modrzejewski, M.; Sokół, A.; Veis, L.; Pernal, K. Efficient Adiabatic Connection Approach for Strongly Correlated Systems: Application to Singlet–Triplet Gaps of Biradicals. _J. Phys. Chem. Lett._ 2022, _13_ , 4570–4578
  * Angeli et al. 2002 Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. $n$-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. _J. Chem. Phys._ 2002, _117_ , 9138
  * Bhaskaran-Nair et al. 2012 Bhaskaran-Nair, K.; Brabec, J.; Aprà, E.; van Dam, H. J. J.; Pittner, J.; Kowalski, K. Implementation of the multireference Brillouin-Wigner and Mukherjee's coupled cluster methods with non-iterative triple excitations utilizing reference-level parallelism. _J. Chem. Phys._ 2012, _137_ , 094112
  * Werner et al. 2012 Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a general-purpose quantum chemistry program package. _Wiley Interdiscip. Rev.: Comput. Mol. Sci._ 2012, _2_ , 242–253
  * Schäfer et al. 1992 Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. _J. Chem. Phys._ 1992, _97_ , 2571–2577
  * Pernal et al. 2023 Pernal, K.; Hapka, M.; Przybytek, M.; Modrzejewski, M.; Sokół, A.; Tucholska, A. GammCor code. https://github.com/pernalk/GAMMCOR, 2023
  * Aprà et al. 2020 Aprà, E. et al. NWChem: Past, present, and future. _J. Chem. Phys._ 2020, _152_ , 184102
  * Neese 2012 Neese, F. _Wiley Interdiscip. Rev.: Comput. Mol. Sci._ 2012, _2_ , 73–78
  * Becke 1988 Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. _Phys. Rev. A_ 1988, _38_ , 3098
  * Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. _Phys. Rev. B_ 1988, _37_ , 785
  * Miehlich et al. 1989 Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Results obtained with the correlation energy density functionals of Becke and Lee, Yang and Parr. _Chem. Phys. Lett._ 1989, _157_ , 200–206
  * Pastorczak and Pernal 2018 Pastorczak, E.; Pernal, K. Correlation Energy from the Adiabatic Connection Formalism for Complete Active Space Wave Functions. _J. Chem. Theory Comput._ 2018, _14_ , 3493–3503
  * Pastorczak and Pernal 2018 Pastorczak, E.; Pernal, K. Electronic Excited States from the Adiabatic-Connection Formalism with Complete Active Space Wave Functions. _J. Phys. Chem. Lett._ 2018, _9_ , 5534–5538
  * Pastorczak et al. 2019 Pastorczak, E.; Hapka, M.; Veis, L.; Pernal, K. Capturing the Dynamic Correlation for Arbitrary Spin-Symmetry CASSCF Reference with Adiabatic Connection Approaches: Insights into the Electronic Structure of the Tetramethyleneethane Diradical. _J. Phys. Chem. Lett._ 2019, _10_ , 4668–4674
  * Beran et al. 2021 Beran, P.; Matoušek, M.; Hapka, M.; Pernal, K.; Veis, L. Density matrix renormalization group with dynamical correlation via adiabatic connection. _J. Chem. Theory Comput._ 2021, _17_ , 7575–7585
  * Matoušek et al. 2023 Matoušek, M.; Hapka, M.; Veis, L.; Pernal, K. Toward more accurate adiabatic connection approach for multireference wavefunctions. _J. Chem. Phys._ 2023, _158_
  * Angeli et al. 2002 Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. $n$-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. _J. Chem. Phys._ 2002, _117_ , 9138