By considering only one electronic state per molecule, charge transport models of molecular solids neglect intramolecular charge transfer. This approximation excludes materials with quasi-degenerate spatially separated frontier orbitals, such as non-fullerene acceptors (NFAs) and symmetric thermally activated delayed fluorescence emitters. By analyzing the electronic structure of room-temperature molecular conformers of a prototypical NFA, ITIC-4F, we conclude that the electron is localized on one of the two acceptor blocks with the mean intramolecular transfer integral of 120 meV, which is comparable with intermolecular couplings. Therefore, the minimal basis for acceptor–donor–acceptor (A–D–A) molecules consists of two molecular orbitals localized on the acceptor blocks. This basis is robust even with respect to geometry distortions in an amorphous solid, in contrast to the basis of two lowest unoccupied canonical molecular orbitals withstanding only thermal fluctuations in a crystal. The charge carrier mobility can be underestimated by a factor of two when using single site approximation for A–D–A molecules in their typical crystalline packings.

## I. INTRODUCTION

Organic semiconductors constitute a broad class of functional materials for optoelectronic and energy applications.^{1} Design of such materials requires an accurate description of their electronic properties,^{2–7} with the most accurate models developed for molecular solids, such as pentacene^{8–10} or naphthalene.^{11,12} Despite these advances, modeling charge transport in molecular solids continues to pose challenges. In molecular crystals, conventional models, such as band transport or hopping transport, are not always appropriate.^{6} The current state-of-the-art approach involves evolving the trajectory of charge carriers using surface hopping or other methods correctly describing electron–phonon interaction in organic semiconductors.^{8} This is only feasible in a *coarse-grained electronic basis*: typically, with one molecular orbital (MO) per molecule, so that the electronic wave-function is approximated by a linear combination of such MOs.^{13} In particular, the highest occupied molecular orbital (HOMO) is used for the hole transport, and the lowest unoccupied molecular orbital (LUMO) is used for the electron transport. Importantly, the coarse-grained basis is *fixed* in the same sense as contracted Gaussian functions are fixed, and only translations and rotations of the entire MO are allowed. The robustness of this description relies on the energetic separation of HOMO/LUMO from the remaining MOs. This separation, typically more than 1 eV, ensures that the solid-state electronic bands composed of different MOs do not overlap. The rigidity of molecules makes the coarse-grained basis robust with respect to fluctuations in molecular geometries.

For large molecules, the energy separation between frontier MOs is comparable to both their vibrational broadening and intermolecular dispersion (a few tenths of eV), making the basis of one MO per molecule incomplete for a description of frontier MOs of molecular solids. If the molecule is rigid, such as polycyclic aromatic hydrocarbons, the simplest reasonable approximation is to include all quasi-degenerate MOs into the coarse-grained basis. In nonrigid molecules, however, geometry fluctuations under ambient conditions can be so large that MOs of one conformation are no longer representable on the basis of canonical MOs of another conformation. This case requires a more sophisticated electronic coarse-graining, which is the subject of this work. Specifically, we will compare the accuracy of different coarse-grained bases in representing frontier orbitals of a nonrigid molecule subject to fluctuations of its geometry under typical conditions for organic semiconductors, mainly in a solid state. It should be noted that inaccurate electronic coarse-graining might result in not only a deficient basis but also a wrong intramolecular charge dynamics in the approximation of instantaneous transitions between degenerate electronic states, which is reasonable for such molecules as fullerenes.^{14}

To perform a detailed investigation, we will consider the case study of the ITIC-4F molecule^{15} whose bulk geometry in both crystalline and smectic phases is known.^{16} ITIC-4F belongs to a class of high-performing non-fullerene acceptors (NFA),^{17} majority of which have an acceptor–donor–acceptor (A–D–A) molecular architecture^{18–21} (Fig. 1) with a large separation between the acceptor blocks, resulting in a quasi-degenerate LUMO.^{22} The most appropriate approach for studying the influence of molecular fluctuations on the electronic structure would be an all-atom non-adiabatic molecular dynamics (MD).^{23} However, it is computationally demanding for a reasonably sized cluster of NFA molecules, not talking about bulk solids.^{24} To investigate structural models with thousands of atoms, such as amorphous ITIC-4F, alternative strategies must be adopted. Here, we sample molecular conformations using a computationally efficient approach: a combination of the classical MD with a force-field and classical/quantum harmonic oscillator models with vibrational modes calculated using density functional theory (DFT). A statistical analysis of electronic states using single-point electronic structure calculations can then be performed based on the sampled structures. At the end, we will determine what is the most appropriate coarse-grained model of NFAs, taking into account a trade off between the accuracy and complexity of the model.

This paper is organized as follows: We start with a description of the methodology, followed by a discussion of the electronic structure of isolated ITIC-4F molecule as a typical representative of the A–D–A class of molecules. Next, we investigate the stability of various coarse-grained bases with respect to different types of geometrical distortions and then analyze fluctuations of tight-binding (TB) elements. Finally, we study how errors in coarse-graining influence calculated charge transport properties.

## II. METHODOLOGY

Localized molecular orbitals (LMO) are obtained by orthogonal projection of canonical MOs onto the basis of the selected block of atoms. In fact, the choice of localization procedure is noncritical for the scope of this work. Nevertheless, a reference to the localization algorithm and code is given in supplementary material. The methodology has been thoroughly tested for organic,^{22} metal-organic,^{25} and inorganic^{26} semiconductors. Multiple LMO models are employed, designated as *N*-LMO, where *N* represents the number of localized molecular orbitals utilized in the model.

One of the aims of this work is to grade various coarse-grained electronic bases. As a measure of basis quality, we take its ability to describe LUMO of a distorted molecule (all bases are naturally defined to give an exact LUMO for the undistorted geometry). Specifically, to evaluate the accuracy of a given basis of molecular orbitals *ψ*_{1}, *ψ*_{2}, …, *ψ*_{n} in approximating a target molecular orbital *ψ*′, we quantify the deviation using a norm of Δ*ψ*′ = *ψ*′ − *∑*_{i}*c*_{i}*ψ*_{i}. The key technical difficulty is how to define the coefficients *c*_{i} and the norm. Indeed, *ψ*_{i} and *ψ*′ are provided in the same atomic orbital (AO) basis set but differently oriented and displaced due to different molecular geometry. Therefore, there is no obvious choice for the norm and for the procedure how coefficients *c*_{i} are determined. Since, in TB models, molecular deformations are encoded in the TB Hamiltonian parametrically, we do not transform the AO basis set (except for molecular block rotations) but consider only a vector space of coefficients itself. However, molecular deformations also change the AO overlap matrix (*S* for *ψ*_{i} and *S*′ for *ψ*′), and this change must be taken into account explicitly. The most straightforward solution is to orthogonalize the AO basis, and we consider this option with the symmetric orthogonalization $\psi i\xb0=S\psi i$, $\psi \u2032\xb0=S\u2032\psi \u2032$. Then, the natural choice for the norm is Euclidean. However, to keep the overall complexity at the level of matrix multiplication, the calculation of the square root of a matrix should be avoided. Although there is no obvious choice for the norm in this case, the generalized Euclidean norm seems to be reasonable: one with the undistorted overlap matrix *S*, which is the default method, and another one with the unit overlap. Since all the mentioned norms are quadratic, we define the *basis deficiency* as ‖Δ*ψ*′‖^{2} with three alternative interpretations of the norm symbol: ‖·‖_{2} for the Euclidean norm, ‖·‖_{S} for the generalized Euclidean norm with the overlap matrix of the undistorted molecule, and ‖·‖_{O} for the Euclidean norm in the orthogonalized basis. If the norm is the generalized Euclidean and *ψ*_{i} are orthonormal in this norm, then coefficients $ci=\psi i+S\psi \u2032$ minimize ‖Δ*ψ*′‖ so that the basis set deficiency becomes $\Vert \Delta \psi \u2032\Vert 2=\psi \u2032+S\psi \u2032\u2212\u2211i|\psi i+S\psi \u2032|2$, where *ψ*^{+} means the Hermitian conjugate. In other cases, the minimization is more complicated, and the resulting expression for *c*_{i} is less meaningful. For this reason, we fix $ci=\psi i+S\psi \u2032$ for all considered norms. If *ψ*_{i} is a LMO, it is rotated together with its biorthogonal vector $\psi i+S$. The resulting set of LMOs becomes non-orthonormal but the LMO overlap is usually close to the unit matrix. Effects of local rotation of AOs due to geometric deformations are not considered.

The crystalline and partially ordered smectic morphologies are taken from classical MD simulations performed in Ref. 16. The MD supercell contains 2000 molecules, 500 of which are used in this analysis.

Technical details are given in Sec. S1 of supplementary material, including parameters of quantum chemistry calculations, comparison of different variations of the localization algorithms (Table S1), treatment of anharmonic vibrational modes (Table S2), and statistical analysis of simulated fluctuations (Table S3).

## III. RESULTS AND DISCUSSION

The electronic structure of ITIC-4F is typical for NFAs with the A–D–A structure:^{22} the HOMO is localized on the donor block and is well separated from the next occupied MO, whereas the LUMO is quasi-degenerate [Fig. 1(a)] with only 0.2 eV separation between frontier MOs (Fig. S3). Further insight can be obtained by block-wise localization of several lowest unoccupied MOs as visualized in Fig. 2 (see also localization of the entire *π*-system in Fig. S5). It shows two weakly coupled lowest unoccupied LMOs localized on the two acceptors (see also Figs. S4 and S5 and Table S4). Consequently, the minimal coarse-grained model should include at least these two LMOs. In the 2-LMO model, the transfer integral is 120 meV, which is comparable with intermolecular acceptor–acceptor couplings in many NFAs, such as Y6, IDTBR, and ITIC.^{22} The intramolecular electron transfer is transmitted through the donor block, but according to Fig. 2, multiple transfer channels exist with at least two lowest donor-block LMOs strongly coupled to the lowest acceptor-block LMOs. Therefore, the 6-LMO coarse-grained model (two LMOs per each block) should provide a more complete description of the low-energy electronic states, while being still separated from the LUMO of the *σ*-symmetry subsystem (see also Figs. S6 and S7). All LMOs in 2-, 4-, and 6-LMO models are visualized in Fig. S4, and the corresponding TB Hamiltonians are given in Table S4.

Inter-block connections are flexible but with a high rotation barrier of almost 0.5 eV, thus suppressing large amplitude librations. In the harmonic approximation, these librations are intermixed with other low-frequency vibrational modes (starting from 1 meV or 5 cm^{−1}), corresponding to the in-plane and out-plane bending of the molecule. The 180°-flip of each dihedral costs only 50 meV, so conformers should be present in the amorphous phase.

Electron polaron is well localized on a single acceptor [Fig. 1(c)] in a highly polarizable environment with large energy gain relative to the symmetric charge distribution (270 meV in the water-like polarizable continuum model). In vacuum, the symmetric distribution is unstable but the relaxed anion natural orbital (NO) differs from the LUMO of the neutral molecule insignificantly (Fig. S16).

To study the influence of molecular fluctuations on the electronic structure, we use four statistical ensembles to mimic the effects of static and dynamic disorder, both classical and quantum in the latter case. The first two ensembles are obtained from a single snapshot of the room temperature classical MD of crystalline and amorphous (“smectic” phase) solid state ITIC-4F. Because the MD uses a classical force field whose accuracy is expected to be lower than DFT, we also generate harmonic vibrations in classical and quantum Boltzmann statistics (the other two ensembles). For this purpose, we use normal modes of the neutral molecule but with renormalized frequencies of several anharmonic modes, see Table S2 for details. Such renormalization does not treat mode mixing, resulting in an excess of higher-energy fluctuations (Fig. S1), as evidenced by the higher value of energy variation (Table S3). The quantum ensemble generates an excess of high-frequency modes via zero-energy vibrations: although these modes interact coherently with electronic degrees of freedom, the “decoupled” description of their influence on the electronic structure seems to be meaningful.

Statistical analysis of the electronic structure sampled in all four ensembles immediately shows that, in contrast to HOMO, the LUMO is unstable with respect to room-temperature fluctuations in the sense that, for the majority of distorted geometries, the LUMO differs substantially from the LUMO of the undistorted molecule (Figs. S8 and S10). Consequently, the combined effect of self-localization and distortion leads to complete localization of extra charge on a single acceptor even in crystal, see Fig. 3(a) [also Fig. 5(a), Fig. S8, Table S5, Fig. S16, and cf. HOMO in Fig. S17].

While the basis of two LUMOs of the unperturbed molecule is capable to describe the LUMO of the majority of molecular conformations in crystal, at least two LMOs are needed to describe the LUMO of molecular conformations in amorphous solid, see Fig. 3(b). In the latter case, the addition of a few other LMOs (6-LMO model) makes the LMO basis as accurate as the entire-band MO basis [Fig. 3(b)] but more robust with respect to large deformations (Fig. S11d).

Detailed analysis of extreme points in the sampled datasets shows three cases. If distortions are small, the LUMO is fully delocalized across the molecule, and any coarse-grained basis is accurate within a few percent, Fig. 4(a). In contrast, if distortions are substantial, the LUMO is localized on one of the acceptors, and the minimal coarse-grained basis should include at least two MOs or LMOs, whereas the use of the more accurate 6-LMO model reduces the basis deficiency to a few percent, Fig. 4(b). In the extreme cases of large-scale molecular deformations such as the 90°-rotation of flexible dihedrals, the use of LMOs remains meaningful albeit with low accuracy, see Fig. S11d.

To make a more rigorous statement about electron localization, we analyze the distribution of coefficients of expansion of electron MO and NO on the 2-LMO basis. We immediately see that the electron NO is localized on one of the acceptors in the majority of molecular conformations in both crystalline and amorphous solids, see Fig. 5. In contrast, LUMO shows no electron localization but only a broad distribution of coefficients in the same range of angles (Fig. S13).

Having determined the optimal coarse-graining of the ITIC-4F molecule, we can analyze fluctuations of the TB Hamiltonian constructed in the 2-LMO and 6-LMO bases. In the first case, the intramolecular Hamiltonian is two-dimensional with three parameters: average onsite energy (*ɛ*_{1} + *ɛ*_{2})/2, energy half-offset *ϵ* = (*ɛ*_{2} − *ɛ*_{1})/2, and transfer integral *t*. Fluctuations of *ϵ* and *t* in amorphous solids are plotted in Fig. 6(a), showing a single-centered distribution without a visual correlation between *ϵ* and *t*. The standard deviations are 50 and 25 meV, respectively, which is much smaller than the absolute value of *t*, see Table S10. This means that the two acceptor blocks are always strongly coupled, and the main disorder effects are in onsite energy fluctuations. Interestingly, both the mean value and standard deviation of *ϵ* and *t* are nearly the same for crystal, amorphous solid, and classical vibrations of the isolated molecule (Table S10), implying that the environmental effects have a minor influence on intramolecular energetics, despite having a large influence on type and amplitude of geometrical distortions. In contrast, quantum effects, such as zero-point vibrations, almost double the standard deviations (Fig. S18), thus emphasizing the importance of quantum treatment of high-frequency modes^{5} (see the caption to Table S10). For the 6-LMO model, fluctuations of acceptor–donor energy offset show a bias toward lower values, see Fig. 6(b), thus explaining why this model provides a noticeably better description of LUMO under large distortions.

To understand how errors in the coarse-graining of individual molecules influence the accuracy of the calculation of electronic properties of a molecular solid, one has to distinguish two cases. If we calculate the electronic structure at fixed positions of atoms, then the only thing important is the deficiency of the coarse-grained basis. Therefore, for crystals, the basis of two LMOs or two MOs gives an accurate description of the bottom of the conduction band, whereas the use of a single MO per molecule misses the intermixing of LUMO and the next LUMO by intermolecular couplings. In contrast, if we calculate charge carrier transport in the hopping regime, we assume that each site orbital remains unchanged upon geometry fluctuations, consequently using MOs for electrons in NFAs is inappropriate at all.

To quantify the error, we consider three topologies of electronic connectivity commonly observed in crystals of A–D–A molecules:^{22,27} one-dimensional slipstack (EH-IDTBR), two-dimensional brickwork (ITIC-4F, Fig. S21), and three-dimensional wiremesh (*o*-IDTBR, Fig. S22). We compare 1-site vs 2-site models in terms of two electronic structure parameters: tensor of inverse effective masses and tensor of squared hopping amplitudes. Accuracy in the evaluation of these quantities contributes to the accuracy of electronic transport modeling because the charge carrier mobility is proportional to these parameters in the two limiting cases: scattering of free charge carriers and hopping of localized ones, see Section S4 in Ref. 22. Results of the comparison depend on multiple parameters including crystal geometry and electronic couplings and are compiled in Table S11. The main conclusion is that under the most typical conditions when the intramolecular transfer integral is substantially larger than intermolecular ones, the mobility is underestimated by a factor of two. This can be explained by the reduction of intermolecular transfer integrals by this factor due to intramolecular delocalization of the wave-function assumed in the 1-site model. In the case of comparable intra- and intermolecular couplings (best performing NFAs), the error of the 1-site model is less systematic but anisotropic by a factor of 2–3, implying that the mobility tensor is substantially distorted upon merge of the acceptor sites. For example, in brickwork geometry, the ratio of mobilities parallel and perpendicular to the long axis of the molecule can be overestimated by a factor of 3. In practice, anisotropy can be directly compared to single-crystal mobility measurements to distinguish between 1-site and 2-site models. In addition, the electronic band structure transforms significantly upon change of transfer integrals: for example, in brickwork geometry, which has a honeycomb topology of electronic connectivity, the band structure consists of two sub-bands connected via Dirac cone if and only if the transfer integrals satisfy the triangle inequality (see Fig. S23).

We have demonstrated that the accurate coarse-grained electronic model of a molecular solid should involve proper coarse-graining of intramolecular electronic degrees of freedom: long flexible *π*-conjugated molecules require at least one electronic site per each electronically active rigid block. A robust coarse-grained basis is provided by localized molecular orbitals, whereas canonical MOs are unstable with respect to fluctuations of intramolecular geometry. In the case of A–D–A type molecules used as acceptor material in solar cells, the electron is localized on one of the acceptor blocks in the majority of molecular conformations in both crystalline and amorphous solids, so that LUMO of the unperturbed molecule is not representative at all. The resulting error in the estimation of charge carrier mobility and its anisotropy can be as large as a factor of 2. Importantly, this is a purely electronic contribution to the error, substantially exceeding errors of state-of-the-art DFT calculations of the electronic structure of organic semiconductors.

## SUPPLEMENTARY MATERIAL

The supplementary material contains details of the computational methodology and additional figures and tables on the electronic structure of undistorted and distorted conjugated backbone, fluctuations of TB Hamiltonian elements, and calculations of charge transport parameters.

## ACKNOWLEDGMENTS

We thank Wenlan Liu, Mukunda Mandal, and Naomi Kinaret for the fruitful discussions. This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science by Los Alamos National Laboratory (Contract No. 89233218CNA000001) and Sandia National Laboratories (Contract No. DE-NA-0003525). D.A. acknowledges the KAUST PSE Division for hosting his sabbatical in the framework of the Division’s Visiting Faculty program. D.A. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for financial support through the collaborative research centers TRR 146, SPP 2196, and Grant No. 460766640.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Andriy Zhugayevych**: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Project administration (lead); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Kun-Han Lin**: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (equal). **Denis Andrienko**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (lead); Software (lead); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

## REFERENCES

*Organic Electronics: Foundations to Applications*

*Ab initio*simulations with mode-specific treatment of molecular vibrations

_{60}fullerene anion: A combined X- and D-band EPR and DFT study