{'fraction_non_alphanumeric': 0.023982340616581863, 'fraction_numerical': 0.01455402574079617, 'mean_word_length': 5.074772727272728, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analyzing dynamical disorder for charge transport in organic semiconductors via machine learning

Organic semiconductors are indispensable for today's display technologies in form of organic light emitting diodes (OLEDs) and further optoelectronic applications. However, organic materials do not reach the same charge carrier mobility as inorganic semiconductors, limiting the efficiency of devices. To find or even design new organic semiconductors with higher charge carrier mobility, computational approaches, in particular multiscale models, are becoming increasingly important. However, such models are computationally very costly, especially when large systems and long time scales are required, which is the case to compute static and dynamic energy disorder, i.e. dominant factor to determine charge transport. Here we overcome this drawback by integrating machine learning models into multiscale simulations. This allows us to obtain unprecedented insight into relevant microscopic materials properties, in particular static and dynamic disorder contributions for a series of application-relevant molecules. We find that static disorder and thus the distribution of shallow traps is highly asymmetrical for many materials, impacting widely considered Gaussian disorder models. We furthermore analyse characteristic energy level fluctuation times and compare them to typical hopping rates to evaluate the importance of dynamic disorder for charge transport. We hope that our findings will significantly improve the accuracy of computational methods used to predict application relevant materials properties of organic semiconductors, and thus make these methods applicable for virtual materials design.


Introduction
Motivation. Organic semiconductors have been successfully integrated into many commercially relevant electronic devices such as in organic light-emitting diodes (OLEDs), 1 organic photovoltaics (OPV), 2 organic lasers, 3 and organic transistors (OFETs) 4 . Their optoelectronic properties and potential flexibility, role-to-role manufacturing 5 or biocompatibility 6 opens many interesting applications. However, concerning charge carrier transport and mobility, organic semiconductors still fall short behind their inorganic counterparts. Electrochemical doping of organic semiconductors does significantly increase conductivity as well as charge carrier injection but can reduce long-term stability and efficiency due to drift and diffusion of dopants. [7][8][9][10] Although organic crystals achieve high mobilities in the range of 1-100 cm 2 /Vs, 11,12 they are usually not used in OLEDs. An OLED made from a single crystalline layer will likely suffer from charge carrier injection and recombination imbalances. 13 To achieve isotropic light emission and homogeneous layer thickness in device stacks with multiple layers, amorphous or semicrystalline materials are required, having significantly lower charge carrier mobilities. Consequently, it is crucial to understand and to improve charge transport, in particular charge mobility, for common OLED and OPV materials, which is determined primarily by the width of the disorder distribution, as shown by Bässler et al. [14][15][16] The enormous suppression of the mobility with increasing disorder is largely responsible for the poor performance of amorphous organic semiconductors with respect to mobility. Multiscale modeling methods had large success in ab-initio mobility prediction for various materials and therefore serve as a screening tool for new or modified molecules. 7,[17][18][19][20][21][22][23][24] In most models, charge transport is described as a sequence of hopping processes, i.e. , thermally assisted tunneling of localized polarons between neighboring molecules. 25 Hopping rates can be estimated based on the Miller-Abrahams model and Marcus theory, 26,27 which require reorganization energies and transition matrix elements from ab-initio DFT calculations, as well as the spatial distribution of molecular energy levels. The hopping rates can then be inserted into master equations, [28][29][30] Monte Carlo simulations 31 or mean field theories 32 to determine the carrier mobility of interest. 33,34 Presently most quantitative models for the mobility assume a static picture with frozen molecules and thus time-independent energy levels, although molecular morphologies are often created by dynamics simulations. 21,35,36 This is partly due to the numerical effort of electronic structure methods, such as DFT, which precludes their use in large systems with full dynamics. For this reason, the role of dynamic disorder is presently not well understood.
Static and dynamic disorder. In amorphous molecular solids, at any given point in time, the distribution of orbital energies across all molecules in the system is typically estimated to be Gaussian with a material-dependent standard deviation, which is referred to as (global) energy disorder in the Gaussian disorder model. 14 The disorder parameter constitutes a sensitive factor for mobility σ prediction: . However, due to thermal vibrations around the individual (and potentially μ ∝ e −C (σ/kT ) 2 different) equilibrium conformations of each molecule in the amorphous matrix, these energy levels are fluctuating in time. The thermally driven conformational disorder can include methyl rotations, 37 pedal motions, 38 resonances in hydrogen bonds 39 or dipole orientational 40 disorder. 41,42 This fluctuation is referred to as dynamic disorder, [43][44][45] which is related to electron-phonon interactions which cause a spread of energy levels. If all molecules had the same equilibrium conformation, all time-averaged energy levels would be equal. However, in most amorphous materials, the molecules are constrained in individual equilibrium structures, leading to a spread of the time-averaged energy levels, which is referred to as static disorder. 46 It will depend on the size and flexibility of the molecules, on the response of energy levels to conformational changes, on the degree of crystallinity and on possible additives such as dopants. 47 The effective energy landscape and energy disorder experienced by the charge carriers thus depends on the typical residence time of charge carriers, 48,49 i.e. , electrons and holes localized on a single molecule. In case the residence time is much smaller than the fluctuation time of energy levels (high electronic coupling between molecules), the dynamic disorder has to be taken into account, and the total disorder becomes an upper bound of the energy disorder. However, in case of lower hopping rates, each charge carrier will experience a fluctuating energy landscape potentially with multiple resonances with neighbouring energy levels, leaving the static disorder as a lower bound. 48 As a consequence, the thermal energy fluctuations and the resulting dynamic disorder can be beneficial or detrimental to charge transport. This will depend on whether the charge transport is limited by either static disorder or electronic couplings between neighbouring molecules, where the transfer integral is influenced by orientation and distance as well. Additionally, electrostatic interaction will modify the energy landscape, as soon as a charge is moving in the system. Depending on the transport model, the energy distribution integrated over in percolation path theory for mobility prediction or used in Monte-Carlo simulations to initialize site energies, is often assumed to be static.
However, the disorder is commonly estimated by the total disorder, or the morphologies are generated by dynamics simulations only, because of the computational cost of energy calculations. 21,33,35,36,50 Therefore, we incorporated faster machine learning (ML) 51 models to upscale energy predictions on long trajectories with thousands of molecules in order to obtain insights into the dynamics of disorder effects and to distinguish static and dynamic disorder contributions. 43,44,[52][53][54] This work. In a recent proof of principle study, we showed that energy level fluctuations in PEDOT:PSS can be evaluated using machine learning models, enabling the calculation of dynamic and static disorder in PEDOT oligomers. 55 In this paper, we present an extended study of static and dynamic disorder, evaluated using ML models for OLED and OPV application relevant materials (depicted in Scheme 1 ). For this study, we focus on conformational disorder and evaluate vacuum single point energies for isolated molecules from DFT reference calculations to train the ML models.
Electrostatic disorder due to the electrostatic potential of neighbouring molecules as well as additional electrons and holes is neglected here, but plays an important role to quantify charge carrier mobility in amorphous organic semiconductors. 21,34 We show that machine learning models can be trained to reliably predict orbital energies with a limited number of DFT calculations and subsequently used to compute static and dynamic disorder on the full set of MD trajectories. Thereby we obtain good statistics for a comprehensive analysis of energy disorder as a function of time and position. This enables us to compare dynamic energy changes between neighboring molecules with typical hopping rates of charge carriers. In this work, we have succeeded in separating the static disorder from the dynamic disorder, which is often larger than the former. The comparison of hopping and energy level fluctuation rates is required to verify the static disorder being used for multiscale models.
Nonetheless, we discuss that for a certain parameter range the dynamic disorder could become relevant to more accurately predict the charge carrier mobility of amorphous organic semiconductors.
We expect that our method will enable accurate prediction of the charge carrier mobility in amorphous organic semiconductors, paving the way to a reliable virtual design of high-mobility amorphous organic semiconductors.

Methods
The workflow used in this study consists of the parameterization of molecule specific classical force

Molecular dynamics.
For each molecule, the geometry is initially optimized with a conformer search using the submodule CREST of semi-empirical program package XTB [56][57][58] to find the minimum energy conformation. The resulting structure is submitted to the online platform ATB 59,60 , which automatically assigns GROMOS force field parameters and partial charges. 61 For the simulation of the melt structures, we use the molecular dynamics (MD) code LAMMPS. 62 Throughout this work we use a timestep of 1 fs, a cutoff of 14 Å for short-range interactions in combination with the PPPM long-range solver for electrostatics and Nosé-Hoover style thermostating and barostating. The MD simulation is initialized by placing randomly rotated molecules on a grid, followed by a cascade of NPT runs to equilibrate the density. During an initial pressure ramp from 10 to 1 atm and at a temperature of 300 K, the loose initial structure is condensed into a solid over a period of 200 ps. Then the long-range solver is activated followed by an additional run of 200 ps at the final pressure of 1 atm. In order to further equilibrate the density, the system is heated from 300 K to 700 K and subsequently cooled back to 300 K, in each case as linear ramps of 200 ps. At the final conditions of 300 K and 1 atm, the system is equilibrated for another 500 ps. Finally, in the production runs, we extract snapshots in an interval of 1 ps over a period of 5 ns and in an interval of 1 fs over a period of 1 ps. A typical example for the development of the mass density during the equilibration and production runs is shown in Figure S13 of the SI .

DFT.
In order to generate a suitable conformer-energy dataset to train machine learning models, we performed density functional theory (DFT) calculations with the Turbomole 7.4 software package on randomly sampled conformers. 63 We used the resolution of identity approximation with a def2-SV(P) basis set on a M3 grid and the B3-LYP functional. 64 We extracted the total energy as well as vacuum HOMO and LUMO energy levels from DFT single point energy calculations. We simply used the last 20-100 snapshots of the MD trajectories with large timesteps, which leads to a set of less correlated conformers compared to the trajectories with a smaller timestep. Each snapshot has 512 molecules (mCP has 1000 molecules), which leads to approximately 10'000 -50'000 conformers and corresponding energy levels of each material for training. Additionally, we randomly picked 6 molecules in each simulation box and calculated DFT-energies for the full trajectories from 0 to 5000 ps with a total of approximately 30'000 conformational energies, which serves as an additional test set and to estimate approximate DFT-based disorder values. Moreover, we calculated energy levels using DFT along 6 trajectories for the small timestep MD trajectories of randomly selected molecules.
We further experimented with different basis set and grid options (see SI for details). There is a small dependence of the DFT results on the basis set size, which likely improves going to def2-TZVP and def2-QZVP, at the cost of a significant increase in computational effort, which makes a DFT reference uncertainty of a few meV seem acceptable.

ML.
For our machine learning model, we chose a fully connected feed-forward neural network (NN) of 3 layers with a leaky RELU activation function ( ) followed by a single linear classification 0.05 α = layer. We trained a separate NN for each material on the respective DFT-energy dataset generated in the previous step. As input features we used rotationally invariant, geometry derived descriptors, such as inverse distances, bond and dihedral angles. Because of the large size of some organic semiconductors in this study, we further restricted the number of inverse distances to a set of close neighbors per atom (including hydrogens). The input features were subsequently concatenated and standardized. The NNs were trained using the Adam optimizer with a learning rate decreasing from 5x10 -4 to 1x10 -5 over 400 epochs. We found a good convergence for a broad set of hyperparameters and selected a NN with 3 hidden layers à 1000 neurons before the final regression layer. To prevent overfitting, we applied L1-regularization for kernel and bias (10 -6 -10 -7 ). Note, an exhaustive hyperparameter optimization was not found to be necessary, which is why only the regularization parameters were tuned. For evaluation we choose a validation-training set split ratio of 0.1. .
Here, represents the variance and denotes the expectation value with respect to time or .
If both distributions can be approximated independently by normal distributions, the total conformational disorder is given by the quadratic sum of dynamic and static ar [U (t)] σ t 2 = V i,t disorder: .
We note that in the following discussion we neglect the electrostatic contribution and focus solely on conformational disorder.

Results and discussion
We

Learning curve for NPD, p-BPD, NPB, Spiro-TAD, and mCP, which shows the mean absolute error (MAE) as a function of training data. d) A zoom into an oscillation of b).
Applying  Table S1 and S2 .  Table S3 and S4 together with the ML validation error and the squared sum of dynamic and static disorder. We find an overall very good agreement of the squared sum with the total conformational disorder value, showing that the disorder is composed of a convolution of independent dynamic and static disorder contributions, which can not be concluded from the six DFT reference trajectories alone due to limited statistics.
In general, we find a notable contribution of the dynamic part to the total conformational disorder. As discussed later, the dynamic disorder has the potential to hinder charge transport. 66 The static disorder seems to be larger for molecules with long rotatable aromatic groups such as TCTA, where the vacuum HOMO level is delocalized and can localize inside groups depending on the conformation (see Figure S5 ). Also, in direct comparison of m-BPD and p-BPD with o-BPD, the more compact and sterically locked ortho-position shows the smallest static disorder. Similarly, with the addition of methyl-groups in the backbone of NPB (yielding a-NPD), static disorder decreases and the ratio between static and dynamic disorder is altered. However, the reduction in static disorder of the methyl-substituted a-NPD is compensated by an increase in dynamic disorder. Figure 3 .

For all n-type molecules the LUMO level and for all p-type molecules the HOMO level is plotted (see Scheme 1 ). All 38.3 million energy levels required for this analysis are generated by our machine learning models. All ML estimates of HOMO and LUMO
disorder are given in Table S3, S4 and Figure S8, S9. The corresponding values obtained from DFT reference trajectories with significantly smaller statistics are given in Table S1, S2 for comparison.
The lower panels show a measure of skewness of the static, dynamic and total energy level distributions. In particular for the static energy level distributions, we find substantial deviations from symmetrical Gaussian distributions.
While the total and dynamic disorder can be well described with a Gaussian distribution, the static disorder is often slightly asymmetric (see Figure 3b , lower panel in Figure 4 , and SI). Electrons and holes with hopping rates that are lower than the typical fluctuation time of the dynamic disorder will thus experience an asymmetrical distribution of energy levels. Depending on the type of asymmetry, this might lead to a decrease or increase in the concentration of tail states / shallow traps. In the example illustrated in Figure 3d , the holes traveling through NPB are in a skewed HOMO distribution that shows more tail states and thus hole traps than what would be expected from a Gaussian distribution of energy levels with the same standard deviation. The opposite behaviour was observed in Friederich et al. 55 for hole transport in PEDOT oligomers, where the distribution of HOMO levels showed an extended tail towards lower HOMO energy levels and a sharp edge towards more shallow HOMO energies. One possible explanation of this behavior might be the correlation between total energy and HOMO energy close to the equilibrium conformation. In case of a positive correlation, the HOMO distribution has a global minimum, leading to a distribution as shown in Figure 3d , whereas in case of a negative correlation, the opposite asymmetry is observed. This can be seen in the example of PEDOT and in various cases in Figure 4 , especially for distributions of static LUMO energy levels and for the distribution of static HOMO energy levels of Spiro-TAD.
Depending on the theoretical model used to describe charge transport, also the structurally resolved energy levels of the surroundings of a molecule can be required. For example in kinetic Monte Carlo simulations, charge transport is modeled as stochastic charge transfer between different sites subject to an external electric field. With our machine learning model, the time-resolved energetic neighborhood to a molecule can be explored.
In the following, we will analyse the time resolved energy difference between neighboring molecules at the examples of a small NPB cluster and compare it to typical hopping rates. The power spectrum of the energy level fluctuations of NBP molecules (see Figure 5d ) shows an approximate 1/f noise-like background, reaching down to a few THz. In addition, we find distinct vibrational modes leading to oscillations with a frequency of~17 THz (59 fs), which corresponds to the more pronounced oscillations plotted in Figure 5b . Furthermore, we find additional frequency components at~37 and~51 THz (59 fs and To compare these characteristic fluctuation and resonance times to typical hopping rates, we have to take into account the electronic coupling between NPB molecules, which quantifies the prefactor of the Marcus hoping rate, and which can be interpreted as a hopping attempt frequency. The electronic coupling strongly depends on the intermolecular distance (see Figure S6 ). The pair correlation function of NPB ( Figure 5a ) shows that a substantial amount of NPB molecules have their nearest neighbour at center-of-mass distances between 3 Å and 8 Å (with a stacking motive in which two NPB molecules are situated in close proximity and rotated approximately 90°relative to each other). The next nearest neighbours which are still in direct contact, thus allowing for non-vanishing electronic coupling matrix elements, are located at distances between 9-10 Å.
Electronic couplings and reorganization energies are typically obtained from quantum chemical J ij calculations. 67 For NPB, we obtain at a distance of~5 Å and a reorganization energy .01 eV J ij = 0 , leading to a maximum hopping rate of 4.2 THz for Marcus theory, assuming that the 0.16 eV λ = resonance condition is fulfilled, i.e. for = 0. This translates to a residence time of E j − E i + λ /W 1 ij 238 fs, which suggests that the oscillations visible in Figure 5b appear averaged for most charge carriers. When taking into account larger hopping distances, the mean electronic coupling is substantially lower, allowing the assumption that the energy levels appear time averaged. This means that the resulting mean energy levels of molecules should be drawn from the static disorder distribution. If we were to remove the distinct oscillations from the frequency spectrum (filter out frequencies higher than 10 THz), the dynamic disorder would shrink by 22 meV from 84 to 62 meV but the average resonance occurrence time increases from 11 fs to 67 fs for 5 neighbours considered.
Nonetheless, we also note that in case of molecules with higher mean electronic couplings than NPB (e.g.~0.1 eV, corresponding to a prefactor of the Marcus rate of 418 THz or 2 fs in resonance J ij condition), even distinct dynamic oscillations become comparable or slower than typical hopping rates, leading to a significantly increased effective disorder seen by charge carriers in a classical description. Although this is a classical picture and the analogue is not strictly true, the fast distinct oscillations in Figure 5b may be regarded as a way to bring neighbouring energy levels into resonance, or as phonons being absorbed to facilitate hopping as in Marcus theory. We would like to point out that with a ML model introduced into atomistic charge transport models, these relative contributions could be captured and quantified by introducing the time dependent energy levels. For example, a larger grid of realistic time-dependent energy sites could in principle be generated for Monte-Carlo simulations. Also in extensive 3D morphologies, a joint energetic and geometric representation will also reveal possible percolation pathways (see Figure S4 ). This could be explored in future work, integrating ML models in charge transport simulations.
To arrive at a static disorder, which can be used in mesoscopic models, e.g., in kinetic Monte-Carlo simulations, the static conformational disorder computed in this work has to be complemented with the electrostatic disorder which can be contributed to intermolecular electrostatic interaction, e.g., σ p molecular dipole moments. The electrostatic disorder can be calculated with the QuantumPatch method 21,34 , which yields meV for NPB. Assuming that electrostatic and static conformational 0 σ p = 7 disorder are not correlated, the total static disorder is equal to . This is close to 8 meV √ σ 2 s + σ 2 p = 8 the disorder parameter (87 meV) used in mesoscopic models of charge transport 68-70 that is used by a top-down approach to reproduce the experimentally observed hole mobility of NPB. Moreover, the hole mobility of NPB computed with the kMC model with a realistic morphology and assuming a disorder parameter 0.93 eV yields a hole mobility of 2x10 -3 cm 2 /Vs 71 close to experimental value (between 2x10 -3 cm 2 /Vs and 3x10 -3 cm 2 /Vs according to Chen et al. 72 and Tse et al. 73 ). This shows the high relevance of the developed method to accurately compute charge carriers mobility in disordered organic semiconductors.

Conclusions
In summary we have investigated static and dynamic conformational disorder for various device-relevant organic semiconductors by using machine learning models integrated in multiscale modeling simulations. The ML models were trained using DFT energies to predict the time-and space resolved energy level distribution of complete MD trajectories of disordered thin films. Thereby, static and dynamic disorder contributions were analysed using~38 million energy level predictions by the ML model, a number that would not have been feasible using quantum mechanical calculations. We firstly find that static disorder distributions of multiple relevant materials are highly asymmetrical.
Consequently the density of states either shows deep tail states in the HOMO-LUMO gap which is unfavourable for charge transport ( e.g. the HOMO distribution of NPB and TPBi), or it shows a steep decline towards the HOMO-LUMO gap which avoids trapping and thus favours charge transport ( e.g. the HOMO distribution of Spiro-TAD). We suspect that the underlying mechanism is related to the HOMO energy surface close to local minima of the potential energy surface. A closer study of this mechanism might lead to new molecular design principles.
We secondly find that characteristic oscillation times of molecular energy levels are below maximally possible Marcus hopping rates in disordered materials like NPB, suggesting that dynamic disorder averages out at typical electronic couplings of disordered materials, leaving static disorder as the key contribution to charge transport. Thirdly, we observed that the amount of static disorder is likely correlated to the size and the amount of rotatable side groups ( e.g. TCTA and BPD), provided that the respective molecular energy level is sufficiently delocalized over the molecule and therefore affected by conformational changes. This effect offers the opportunity to be used as a design rule for future organic semiconductors.
In summary, we showed that the integration of ML models as fast and accurate predictors in multiscale modeling workflows offers the opportunity to study a wide range of interesting and application relevant phenomena, e.g. time resolved energy level fluctuations, dynamic disorder and asymmetries in static energy level distributions, which were not accessible before due to prohibitively high computational costs.