Molecular library of OLED host materials—Evaluating the multiscale simulation workflow

In recent years, display and lighting technologies based on organic light emitting diodes (OLEDs) have steadily increased in popularity, from the automotive industry to smart phones and televisions. Not only do they offer improved image quality via a greater contrast ratio when compared with the liquid crystal display (LCD) technology, but the flexible and transparent possibilities of OLEDs make them an innovative choice in a competitive market. OLEDs usually comprise multiple layers of small organic molecules, sandwiched between two electrodes. Each layer is tailored for a specific functionality to facilitate balanced charge carrier injection and transport to the emissive layer, where the charge carriers recombine to form excitons and consequently photons, achieving a desirable wavelength of emission.^1–3 

One of the major advantages of using organic molecules is the possibility of manipulating chemical composition to target specific properties, such as emission color and charge carrier transport abilities. The aim of OLED material design is to achieve an optimal balance between stability, efficiency, operational driving voltage, and color coordinate of the device. However, molecular design often relies on chemical intuition and with a vast number of potential candidates for each of the layers, this approach is inherently flawed. Therefore, in order to optimize design a systematic approach which links chemical structure to macroscopic properties would be greatly beneficial.^4,5 Ideally, this would focus solely on in silico prescreening, prior to synthesis. Predictive computational modeling is, however, not yet sufficiently accurate for this task exclusively⁶ and with experimental prescreening being an extensive and costly procedure, a collective experimental and in silico approach presents a favorable alternative.

The in silico evaluation of organic materials involves accurate prediction of device characteristics from the corresponding molecular building blocks, which requires simulations over a broad range of length and timescales. The multiscale simulation techniques, as depicted in Fig. 1, help us to establish links between microscopic and macroscopic material properties. Here, quantum chemical calculations are used to obtain gas-phase geometries of a neutral molecule, its cation and anion, as well as ionization energies and electron affinities. The material morphology is then simulated by molecular dynamics (MD), with classical force fields describing atomistic interactions. Polarizable force fields are then used to account for electrostatic effects upon charge/exciton transfer. Kinetic Monte Carlo (KMC) is then used to perform charge transport simulations, whereby it is possible to extract macroscopic observables, such as charge carrier mobilities, by solving the master equation.^5,7–13

It is clear that the morphology is an integral part of the simulations, and it is often a challenge to generate morphologies representative of experimental systems. For this task, both the classical force fields and the more computationally demanding polarizable force fields have to be parameterized. Experimental input at this stage ensures the accuracy of the structural predictions using the classical force fields, while polarizable force fields can be parametrized from the first principles. The parametrization of these force fields is a tedious task and impossible for the vast number of organic compounds required for prescreening.^14–16 To overcome this, it is possible to use the similarities among the molecules most likely to be experimentally investigated, in order to create molecular fragments or building blocks, including the force field parameters. The concept of an extendable molecular library containing these well-defined building blocks required to generate realistic morphologies would then permit the swift characterization of new systems.

For this concept to be brought to a realization, a well-defined simulation workflow capable of predicting relevant system properties from the chemical structure, must first be approved. Therefore, the accuracy and reliability of the various simulation techniques of this workflow, is the subject of the present work. To establish a starting point, simulation results for 12 small organic molecules are summarized, the molecular structures of which are shown in Fig. 2. The individual steps of the simulation workflow are scrutinized and directly compared to experimental results for the glass transition temperature, energetic disorder, ionization energy in the solid-state, and charge carrier mobilities. By doing so, the outlined simulation workflow and the force fields used can be validated, allowing for the expansion of the library and further structures to be investigated.

While this review focuses on one multiscale simulation procedure, it is worth emphasizing that there are several different computational approaches for the consideration of organic semiconductors. One such method highlights the importance of using a fully quantum mechanical approach for charge dynamics to improve on semiclassical Marcus rates.¹⁷ Quantum chemical methods can also be used for evaluating ground and excited state properties of organic electronic materials. However, standard DFT functionals can lead to errors surrounding localization and delocalization of electron and hole densities, prompting the use of tuned long-range corrected hybrid functionals.¹⁸ Excited state properties have also been evaluated using many-body Green's function theory and GW approximation with the Bethe–Salpeter equation (GW-BSE).^19–24 In terms of charge transport, KMC has been used to study degradation and the sensitivity of OLED device lifetime and efficiency, in relation to material-specific parameters, to aid with design strategies with regard to energy levels.²⁵ Additionally, KMC models with molecular-level resolution have been developed for investigating organic field-effect transistors (OFETs).²⁶ Computational methods have also been used to investigate chiral composition–dependent charge mobilities, helping to link variations of molecular chirality to charge transport properties.²⁷ Multiscale simulations have also aided with the study of photoluminescence quenching mechanisms in phosphorescent OLEDs.²⁸ Therefore, it is clear that predictive multiscale protocols are essential in the search and the design of novel materials for organic electronic devices.⁴ 

In addition to electronic properties, the thermal stability of an organic material is essential in determining its suitability to be used in OLED devices. To enhance the device's lifetime, materials which are less susceptible to thermal degradation are targeted, particularly making use of materials with a high glass transition temperature, $Tg$⁠. When the device is operational, the flow of current results in Ohmic heating within the organic layers and due to this local heating, materials which have a low $Tg$ experience a higher degree of molecular vibrations²⁹ and chemical decomposition, resulting in lower stability. Furthermore, high $Tg$ also prevents morphology changes during post-processing of devices, like unwanted partial re-crystallization. Therefore, accurate predictions of the $Tg$ value for OLED materials can serve as a method of screening thermally stable candidates. The $Tg$ for the 12 organic materials was obtained by MD simulations, as outlined in the Methods section, the values are listed in Table I, and the comparison to experimental values is shown in Fig. 3. The experimental values are listed as referenced values from previous studies, or as new experimental values obtained by differential scanning calorimetry (DSC), as outlined in the Methods section.

The simulation results show a systematic overestimation of T_g for all systems, in comparison to the experimental values, except for TMBT. This overestimation is expected as the process of obtaining the simulated morphologies, outlined in the Methods section, involves thermal annealing above T_g, followed by a fast-cooling process. The simulated morphologies often disagree with experiment, as the cooling rates are much faster in simulation compared to experiment and reproducing the realistic molecular packing requires long simulation times. Additionally, the relaxation times of typical OLED depositions cannot be achieved by the slow molecular dynamics close to T_g. Furthermore, the deposition conditions can have a significant impact on molecular orientation,^39,40 such that experimental T_g values may also show variation, particularly for systems with low T_g, for example TMBT or also BCP, where T_g was not detectable (as there was no significant kink in calorimetric measurements, described in the Methods section) despite a previously reported value in the literature.³⁰ For some of the compounds, such as TMBT, it is known that vacuum deposition leads to molecular alignment, i.e., molecular orientations are not random. This effect is not accounted for in our outlined protocol, as it is difficult to reproduce in simulations; large timescales are necessary to cover the diffusion process of evaporated materials, mimicking experimental film-growth rates of 1 Å/s.^39,41,42

In order to accurately predict the glass transition temperature from the morphology, it is clear that slower deposition rates must be employed, for example with simulation methods such as coarse graining.^39,43 An alternative to predict T_g for OLED materials is a quantitative structure–property relationship approach,⁴⁴ with the use of topological indices⁴⁵ and various descriptors. However, this only allows to interpolate within known chemical space. Despite the inaccurate prediction of the T_g values, the proceeding simulation results will show that this does not have a significant impact on the charge transport simulations and mobilities for the materials in this study.

In disordered organic materials, due to weak intermolecular coupling⁴⁶ and relatively large energetic disorder, the charges are localized and propagate by a successive hopping process from one molecule in the system to the next. The hopping rate between two sites in a material is characterized as a thermally activated transport process, where the rate can be expressed in terms of the Marcus rate equation,^47–49 

Here $T$ is the temperature, $kB$ is Boltzmann's constant, $λij$ is the reorganization energy, $Jij$ is the electronic coupling matrix element, and $ΔEij=Ei−Ej$ is the driving force or site energy difference between two neighboring sites, where $Ei$ is the site energy of molecule i.⁵⁰ The site energy includes the internal molecular energy (⁠$Eiint$⁠) and interaction with the environment, including the electrostatic (⁠$Eielec$⁠) and induction (⁠$Eiind$⁠) contributions. It is calculated as $Ei=Eiint+Eielec+ Eiind+qF · ri$⁠, where q is the hopping carrier charge, F is an applied external field, and $ri$ is the center-of-mass of molecule i.⁵ For computational efficiency, here we use the Marcus rate expression; the entire scheme can be also adapted to the quantum treatment of vibrational modes.¹⁷ The microscopic quantities, such as reorganization energy, electronic coupling, etc., are computed from first principles and serve as an input to the master equation, the solution of which gives charge carrier mobilities.^51,52

The methodology of obtaining the reorganization energies and electronic coupling elements is described in the Methods section. The reorganization energies for each material are listed in Table S1 and electronic coupling elements are shown in Fig. S1, of the supplementary material. Even though the variations of the reorganization energies and electronic coupling elements lead to variations in the simulated mobility, μ, the most significant parameter is the distribution of site energies $Ei$ within the system, characterized by the energetic disorder σ. To a certain extent, this is anticipated, as the mobility is exceptionally sensitive to changes in the width of the disorder distribution, for example, $μ ∝exp[−C (σkBT)2]$⁠.^53–55 The energetic disorder stems from the disorder on the local electronic states which, as we will see in Fig. 4, is Gaussian-distributed.

To evaluate the site energies of holes and electrons in the 12 systems, a perturbative scheme was used, as outlined in the Methods section. The distributions of the on-site energy differences, i.e., the differences between the energies of the system when a selected molecule is in the anionic/cationic or neutral state, including the constant internal contribution due to the gas-phase electron affinity/ionization potential, are displayed in Figs. 4(a) and 4(b) for electrons and holes, respectively.

The corresponding energetic disorder (widths of site energy distributions) are summarized in Table II. In amorphous organic materials, the energetic disorder is predominantly electrostatic; such electrostatic interaction originates from the potential exerted on a molecule from its specific environment. Therefore, the disorder is governed by the molecular static multipoles, as well as the positional and conformational order in a given material. On the other hand, the induction contribution stemming from the interaction of microscopic dipoles (the distributions of the electrostatic potential in the ground state are shown in Fig. S2 of the supplementary material) with the localized charge carrier, reduces the energetic disorder.^13,56 The electrostatic and induction contributions for each system, for both electrons and holes, are listed in Table II.

Apart from energy distributions, the long-range electrostatic interactions of a charge with molecular dipoles can lead to spatial correlations of site energies.^55,57 To unveil such spatial correlations, we computed these correlations for holes and electrons, which are shown in Fig. S3 of the supplementary material. We found that the hole and electron site energy correlations extend up to 2–3 nm. The correlations exhibit a decay profile with distance, approximately following the r⁻¹ signature expected for the random dipolar disorder.⁵⁵ Therefore, all of these materials possess a correlated disorder even though their ground state dipole moments are small for a large set of molecules. We now analyze in more detail ionization energies, electron affinities, and the widths of the DOSs, which we consider to be the most important parameters for charge transport.

The gas-phase electron affinities (EA₀) and ionization energies (IE₀) were calculated using density functional theory at the M062X/6–311g(d,p) as well as IP/EA-EOM-DLPNO-CCSD/aug-cc-pVTZ level of theories, as described in the Methods section and in the supplementary material. The comparison of EA₀ and IE₀ obtained at various levels of theory is shown in Figs. S4–S6 (see also Table S2) of the supplementary material. In inverse photoemission spectroscopy (IPS) and ultraviolet photoelectron spectroscopy (UPS), the electron affinity and ionization energy are measured as an onset of the spectra. In simulations, to account for the finite width of the density of states (σ), the solid-state electron affinities (EA_tot) and ionization energies (IE_tot) were shifted by 2σ, such that EA = α + 2σ or IE = α − 2σ, where α represents the mean of the DOS for electrons or holes, respectively. The EA_tot and IE_tot values are summarized in Table II for each material. Experimental ionization energies are also listed in the table and shown in Fig. 4(b), for comparison. Further to this, the comparison of solid-state IE values to the corresponding experimental values are shown in Fig. 5. The method of obtaining IE values by photoelectron yield spectroscopy in air (PESA) is described in the Methods section, the UPS values are taken from the literature, with references included in Table II. Both gas-phase and solid-state ionization energies are compared to PESA and UPS values in the supplementary material, shown in Fig. S7. Additionally, the gas-phase and solid-state ionization energies and electron affinities are compared to experimental cyclic voltammetry (CV) data, shown in Fig. S8 of the supplementary material, where the experimental method is also described.

The IE values obtained from the DOS of the various materials are in good agreement with the experimental IE values, as depicted in Fig. 4(b) for the individual systems and Fig. 5 as the total correlation. The computed IE_tot values are a combination of the gas-phase IE₀ as well as the electrostatic and induction contributions. This is necessary as the gas-phase IE is for a single isolated molecule and does not account for the solid-state effects required to accurately determine electronic properties. Interestingly, as shown in Fig. S7 of the supplementary material, the linear correlation to experimental solid-state values is already reproduced from the gas-phase simulations, however, with a proportionality coefficient different from unity. The agreement of the computed (IE_tot) and experimental IE values signifies a reasonable degree of accuracy of the polarizable force fields.

The method deployed in the present work also enables the distinction of different contributions to the EA and IE, e.g., gas-phase, electrostatic, and induction, shown in Table II. For all compounds, the stabilization due to induction (E_ind) contributes around 0.5–0.8 eV and 0.6–1.4 eV for holes and electrons, respectively. On the other hand, the electrostatic contribution (E_elec) varies from one system to another (0.01–0.4 eV for both holes and electrons) governed by the distribution of molecular dipoles in a given system. A closer scan of the numbers in electrostatic contribution (E_elec in Table II) reflects a salient correlation between energetic disorder and electrostatic interactions for molecules such as NBPhen, TMBT, and BCP possessing large charge-dipole interaction. Here, the magnitudes of E_elec are almost two times higher than the rest of the materials, as well as the higher σ values, shown for both holes and electrons.

The energetic disorder for hole transport for each material is taken as the width of the DOS in simulations. These values are directly compared in Fig. 6 to experimental values from previously reported space charge–limited current (SCLC) measurements,⁶¹ and/or newly carried out thermally stimulated luminescence (TSL) measurements. Details of the SCLC measurements can be found in Ref. 61, and the TSL measurements are described in the Methods section. The simulations predict a significant variation of the energetic disorder among these materials which features an almost twofold change spanning from $σ=0.09 eV$ observed for a weakly disordered NPB, to $σ=0.19 eV$ obtained for highly disordered NBPhen. The experimental results demonstrate a remarkably similar trend (Fig. 6), which testifies to the intrinsic nature of the large width of the DOS (σ > 0.1 eV) observed experimentally in most of these OLED materials, implying the DOS of a chemically pure disordered material rather than being significantly affected by impurity-related traps. The energetic disorder inferred from TSL normally provides the DOS of shallower trapped carriers, which are thermally released and recombine with the more deeply trapped counter charges. This is because TSL is an electrode-free technique, so the number of electrons and holes in the film is always the same to maintain the material neutrality at any time after the carriers were photogenerated. Once a shallower trapped carrier is released from its trap, it will recombine with its deeper trapped counterpart. Therefore, the latter sort of carriers is expected to already be completely annihilated, and there are simply no carriers left at a temperature relevant to their anticipated release from their deeper states.

The ionization energies of most of the materials measured in this study are less than, or around 6 eV (Table II). Therefore, their hole DOS should not be affected much by extrinsic traps, as the IE values fall within an energy window, inside which organic semiconductors normally do not experience charge trapping.⁶⁷ Materials with an ionization energy above 6 eV will exhibit trap-limited hole transport, similarly, an electron affinity lower than 3.6 eV will cause electron trapping to limit electron transport.⁶⁷ As the materials in this study have electron affinities significantly less than 3.6 eV (Table II) shown by simulations, the electron transport in these systems is most likely hindered by electron trapping (as oxygen-related traps). Due to this reasoning and taking into account that all samples were exposed to air prior to TSL measurements, which makes relatively deep electron traps unavoidable, our TSL measurements probe the hole DOS in these materials. Although the simulated σ-parameters for holes and electrons (Table II) are not appreciably different, it should be mentioned that TSL measures the “effective DOS” and even a small concentration of shallow extrinsic traps (depending on the trap depth it can even be at trap concentration level of ∼0.1%) can give rise to a notable DOS broadening. This has been previously demonstrated by deliberate doping of photoconductive polymers with traps.⁶⁸ Since the IE of CBP, mCBP, and TPBi slightly exceeds the trap-free window threshold value of 6 eV, a certain extrinsic hole trapping cannot be excluded here. This can explain why TSL measurements for these materials yield slightly overestimated energetic disorder parameters, compared to the simulated values (Fig. 6).

Finally, we comment on some differences in $σ$-parameters determined by SCLC and TSL techniques observed for some materials, e.g., for CBP. This may be partially due to different film morphologies related to different film deposition techniques used: while films for SCLC measurements were vacuum-deposited, most of the materials for TSL measurements were spin-coated from a solution, except NPB, TCTA, and Spiro-TAD, which were vacuum-deposited. Another likely reason might be due to the fact that the SCLC probes the DOS in a large-carrier-concentration transport regime, for which deep tail states of the DOS (and/or traps) might be filled with carriers during the measurements and therefore, $JV$-characteristics might be governed mostly by a shallower portion of the DOS. On the other hand, carrier concentration in TSL experiments is significantly smaller (comparable to time-of-flight experiments) and deep tail states can be well probed by this technique.

One can also attempt to relate the observed energetic disorder to molecular dipoles, since lattice models with randomly oriented dipoles of equal magnitude $d$ and lattice spacing $a$⁠, yield disorder $σ∼d/a$⁠. However, in a realistic morphology, both distances and dipoles can vary from molecule to molecule, which then increases the disorder further. It is evident from the distributions in Fig. S9 that the dipole moments of individual molecules in an amorphous morphology vary significantly from their equilibrium values (see Table S3 of the supplementary material). Such broad distributions of dipole moments are attributed to the presence of one or more soft dihedrals in the organic materials. At finite temperature, rotation around such soft degrees of freedom leads to multiple conformers. In general, systems with a narrow distribution of molecular dipoles in their amorphous morphologies possess small fluctuations of electrostatic multipoles in the amorphous matrix which results in smaller energetic disorders. This conclusion is, however, not universal since large local electrostatic potential variations of TMBT, for example, can result in a sizable disorder despite a zero molecular dipole moment. The correlation of average dipole moment and energetic disorder for the 12 systems are shown in Fig. S10 of the supplementary material, for both holes and electrons. Although there is a slightly better correlation for holes compared to electrons, a number of outliers remain for both, implying that it is difficult to predict the trend in sigma, solely from the molecular dipole moment.

Charge transport rates were computed using the high temperature limit of classical charge transport theory^47–49 as given by the Marcus rate equation. The master equation can then be solved with kinetic Monte Carlo (KMC), providing the time evolution of the system, giving a randomly generated trajectory of charge carrier movement. This was carried out for one charge carrier (hole or electron) in the presence of an applied electric field (⁠$F=1×104$ V/cm), using a periodic simulation box. Mobilities were extracted as outlined in the Methods section. The extrapolated dependence of mobility on temperature is shown in Fig. S11 and Fig. S12 of the supplementary material, for holes and electrons, respectively. The single-carrier mobilities at room temperature for holes and electrons are summarized in Table III. The experimentally measured mobility and the corresponding experimental techniques used, are also listed for comparison.

Table III and Fig. 7 show the correlation between simulated and experiment mobility. The remarkable agreement of simulation and experiment is particularly evident for hole mobilities. On the other hand, for electron mobilities, a larger deviation is observed between experiment and simulation, where simulated results indicate a systematic underestimation of experimental measurements. There are several possible explanations to account for these discrepancies. First, due to the much larger energetic disorder for electrons in certain materials (Table II) when compared to holes, the electron mobility will be inherently lower. As a result of the large disorder found in 2-TNATA, the simulated and experimental mobilities show large variation. This may stem from energetic traps in the simulated morphology, which can lead to lower mobility values. To highlight the significant role of energetic disorder on the simulated mobility, a correlation plot is shown in Fig. 8, for both holes and electrons.

Due to different morphologies, the structural and energetic disorder can differ significantly between simulation and experiment. Despite the reasonable agreement for energetic disorder for hole transport, shown in Fig. 6, the experimental systems used for the mobility correlation are a collection of referenced values from various studies, with potentially significant variations in disorder.

It should also be noted that the simulated morphologies for the 12 systems do not account for the presence of carrier traps formed by structural defects or impurities such as water, which are typically unavoidable in reality. However, the inclusion of carrier traps in the simulated morphology would, in fact, lead to larger deviation between simulated and experimental mobilities. As previously stated, hole or electron transport have been shown to become trap-limited in materials with an IE greater than 6 eV or an EA less than 3.6 eV, respectively.⁶⁷ Therefore, direct comparison of simulation and experiment mobilities may be difficult when considering low EA and high IE materials.

Finally, the takeaway message here is that the width of the density of states is the key property in determining the mobility (Fig. 8). Hence, accurate predictions of the DOS should be given priority when prescreening OLED hosts. DOS width correlates with some extent with the molecular dipole, but this correlation has too many outliers, e.g., due to conformational freedom or higher order multipoles and therefore, cannot be used as a reliable descriptor for mobility predictions.

It is clear that a molecular library of OLED hosts would be invaluable, permitting the swift evaluation of new materials. The key question is how accurately and reliable a combination of various simulation techniques can predict relevant material properties, to bring prescreening a step closer.

For the simulated morphological properties, the T_g comparison showed large variation to experimental values, suggesting that the current approach needs to be improved for accurate evaluation of thermally stable materials. Reliable methods for predicting T_g and accurate atomistic force fields are the key improvements required for more accurate simulation results.

The morphologies used for charge transport simulations revealed that the inaccuracies in T_g predictions had no significant impact on simulated ionization energies. In this respect, PESA measurements and UPS data taken from the literature, as well as cyclic voltammetry measurements, are in an excellent agreement with simulation results, supporting our confidence in the polarizable force fields used for evaluation of the solid-state electrostatic contribution. The situation is somewhat different for electron affinities, where using different computational techniques led to large variation of the gas-phase electron affinity values, even at a computationally affordable variation of the coupled cluster implementation. Moreover, there is no clear benchmark possible for the solid-state because of the sparse availability of the inverse photoemission spectroscopy measurements.

Accurate solid-state energetics allowed us to predict the density of states which, when compared to the thermally stimulated luminescence measurements, showed a similar trend. The energetic disorder can be potentially correlated with the distribution of molecular dipoles, but the extent of this will require further investigation to be conclusive.

Finally, the simulated charge carrier mobility showed a remarkable agreement with experimental values, particularly for hole transport where energetic disorder is typically lower. The accurate prediction of energetic disorder is therefore vital, as it has a significant impact on mobility.

Overall, the correlation of simulation and experimental results has been used to validate the accuracy of the force fields and the simulation methods, as an initial step toward building a larger molecular library. The agreement of simulation and experiment for the various parameters, particularly mobility, highlights the predictive capability of the outlined methods and the simulation workflow. The next step is to expand this library with further materials, in an effort to draw structure-property conclusions for effective prescreening.

For the isolated molecules, density functional theory (DFT)–based electronic structure methods were used to compute gas-phase ionization energy (IE₀) and electron affinity (EA₀) using the Gaussian09 program.⁷⁸ For this, the neutral molecule in the neutral geometry (E_nN), as well as the charged molecule in the charged geometry (E_cC), are computed. The IE₀ and EA₀ values are then calculated as E_cC – E_nN, where E_cC represents the cationic and anionic state, respectively.

The prediction of accurate IEs and EAs remain a challenge in electronic structure theory, primarily due to the self-interaction error⁷⁹ (SIE) or localization/delocalization⁸⁰ error, inherent to commonly used DFT functionals. In a series of recent reports,^18,81–84 the importance of considering long-range corrected hybrid functionals has been demonstrated, such that the SIEs in DFT description of molecules can be reduced. Therefore, adiabatic IE₀ and EA₀ calculations were performed by employing a range of DFT functionals: PBEPBE, B3LYP, CAM-B3LYP, ωB97X-D, M062X, and LC-ωPBE, combined with the basis set 6–311+g(d,p). The IEs and EAs are compared for each of the 12 organic molecules and each DFT functional, as shown in Fig. S4 of the supplementary material. Different levels of theory are also compared in Fig. S5 and Fig. S6.

It is evident that calculations performed using PBEPBE and B3LYP underestimate IE₀ compared to other functionals, an observation which is attributed to an underestimation of Kohn–Sham eigenvalue by hybrid functionals, leading to significant overscreening of the Coulomb interaction. Overall, the IE₀ is better predicted than the EA₀ for the given molecules. This is clear when comparing the Cam-B3LYP, ωB97X-D, and M062X functionals, as there is greater variation in the EA₀ prediction. The small deviation among these functions, with regard to the IE₀, makes any of the three a suitable choice (M062X is the chosen functional for comparison to experimental IE and EA values in this study).

DFT methods were then used to accurately parameterize the empirical OPLS-AA force field,^14–16 for the 12 chemically diverse molecules. All Lennard–Jones parameters were taken from this force field in combination with the fudge-factor of 0.5 for 1–4 interactions. Atomic partial charges were computed using the ChelpG⁸⁵ scheme for electrostatic potential fitting as implemented in Gaussian09,⁷⁸ employing the ground state electrostatic potential determined at the B3LYP/6–311+g(d,p) level of theory.

In order to generate amorphous morphologies, molecular dynamics (MD) simulations were carried out using the GROMACS simulation package.^86,87 The amorphous state was generated by an annealing step, followed by a rapid quenching to lock the molecules in a local energy minimum. This procedure has been previously applied for the preparation of amorphous structures of OLED materials.^61,88 The starting configurations used in the MD simulations were prepared by randomly arranging 3000 molecules in a simulation box using the Packmol program.⁸⁹ These initial structures were energy-minimized using the steepest-descent method and annealed from 300–800 K, followed by fast quenching to 300 K. Further equilibration for 2 ns and 1 ns production runs were performed at 300 K. All simulations were performed in the NPT ensemble using a canonical velocity rescaling thermostat,⁹⁰ a Berendsen barostat for pressure coupling,⁹¹ and the smooth particle mesh Ewald technique for long-range electrostatic interactions. A time step of 0.005 ps was used to integrate the equations of motion. Non-bonded interactions were computed with a real-space cutoff of 1.3 nm.

To obtain the glass transition temperature of each material, a similar procedure was carried out with a time step of 0.001ps. The material was first equilibrated at 800 K, followed by gradual cooling to 0 K at a rate of 0.1 K/ps. The change in density as the material was cooled was used, in combination with two linear fits (the linear fitting procedure is described in the supplementary material) to extract the intersection point, giving the T_g value.

We used MD simulation trajectories to evaluate the site energies of holes and electrons by employing a perturbative scheme. In this approach, the electrostatic and induction energies are added to the gas-phase energies (IE₀ or EA₀), to obtain the total site energy. The electrostatic contribution is calculated with the use of Coulomb sums based on distributed multipoles (obtained from the GDMA program⁹²) for neutral and charged molecules in their respective ground states. The polarization contribution is computed using a polarizable force field based on the Thole model^93,94 with isotropic atomic polarizabilities (α_ai) on atoms a in molecules i. Aperiodic embedding of a charge method⁹⁵ as implemented in the VOTCA^13,43 package, was used for these calculations.

The transfer integral or coupling elements, $Jij=ϕi Ĥ ϕj$⁠, represent the strength of the coupling of the two frontier orbitals |$ϕi$⟩ and |$ϕj$⟩ localized on each molecule in the charge transfer complex. It is highly sensitive to the characteristic features of the frontier orbitals as well as the mutual orientations of the two molecules and follows an exponential decay with distance. The electronic coupling elements shown Fig. S1 of the supplementary material, were computed for each neighboring molecular pair (ij) using a projection method.^96,97 Molecular pairs were added to the neighbor list, with a center-of-mass distance cutoff (between rigid fragments) of 0.7 nm. These calculations were performed at PBEPBE/6–311+g(d,p) level of theory using the Gaussian09⁷⁸ and VOTCA^13,43 packages. The frozen core approximation was used with the highest occupied molecular orbitals providing a major contribution to the diabatic states of the dimer.

The reorganization energy of the system takes into account the charging and discharging of a molecule. When a charge moves from molecule i to molecule j, there is an intramolecular contribution (⁠$λijint$⁠), due to the internal reorganization of the two molecules and an intermolecular, known as an outersphere contribution (⁠$λijout$⁠), due to the relaxation of the surrounding environment.^13,50 The internal reorganization energy is calculated as $λijint=(UinC−UinN)+(UjcN−UjcC)$⁠, for molecule i and j, where the lowercase represents the neutral (n) or charged (c) molecule and the uppercase represents the neutral (N) or charged (C) geometry. The reorganization energies are summarized for holes and electrons in Table S1 of the supplementary material, for each of the systems. The individual contributions were calculated by DFT using the B3LYP/6–311+g(d,p) level of theory. In the current work, we ignore the outersphere contribution to the reorganization energy, since in the amorphous solids considered here, the Pekar factor is on the order of 0.01,¹³ leading to a relatively small contribution to the total reorganization energy.

Mobilities were extracted as μ = ⟨v⟩/F, where ⟨v⟩ is the average projection of the carrier velocity in the direction of the field (⁠$F=1×104$ V/cm). The convergence of simulated mobilities with respect to the system size (i.e., a sufficient number of sites for the simulated transport to be nondispersive) due to energetic disorder must be ensured. For this purpose, the critical temperature (within Gaussian disorder model), $Tc$⁠, at which the transition from dispersive to nondispersive regime takes place, was estimated and is shown in Fig. S11 and Fig. S12 of the supplementary material. The mobilities obtained in the nondispersive regimes are then fitted with an empirical temperature dependence, allowing for extrapolation of the nondispersive charge carrier mobility at room temperature. Details of such calculations can be found elsewhere in Refs. 98 and 99. All charge transport calculations were performed using the VOTCA package.^13,43

For determination of the glass transition temperature (T_g) at Merck KGaA, Darmstadt, Germany, we used differential scanning calorimetry (DSC) analyzing samples of 10–15 mg in DSC 204/1/G Phönix from Netsch. Samples were heated by 5 °C/min up to 370 °C then cooled by 20 °C/min to 0 °C and finally heated again by 20 °C/min to 370 °C, where T_g was determined by the kink in heat flow vs temperature using the temperature corresponding to half the drop in heat flux. Only for BCP and TMBT, this protocol did not yield a significant kink. TMBT was expected to be the lowest T_g material from simulation, so we tried other protocols to measure T_g. We finally used a 5 mg TMBT sample in DSC Discovery from TA Instruments in nitrogen atmosphere and first heated by 20 °C/min up to 320 °C, then for cooling, quenched the sample by liquid nitrogen and finally heated by 20 °C/min up to 320 °C, where the Tg was observed. Other protocols we tried for TMBT without the cooling quench did not lead to observation of T_g.

Photoelectron yield spectroscopy in air (PESA) was performed on 50 nm thick thermally evaporated films at Merck KGaA, Darmstadt, Germany, using the surface analyzer AC-3 from RIKEN KEIKI Co., Ltd. The film is exposed to monochromatic light from a deuterium-lamp, with incident photon energy between 4 and 7 eV (increased in steps of 0.05 eV),¹⁰⁰ while an open counter analyzes the photoelectron yield.¹⁰¹ The square root of this photoelectron yield is plotted vs energy of the incident photons. The underground is taken as the average horizontal line through the measurement points for low photon energies. The ionization potential is calculated as intersection of the underground and a fit to the onset of the square root of the photoelectron yield.

Due to the energy step size and this fitting procedure, which is done manually, an error bar of roughly $±0.06$ eV must be associated with the obtained ionization energies as we have verified by repeated measurement and evaluation. As compared to UPS measurements of films, PESA has the advantage that the penetration depth of photons of 7 eV is roughly 10 nm, while for UPS, only the top 0.5 nm of the film can contribute, so PESA is measuring bulk properties of the film.

TSL is the phenomenon of luminescent emission after removal of excitation (UV light in our case) under conditions of increasing temperature. Our TSL measurements were carried out over a wide temperature range from 5 to 330 K using an optical temperature-regulating helium cryostat. When exciting a sample optically with 313 nm light for typically 3 min at 4.2 K, the charge carriers are generated and populate trapping states. Once the sample is heated up, trapped charge carriers are released and then recombine, producing a luminescence emission. TSL measurements were performed in two different regimes: upon uniform heating at the constant heating rate 0.15 K/s and under the so-called fractional heating regime.¹⁰² In the latter regime, we apply a temperature-cycling program in which a large number of small temperature oscillations are superimposed on a constant heating ramp that allows determining the trap depth with high accuracy. This applies when different groups of traps are not well separated in energy or are continuously distributed, which is of special relevance for disordered organic solids where the intrinsic tail states can act as traps at low temperatures. The details of the TSL measurements have been described elsewhere in Refs. 103–105.

Technically, TSL is a relatively simple technique with a straightforward data analysis. However, it should be noted that the mechanism of TSL in amorphous organic semiconductors with a broad energy distribution of strongly localized states, differs from the mechanism commonly accepted for crystalline materials with band-type transport where a discrete trapping level model is applicable. A specific feature of amorphous solids is that the localized states within the lower energy part of the intrinsic DOS distribution can give rise to shallow charge trapping at very low temperatures, and as a consequence, TSL can be observed even in materials where the “trap-free limit” has been postulated. Since TSL measurements are normally performed after a long dwell time after photoexcitation, the initial energy distribution of localized carriers is formed after low-temperature energy relaxation of photogenerated carriers within a Gaussian distribution of DOS. Theoretical background for application of TSL for probing the DOS distribution in disordered organic systems has been developed,^105,106 using a variable-range hopping formalism and the concept of thermally stimulated carrier random walk within a positionally and energetically random system of hopping sites. The theory proves that the high temperature wing of the TSL curve in such amorphous materials should be an exact replica of the deeper portion of the DOS distribution^105,106 and yields the effective DOS width. Kadashchuk and co-workers have applied low-temperature fractional TSL to investigate the intrinsic energetic disorder in a variety of important semiconducting polymers, oligomers, and hybrid metalorganic perovskites (see, for instance Refs. 103–109). A clear advantage of TSL is that it is a purely optical and electrode-free technique. This helps to eliminate interface/contact effects and, most importantly, it allows DOS measurements even in materials with large energy disorder where the charge transport is very dispersive.

All simulations were performed with the open-source software package VOTCA.¹¹⁰ Atomistic and polarizable force fields, VOTCA input files, and analysis scripts are available in the Materials Library git repository, version 1.0.¹¹¹ 

See the supplementary material for additional material properties, including: reorganisation energies, transfer integrals, comparison of gas phase ionisation energies, and dipole moments. Experimental TSL curves are also shown.

A.M. and L.P. contributed equally to the work.

D.A. acknowledges the BMBF Grant InterPhase (Grant No. FKZ13N13661) and the European Union Horizon 2020 Research and Innovation Programme ‘Widening Materials Models’ under Grant No. 646259 (MOSTOPHOS). This research has been supported by the King Abdullah University of Science and Technology (KAUST), via the Competitive Research Grants (CRG) Program. D.A. acknowledges KAUST for hosting his sabbatical. The Ukrainian team acknowledges funding through the EU Marie Skłodowska-Curie ITN TADFlife grant (Grant No. 812872). This research was also supported by the European Research Council under the ERC Grant No. 835133 (ULTRA-LUX), VW Foundation, and by the National Academy of Science of Ukraine (Project No. VC/205) and NRFU 2020.01/0144. K.-H.L. acknowledges the financial support from the Swiss NSF Early Postdoc Mobility fellowship (Grant No. P2ELP2_195156).

The data that support the findings of this study are openly available in https://gitlab.mpcdf.mpg.de/materials, version number 1.0.