Read between the Molecules: Computational Insights into Organic Semiconductors.

The performance and key electronic properties of molecular organic semiconductors are dictated by the interplay between the chemistry of the molecular core and the intermolecular factors of which manipulation has inspired both experimentalists and theorists. This Perspective presents major computational challenges and modern methodological strategies to advance the field. The discussion ranges from insights and design principles at the quantum chemical level, in-depth atomistic modeling based on multiscale protocols, morphological prediction and characterization as well as energy-property maps involving data-driven analysis. A personal overview of the past achievements and future direction is also provided.

Introduction

Since the early experiments in the 1940s, immense progress in fabrication, particularly solution processing, and characterization of small-molecule and polymer organic semiconductors transported these systems from the laboratory bench to billions of hands and households. From phone and television displays to nanoscale memory and sensing devices, organic field effect transistors (OFETs) and light emitting diods (OLEDs) play an increasingly important role in modern technology. π-Conjugated cores are the typical buildings blocks in such molecular semiconductor materials, offering facile tuning of the key electronic properties via diverse chemical modifications. Further advantages of organic semiconductors compared to conventional silicon-based materials include mechanical flexibility, lightweight and easy and inexpensive solution processability.[1−3] Though these cover a wide range of materials, from molecules to polymers and from single- to multicomponent blends, in this Perspective we focus on molecular, primarily crystalline single-component organic semiconductors. Continuous refinement of design strategies and fabrication techniques allowed them to now routinely display impressive charge carrier mobilities of over 1 cm2 V–1 s–1, reaching as high as 20–40 cm2 V–1 s–1 in single crystals and even hundreds of cm2 V–1 s–1 in ultrapure samples at low temperatures.[1,4,5]

Advanced understanding of the chemical and physical factors determining the properties and performance of molecular semiconductors is an obvious prerequisite toward systematic improvement of their mobility and stability. In this context, two fundamental challenges are (1) the lack of a universal theory of charge transport and (2) the complex dependence of transport characteristics on the material’s morphology. These challenges have been discussed in a number of comprehensive reviews.[6−13] Here, we focus on the ingenious and, at times, only available solutions to these problems offered by computational chemistry. Above that, we select illustrative examples that yield conceptual insights into the structure–property relationships of molecular semiconductors and to allow for better understanding of their charge transport properties, ultimately aiming at better performing systems. We especially discuss how diverse computational approaches, from quantum chemistry-based understanding of intermolecular interactions to multiscale modeling of (disordered) crystalline morphologies and high-throughput computations of charge mobility, advance the field of organic semiconductors (Chart ). While the chemical nature of the molecular building blocks is unquestionably important on its own, here we offer our perspective on the intermolecular context of organic semiconductors, crucial both in terms of the interactions within the material and its structural organization (packing). Overall, this Perspective assembles a relevant selection of recent studies demonstrating that the field of organic semiconductors represents an ideal playground for further developing and exploiting all types of atomic-scale modeling methods ranging from advanced quantum chemical tools to the latest innovations such as machine learning (ML).[14−17] Of interest to the experimental community, it also highlights strategies to acquire in-depth chemistry-based information, necessary to advance the field.

Chart 1

Focus Areas of this Perspective

Fundamental Aspects of Atomic-Scale Modeling of Semiconductors

Before discussing specific computational approaches to organic semiconductors, we briefly outline two fundamental challenges in this area, i.e., the choice of the theoretical model underlying the transport computation and the treatment of the material’s structure in it.

Charge Transport Theories

A number of theoretical descriptions, corresponding to distinctly different regimes of charge transport, are continuously debated in the literature.[6−13] Band transport describes conduction of delocalized charges in a low temperature limit. Within this framework, charge mobility depends upon the effective mass of the charge carrier and the relaxation time of the band (i.e., the average time between collisions). Band-structure computations provide access to these and related parameters in order to quantify the rate of band transport.[6,18,19] As temperature rises, the band narrowing[20] leads to charge localization, switching on the thermally activated transport. Traditionally described as a charge hopping between the neighbor molecules, this mechanism has been challenged in the recent years, with the observed mobility being instead attributed to the presence of structural defects and trap sites. Nonetheless, the hopping model in the framework of the Holstein small polaron theory is still very informative and routinely employed to rationalize the transport properties of organic semiconductors due to its simplicity and intuitiveness. Depending on the treatment of the molecular vibrations, several flavors of the thermally activated hopping exist: the popular Marcus theory,[21−23] the Levich–Jortner formalism,[10,24,25] and the spectral overlap method.[26−28] Furthermore, the band theory can be extended to include not only the local electron–phonon coupling (also called Holstein coupling, diagonal dynamic disorder, and corresponds to the reorganization energy in Marcus theory)[29−31] but also the nonlocal electron–phonon coupling (Peierls coupling, off-diagonal dynamic disorder).[32,33] The latter is closely associated with the static and dynamic disorders, which are often inseparable in organic semiconductors, further complicating the modeling. An ultimate goal, to cater to various transport regimes for a broad range of temperatures, becomes achievable with the development of methods capable of describing the band-to-hopping crossover.[34−36] Importantly, regardless of the actual equation used to compute the charge mobility μ, it is generally greater in systems with smaller local electron–phonon coupling (reorganization energy, λ) and larger electronic coupling (also often called “transfer integral”, V).[1,4] Methods for computing[13,37,38] and chemical means for tuning the reorganization energy[39−42] are discussed in the cited literature. Electronic coupling, the measure of “communication” between the neighbor cores within the crystal structure, can be evaluated using one of several available theoretical models.[43−52] Obviously, V is strongly dependent on the morphology of the material (the nature and diversity of the nearest-neighbor dimers within the crystal), a topic that is extensively addressed in the next sections. Ultimately, any study modeling organic semiconductors is bound to choose the transport theory/regime, but this Perspective essentially covers examples based on the hopping transport model in the framework of Marcus theory.

Morphology of Molecular Semiconductors

The other crucial step in modeling organic semiconductor involves choosing a way to approximate and describe the system’s structural arrangement. Assemblies of π-conjugated cores feature such diverse morphologies as the disordered amorphous and the ordered crystalline phases. The latter exist in an array of packing motifs, from herringbone, typically associated with suppressed orbital overlap between the adjacent cores but often featuring large charge mobilities due to polarization effects (e.g., OFET mobility of >3.0 cm2 V–1 s–1 in dinaphtho[2,3-b:2′,3′-f]thieno[3,2-b]thiophene, DNTT),[48,53,54] to π-stacking involving strong electronic couplings (Figure ).[1] In some systems, the existence of several polymorphs,[55] whose stability in molecular materials has been recently shown to correlate with the number of short intermolecular contacts,[56] further inflates the structural diversity and complicates computational analysis. Variation in the electronic couplings between different pairs of neighbor cores (i.e., in different directions) results in the anisotropy of charge mobility in many crystalline organic semiconductors. In cases where the experimental crystal structure is available, it can be used to extract the representative nearest-neighbor dimers, compute their electronic couplings and construct the angular resolution mobility anisotropy curve.[57−59] A more sophisticated alternative, often referred to as the multiscale approach, takes into account the entire crystal structure and involves computing electronic and local electron–phonon couplings, running molecular dynamics simulations to access the energetic, configurational and dynamic disorder, and finally performing the diffusive charge dynamics simulation to evaluate the bulk mobility.[60−62] In the absence of experimental crystal structures, which is often the case for newly designed systems, the morphology of ordered systems is often approximated by scanning various dimer geometries, generated either manually by systematically varying their structural parameters[43,63,64] or in an automatized manner using, for instance, random search algorithms,[65] and constructing two-dimensional maps of their electronic couplings; for disordered systems, molecular dynamics (MD) simulations are often employed to probe the morphologies.[66] However, neither of these methods affords insights beyond the dimer level. Thus, crystal structure prediction (CSP) becomes necessary to achieve a complete representation of the complex material morphology.[67] In this Perspective, we discuss these different approaches to miscellaneous packing scenarios.

Figure 1

Typical aggregates and crystal packing motives of the π-conjugated cores: (A) lamellar π–π stacking motifs with one-dimensional charge carrier channel, (B) brick-stone or brick-wall (also called β-sheet) arrangement with two-dimensional π–π stacking; (C) γ-packing with slipped face-to-face π–π stacking; (D) herringbone face-to-edge packing without face-to-face π–π overlap.

Dimer Model

Focus on Energy

The most straightforward, albeit limited, approach to understand the relationship between charge transport properties and geometrical structures involves computing and comparing the interaction energies, electronic couplings (in the framework of properly orthogonalized monomer orbitals[46−48]) and reorganization energies of the representative dimers of diverse molecular cores.[68] Aside from this routine, rationalizing the physical nature of the interactions behind the molecular packing and their connection to the resulting charge transport properties is crucial for modeling and designing crystalline organic semiconductors.[69−71] In this regard, quantum chemistry offers diverse tools to analyze and classify typical noncovalent patterns (e.g., the π–π stacking, common in these systems) or to quantify and identify the nature of the dominant interactions.[72−76]

In the present context, schemes that allow for the decomposition of the intermolecular interactions between two molecules into physically meaningful energy contributions (e.g., London dispersion, exchange, induction, electrostatics, etc.) are rapidly becoming a popular means of rationalizing the morphological and transport properties of organic semiconductors. Numerous flavors of these so-called energy decomposition analysis (EDA) schemes exist, each relying on simple approximations or on a different central quantum chemical concept (i.e., molecular orbitals, electron density, the Hamiltonian) to proceed with the energy decomposition. An early work by Azumi et al.[77] demonstrated the utility of comparing low and high-level interaction energies (Hartree–Fock, MP2, extrapolated CCSD(T) limit) of 17 model thiophene dimers to rationalize the crystal packing in substituted quaterthiophene crystals. In particular, they exploited correlation energies together with distributed multipole/polarizability analysis[78] to assess the relative importance of London dispersion and electrostatics in the model systems. They concluded that dispersion interactions are the major source of attraction although electrostatic contributions further stabilize the perpendicular thiophene dimers considerably, thus explaining the preference for the herringbone structures in the crystals of nonsubstituted oligothiophenes. Symmetry-adapted perturbation theory[72,74,79−81] (generally at the SAPT0 level of approximation) is the most frequently employed scheme to offer a nonempirical definition of the individual energy contributions and to provide insight into the driving forces behind the resulting material morphology. In an extensive study of naphthodithiophenediimide (NDTI) thiophene α-substituted derivatives,[82] SAPT along with an exploration of the interactions in terms of intermolecular contact type (i.e., Hirshfeld surface analysis, vide infra)[83] were used to rationalize the change from a less conductive herringbone packing in the nonsubstituted NDTI to a more efficient 2D brickwork in its chlorinated derivative. In the substituted NDTI-Cl, the electrostatic contribution to the total interaction energies originating from halogen or hydrogen bonds was shown to play a key role in the change of packing patterns. In the crystals of another common organic semiconductor, rubrene, a desired π-stacked arrangement is conditional upon the planarity of the central tetracene core, which in turn depends on its substitution pattern. Risko, Brédas et al.[84] employed SAPT to explain the energetic origins of this effect. They showed that in isolated rubrene, the central tetracene core is twisted such that the Pauli repulsion between the neighbor phenyl moieties is minimized. However, its planarity, associated with more favorable π–π packing, can be restored in the crystal bulk by means of mitigating the increasing Pauli repulsion through enhancing the stabilizing contributions (dispersion, electrostatic and induction terms), tunable via chemical modification of the substituents.

Side-chain engineering to modulate the crystal packing often entails the introduction of long alkyl chains. Density functional theory (DFT) computations (at the B97-D level) of the intermolecular interaction energies in alkyl-substituted benzothieno[3,2-b][1]-benzothiophene (BTBT) crystals indicated that long alkyl chains enhance the stability of the layered-herringbone packing, which is known to afford high-performance organic thin-film transistors.[85] The interaction energies were decomposed in a fairly crude way into the correlation (considered there to be mainly composed of dispersion) and Hartree–Fock contributions, the latter was further broken down into the orbital–orbital, electrostatic and induction (computed from distributed multipole and atomic polarizabilities) energies. This allowed identification of dispersion as the main stabilizing component and deduced a rule-of-thumb for preferentially stabilizing the layered-herringbone packing: the ratio of the total intermolecular attractive forces between the T-shaped and slipped parallel contacts should be ca. 3:2. In the absence of side-chain groups, the symmetry of the molecular cores itself can influence the morphology and transport properties. In this context, the term “disordermer” has been introduced to describe the packing isomerism in molecular crystals. For the case of benzodithiophene, the interaction energies in various disordermers were compared using SAPT, which indicated that the main difference between isomeric dimers arises from the exchange contribution and that arrangements with direct S–S contacts and correspondingly high electronic couplings are generally disfavored by this term.[86] Polymorphism is yet another facet of the morphological diversity of organic semiconductors. Risko et al. employed SAPT to compare the intermolecular interactions in triisopropylsilylethynyl (TIPS) pentacene and its triethylsilylethynyl (TES) analog, for which they also generated different crystal packing configurations in silico and mapped their dimer potential energy surface.[87] The functionalization of pentacene typically leads to polymorphic situations, where the chromophore π-cores spatially overlap with TIPS and TES adopting brickwork and slipped-stack arrangement (see Figure ) respectively. In fact, dispersion interactions between the trialkylsilylethynyl groups (as opposed to those between the cores) were found to be the driving force between the tighter brickwork packing in TIPS-pentacene and the less compressible slipped-stack-in TES-pentacene. Processing strategies to impose more favorable packing in TES-pentacene were proposed. Overall, these studies generally highlight the crucial role of dispersion as the key stabilizing force counterbalancing the repulsive exchange and imposing specific packing in the engineered dimer-based structures.

However, alternative chemical strategies can be employed to amplify a different stabilizing component–the charge penetration contribution to the electrostatic interactions.[88−93] For example, combined SAPT and distributed multipole analyses (DMA)[78] of the interaction energetics in the dimers of heteroaromatic π-conjugated cores demonstrated that, in contrast to the dispersion-driven acene dimers, systems featuring heavier and more diffuse heteroatoms (e.g., sulfur, selenium and phosphorus) are driven by charge penetration.[94] Such systems also typically feature greater electronic couplings, suggesting that charge penetration can be utilized to enhance stability and mobility in molecular semiconductors. Extending this concept from organic semiconductors to noncovalent (dimer) molecular junctions, we have shown that similar chemical patterns determine the transport in these two types of molecular electronic assemblies. The interplay between the molecular and intermolecular parameters combined with the role of the packing motif generates a unified picture of noncovalent molecular electronics (Chart ).[95]

Chart 2

Spectrum of Relationships between the Molecular and Intermolecular Factors, which Dominate the Transport Properties of Different Noncovalent Molecular Electronic Architecturesa

Reprinted with permission from ref (95). Copyright 2018 American Chemical Society.

Despite all of the insights obtained from dimer-level studies, the shortcomings of this approach, such as neglecting bulk morphology, are rather obvious. Thus, it is necessary to continue developing EDA schemes that permit the treatment of many-fragment systems.[96−98] Moving toward realistic morphologies, EDA can also serve to parametrize and/or benchmark force fields.[99−101] For example, to assess the performance of commonly used force fields for describing the intermolecular interactions in molecular semiconductors, Engels et al. compared the results of different energy decomposition schemes (i.e., SAPT and localized molecular orbital energy decomposition analysis LMO-EDA) with those obtained using the MM3, OPLS-AA and AMOEBA force fields.[102] They found that in the model dimers of both the apolar (e.g., acenes) and highly polarized π-systems (merocyanines and squaraines), the shape of the potential energy surface (PES) is determined by variations in the highly specific short-range (exchange) repulsion forces. Several approaches were introduced to mimic the intermolecular potentials with electrostatics and exchange using distributed quadrupoles,[103] van der Waals potentials[104] or a simple π-orbital overlap-based force-field ansatz.[105]

Focus on the Electronic Structure

In addition to analyzing the dominating energetic contributions between the building blocks of organic semiconductors, insight can be gained by considering the features of the spatial distribution of their charge, multipole or intermolecular contact. For example, a combination of contact distance mapping (Hirshfeld surface analysis) with molecular electrostatic potential plots, polarizabilities and quadrupole moments was employed to compare a range of crystalline phases of nine heteroaromatic pentacene cores featuring diverse packing motifs.[106] Hirshfeld surfaces served to analyze the prevalence of contact points and their spatial distribution in the molecular crystals, which were then sorted into three groups according to their resemblance. The electrostatic nature of the underlying interactions (contacts) was then analyzed to identify that the herringbone packing is favored for systems with uniform electrostatic potentials. Cores with peripherally perturbed ESPs (due to electronegative substituents) instead crystallize in columnar and brickwork morphologies.

Other qualitative classes of approaches, which can reveal the presence of noncovalent interactions in real space, have also made their way into the area of organic semiconductors. Relevant examples are the noncovalent interactions index (NCI),[107] which exploits the reduced density gradient at low-density values, and the Density Overlap Regions Indicator (DORI),[108] which reveals regions with pronounced density overlap. In essence, the noncovalent domains identified by these approaches should not be used to rationalize the trends in interaction energies. Instead, these approaches can serve to visualize dominant nonbonded contacts and help understanding the influence of the density-based features of a given packing motif on the charge transport properties. A recent application of NCI in the field of organic semiconducting molecular crystals includes the computations on all nearest-neighbor dimers of dihydroindolocarbazole isomers extracted from their crystals. The analysis revealed that multiple intermolecular NH···π and CH···π interactions with energies close to common NH···N hydrogen bonds are associated with higher electronic couplings for hole transport and, correspondingly, better charge mobility in the crystals.[109] In this latter study as well as in ref (82), topological analysis of the density using Bader’s Quantum Theory of Atoms in Molecules (QTAIM)[110] has also served to find bond critical points between nearest-neighbor dimers and identify the nature of the dominant nonbonded contacts. Relying upon DORI, some of us examined the effectiveness of H-bonding and dispersion-driven side-chain aggregators to impose the tighter π-stacked arrangement in one-dimensional (1D) quaterthiophene nanowires.[111,112] It was also demonstrated that the electronic compactness, as probed by DORI, largely mirrors the computed charge transport properties in the derivatives of quaterthiophene crystal, 1D nanofibrils of quaterthiophene- and oligothienoacene-type cores.[113]

More controversial connections have been made between the “aromatic” character of dimerized trithiophene units and their hole mobility. However, these relationships, which suggest that decreased aromaticity (based on nucleus independent chemical shifts (NICS(1)) and harmonic oscillator model of aromaticity (HOMA) indices) is beneficial for transport,[114] suffer from the lack of unique quantitative assessments of fuzzy chemical concepts such as aromaticity.[115,116]

Despite extensive use of these quantum chemistry-based analytical tools and the valuable insights they provide, their main impact has so far been restricted to the rationalization of structure–energy–property trends.[86,87,109,113] Some have led to the formulation of concrete design principles,[84,85,94,95,111] albeit they remain yet to be fully exploited experimentally. Though design strategies based on the static dimer models are certainly insufficient to predict systems with better mobilities, some of the more rigorous and creative studies, discussed above, offer reliable opportunities to go beyond standard trial and error approaches. What is largely lacking at present is an efficient transfer of this unquestionably useful knowledge into an experimental domain.

Beyond dimers

Multiscale Approaches

The dimer approach became the workhorse of organic semiconductor quantum chemical modeling, habitually employed to either predict mobility or rationalize experimentally measured transport properties. However, it is inherently limited to only capture the elementary event of a charge hop between the neighbor molecules and thus neglects other critical factors, including the crystal packing, the presence of defects, the various disorder effects, etc. Multiscale approaches to charge transport deliver the remedy to these limitations.[60−62] They generally combine classical molecular dynamics of the semiconductor morphology with an estimation of the drift mobility within the high-temperature Marcus theory limit using, e.g., kinetic Monte Carlo (kMC) or master equation techniques. Flavors of multiscale approaches exist due to different methods and models used for computing the reorganization energies and electronic couplings and evaluating the disorder effects. Some studies also include crystal structure prediction (vide infra) or free energy estimation of the crystalline assemblies in their multiscale protocol. In any case, the multistage nature of these computations implies the need to use and interface several codes and methods. In this context, integrated computational workflows,[62,117,118] such as the Versatile Object-oriented Toolkit for Coarse-graining Applications (VOTCA),[119] significantly facilitate the multiscale modeling of microscopic transport (Chart ).

Chart 3

Workflow for Microscopic Simulations of Charge Transport Using VOTCAa

Reprinted with permission from ref (60). Copyright 2011 American Chemical Society.

Disorder Effects

When a molecular solid loses translational symmetry, both site energies (ε) and electronic couplings show a widespread static energetic and static positional disorder. Static disorder is naturally present in amorphous organic semiconductors and is detrimental to charge transport. In a perfectly ordered organic crystal (without chemical and structural defects), on the other hand, the only type of disorder is dynamic, arising from the thermal vibrations of the molecules.[120] Local and nonlocal electron–phonon couplings in the bulk structure, which can be extracted from multiscale simulations, allow for assessment of the extent of thermal fluctuations of the site energy ε and electronic coupling V.[121] Significant efforts have been carried out to develop fast and accurate methods for evaluating site energies.[122] In VOTCA, they are computed based on the microelectrostatic approach, in which correction terms (i.e., electrostatic and induction) are added to the isolated molecule’s ε.[60] Cornil et al. exploited constrained-DFT to compute ε from the energy difference between charged and neutral molecular clusters.[123] Alternatively, site energies can be computed with quantum mechanical/molecular mechanical (QM/MM) frameworks,[124] intramolecular charge redistribution,[125] quantum patch approach[126] or a mixed valence bond/Hartree–Fock (VB/HF) model.[127]

The important role of the dynamic energetic disorder in crystals has been illustrated by simulating hole mobilities of dicyanovinyl-substituted oligothiophene (DCVnT) crystals utilizing VOTCA.[128] The conformational disorder arising from thermal fluctuations of DCV-thiophene dihedral angles leads to substantial energetic disorder σE, resulting in a 1–2 orders of magnitude decrease in hole mobility. This reduction originates from interactions between the charge carriers and local electric fields induced by the fluctuating multipoles (conformational disorder). In addition, the presence of the energetic disorder can further lower the mobility by introducing energetic traps in the original (disorder-free) percolating network, as seen in the DCV3T and DCV4T cases. Computationally, the contributions to σE can be partitioned into unscreened Coulomb interactions and polarization with the latter causing a significant reduction. This, in turn, sheds light on the future rational design of organic semiconductors that feature large molecular polarizabilities to reduce σE.

Comparisons between computed and experimental mobilities for 22 π-conjugated cores (as crystalline and thin-film semiconductors) were performed based on a similar protocol.[129] Three different models were used to account for different degrees of dynamic disorder: (1) perfect crystal based on experimentally derived structures (zero disorder), (2) an MD-equilibrated structure with only dynamic positional disorder and (3) an MD-equilibrated structure with both disorders. A significant decrease in hole mobilities was observed after taking into account the energetic disorder (model (3) vs model (2)). Interestingly, the best agreement for experimental organic crystal mobilities was achieved using the perfect crystal model (rather than the MD-equilibrated morphology) and neglecting the σE. However, mobilities computed based on a more realistic model (3) are much lower than the experimental OFET ones. This fictitious agreement between experiment and the simplified perfect crystal model may be due to the discrepancy in charge carrier concentration between simulations and experiments. The experimentally determined OFET mobilities are usually measured at high carrier concentration, which is less sensitive to energetic disorder due to a trap-filling effect.[130] Yet, the kMC simulations are instead typically performed with only one charge carrier, coinciding with lower mobilities upon incorporation of σE.[128]

Given the detrimental effect the energetic disorder has on charge transport, it is crucial to account for it when optimizing a given semiconductor. Failing to do so may cause the desired improvements in material’s other properties to be negated by its amplified σE. For example, adding solubilizing groups, which is a common practice for improving the solution processability of the semiconductor, will also likely alter molecular packing and σE. Brédas, Coropceanu and co-workers demonstrated that, among various fullerene solubilizing adducts, indene-based adducts in general lead to smaller dynamic energetic disorders as compared to butyric-acid-methyl-ester ones.[131]

In addition to energetic disorder, the thermal fluctuations of the electronic couplings (dynamic positional disorder) have a large impact on the charge transport. Combining transmission electron microscopy (TEM) and MD simulations, Illig et al. showed that substitution of the side chains along the long axis of a conjugated core can be an effective strategy for reducing intermolecular thermal vibrations, leading to a decreased dynamic positional disorder. The small positional disorder present in the crystal can explain the outstanding charge mobility of C8-BTBT and C10-DNTT.[132] Dynamic positional disorder can sometimes enhance mobility, as seen in fluorinated perylene bisimide organic semiconductors.[133] In this case, the optimal molecular packing in the crystal, induced by intermolecular interactions, leads to minimal electronic coupling (Figure ). Therefore, deviations from the equilibrated packing configuration result in larger electronic couplings and thereby higher charge mobility. The importance of the dynamic positional disorder recently prompted Landi and Troisi to develop a computational methodology for the fast evaluation of nonlocal electron–phonon couplings. Their method affords excellent agreement with the results of more accurate computations and allows the screening of large databases for promising organic semiconductors with optimal molecular arrangements.[134]

Figure 2

Dependence of the electronic coupling on the relative transversal shift between the cores of alkyl-substituted perylene bisimide starting from the perfectly π-stacked (red) and crystal (blue) dimer geometries; intermolecular distance is kept equal to that in a crystal. Reprinted with permission from ref (133). Copyright 2010 American Chemical Society.

Although useful insights can be obtained from a multiscale simulation based on pure single crystal morphology, the computed charge mobility is usually much higher than that of a polycrystalline organic thin film. Charge hopping processes across grain boundaries serve as limiting steps of the overall charge transport process despite a relatively fast intragrain charge transport. Practically, a calibration can be applied to the simulated single-crystal mobility to account for the difference between the single crystal and the realistic polycrystalline morphology. On the basis of a one-dimensional microstructure model, Shin et al. derived a calibration equation with the unknown coefficients fitted by μcal/μexp – grain size plot for 4 different organic semiconductors.[135] The transferability of the calibration equation is further examined by comparing experimental charge mobilities with the calibrated single-crystal ones for 17 organic semiconductors. A large mismatch (1–2 orders of magnitude) is observed between experimental and calibrated mobilities for several molecules, which implies that the detailed microscopic information is critical, but is missing in the simple calibration model. Nelson et al. further investigated the energetic profile across the grain boundary and the dependence of the charge transport on the relative orientations of the grains.[136] The disruption of the crystalline packing at the grain boundary results in energetic barriers for holes and potential wells for electrons, which lower the charge mobility by ∼2 orders of magnitude. As the contact length between the grains shortens, the charge mobility becomes correlated to the relative orientation of the two grains, where the charge transport is confined to a high-energy pathway. Apart from the grain boundaries, crystal step edges can hamper electron transport in n-type organic semiconductors.[137] They have a relatively positive surface potential as compared to the flat surface, which can form in-gap states that trap electrons. To our knowledge, the effect of the crystal step edges has not yet been investigated in any previous work utilizing a multiscale approach, presumably due to the limitation of the simulated system size.

For amorphous organic semiconductors, the static disorder becomes dominant over the dynamic one,[138] thus the simulated mobility, based on a crystal structure, overestimates the real charge mobility. Several computational protocols have been utilized to construct organic amorphous morphologies relying upon classical MD or Monte Carlo simulations.[139] Owing to the enormous conformational disorder characteristic of the amorphous morphology, it is necessary to account for the intrinsic energetic disorder (σi) resulting from a spread of HOMO energy levels, induced by various molecular conformations, in addition to the contributions from the Coulomb interactions and polarization (σp).[140] The former can be evaluated by the HOMO energy variation with respect to the conformational change, while the latter shows a correlation with the molecular dipole moment and intermolecular distance.[138] From an a priori molecular design perspective, it is essential to estimate the extent of energetic disorder on the charge mobility.[141−143] Thus, single molecule properties (HOMO energy, dipole moment, etc.) are likely to represent a helpful screening criterion to rule out unfavorable candidates.

Packing Effects

The charge mobility of an organic crystal is usually anisotropic. For this reason, accessing the full mobility tensor is crucial to align the fastest transport pathway with the charge transport direction (electric field direction). The anisotropic mobility can be reflected in a polar plot (Figure ),[144] where kMC simulations are performed with the electric fields applied in various directions. Even at the pairwise (dimer) level, this type of analysis is useful to rationalize the transport properties of different polymorphs. For example, among the three polymorphs of a dipyrrolyldiketone difluoroboron complex, two were shown to feature a one-dimensional transport along the π-stacks (hampered by the localized trap sites in one of them), while the third instead involves a three-dimensional transport.[145] Multiscale simulations also provide an insight into the relationship between molecular structure, morphology, percolation network and charge carrier mobility.[146] This protocol involves computing the electronic coupling between all neighbor molecules within a crystal structure snapshot. But in addition to simply evaluating the charge transfer rates using the Marcus theory, this information can also serve to visualize the charge transport dimensionality and directionality and the impact of the disorder on the transport topology via the so-called connectivity graphs (Figure A). In these graphs, the electronic couplings above certain threshold value are shown as “bonds”. This work illustrated that, despite the common perception, a perfectly π-stacked arrangement is not necessarily the one with the highest mobility. Instead, a shifted cofacial alignment is a better alternative as it allows tighter packing (decreasing the hop distance) and a two-dimensional transport that reduces the influence of defects. This beneficial packing can be enforced by the attachment of side chains perpendicular to the conjugated core, as in the case of rubrene.[84] Despite the usefulness of connectivity graphs based on electronic couplings, their direct comparison for different organic crystals does not provide information concerning relative charge transport performance. Since the threshold of the electronic coupling V used in the graph is arbitrary and high V alone does not guarantee fast transport, replacing this parameter with a more representative quantity would be very advantageous. Recently, we constructed connectivity graphs based on a single-hop mobility μsh between each dimer (Figure B).[147] To construct the graphs, it is necessary to identify the maximum single-hop mobility μsh,max, then plot all μsh within a certain factor of μsh,max as ‘bonds’. The resulting graphs based on this quantitative reference provide in-depth topological transport behavior and are particularly beneficial for systematically investigating and comparing systems within a database. Besides their application to crystals, topological connectivity plots based on electronic couplings have been constructed for DCVnT of amorphous and smectic mesophases.[148] In addition to the electronic coupling, spatial information on the site energy helps rationalizing charge transport in disordered phases. This information can be represented by an energetic color map (Figure ) showing the energetic spatial correlation induced by large molecular dipole moments. Important molecular pairs that contribute significantly to the total current can be further identified by the edge current analysis. The significantly lower charge mobility of a more ordered smectic phase of DCV6T compared to that of its amorphous counterpart can only be explained using this three-dimensional information.

Figure 3

Polar plots of hole (left, filled circles) and electron (right, empty circles) anisotropic mobility, computed for different disorder conditions (black and blue). Reprinted with permission from ref (144). Copyright 2018 American Chemical Society.

Figure 4

(A) Nearest-neighbor alignment and corresponding percolation (connectivity) network of electronic couplings in representative semiconductor crystals. Reprinted with permission from ref (146). Copyright 2010 American Chemical Society. (B) Transport patterns of representative molecular semiconductors, in which the topological connectivity is based on single-hop mobilities within a factor of 2 (TC-2) and 10 (TC-10) of μsh,max. Adapted from ref (147).

Figure 5

Equilibrated simulation boxes (a) and topology of hole hopping in amorphous and smectic DCV6T, depicted via connectivity plots (b) and color maps (c). Reproduced with permission from ref (148). Copyright 2012 Royal Society of Chemistry.

Overall, multiscale approaches have been used to (i) reproduce experimental charge mobilities, (ii) explain their trends from a molecular perspective and (iii) guide molecular design.[128,129,149] Akin to the preceding section, the discovery of new high-performance molecular semiconductors based on these approaches is still scarce. One major obstacle is the inaccessibility of reasonable crystal morphologies of new organic molecules, highlighting the crucial role of crystal structure prediction, discussed in next section. In addition to bulk morphology, various interfacial effects are likely to impact the mobility. On one hand, within any electronic device the charge is transported across numerous interfaces and the resulting device mobility is a combination of transport within (including the organic semiconducting ones) and between the layers. This is addressed by complete device simulations[150−154] and is outside of the direct scope of this Perspective. On the other hand, the interface can also affect the intrinsic charge transport characteristics of a semiconducting layer itself in a number of ways. For example, in organic thin film transistors (OTFTs), charge transport occurs primarily in only the several molecular layers neighboring the organic semiconductor/dielectric interface.[155] Therefore, molecular alignment near the interface is crucial to the measured charge mobility of OTFTs. Different experimental techniques have been developed to modify the interfacial properties in order to alter the molecular packing and achieve optimal charge transport properties.[156,157] Apart from the interfacial molecular alignment, it has been found that the energetic disorder of the organic semiconductor layer in an OTFT device can be enhanced due to static dipolar disorder in the adjacent dielectric layer to an extent, proportional to the dielectric constant of the latter.[158,159] Thus, evaluating energetic disorder on the basis of bulk morphology only will likely lead to its underestimation. In silico prediction of such interfacial effects on site energy can, in principle, be achieved if the constructed atomistic model is a good estimate of a real morphology of the interface, which is still challenging for amorphous systems.[160]

Crystal Structure Prediction

The design of novel organic semiconductors often stumbles upon the lack of experimental crystal structures. Molecular crystal structure prediction has been, for a long time, an insurmountable challenge. However, enormous research efforts to resolve it have been slowly but surely coming to fruition in the last 10 years, as evidenced by the growing successes in the increasingly more complex CSP Blind Tests, set up by the Cambridge Crystallographic Data Centre.[67,161,162] Though the field has been historically driven by the intellectual and financial involvement from the pharmaceutical industry,[163] it is nowadays reaching the organic semiconductors domain. An early CSP effort in this area included using experimental crystal structures of reasonably similar molecular cores as starting points for force field optimizations of the candidate systems. For example, Aspuru-Guzik et al. employed this approach to predict crystal structures of the derivatives of dinaphtho[2,3-b:2′,3′-f]thieno[3,2-b]thiophene (DNTT),[164] which were identified as promising candidates due to their low reorganization energies.[165] The compound with the best predicted charge transport characteristics was then tested experimentally and found to afford an impressive mobility of up to 16 cm–2 V–1 s–1. A comparison between the predicted and experimental crystal structures and resulting mobilities indicates that despite the quantitative disagreements (e.g., mobility of 3.3 cm–2 V–1 s–1 was predicted for the aforementioned compound), the computational modeling does capture the relative charge transport trends. A more sophisticated approach to crystal structure prediction of organic semiconductors involving global exploration of the lattice energy surface (i.e., the potential lattice energy landscape depending on the crystal density) was employed to evaluate the effect of small chemical changes on crystal packing and charge mobility in a set of azapentacenes.[166] A total of 212 000 trial crystal structures were generated using the Global Lattice Energy Explorer software and their lattice energies were minimized using a model potential (i.e., W99 exp-6). The resulting crystal structures were categorized according to their packing types and used to compute the charge mobility within the Marcus theory framework. The results were then combined into the energy-structure–function (ESF) maps to identify which systems offer the best combination of structural stability and charge mobility. A large spread of electron mobility is observed among the low energy crystal structures for all studied molecules due to their diverse packing motifs. Crucially, for none of the studied molecules does the crystal structure with the highest charge mobility correspond to the global lattice energy minimum. This dramatic relationship (more precisely, a lack thereof) is due to the fact that both the destabilizing (exchange) contribution to the interaction energy and the high electronic coupling are greater in systems with better spatial overlap between the cores,[69,94,167] e.g., the perfectly π-stacked arrangement. However, when a similar computational protocol combining crystal structure prediction and energy–structure–function map was applied to a potential chiral organic semiconductor, [6]helicene (Figure ), many low-energy structures were found to have some of the highest mobilities.[6]Helicene Molecule. Nanoscale. 2018 ">168] Though none of these structures is the lowest-energy polymorph, they are all within the thermodynamically accessible range of 1.1–4.6 kJ mol–1 from it. Interestingly, neither the energetically favorable nor the high charge-carrier mobility packing motifs of [6]helicene were intuitively obvious, illustrating the potential of such screening techniques for the development of novel molecular semiconductors. Prototypical n-type semiconductor, N,N′-ditridecylperylene-3,4,9,10-tetracarboxylic diimide, exemplifies yet another type of energy-mobility relationship. For this system, a combination of MD and metadynamics served to map the free energy profile of various morphologies (z-component of the distance between the terminal carbon atoms of contiguous alkyl chains was chosen as the collective variable) and access kinetic information.[169] It was demonstrated that the global energetic minimum configuration also affords some of the highest charge mobility among all considered structures.

Figure 6

Lattice energy landscape (A), labeled by the corresponding molecular substructures of [6]helicene in the order of their decreasing frequency (B). Reproduced with permisison from ref (168). Copyright 2018 Royal Society of Chemistry.

When advantageous packing motifs are unachievable under typical fabrication conditions, external strain can be applied to induce molecular reorganization and a shift into nonequilibrium morphology. In this context, multiscale simulations were employed to predict the nonequilibrium structures of TIPS-pentacene thin films, forming under strain, and evaluate the associated mobility.[170] A combination of the shear and tensile lattice strain was found not only to significantly improve the charge mobility, but also to decrease its anisotropy.

The exploration of the crystal energy landscapes, mentioned above, generally involved simple force-field models for describing the intermolecular forces (model potentials, distributed multipole) and, eventually, biased MD simulations. Considering the tremendous efforts placed in developing the next generation of “physics-based” machine learning models for reproducing both molecular energies[14,15] and noncovalent interactions[171,172] or lattice energies,[171,173] we can anticipate that the prediction of entire crystal energy landscapes will be dramatically accelerated in the near future.

Data-Driven Searches

Not only the morphology can be predicted by atomic-scale modeling or increasingly popular machine learning algorithms. Large scale screening strategies can be employed to predict the semiconducting properties of materials from various easily computable descriptors. In a 2013 study, Manuel et al. constructed a database made of organic compounds (77 molecules in the training set and 19 in the test set) and exploited four ML methods that revealed the conducting properties of Schiff base molecules.[174] About 1500 chemoinformatics descriptors were used to train the ML models to distinguish semiconductors and nonsemiconductors.

Using a “quantum” machine learning perspective (as opposed to chemoinformatics),[14] Ceriotti, Day and co-workers elegantly combined the accelerated, albeit accurate predictions, of lattice energies and electronic coupling values of polymorphs with a data-driven classification to explore the structure–energy–property relationships of organic semiconductors.[173] The predicted crystal structures of pentacene and azapentacenes (including crystal structures up to an energetic cutoff above the global minimum) were taken from ref (166) and classified using the SOAP structural similarity kernel[175] and clustering techniques.[176] Gaussian Processes Regression machine learning models based on the same kernel were employed to predict stability- and mobility-related properties. This modern and efficient framework offered an appealing and unbiased classification of the predicted crystal structures combined with dramatically accelerated computations of the electronic couplings that enter the evaluation of charge mobility. Ultimately, it yielded intuitive insights into relationships between the packing motifs in polymorphs, types of intermolecular interactions and resulting charge mobility in molecular semiconductors. This approach was also used to screen a set of 28 structural isomers of a recently reported promising candidate for molecular electronics, pyrido[2,3-b]pyrido-[3′,2′:4,5]pyrrolo[3,2-g]indole.[177] As before, ML allowed rapid classification of large numbers of predicted crystal structures by their packing motifs, construction of a multilandscape sketch-map clustering structural patterns with predicted energetic and transport characteristics, and identification of the promising candidate molecules (Figure ).

Figure 7

Sketch-map representations of the low-energy crystal structure landscapes for the 28 studied isomers of pyrido[2,3-b]pyrido-[3′,2′:4,5]pyrrolo[3,2-g]indole. The sketch-maps are color-coded according to (a) clusters detected by the clustering method; (b) conventional packing motifs; (c) different isomers; (d) calculated lattice energies (kJ/mol) and (e) predicted electron mobilities (in cm2 V–1s–1). Reprinted with permission from ref (177). Copyright 2018 American Chemical Society.

Automated screening workflows are alternative high-throughput strategies, which do not rely on statistical learning. Oberhofer et al. developed an extensive high-throughput workflow to compute the charge transport parameters of 95 445 molecular crystals, extracted from the Cambridge Structural Database (CSD).[178] Electronic couplings were estimated using the fragment molecular orbital approach (FO-DFT), whereas intramolecular reorganization energies were obtained from nonperiodic QM/MM computations (Figure ). Such a large database of descriptors was aimed at future materials discovery through enabling in-depth exploration of the relationships between the chemical nature of molecular building blocks, the crystal structures and percolation pathways. When the charge transport parameters are available, mobility can be computed via kinetic Monte Carlo simulations, for which van der Kaap and Koster have also developed a high-throughput approach.[179] Their massively parallelized lattice-based kMC method includes Coulombic particle–particle interactions and runs on general-purpose graphic processing units. Comparison with the mobilities, obtained by numerically solving the corresponding master equation, validated the accuracy of this new method.

Figure 8

High-throughput workflow, developed in ref (178). Reprinted with permission from ref (178). Copyright 2016 American Chemical Society.

On a much smaller scale, DFT (i.e., PBE-TS) screening performed on the crystal structures of all 91 polyaromatic hydrocarbons (PAH), available in CSD,[180] demonstrated the existence of a limit between the maximum optical gap and the intermolecular cohesive energy or the C···C nonbonded contacts. Effort was spent in bypassing the demanding optical gap computations in gas, solution and crystalline phase via the prediction of the Kohn–Sham gaps.

There is no doubt that the aforementioned advancements in generating and characterizing numerous (up to the order of thousands) crystal structures, as well as a rapidly increasing number of available experimental morphologies (e.g., the Cambridge Structural Database, CSD, and the recent Organic Materials Database, OMDB[181]) will coincide with the development of additional cost-efficient electronic structure approaches or composite electronic structure-ML methods,[182,183] offering dramatic speed-up in computations of the related charge transport parameters[173,184] and energy landscapes.[15,171,185]

Conclusions and Outlook

Organic semiconductors are certainly among the key players in the arena of functional materials, performing as well, and sometimes even better, than their conventional inorganic counterparts. These successes, however, are conditional upon resolving diverse challenges at the atomic, molecular and morphological length scales. This poses a formidable yet exciting task for the computational community at large, including chemists, physicists, materials and data scientists. In response, various methods and models have been developed and applied to quantify and rationalize the properties of organic electronic materials.

Computational assessment of the key charge transport parameters, electronic and local electron–phonon coupling (e.g., electronic coupling and reorganization energy), became a routine exercise, nowadays accompanying many, if not most experimental papers reporting new molecular semiconductors. However, the output from computations is seemingly hampered in its move beyond this useful, but rather limited end result. Elegant and insightful quantum chemical approaches toward the intricate details of charge transport on a molecular/pairwise level have so far succeeded in explaining the known, but not predicting unknown organic transport materials.

An elephant in their room is, of course, a complete lack of recognition of the semiconductor morphology beyond several representative dimers. Dreaming big, we can hypothesize that some of the bulk features, such as disorder or lattice energy, might somehow be traced back to the basic molecular and intermolecular features, which would certainly facilitate and conceptualize the design and screening of organic semiconductors. This has not been achieved so far, necessitating the development and use of more sophisticated multiscale approaches, which combine simulations of the bulk morphology with the modeling of the charge carrier dynamics and, in this way, provide a considerably more realistic description of charge transport.

Still, in their core, these approaches rely upon the availability of materials morphology, which might arguably be the true stepping stone for computational modeling of organic semiconductors. The validity and significance of the most general proposed design principles would take on another dimension if they could be readily tested on realistic (yet unavailable for newly designed systems) morphologies. Thus, crystal structure (and, more generally, morphology) prediction remains the Holy Grail for the in silico organic semiconductors community, which can hopefully build upon the successes of the pharmaceuticals-oriented research in this domain. This road can now be traveled quicker with the help of rapidly emerging data-driven approaches, e.g., modern machine learning techniques, which can dramatically expedite the evaluation of the energetic and transport characteristics of the plentiful, both computer-generated or available experimentally, crystal structures.

Computations aside, for experimentalists their Holy Grail lies in the precise control with on-demand manipulation of the morphology of both the bulk layers and their interfaces within the electronic device. Ultimately, the bright future of organic semiconductors hinges upon the ability of the computational community to move from rationalizing to predicting materials and the willingness of experimental colleagues to take the resulting design principles on board.