Probing disorder in 2CzPN using core and valence states

Abstract

Molecules which exhibit thermally activated delayed fluorescence (TADF) show great promise for use in efficient, environmentally-friendly OLEDs, and thus the design of new TADF emitters is an active area of research. However, when used in devices, they are typically in the form of disordered thin films, where both the external molecular environment and thermally-induced internal variations in parameters such as the torsion angle can strongly influence their electronic structure. In this work, we use density functional theory and X-ray photoelectron spectroscopy to investigate the impact of disorder on both core and valence states in the TADF emitter 2CzPN (1,2-bis(carbazol-9-yl)-4,5-dicyanobenzene). By simulating gas phase molecules displaying varying levels of disorder, we assess the relative sensitivity of the different states to factors such as varying torsion angle. The theoretical results for both core and valence states show good agreement with experiment, thereby also highlighting the advantages of our approach for interpreting experimental spectra of large aromatic molecules, which are too complex to interpret based solely on experimental data.

---

1 Introduction

Organic light-emitting diodes (OLEDs) have established themselves as a leading technology in a variety of commercial devices, ranging from lighting to full-color displays.^1–3 However, their internal quantum efficiency (IQE) is not always optimal due to the unfavourable spin statistics observed when charges combine to create excitons,⁴ a situation which is typically remedied by employing emitters containing heavy metals such as Ir or Pt. A promising alternative to such potentially costly and environmentally harmful compounds is provided by a novel class of OLEDs which can be fully organic, and rely on the mechanism of thermally activated delayed fluorescence (TADF). TADF emission is characterised by reverse intersystem crossing taking place between non-radiative triplet states, which act as a reservoir to the first excited radiative singlet state, giving a maximum potential IQE of 100%.⁵ Such a mechanism is facilitated by a small singlet–triplet energy splitting, ΔE_ST. Being able to predict which molecules display this mechanism, as well as understand which factors influence ΔE_ST, is therefore crucial for determining the structural features that facilitate TADF.

Established computational tools such as density functional theory (DFT)^6,7 are commonly used to model the excited states of ideal isolated gas phase TADF molecules. However, these molecules interact with their neighbouring molecules in the thin films used in devices, and are also susceptible to conformational fluctuations due to thermal effects. Therefore, in order to gain a more realistic description of TADF materials it is crucial to also properly account for thermal fluctuations of the molecular conformations and microscopic electronic polarization effects in amorphous materials.^8,9 For example, Olivier et al.¹⁰ combined classical molecular dynamics (MD) with time-dependent DFT (TD-DFT) to study two TADF molecules, including 2CzPN (1,2-bis(carbazol-9-yl)-4,5-dicyanobenzene). They demonstrated the influence of conformational fluctuations on emission properties, highlighting the effectiveness of coupling experimental observations with theoretical insights for investigating disorder in TADF materials.

While understanding the effects of disorder on excited states is crucial for applications, excited state methods are challenging in terms of both cost and accuracy. Nonetheless, meaningful information can also be gained by looking at the influence of disorder on the ground state electronic structure, an area which is relatively unexplored for many TADF emitters. To this end, we present a combined computational and experimental investigation into the ground state electronic structure of the prototypical TADF emitter 2CzPN, including both core and valence states. While valence states can be used to gain useful insights into the effects of disorder, core states have the additional advantage that they can be directly mapped to specific atomic environments, and thus experimental approaches such as core X-ray photoelectron spectroscopy (XPS) can be a useful tool for acting as a local probe of disorder. However, the interpretation of core XPS typically relies on comparison to reference spectra, which is challenging in the case of large molecules containing many different local chemical environments. Despite having received relatively little attention compared to valence state calculations, the theoretical prediction of core binding energies (BEs) is invaluable for the interpretation of core XPS. While calculations based on more expensive quantum chemical approaches or methods such as GW^11,12 show great promise, DFT has the advantage of being less computationally expensive, thereby allowing the treatment of much larger systems.

One promising DFT-based method for core BE calculations is ΔSCF, in which an electron is removed from a core state of interest, with the corresponding BE being defined as the difference between the energy of the core-excited state and that of the ground state. It thereby enables final state effects to be taken into account with modest computational cost, and has been successfully employed in the context of both Hartree Fock¹³ and DFT.¹⁴ ΔSCF has been successfully employed for benchmark calculations of a range of molecules, however, it is not yet widely used in applications, and, with relatively few exceptions (e.g. ref. 15–20), most calculations have also been restricted to small molecules. One challenge for ΔSCF in Gaussian basis codes is the need to either use very large or specifically constructed basis sets which are adapted for core-excited calculations, e.g. by including basis functions for the Z + 1 element.²¹ We recently introduced an alternative approach, which uses an adaptive multi-wavelet basis, as implemented in the MADNESS code.²² The multi-wavelet basis leads to a guaranteed, high precision without requiring any user intervention or specialist knowledge. In combination with a mixed all-electron/pseudopotential scheme, this facilitates a highly accurate and computationally efficient approach,²³ which also avoids the common problem of core-hole hopping and associated convergence problems. We have previously applied this approach to both simple²³ and aromatic amino acids,²⁴ with final state effects playing a significant role with the latter, so that experimental results cannot be interpreted using chemical intuition alone. In this work we investigate the effects of disorder on the core BEs of 2CzPN, demonstrating how our approach may be combined with experiment to provide a powerful tool for investigating core states in large molecules.

2 Methods

2.1 Atomic structure

The atomic structure of 2CzPN is depicted in Fig. 1a, and is characterised by three ‘subfragments’ (highlighted in Fig. 1b): a central dicyano-subsituted phenylene core, which is linked to two carbazole units, whose relative orientations are described by torsion angles θ₁ and θ₂. The energy barrier is small enough such that thermal effects may lead to disorder in the form of diverse torsion angles being present in thin films of 2CzPN. This variation can in turn influence both the ground and excited state electronic structure. Starting from an amorphous MD snapshot extracted from an equilibrated simulation, which was previously used to explore such effects on excited states,¹⁰ in the following we assess the influence of various levels of disorder on the ground state electronic structure of 2CzPN. In order to reduce the computational cost, while still sampling enough of the conformational space, a cluster of 500 molecules was extracted out of a given 1000 molecule MD snapshot. The focus of this work is the internal disorder present in 2CzPN, and thus each of the 500 molecules were calculated in separate gas phase calculations without further relaxation, while also being compared with the relaxed gas phase molecule. Free boundary conditions are used throughout.

Fig. 1 Depiction of the atomic structure of 2CzPN, showing the torsion angles in (a). The subfragments depicted in (b) are used to perform the rigidification in (c). (d) Shows the atoms grouped according to similar chemical environments and thus core binding energies, and their corresponding labels.

Since both torsion angles may vary independently, using the angle alone does not provide a single parameter which can characterise the structure of a molecule. Although one could use the average value, some information is lost when comparing how similar two structures are. We therefore propose an alternative parameter, namely the cost function J, which is defined as

where [scr R, script letter R]^β→α_ab is a rotation matrix between two instances of a molecule characterised by their atomic coordinates R^α and R^β, and N_at is the number of atoms in the molecule. The optimal rotation, [scr R, script letter R]^β→α_ab, may be efficiently found by minimising J,^25–27 where J = 0 implies a rigid transformation. J has previously been employed in BigDFT as a measure of similarity between fragment geometries for the purposes of reusing basis functions between calculations,^28–30 and thus provides a computationally inexpensive means of comparing extracted molecules. At present, inversion symmetry is not taken into account, this could, however, be included in the future.

A small value of J_αβ implies that molecules α and β have similar torsion angles, while also not showing any other significant differences in geometry. Indeed, beyond variations in the torsion angles, thermal fluctuations also result in other distortions, which may be captured by the use of J. This is evidenced by the fact that in the relaxed molecule, the three subfragments shown in Fig. 1b are close to planar, while the equivalent subfragments in the molecules extracted from the MD snapshot show noticeable distortions, and are no longer planar. In order to assess the impact of this additional disorder, for each of the 500 extracted molecules, we also generate a ‘rigidified’ molecule. This is achieved using the relaxed molecule as a template. First, the rototranslations between each of the subfragments in the template and those of the target molecule are identified. Then the rototranslations are applied to each subfragment in the template molecule, thereby constructing a new molecule which preserves the internal structure of the subfragments of the relaxed molecule. An example molecule to which this process has been applied is depicted in Fig. 1c.

Following the rigidification, there are thus two sets of unaltered and rigidified molecules, containing 500 molecules each. In order to reduce the number of computationally intensive calculations to be performed, while nonetheless maintaining a balanced coverage of configuration space, the cost function may be used to identify a subset of distinct molecules for each set. This is achieved by means of a hierarchical clustering approach, as implemented in SciPy,³¹ using the matrix J_αβ as the linkage matrix, where α and β run over all molecules in the set. Based on a threshold value of J, the molecules are grouped into clusters, with a lower value of J giving rise to more clusters. One molecule from each cluster is then selected as a representative. In this work, the first molecule from each cluster is selected, but alternative choices could be used in future. The representative molecules are then used to perform both the core calculations and hybrid functional valence state calculations, where the computational cost of calculating all molecules would be very high. The clustering is performed for the unmodified extracted molecules only; the corresponding rigidified molecules are directly employed.

Finally, as an independent comparison to the extracted molecules, a set of explicitly rotated molecules were also generated. Initially these were generated from the DFT-relaxed molecule, keeping the central subfragment fixed, and rotating the two carbazole units in tandem. However, the interaction between the two carbazole units in the relaxed molecule is such that a small degree of asymmetry is introduced, resulting in an asymmetric potential energy surface (see Fig. S3 in the ESI†). Therefore, the molecule was also relaxed with constraints on the symmetry and torsion angle, with θ₁ ≡ θ₂ ≡ 90°. Twenty additional rotated molecules, with angles ranging from 40° to 140° were then generated starting from the constrained molecule, keeping θ₁ ≡ θ₂. Due to the imposed symmetry, only molecules with θ_Av. ≤ 90° were explicitly calculated. The symmetry-constrained molecule was also used as an alternative template for performing the rigidification, however the differences were found to be small (less than 0.1 eV average difference in the PBE-calculated frontier orbital energies), and so the unconstrained relaxed molecule was retained as the template for the rigidification.

To summarise, three sets of molecules with varying levels of disorder were therefore generated:

(1) Rotated molecules: the torsion angle may vary, but θ₁ ≡ θ₂, while the subfragments are kept rigid;

(2) Rigid molecules: both torsion angles may vary independently, while rigidity is enforced in the subfragments;

(3) Extracted molecules: no constraints are present, with all atoms able to move freely.

By calculating core and valence quantities for all three sets of molecules, we are therefore able to assess the relative importance of the different types of disorder present in 2CzPN.

2.2 Density functional theory

Valence state calculations were performed using the BigDFT code,³⁰ which uses a Daubechies wavelet basis set.³² Calculations employed the cubic scaling version of BigDFT,³³ using HGH-GTH pseudopotentials (PSPs)^34,35 with non-linear core corrections.³⁶ A wavelet grid spacing of 0.5 bohr was used, with coarse (fine) multipliers of 5 (7). The PBE functional was used for all molecules,³⁷ while only the representative molecules described above were calculated using the hybrid PBE0 functional.³⁸ The molecule was relaxed using a smaller wavelet grid spacing of 0.45 bohr, using a maximum force threshold of 0.02 eV Å⁻¹. Projected densities of states (PDOS) were calculated using a Mulliken-style population analysis, with Gaussian smearing of 0.44 eV applied to reflect the experimental broadening. For comparison with the valence XPS measurements, Scofield photoionisation cross sections for 1.48667 keV^39,40 were applied to the calculated PDOS using the Galore package.⁴¹ Theoretical spectra were aligned and normalised to the lowest binding energy peak in the experimental spectrum.

Core binding energies (BEs) were calculated using the MADNESS code,²² which uses a multi-wavelet approach to achieve guaranteed high precision and can be used for all-electron (AE), PSP, or mixed AE/PSP calculations.⁴² Calculations were performed at both the Koopmans' and ΔSCF level of theory, using a recently developed approach in which the core-excited atom is treated at the AE level, and the remainder of the atoms are treated at the PSP level, as described in detail in ref. 23. This approach allows explicit access to the core states of interest, while forcing the core hole to remain localised, and at the same time reducing the computational cost compared to a pure AE approach. Following ref. 23 and 24, ground state calculations used a wavelet threshold of 10⁻⁴ followed by 10⁻⁶ (wavelet orders k = 6 and k = 8 respectively), with core hole calculations using a wavelet threshold of 10⁻⁶ (k = 8). HGH-GTH PSPs were used, while all calculations were spin restricted. Except where otherwise stated, calculations employed PBE. Since we are only interested in relative BEs, relativistic effects were neglected. All other parameters are the same as in ref. 23 and 24. A combination of 0.44 eV Gaussian smearing based on the experimental resolution and a 0.2 eV Lorentzian smearing to reflect lifetime contributions was applied to the calculated BEs. C 1s theoretical spectra were aligned and normalised to the main experimental peak, while N 1s spectra were aligned to the higher BE peak.

2.3 Experimental approach

X-ray photoelectron spectra were collected on a Thermo Scientific K-Alpha+ spectrometer at Imperial College London, London, UK. The system uses a microfocused Al Kα X-ray source with a photon energy of 1486.7 eV. The X-ray source is operated at 12 kV and 6 mA, with a base pressure in the analysis chamber of 2 × 10⁻⁹ mbar. The X-ray spot size used for the measurements was 400 μm. The 2CzPN powder, acquired from Ossila (Sublimed >99% (HPLC)), was mounted on the sample plate using conductive carbon tape and the spectrometer flood gun was used to compensate for possible sample charging.

In order to avoid any beam induced sample damage, the sample was rastered and the main core level spectra (C 1s and N 1s) and valence spectra measured at equidistant spots across a defined rectangular region on the sample (6 × 3 for the core and 5 × 3 for the valence spectra). This method enabled each spot to be irradiated for a significantly shorter time than would be necessary in a single spot measurement, achieving equally good signal to noise. The core level spectra were collected at a pass energy of 20 eV and a dwell time of 50 ms at each measurement point, whilst the valence spectra were measured at 40 eV pass energy and 75 ms dwell time. The survey spectrum, obtained from a single spot, was measured at 200 eV pass energy. The Multiplex function on the data acquisition software (Thermo Scientific Avantage) enabled the measurement of single spectra of each core level or valence band until the predefined total number of scans had been collected, in order to equally distribute the X-ray dose received. The Thermo Scientific Avantage software package was also used for peak fit analysis of the collected spectra. Details of the fit procedure can be found in the ESI.†

3 Results

3.1 Atomic structure

Before discussing the electronic structure of 2CzPN, we first consider the atomic structure, including the variation in torsion angle. Although fewer molecules are considered in this work compared to the full trajectory of ref. 10, as shown in Fig. 2a, 500 molecules are nonetheless sufficient to see the observed bimodal distribution in the average torsion angle, with peaks at around 60° and 120°. In the PBE-relaxed molecule θ₁ = θ₂ = 51.9°; although there has been shown to be some sensitivity with respect to functional choice, this value nonetheless falls within the range of values calculated using different functionals.⁴³ In contrast to the relaxed molecule, Fig. 2b shows significant differences between θ₁ and θ₂ for a large number of the molecules, by more than 50° in some cases. At the same time, there is also disorder present in the bond lengths, with C–N and C–C bonds showing a wider distribution of ranges compared to the relaxed molecule (see Fig. S1 in the ESI†), further supporting the need for a metric for comparing molecules which goes beyond the average torsion angle. To this end, Fig. 2b also shows the relationship between the torsion angles and the cost function J between each molecule and the relaxed structure. As expected, molecules with smaller torsion angles, close to the relaxed value, have a smaller cost function value, in other words J very loosely follows the difference in average torsion angle, while also taking into account other differences in structure.

Fig. 2 (a) Distribution of average torsion angle, θ_Av., across the 500 unmodified extracted molecules; and (b) comparison between the two torsion angles, θ₁ and θ₂, with the corresponding cost function J_i,j=relaxed between each molecule and the relaxed gas phase structure. Also depicted on both plots are the corresponding angles for the relaxed molecule.

The full cost function matrix between the 500 extracted molecules is depicted in Fig. 3a, with the corresponding dendrogram showing the hierarchical clusters in Fig. 3b. By varying the cost function threshold, the number of clusters can also be varied. Using a cost function threshold of J = 1.2 a₀ results in 42 clusters of similar molecules; the first molecule of each cluster is depicted in Fig. 3c. These molecules show a relatively well distributed range of torsion angles, while there is also a mix of molecules with both similar torsion angles and larger differences between the two torsion angles (see also Fig. S2 in the ESI†). Similarly, the molecules are well distributed across the total energy landscape, with sufficient sampling to give a reasonable approximation of both the probability density and average energy for both the unmodified extracted and rigidified molecules (Fig. S4 in the ESI†). This gives confidence that the chosen threshold is sufficient to capture the variety of torsion angles and other disorder present in the full dataset, an assumption which is further validated in the following investigation of the valence states.

Fig. 3 Depiction of (a) the cost function matrix, (b) the resulting hierarchical relationship, and (c) the representative molecules.

3.2 Valence states

Before investigating the core states, we first consider the valence states as well as the HOMO–LUMO gap. The relatively lower computational cost of calculating valence states allows the choice of selecting molecules to be validated, while it is also interesting to compare the relative sensitivity to disorder between core and valence states.

3.2.1 Frontier orbitals. First, in order to assess the effectiveness of our approach for selecting representative molecules, including the chosen threshold for J, PBE calculations were performed for all 500 extracted molecules, as well as the 500 corresponding rigidified molecules. The frontier orbital energies and HOMO–LUMO gaps are depicted in Fig. 4a–c alongside the results for the explicitly rotated molecules. Fig. 4d–f show the same quantities, in the case where only the 42 select molecules are considered (both the unmodified and rigidified versions). As with the total energy landscape, the representative molecules are able to capture the general trends well, including both the distribution and the average values. In the following we therefore consider the results for the representative molecules only.

Fig. 4 Variation in the highest occupied molecular orbital energy, ε_HOMO [left], the lowest unoccupied molecular orbital energy, ε_LUMO [centre], and the HOMO–LUMO gap, E_g [right], with respect to the average torsion angle θ_Av., as well as the probability density. Results are shown for PBE calculations of all the molecules considered [top], and for PBE and PBE0 for select representative molecules [middle and bottom respectively]. Horizontal lines denote average values, except in the case of the rotated molecules, where they correspond to the relaxed molecule. The data include the unmodified extracted molecules, rigidified molecules, and explicitly rotated molecules.

Considering first the frontier orbital energies, the highest occupied molecular orbital (HOMO) shows a weaker sensitivity to the torsion angle than the lowest unoccupied molecular orbital (LUMO), so that the trend in HOMO–LUMO gap is more similar to that of the LUMO. As shown in Fig. 5, the HOMO is more concentrated on the carbazole units, while the LUMO is centred on the phenylene core, suggesting that this central core is more sensitive to changes in torsion angle, as will be discussed in more detail in the context of core states. The HOMO and LUMO also follow different trends with respect to the change in torsion angle, with the HOMO showing a local maximum at 90° and two local minima around 50° and 130°. In contrast, both the LUMO and the HOMO–LUMO gap have a minimum at 90°.

Fig. 5 Depiction of the PBE-calculated HOMO and LUMO of the relaxed molecule, visualized in VESTA⁴⁴ using an isosurface value of 0.002 a₀^−3/2.

Considering next both the unmodified and rigidified extracted molecules, the trend with respect to average torsion angle is less evident for all three quantities. In other words, the effects of independently varying torsion angles, and to a lesser extent, variation in bond length between subfragments, lead to a larger spread in values. Indeed, the HOMO–LUMO gap of the unmodified extracted values varies over a range of 0.8 eV compared to 0.6 eV for the rotated molecules, a larger range than that over which the singlet–triplet splitting was previously found to vary.¹⁰ In all cases the variation is large enough such that the effects of thermal fluctuations cannot be neglected. In addition, there is a net shift in the average value away from the values for the relaxed molecule, which reflects the distribution of angles present in the original structure. On the other hand, when going from the rigidified to the unmodified extracted molecules, the changes are much smaller, indicating that the effect of non-rigidity is relatively unimportant.

The above results are with PBE, however semi-local functionals such as PBE are expected to significantly underestimate the gap compared to experiment. Therefore, the frontier orbital energies of the representative molecules were also calculated using PBE0, the results for which are depicted in Fig. 4g–i. Unsurprisingly, there is a net shift in values of around −0.8 eV for the HOMO and 0.8 eV for the LUMO, leading to a shift in the PBE0 gaps of around 1.6 eV relative to the PBE values. However, the overall trends are similar, as is further evidenced in Fig. S5 in the ESI.† The external environment can have an important effect on the measured values of the HOMO, LUMO and HOMO–LUMO gap, contributing to the wide range of experimental values found in the literature. In addition, even in the case where the same solvent (toluene) is used, significant differences exist.^45,46 Nonetheless, the reported HOMO–LUMO gaps lie between the PBE and PBE0 values for the relaxed molecule (2.1 and 3.7 eV, respectively). This is also the case for the measured HOMO and LUMO values. Although the differences between experimental and theoretical setups limit the extent to which agreement should be expected, it is nonetheless interesting to note that the average PBE0 gap over the extracted molecules (3.5 eV) lies closer to the experimental values than the PBE0 value for the relaxed molecule (3.7 eV), as do the average HOMO and LUMO values.

3.2.2 Valence spectra. Having discussed the frontier orbital energies, we now consider the full range of occupied states. These are also compared to measured valence XPS data, by applying appropriate photoionisation cross-sections to the PDOS, as described in Section 2.2. In order to verify that the use of representative molecules remains a reasonable approximation, the DOS was calculated using PBE for all molecules. The averaged DOS for the full set of molecules is indistinguishable from the average over the representative molecules, for both the rigidified and unmodified extracted molecules (see Fig. S6 in the ESI†). Therefore, as in the previous section we only consider the 42 representative molecules in the following.

Before investigating the effects of disorder, we first consider the unweighted PDOS for the relaxed molecule, which is depicted in Fig. 6a for PBE0 (see Fig. S7a in the ESI† for PBE). The PBE0 and PBE PDOS exhibit similar features and relative peak heights, with the key difference being that the PBE0 DOS has a larger overall bandwidth. This also affects the DOS close to the valence band edge, where a shoulder in the PBE DOS becomes a separate peak in the PBE0 DOS. On the other hand, the projection onto the different atomic orbitals is largely unaffected by the functional. The cross-section weighted PDOS for PBE0 as well as the experimental valence spectrum are shown in Fig. 6b (see Fig. S7b in the ESI† for PBE). The experiment shows a number of clearly discernible features and the PBE0 simulated spectra match experiment much more closely than PBE regarding both the total width of the valence states as well as their relative positions. Indeed, the PBE0 spectrum shows very good agreement with the experimental binding energy positions of the valence features. The relative theoretical peak heights are in less good agreement with experiment, in large part due to an overestimation of the N s contribution. This is most likely a limitation of the employed projection scheme; in future work it could therefore be interesting to test alternative projection schemes.

Fig. 6 PBE0-calculated occupied projected density of states for the relaxed molecule, both (a) unweighted and (b) photoionisation cross-section weighted. Also shown are the sum over contributions for the relaxed molecule, the average over select unmodified extracted molecules, and the average over rigidified molecules, for both (c) unweighted and (d) photoionisation cross-section weighted. The weighted spectra have been aligned to experiment.

We now consider the effects of disorder, for which the average extracted and rigid results are compared to the relaxed molecule in Fig. 6c and d for PBE0 (and Fig. S7c and d in the ESI† for PBE). Although the unweighted PDOS shows some sensitivity to average torsion angle across the full range of energies (see also Fig. S8 in the ESI†), for the weighted PDOS significant deviations between the averaged disordered spectra are only present for higher binding energies (above 15 eV). These changes are not attributable to one single factor, with the extracted, rigid and relaxed spectra deviating to different extents for different peaks. In other words, as well as showing differences in sensitivity to average torsion angle, some peaks are more sensitive to the effects of non-rigidity than others. This also proves to be the case for the core states, as will be demonstrated in the following. The spectrum was measured for crystalline 2CzPN, where it is not clear to what extent the range of accessible torsion angles might be limited by intermolecular interactions. Nonetheless, the average spectra over extracted molecules is in slightly better agreement with respect to experiment in terms of peak heights, compared to the relaxed molecule. Although subtle, this suggests that the presence of disorder does indeed have an influence on the spectrum, an assumption which is also examined for the core states.

3.3 Core states

3.3.1 Level of theory. The computational cost of calculating core BEs using ΔSCF is significantly higher than for valence quantities, since each C and N atom requires a separate core-excited calculation for each molecule. Furthermore, due to the presence of disorder, the number of atoms to be treated cannot be reduced by relying on symmetry, as could be done for example in an ideal crystal, since each atom may have its own distinctive environment. In order to keep the computational cost tractable for core BE calculations we therefore do not consider the full set of molecules, only the representative molecules. Based on the results of Section 3.2, we expect that this should give sufficient coverage of the potential energy landscape.

Given the high computational cost, it is important to balance cost and accuracy when deciding the level of theory. Koopmans' approach enables the calculation of core BEs at much lower cost than ΔSCF, since unlike ΔSCF, Koopmans' BEs can be calculated using a single (all electron) calculation per molecule. Although there is often a large shift in absolute BEs between the two approaches, in this work we focus on relative BEs, and thus this shift can be ignored. Indeed, the initial state effects which are taken into account using Koopmans' approach¹⁴ have been shown to provide a good estimate of relative BEs in many cases.⁴⁷ On the other hand, ΔSCF is able to take into account final state effects, with some limitations, e.g. multiplet splittings. In some systems, as we have previously shown for example in the aromatic amino acids,²⁴ these final state effects, which are not captured by Koopmans', can sometimes lead to a quantitatively or even qualitatively wrong BE order. It is therefore important to determine how significant final state effects are for 2CzPN, and thus whether Koopmans' approach is sufficient, or if ΔSCF needs to be employed.

A comparison between the two approaches is shown in Fig. S9 in the ESI.† Aside from the expected large systematic offset in absolute energies, the Koopmans' and ΔSCF BEs show similar trends for most atoms, although the Koopmans' BEs are more spread out compared to ΔSCF. In other words, for the majority of environments, the initial state effects capture the relative BEs quite well, for both C and N 1s BEs. However, the C atoms in the cyano group follow a notably different trend, with their BEs being significantly underestimated using the Koopmans' approach relative to the BEs associated with other C atoms. Given the apparent importance of final state effects for the cyano C atoms, we therefore employ ΔSCF for all results in the following.

The choice of exchange correlation functional is known to have a significant effect on calculated core BEs (see e.g. ref. 48), indeed we have previously observed an improvement in agreement with experiment when using PBE0 compared to PBE, in the case of amino acids.^23,24 Although prohibitively expensive to calculate ΔSCF BEs for all molecules using PBE0, we therefore nonetheless performed a comparison between PBE and PBE0 for the relaxed molecule, to determine to what extent the results for 2CzPN are sensitive to the functional. In addition to a systematic offset, the C 1s BEs are more spread out for PBE0 than PBE, however, the PBE and PBE0 results are well correlated, with the cyano C atoms being much less sensitive to the choice of functional than the inclusion of final state effects (see Fig. S10 in the ESI†). In other words, PBE represents a reasonable approximation for C 1s, and so except where otherwise stated, all results in the following are for PBE. For N 1s the BEs are more strongly affected by the choice of functional than C 1s, with the separation between the two pairs of BEs being 1.7 eV for PBE compared to 1.4 eV for PBE0. In addition to the relaxed molecules, N 1s BEs were therefore also calculated for three additional rotated molecules, as shown in Table 1. Importantly, the difference between the PBE and PBE0 separations is consistently 0.3 eV, so that although the relative BEs are sensitive to the choice of functional, the trends for PBE and PBE0 are nonetheless expected to be similar.

Table 1 N 1s binding energies and binding energy differences, Δ, in eV, calculated for select explicitly rotated molecules using PBE and PBE0, and compared to experiment. The experimental values have an estimated error of ±0.1 eV

	PBE			PBE0
	N1	N2	Δ	N1	N2	Δ
40.0°	413.0	411.2	1.8	418.6	417.1	1.5
51.9° (relaxed)	413.0	411.3	1.7	418.5	417.2	1.4
65.0°	412.9	411.4	1.5	418.5	417.3	1.1
90.0°	412.8	411.8	1.0	418.4	417.7	0.7
Experiment	400.5	399.5	1.0

---

3.3.2 Experimental spectra. X-ray photoelectron spectra of the C 1s and N 1s core levels of 2CzPN are presented in Fig. 7. Fig. 7e and f show the main chemical state contributions identifiable from peak fit analysis. The C 1s spectrum exhibits three features at 284.5, 285.6 and 286.6 eV. Based on the theoretical BE values, these three peaks can be assigned to specific C atoms within the 2CzPN molecule as follows. The lowest BE contribution corresponds to the C₁ and C₂ atoms of the two carbazole units that are only bound to other C and H atoms, commensurate with basic arguments of chemical shift in carbon spectra. In contrast, the peak at 285.6 eV, assigned to C₃ and C₄ atoms, and the peak at 286.6 eV, associated with the C₅, C₆ and C₇ atoms, show BE positions not fully expected from basic arguments and experimental experience. The BE position of the C₃ and C₅ atoms are considerably higher than what is expected from atoms bound only to other C and H atoms. This observation is commensurate with our previous work on aromatic amino acids, where next-nearest neighbour effects as well as the influence of aromatic rings on BEs are in effect.²⁴ The relative peak intensities in the experiment match well with the expected atomic ratios for the C₁ and C₂ peak (61 rel. at% from experiment and 62.5% from molecular structure) as well as with the sum of the two higher BE contributions (39 rel. at% from experiment and 37.5% from molecular structure). It is worth noting, that whilst the latter two features would be expected to appear in a 1 : 1 intensity ratio from the theoretical calculations, the experiment shows relative contributions of 24 and 15 rel. at%.

Fig. 7 Comparison between experimental (dots) and theoretical (lines) core spectra, including the relaxed molecule, the average over select unmodified extracted molecules, and the average over rigidified molecules, calculated using PBE for (a) and (c) C 1s and (b) and (d) N 1s. The colours correspond to different atomic sites, as depicted in Fig. 1d. Also shown are vertical lines denoting the calculated binding energies for the relaxed molecule in (a) and (b) and the unmodified extracted molecules in (c) and (d). Theoretical results have been aligned to experiment as described in Section 2.2. Also shown are the experimental peak fits for (e) C 1s and (f) N 1s.

The N 1 s core level spectrum has two major contributions at 399.5 and 400.5 eV assigned to the N₁ and N₂ atoms of the molecular structure, respectively. The relative contributions of the two N environments from peak fit analysis is 44 and 56 rel. at%, respectively. In addition to the discussed expected chemical environments of both C and N, the core level spectra show small additional contributions from surface species, expected for samples handled in air.

3.3.3 Influence of disorder. Considering now the effect of the torsion angle on the core state BEs, Fig. 8c and f show the dependence on average torsion angle for the rotated molecules. Although there are variations between both the degree of sensitivity to the torsion angle and the exact shape of the trend, there are two distinct scenarios: the atoms associated with the phenylene core (C₃, C_5–7 and N₁) show a maximum at 90° and a higher sensitivity to the torsion angle, while the atoms associated with the carbazole units (C_1–2, C₄, N₂) show a much shallower minimum at 90°. This confirms earlier observations that the phenylene core is more sensitive to changes in torsion angle (see Section 3.2.1). The strong sensitivity of the N₁ BE to torsion angle leads to a large variation in BE separation between N₁ and N₂ – ranging from 1.1 to 1.8 eV (Table 1), a variation which is clearly visible in the calculated spectra (see Fig. S11 in the ESI†). In contrast, the C 1s theoretical spectra show much more subtle spectral changes with varying torsion angle.

Fig. 8 Variation in core binding energies with respect to the average torsion angle θ_Av., as well as the probability density, for C 1s [top] and N 1s [bottom], where the colours correspond to different atomic sites, as depicted in Fig. 1d. Horizontal lines denote average values, except in the case of the rotated molecules, where they correspond to the relaxed molecule. The data include (a) and (d) the unmodified extracted molecules, (b) and (e) the rigidified molecules, and (c) and (f) the explicitly rotated molecules.

Going beyond the explicitly rotated molecules, Fig. 8 also shows the variation in core BEs for the extracted and rigidified molecules. As with the frontier orbitals, the additional disorder dilutes the strong trend with respect to torsion angle in both cases, leading to wider variations in BEs. The extent to which particular BEs are sensitive to the different levels of disorder varies. For example the N₁ BEs show similar spreads for the rotated, rigidified and extracted molecules (0.6, 0.5 and 0.7 eV respectively), whereas N₂ BEs have a spread of 0.2 eV for the rotated molecules, increasing to 0.5 eV for both the rigidified and extracted molecules (see also Fig. S12 and Table S1 in the ESI†). The C 1s BEs also show different patterns of spreads, however, there is no obvious trend concerning which atoms are more sensitive to rigidity than others. Despite this variation in behaviour, both the ranges of values and shifts in average relative to the relaxed molecule are the same order of magnitude as for the frontier orbitals.

As shown by comparing Fig. 7a and c, the large spreads in C 1s BEs of the extracted molecules (up to 0.8 eV) nonetheless have only a relatively subtle impact on the overall spectrum. Indeed, while the C 1s spectrum for the relaxed molecule already agrees well with experiment, the higher BE shoulder reduces slightly in intensity going from the relaxed to rigid to extracted molecules, leading to a slightly improved experimental agreement. This is in part due to the spreading out and resulting decreased separation between C₃ and C_5/6 BEs, which also helps explain the discrepancy between the experimental and theoretical relative contributions of the two higher BE peaks.

On the other hand, the variation in N 1s BEs (up to 0.7 eV) lead to larger variations in the theoretical spectra. Indeed, as shown in Fig. 7b and d, the extracted and rigid molecules give rise to a better match to experiment than the relaxed molecule, although the BE separation is still considerably overestimated. However, as previously discussed this separation is sensitive to both the functional and the torsion angle. Indeed, as shown in Table 1, the PBE0 separation for the relaxed molecule is 0.4 eV higher than experiment, whereas the PBE value is 0.7 eV higher (see also Fig. S13 in ESI†). At the same time, the difference between average N₁ and N₂ PBE-calculated BEs for the extracted molecules reduces to 1.4 eV, compared to 1.7 eV in the relaxed molecule.

The high computational cost prohibits calculating the BEs for all of the extracted molecules with PBE0. However, we may instead employ a simplified variation of the so-called spectral warping approach, also known as colour morphing, which has been employed in the context of optical excitations calculated using TDDFT.^49–51 In this approach, a correction is imposed on semi-local calculations by comparing with representative hybrid functional calculations. In this case, we may simply impose a shift of 0.3 eV on the N₁ BEs. The smaller separation of the two peaks in the warped spectrum does indeed lead to a closer match between the theory and experiment, albeit with the overall bandwidth still being slightly overestimated in the simulated spectrum (see Fig. S14 in the ESI†). As a result, the agreement remains slightly worse than for C 1s. However, such a strong sensitivity of the N 1s BEs suggests that they may also be affected by other factors, such as intermolecular interactions with neighbouring molecules. Such effects may also influence the calculated valence spectra, therefore, it would be interesting in future work to assess the impact of going beyond gas phase by including neighbouring molecules. This could be achieved either by extracting a small cluster from the MD trajectory, or by simulating 2CzPN in its crystalline form.

4 Conclusion

In this work, theory and experiment are used to investigate the influence of disorder on the core and valence states of the prototypical TADF emitter 2CzPN, by simulating gas phase molecules with varying levels of disorder. In order to make such calculations tractable, we employ a novel approach for identifying representative molecules drawn from an MD snapshot, which is validated by comparing valence quantities calculated for both the representative molecules and the full data set. The results show that the frontier orbitals exhibit differing trends with respect to torsion angle, with the LUMO also being more sensitive than the HOMO. However, the influence on the simulated X-ray photoelectron spectra is subtle, with the theoretical spectrum already showing good agreement with experiment without taking into account the disorder.

While excited state properties are more directly applicable for understanding the performance of TADF-based OLEDs in devices, core XPS has the advantage that the influence of disorder can be mapped at a more local level. However, the interpretation of such spectra for large molecules based purely on chemical intuition or reference spectra is challenging. Indeed, in the case of aromatic molecules such an interpretation may even be qualitatively wrong, due to the importance of both non-nearest neighbour and final state effects. As such, the theoretical prediction of core binding energies can be invaluable both for interpreting experiments, and for assessing the influence of disorder. Nonetheless, core BE calculations based on for example DFT are relatively uncommon, and have typically been restricted to small molecules, with a few recent exceptions. In this work, we employ the multi-wavelet-based MADNESS code to calculate core BEs of 2CzPN, in which the ΔSCF approach has recently been implemented, and which is both precise and computationally efficient enough to simulate large molecules. Although we find that Koopmans' approach captures the trends in relative BEs well for most environments in 2CzPN, final state effects prove to be important for the cyano C atoms, motivating the use of ΔSCF. The results show that both C 1s and N 1s BEs are sensitive to disorder, with some sites more strongly influenced by factors such as varying torsion angle than others. The theoretical C 1s spectrum shows excellent agreement with experiment, and, like the valence XPS, is relatively insensitive to disorder. The experimental spectrum is less well predicted for N 1s, while it is also found to be more sensitive to both the employed exchange correlation functional and the disorder. Our approach could be extended in future to also consider environmental effects, by going beyond gas phase to the treatment of e.g. clusters of molecules.

This work demonstrates the importance of taking into account disorder in TADF materials, and the resulting implications for providing an accurate picture of their electronic structure. A similar approach could be used in future to investigate other TADF emitters. Furthermore, beyond the implications for TADF-based OLEDs, this work also demonstrates the benefits of using theory to aid in the interpretation of experimental spectra, especially for large aromatic molecules where a reliable interpretation cannot be achieved purely by relying on comparison to reference spectra.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

LER and MS acknowledge support from an EPSRC Early Career Research Fellowship (EP/P033253/1). NKF acknowledges support from the Engineering and Physical Sciences Research Council (EP/L015277/1). AR acknowledges support from the Analytical Chemistry Trust Fund for her CAMS-UK Fellowship. Calculations were performed on the Imperial College High Performance Computing Service and the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk). LER and MS would like to thank Yoann Olivier for providing the MD snapshot and for useful discussions. Data associated with this paper, including a Jupyter notebook demonstrating the workflow for identifying representative molecules, are available at https://gitlab.com/lratcliff/2czpn.

References

1. F. De Angelis, Acc. Chem. Res., 2014, 47, 3349–3360.
2. Y. Tao, K. Yuan, T. Chen, P. Xu, H. Li, R. Chen, C. Zheng, L. Zhang and W. Huang, Adv. Mater., 2014, 26, 7931–7958.
3. M. Y. Wong and E. Zysman-Colman, Adv. Mater., 2017, 29, 1605444.
4. U. Hiroki, G. Kenichi, S. Katsuyuki, N. Hiroko and A. Chihaya, Nature, 2012, 492, 234–238.
5. Q. Zhang, D. Tsang, H. Kuwabara, Y. Hatae, B. Li, T. Takahashi, S. Y. Lee, T. Yasuda and C. Adachi, Adv. Mater., 2015, 27, 2096–2100.
6. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
7. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
8. P. de Silva, C. A. Kim, T. Zhu and T. Van Voorhis, Chem. Mater., 2019, 31, 6995–7006.
9. C. A. Kim and T. Van Voorhis, J. Phys. Chem. A, 2021, 125, 7644–7654.
10. Y. Olivier, B. Yurash, L. Muccioli, G. D'Avino, O. Mikhnenko, J. C. Sancho-García, C. Adachi, T.-Q. Nguyen and D. Beljonne, Phys. Rev. Mater., 2017, 1, 075602.
11. M. J. van Setten, R. Costa, F. Viñes and F. Illas, J. Chem. Theory Comput., 2018, 14, 877–883.
12. D. Golze, L. Keller and P. Rinke, J. Phys. Chem. Lett., 2020, 11, 1840–1847.
13. P. S. Bagus, Phys. Rev., 1965, 139, A619–A634.
14. N. Pueyo Bellafont, P. S. Bagus and F. Illas, J. Chem. Phys., 2015, 142, 214102.
15. T. Susi, D. J. Mowbray, M. P. Ljungberg and P. Ayala, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 081401.
16. J. M. Kahk and J. Lischner, Phys. Chem. Chem. Phys., 2018, 20, 30403–30411.
17. J. M. Kahk and J. Lischner, Phys. Rev. Mater., 2019, 3, 100801.
18. M. Figueras, I. J. Villar-Garcia, F. Viñes, C. Sousa, V. A. de la Peña OShea and F. Illas, J. Phys. Chem. C, 2019, 123, 11319–11327.
19. R. Costa, A. Morales-García, M. Figueras and F. Illas, Phys. Chem. Chem. Phys., 2021, 23, 1558–1565.
20. J. M. Kahk, G. S. Michelitsch, R. J. Maurer, K. Reuter and J. Lischner, J. Phys. Chem. Lett., 2021, 12, 9353–9359.
21. M. W. Hanson-Heine, M. W. George and N. A. Besley, Chem. Phys. Lett., 2018, 699, 279–285.
22. R. J. Harrison, G. Beylkin, F. A. Bischoff, J. A. Calvin, G. I. Fann, J. Fosso-Tande, D. Galindo, J. R. Hammond, R. Hartman-Baker, J. C. Hill, J. Jia, J. S. Kottmann, M.-J. Y. Ou, J. Pei, L. E. Ratcliff, M. G. Reuter, A. C. Richie-Halford, N. A. Romero, H. Sekino, W. A. Shelton, B. E. Sundahl, W. S. Thornton, E. F. Valeev, Á. Vázquez-Mayagoitia, N. Vence, T. Yanai and Y. Yokoi, SIAM J. Sci. Comput., 2016, 38, S123–S142.
23. J. M. Pi, M. Stella, N. K. Fernando, A. Y. Lam, A. Regoutz and L. E. Ratcliff, J. Phys. Chem. Lett., 2020, 11, 2256–2262.
24. A. Regoutz, M. S. Wolinska, N. K. Fernando and L. E. Ratcliff, Electron. Struct., 2020, 2, 044005.
25. G. Wahba, SIAM Rev., 1965, 7, 409.
26. W. Kabsch, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1978, 34, 827–828.
27. F. L. Markley, J. Astronaut. Sci., 1988, 36, 245–258.
28. L. E. Ratcliff, L. Genovese, S. Mohr and T. Deutsch, J. Chem. Phys., 2015, 142, 234105.
29. L. E. Ratcliff and L. Genovese, J. Phys.: Condens. Matter, 2019, 31, 285901.
30. L. E. Ratcliff, W. Dawson, G. Fisicaro, D. Caliste, S. Mohr, A. Degomme, B. Videau, V. Cristiglio, M. Stella, M. D'Alessandro, S. Goedecker, T. Nakajima, T. Deutsch and L. Genovese, J. Chem. Phys., 2020, 152, 194110.
31. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt and SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261–272.
32. I. Daubechies, Ten Lectures on Wavelets, SIAM, 1992.
33. L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Willand, D. Caliste, O. Zilberberg, M. Rayson, A. Bergman and R. Schneider, J. Chem. Phys., 2008, 129, 014109.
34. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
35. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
36. A. Willand, Y. O. Kvashnin, L. Genovese, A. Vázquez-Mayagoitia, A. K. Deb, A. Sadeghi, T. Deutsch and S. Goedecker, J. Chem. Phys., 2013, 138, 104109.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
39. J. H. Scofield, Theoretical Photoionization Cross Sections from 1 to 1500 keV, U.s. atomic energy commission, divison of technical information extension technical report, 1973.
40. C. Kalha, N. K. Fernando and A. Regoutz, Digitisation of Scofield Photoionisation Cross Section Tabulated Data, 2020.
41. A. J. Jackson, A. M. Ganose, A. Regoutz, R. G. Egdell and D. O. Scanlon, J. Open Source Softw., 2018, 3, 773.
42. L. E. Ratcliff, W. S. Thornton, A. Vázquez Mayagoitia and N. A. Romero, J. Phys. Chem. A, 2019, 123, 4465–4474.
43. L. Lin, Z. Wang, J. Fan and C. Wang, Org. Electron., 2017, 41, 17–25.
44. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
45. G. H. Kim, R. Lampande, J. B. Im, J. M. Lee, J. Y. Lee and J. H. Kwon, Mater. Horiz., 2017, 4, 619–624.
46. Y. Zhang, D. Zhang, T. Tsuboi, Y. Qiu and L. Duan, Sci. China: Chem., 2019, 62, 393–402.
47. N. Pueyo Bellafont, F. Illas and P. S. Bagus, Phys. Chem. Chem. Phys., 2015, 17, 4015–4019.
48. N. Pueyo Bellafont, F. Viñes and F. Illas, J. Chem. Theory Comput., 2016, 12, 324–331.
49. X. Ge, I. Timrov, S. Binnie, A. Biancardi, A. Calzolari and S. Baroni, J. Phys. Chem. A, 2015, 119, 3816–3822.
50. T. J. Zuehlsdorff, P. D. Haynes, M. C. Payne and N. D. M. Hine, J. Chem. Phys., 2017, 146, 124504.
51. J. C. A. Prentice and A. A. Mostofi, J. Chem. Theory Comput., 2021, 17, 5214–5224.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp02638d

---