Microscopic Simulations of Charge Transport in Disordered Organic Semiconductors.

Charge carrier dynamics in an organic semiconductor can often be described in terms of charge hopping between localized states. The hopping rates depend on electronic coupling elements, reorganization energies, and driving forces, which vary as a function of position and orientation of the molecules. The exact evaluation of these contributions in a molecular assembly is computationally prohibitive. Various, often semiempirical, approximations are employed instead. In this work, we review some of these approaches and introduce a software toolkit which implements them. The purpose of the toolkit is to simplify the workflow for charge transport simulations, provide a uniform error control for the methods and a flexible platform for their development, and eventually allow in silico prescreening of organic semiconductors for specific applications. All implemented methods are illustrated by studying charge transport in amorphous films of tris-(8-hydroxyquinoline)aluminum, a common organic semiconductor.

Introduction

The progress currently observed in the field of organic electronics is a result of a combined effort of several communities. Synthetic chemists have identified classes of promising compounds, ranging from small conjugated molecules to self-assembling oligomers and conjugated polymers, and developed new synthetic routes, improving both stability and processability of the materials.[1−7] At the same time, material processing, such as doping, annealing, use of a secondary solvent, and composition tuning, has been adjusted to the demands of the field.[8−13] In parallel, increased device efficiencies could be achieved, e.g., by optimizing light in- and out-coupling and introducing tandem concepts.[14,15]

Compound design requires an in-depth understanding of elementary processes occurring in organic semiconductors.(16) In particular, linking the chemical structure to charge dynamics is a nontrivial task, since several factors can influence charge carrier mobility: the molecular electronic structure, the relative positions and orientations of neighboring molecules, and spatial inhomogeneities in the morphology, which determine charge carrier pathways on a macroscopic scale.(17)

Furthermore, the choice of the model Hamiltonian depends on the specific situation.(18) For perfectly ordered defect-free crystals at low temperatures, the Drude model based on band theory(19) or its extensions, which account for local electron–phonon coupling,[20−22] are often used. At ambient conditions, however, the thermal fluctuations of the transfer integral, i.e., the nonlocal electron–phonon coupling, are on the same order of magnitude as its average value, and charge transport should be treated as diffusion limited by thermal disorder. This can be achieved using semiclassical dynamics based on a Hamiltonian with interacting electronic and nuclear degrees of freedom.[23−25] If nuclear dynamics are much slower than the dynamics of charge carriers and electronic coupling is weak, charge transport can be described by a Hamiltonian with static disorder, based on simple assumptions on the electronic density of states and on the hopping rates between localized states.

The latter approach is by now routinely used to study charge transport in amorphous and partially disordered small-molecule-based organic semiconductors.[26−40] Its key ingredients are material morphology and intermolecular charge transfer (hopping) rates.(41) The rates not only depend on the molecular electronic structure but are also sensitive to the mutual positions and orientations of molecules. Hence, in order to evaluate the rates, the material morphology must be known at an atomistic resolution. This can be achieved by performing molecular dynamics simulations and thus relies on force-field development for new compounds. If the required time and length scales exceed the range available to atomistic molecular dynamics, coarse-graining techniques can be used.(42) These techniques need to be capable of back-mapping the coarse-grained representation to an atomistic resolution.

Charge transfer rates can be postulated on the basis of intuitive physical considerations, as is done in the Gaussian disorder models.[43−46] Alternatively, charge transfer theories can be used to evaluate rates from quantum chemical calculations.[28,47−51] In spite of being significantly more computationally demanding, the latter approach allows one to link the chemical and electronic structure, as well as the morphology, to charge dynamics.

The high temperature limit of classical charge transfer theory[52,53] is often used as a trade-off between theoretical rigor and computational complexity. It captures key parameters which influence charge transport while at the same time providing an analytical expression for the rate. Within this limit, the transfer rate for a charge to hop from a site i to a site j reads:where T is the temperature, λ = λint + λout is the reorganization energy, which is a sum of intra- and intermolecular (outer-sphere) contributions, ΔE is the site-energy difference, or driving force, and J is the electronic coupling element, or transfer integral.(54) A more general, quantum-classical expression for a bimolecular multichannel rate is derived in the Supporting Information.

All of the ingredients entering eq 1 can be calculated using electronic structure techniques, classical simulation methods, or their combination. With the rates at hand, one can study charge transport by solving the differential (master) equation, e.g., by using the kinetic Monte Carlo method, which is capable of simulating charge dynamics of non-steady-state systems.

Altogether, the task of charge transport characterization is rather tedious and time-consuming to perform, even for a single compound. The main aim of this work is to introduce a software package which implements a set of techniques for charge transport simulations as well as provides a flexible modular platform for their further development.

The paper is organized as follows. In section , we describe the workflow and the basic ideas behind each method. As an illustration, we study charge transport in the bulk of amorphous tris-(8-hydroxyquinoline)aluminum (Alq3). Section deals with the analysis and visualization of charge dynamics. A brief summary of implementation is given in section .

Methods

A workflow of charge transport simulations is depicted in Figure 1. The first step is the simulation of an atomistic morphology (section ), which is then partitioned into hopping sites (section ). The coordinates of the hopping sites are used to construct a list of pairs of molecules (neighbor list). For each pair, an electronic coupling element (section ), a reorganization energy (section ), a driving force (section ), and eventually the hopping rate are evaluated. The neighbor list and hopping rates define a directed graph. The corresponding master equation is solved using the kinetic Monte Carlo method (section ), which allows one to explicitly monitor the charge dynamics in the system as well as calculate time or ensemble averages of occupation probabilities, charge fluxes, correlation functions, and field-dependent mobilities (section ).

Figure 1

Workflow for microscopic simulations of charge transport. Ground state geometries, partial charges, and a refined force field are used to simulate atomistically resolved morphologies (section ). After partitioning into conjugated segments and rigid fragments (section ), a list of pairs of molecules (neighbor list) is constructed. Transfer integrals (section ), reorganization (section ) and site energies (section ), and eventually hopping rates are calculated for all pairs from this list. A directed graph is then generated, and the corresponding master equation is solved using the kinetic Monte Carlo method (section ).

Material Morphology

There is no generic recipe on how to predict a large-scale atomistically resolved morphology of an organic semiconductor. The required methods are system-specific: for ultrapure crystals, for example, density-functional methods can be used provided the crystal structure is known from experimental results. For partially disordered organic semiconductors, however, system sizes much larger than a unit cell are required. Classical molecular dynamics or Monte Carlo techniques are then the methods of choice.

In molecular dynamics, atoms are represented by point masses which interact via empirical potentials prescribed by a force field. Force fields are parametrized for a limited set of compounds, and their refinement is often required for new molecules. In particular, special attention shall be paid to torsion potentials between successive repeat units of conjugated polymers or between functional groups and the π-conjugated system. First-principles methods can be used to characterize the missing terms of the potential energy function. The parametrization must take into account existing force-field contributions, e.g., due to nonbonded interactions. If q is the degree of freedom of interest, constrained geometry optimizations must be performed using both first principles and the force-field levels, yielding the total energies Ufp(q) and Uff(q), respectively. Then, the missing force-field terms are fitted to their difference, Ufp(q) – Uff(q), using a prescribed functional form. For several identical or coupled degrees of freedom, a multidimensional fit can be used. Note that force-field validation is as important as its refinement. For instance, X-ray scattering and solid-state NMR provide information about averaged molecular arrangements to which simulation results can be compared.

As an example, the refinement of the OPLS force field for Alq3 is described in the Supporting Information. In total, 16 bonded interactions were parametrized. To validate the force field, the glass transition temperature and bulk density of amorphous Alq3 were compared to the experimental values. An amorphous morphology of Alq3 was then obtained by quenching the system after equilibrating it above the glass transition temperature.

Self-assembling materials, such as soluble oligomers, discotic liquid crystals, block copolymers, partially crystalline polymers, etc., are the most complicated to study. The morphology of such systems often has several characteristic length scales and can be kinetically arrested in a thermodynamically nonequilibrium state. For such systems, the time and length scales of atomistic simulations might be insufficient to equilibrate or sample desired morphologies. In this case, systematic coarse-graining can be used to enhance sampling.(42) Note that the coarse-grained representation must reflect the structure of the atomistic system and allow for back-mapping to the atomistic resolution.

Conjugated Segments and Rigid Fragments

With the morphology at hand, the next step is the construction of the effective electronic Hamiltonian of the system. In a static disorder approximation, this is equivalent to partitioning the system into hopping sites, or conjugated segments, and calculating charge transfer rates between them. Physically intuitive arguments can be used for the partitioning, which reflects the localization of the wave function of a charge. For most organic semiconductors, the molecular architecture includes relatively rigid, planar π-conjugated systems, which we will refer to as rigid fragments. A conjugated segment can contain one or more such rigid fragments, which are linked by bonded degrees of freedom. The dynamics of these degrees of freedom evolves on time scales much slower than the frequency of the internal promoting mode. In some cases, e.g., glasses, it can be “frozen” due to nonbonded interactions with the surrounding molecules.

To illustrate the concept of conjugated segments and rigid fragments, three representative molecular architectures are shown in Figure 2. The first one is a typical discotic liquid crystal, hexabenzocoronene. It consists of a conjugated core to which side chains are attached to aid self-assembly and solution processing. In this case, the orbitals localized on side chains do not participate in charge transport, and the π-conjugated system is both a rigid fragment and a conjugated segment.

Figure 2

The concept of conjugated segments and rigid fragments. Dashed lines indicate conjugated segments, while colors denote rigid fragments. (a) Hexabenzocoronene: the π-conjugated system is both a rigid fragment and a conjugated segment. (b) Alq3: the Al atom and each ligand are rigid fragments, while the whole molecule is a conjugated segment. (c) Polythiophene: each repeat unit is a rigid fragment. A conjugated segment consists of one or more rigid fragments. One molecule can have several conjugated segments.

In Alq3, a metal-coordinated compound, a charge carrier is delocalized over all three ligands. Hence, the whole molecule is one conjugated segment. Individual ligands are relatively rigid, while energies on the order of kBT are sufficient to reorient them with respect to each other. Thus, the Al atom and the three ligands are rigid fragments.

In the case of a conjugated polymer, one molecule can consist of several conjugated segments, while each backbone repeat unit is a rigid fragment. Since the conjugation along the backbone can be broken due to large out-of-plane twists between two repeat units, an empirical criterion, based on the dihedral angle, can be used to partition the backbone on conjugated segments.(36) However, such intuitive partitioning is, to some extent, arbitrary and shall be validated by other methods.[55−57]

After partitioning, an additional step is often required to remove bond length fluctuations introduced by molecular dynamics simulations, since they are already integrated out in the derivation of the rate expression. This is achieved by substituting respective molecular fragments with rigid, planar π systems optimized using first-principles methods. Centers of mass and gyration tensors are used to align rigid fragments, though a custom definition of local axes is also possible. Such a procedure also minimizes discrepancies between the force-field and first-principles-based ground state geometries of conjugated segments, which might be important for calculations of electronic couplings, reorganization energies, and intramolecular driving forces.

Finally, a list of neighboring conjugated segments is constructed. Two segments are added to this list if the distance between centers of mass of any of their rigid fragments is below a certain cutoff. This allows neighbors to be selected on a criterion of minimum distance of approach rather than center of mass distance, which is useful for molecules with anisotropic shapes.

Transfer Integrals

The electronic transfer integral J entering the Marcus rates in eq 1 is defined aswhere ϕ and ϕ are diabatic wave functions, localized on molecules i and j, respectively, participating in the charge transfer, and Ĥ is the Hamiltonian of the formed dimer. Within the frozen-core approximation, the usual choice for the diabatic wave functions ϕ is the highest occupied molecular orbital (HOMO) in the case of hole transfer and the lowest unoccupied molecular orbital (LUMO) in the case of electron transfer, while Ĥ is an effective single particle Hamiltonian, e.g., a Fock or Kohn–Sham operator of the dimer. As such, J is a measure of the strength of the electronic coupling of the frontier orbitals of monomers mediated by the dimer interactions. Intrinsically, the transfer integral is very sensitive to the molecular arrangement, i.e., the distance and the mutual orientation of the molecules participating in charge transport. Since this arrangement can also be significantly influenced by static and/or dynamic disorder,[24,31,34,35,53] it is essential to calculate J explicitly for each hopping pair within a realistic morphology. Considering that the number of dimers for which eq 2 has to be evaluated is proportional to the number of molecules times their coordination number, computationally efficient and at the same time quantitatively reliable schemes are required.

In general, information about three objects is needed: the two monomer wave functions ϕ and ϕ and the dimer interaction Hamiltonian Ĥ. An approximate method based on Zerner’s independent neglect of differential overlap (ZINDO) has been described in ref (51). This semiempirical method is substantially faster than first-principles approaches, since it avoids the self-consistent calculations on each individual monomer and dimer. This allows one to construct the matrix elements of the ZINDO Hamiltonian of the dimer from the weighted overlap of molecular orbitals of the two monomers. Together with the introduction of rigid segments, only a single self-consistent calculation on one isolated conjugated segment is required. All relevant molecular overlaps can then be constructed from the obtained molecular orbitals. This molecular orbital overlap (MOO) method has been applied successfully to study charge transport, for instance, in discotic liquid crystals,[31,37,38] polymers,(36) or partially disordered organic crystals.[33−35]

While the use of the semiempirical ZINDO method provides an efficient on-the-fly technique to determine electronic coupling elements, it is not generally applicable to all systems. For instance, its predictive capacity with regards to atomic composition and localization behavior of orbitals within more complex structures is reduced. Moreover, transition or semimetals are often not even parametrized. In this case, ab initio based approaches, e.g., density-functional theory, can remedy the situation.[50,58−62] Within the dimer projection method described in detail in ref (50), explicit quantum-chemical calculations are required for every molecule and every hopping pair in the morphology. As a consequence, this procedure is significantly more computationally demanding. The code currently contains scripts which support an evaluation of transfer integrals from quantum-chemical calculations performed with the GAUSSIAN and TURBOMOLE packages.

As an example, distributions of transfer integrals calculated using ZINDO and DFT (with the gradient-corrected B–P functional[63,64] and a TZVP basis set(65)) methods are shown in Figure 3. While both distributions are similar, ZINDO integrals are, on average, smaller than DFT ones.

Figure 3

Distributions and correlation of transfer integrals calculated using ZINDO and DFT methods.

Reorganization Energy

The reorganization energy λ takes into account the change in nuclear (and dielectric) degrees of freedom as the charge moves from donor i to acceptor j. It has two contributions: intramolecular, λint, which is due to a reorganization of nuclear coordinates of the two molecules forming the charge transfer complex, and intermolecular (outersphere), λout, which is due to the relaxation of the environment. In what follows, we discuss how these contributions can be calculated.

Intramolecular Reorganization Energy

If intramolecular vibrational modes of the two molecules are treated classically, the rearrangement of their nuclear coordinates after charge transfer results in the dissipation of the internal reorganization energy, λint. It can be computed from four points on the potential energy surfaces (PES) of both molecules in neutral and charged states, as indicated in Figure 4. Adding the contributions due to discharging of molecule i and charging of molecule j yields(49)Here, U is the internal energy of the neutral molecule i in the geometry of its charged state (small n denotes the state and capital C the geometry). Similarly, U is the energy of the charged molecule j in the geometry of its neutral state.(66) Note that the PESs of the donor and acceptor are not identical for chemically different compounds or for conformers of the same molecule. In this case, λ ≠ λ and λ ≠ λ. Thus, λint is a property of the charge transfer complex and not of a single molecule.

Figure 4

Potential energy surfaces of (a) donor and (b) acceptor in charged and neutral states. After the change of the charge state, both molecules relax their nuclear coordinates. If all vibrational modes are treated classically, the total internal reorganization energy and the internal energy difference of the electron transfer reaction are λint = λ + λ and ΔEint = ΔU – ΔU, respectively.

In Alq3, the three ligands can easily change their mutual orientations. Molecular conformations are then “frozen” due to nonbonded interactions in an amorphous glass. The internal energies entering eq 3 were calculated after optimizing molecular geometries of all 512 molecules in charged and neutral states with the soft degrees of freedom constrained to their average values (see the Supporting Information for details). The distribution of λint for holes, which is shown in the Supporting Information, is sharply peaked with a maximum at 0.21 eV and variance of 0.03 eV. Computing λint from the PES of two unconstrained molecules leads to a similar value of 0.23 eV. Since Alq3 has high energetic disorder arising from its large dipole moment, this small variance in reorganization energy does not affect charge carrier mobility or Poole–Frenkel behavior.

Outersphere Reorganization Energy

During the charge transfer reaction, also the molecules outside the charge transfer complex reorient and polarize in order to adjust for changes in electric potential, resulting in the outersphere contribution to the reorganization energy. λout is particularly important if charge transfer occurs in a polarizable environment. Assuming that charge transfer is much slower than electronic polarization but much faster than nuclear rearrangement of the environment, λout can be calculated from the electric displacement fields created by the charge transfer complex:(67)where ε0 is the permittivity of free space, I,F() are the electric displacement fields created by the charge transfer complex in the initial (charge on molecule i) and final (charge transferred to molecule j) states, Vout is the volume outside the complex, and cp = 1/εopt – 1/εs is the Pekar factor, which is determined by the low (εs) and high (εopt) frequency dielectric permittivities.

Equation 4 can be simplified by assuming spherically symmetric charge distributions on molecules i and j with total charge e. Integration over the volume Vout outside of the two spheres of radii R and R centered on molecules i and j leads to the classical Marcus expression for the outersphere reorganization energy:where r is the molecular separation. While eq 5 captures the main physics, e.g., predicts smaller outersphere reorganization energies (higher rates) for molecules at smaller separations, it often cannot provide quantitative estimates, since charge distributions are rarely spherically symmetric.

Alternatively, the displacement fields can be constructed using the atomic partial charges. The difference of the displacement fields at the position of an atom b outside the charge transfer complex (molecule k ≠ i, j) can be expressed aswhere q (q) is the partial charge of atom a of the neutral (charged) molecule i in a vacuum. The partial charges of neutral and charged molecules are obtained by fitting their values to reproduce the electrostatic potential of a single molecule (charged or neutral) in a vacuum. Assuming a uniform density of atoms, the integration in eq 4 can be rewritten as a density-weighted sum over all atoms excluding those of the charge transfer complex.

Using eq 6, λout/cp was calculated for all pairs from the neighbor list of a system of 512 Alq3 molecules. The neighbor list was constructed using a cutoff of 0.8 nm for the centers of mass of the three Alq3 ligands, which results in an average of 12 neighbors in the first coordination shell. The electrostatic potential of a single molecule in a vacuum was calculated using the B3LYP functional and a 6-311G(d,p) basis set. The CHELPG method(68) was used to obtain atomic partial charges. The resulting distribution of λout/cp is shown in Figure 5, together with a fit to eq 5. The fit yields RAlq3 = 0.57 nm and predicts negative λout for separations smaller than this radius, which is unphysical.

Figure 5

Outersphere reorganization energy divided by the Pekar factor as a function of the distance between two molecules and its fit to eq 5.

The remaining unknown needed to calculate λout is the Pekar factor, cp. In polar solvents εs ≫ εopt ∼ 1, and cp is on the order of 1. In most organic semiconductors, however, molecular orientations are fixed, and therefore the low frequency dielectric permittivity is on the same order of magnitude as εopt. Hence, cp is small, and its value is very sensitive to differences in the permittivities.

For Alq3, εs = 3.0 ± 0.3 is the experimentally measured dielectric constant at low frequencies,(69) while at optical frequencies below electronic absorption, εopt = 2.9 ± 0.1.(70) Thus, cp = 0.01 ± 0.04, yielding outersphere reorganization energies of λout < 0.08 eV, which are small compared to λint. Similar results have been reported for other organic semiconductors and different methods for computing λout.[71−73]

Site Energy Difference

A charge transfer reaction between molecules i and j is driven by the site energy difference, ΔE = E – E. Since the transfer rate, ω, depends exponentially on ΔE (see eq 1), it is important to compute its distribution as accurately as possible. The total site energy difference has contributions due to an externally applied electric field, electrostatic interactions, polarization effects, and internal energy differences. In what follows, we discuss how to estimate these contributions by making use of first-principles calculations and polarizable force fields.

Externally Applied Electric Field

The contribution to the total site energy difference due to an external electric field is given by ΔEext = q(·), where q = ±e is the charge and = – is a vector connecting molecules i and j. For typical distances between small molecules, which are on the order of 1 nm, and moderate fields of F < 108 V/m, this term is always smaller than 0.1 eV.

Electrostatic Energy

Variations of the local electric field can result in large electrostatic contributions to the energetic disorder.(74) When the atomic partial charges of charged and neutral molecules are used, as introduced in section , ΔEel can be computed from the site energies(31)where r = | – | is the distance between atoms a and b and εs is the static relative dielectric constant. The first sum extends over all atoms of molecule i, for which the site energy is calculated. The second sum reflects interactions with all atoms of neutral molecules k ≠ i. By using eq 7, one assumes that the influence of conformational changes on partial charges and changes of the molecular geometry upon charging are small.

In order to minimize finite size effects, we do not use a spherical cutoff but apply the nearest image convention, that is, sum over all neutral molecules in the box after centering the box around the charged molecule. For Alq3, with long-range interactions due to its large dipole moment, this procedure converges already for a few hundred molecules.

The resulting distribution of the site energy differences without screening (εs = 1), shown in Figure 6, is Gaussian, with variance of σ = 0.30 eV. Note that ΔEel is constructed on the basis of the neighbor list as described in section .

Figure 6

Distributions of site energy differences without (blue) and with (red) polarization effects for pairs from the neighbor list. Solid lines are fits to Gaussian distributions. The dashed line corresponds to a parametrized distance-dependent ε(r) according to eq 8.

Polarization Effects

The influence of polarization effects on the Coulomb interactions can be taken into account by using a relative dielectric constant in eq 7. Bulk values of εs = 2–5 for typical organic semiconductors uniformly scale all site energies but are not capable of describing polarization effects on a microscopic level. The contribution to Eel from the first coordination shell is then underestimated due to overscreening, and as a result, the site-energy differences become artificially small. Alternatively, one can introduce a phenomenological distance-dependent screening function ε(r) in eq 7(30)where the parameter s is the inverse screening length. For a monovalent ion in water, for example, εs = 80 and s = 3 nm–1.(75) This screening function ensures that neighboring atoms interact via an unscreened Coulomb potential (ε ∼ 1), while the electrostatic interaction between atoms at large separations is screened as in the bulk.

While phenomenological distance-dependent screening is computationally efficient, it cannot be used for inhomogeneous systems or systems with anisotropic molecular polarizabilities. Moreover, εs and s are not known for newly synthesized compounds. A more general approach relies on self-consistent methods to obtain polarization fields.(76) Here, we use a polarizable force field based on the Thole model(77) as implemented in the TINKER package.(78)

The polarization contribution is refined iteratively. After evaluating the electric field at atom a in molecule i, (0), created by all atomic partial charges (εs = 1, nearest image convention), the induced dipole moments, μ(0), are computed. During this first step, intramolecular interactions are excluded. Induced dipole moments are then iteratively refined as μ( = ω(α +(1−ω)μ(, where α is the isotropic atomic polarizability and ω = 0.5 is a damping constant for successive over-relaxation. The new electric fields are computed using the induced dipole moments, which now interact with each other also within molecules, allowing for anisotropic molecular polarizabilities. The procedure is repeated until ∑|μ( – μ(| < 10–6 Debye.

Such self-consistent calculations can, however, become computationally prohibitive for large systems. For homogeneous systems and isotropic molecular polarizabilties, one can perform self-consistent calculations for small systems, parametrize eq 8 accordingly, and use ε(r) to study larger systems. To this end, the bulk dielectric constant is obtained from the Clausius–Mosotti relation(79)where α is the molecular polarizability volume and N/V is the number density. Using this value of εs, the parameter s is then fitted to reproduce the distribution of site-energy differences for molecules from the neighbor list.

For a neutral Alq3 molecule, the Thole model (using atomic polarizabilities αH,C,N,O,Al = 0.696, 1.75, 1.073, 0.837, 5.5 Å3, respectively, and a damping factor of a = 0.39 for interactions with induced moments(78)) gives a practically isotropic polarizability volume tensor with α = 54.9 Å3. This agrees with DFT calculations (B3LYP functional and 6-311G(d,p) basis set), yielding α = 55.2 Å3.(80) Using N = 512 molecules in a cubic box of length L = 67.8 Å, we obtain εs = 2.84, which reproduces the experimental value of 3.0 ± 0.3.(69)

The corresponding distributions of site energy differences, shown in Figure 6, are practically Gaussian. For the Thole model, the variance of 0.21 eV is obviously larger than the 0.10 eV obtained using bulk screening with εs = 3.0 (not shown) and is smaller than 0.30 eV for εs = 1. A fit to eq 8 gives s = 1.3 nm–1, which is significantly smaller than the inverse screening length of water, 3 nm–1.

Internal Energy Difference

The contribution to the site energy difference due to different internal energies (see Figure 4) can be written aswhere U is the total energy of molecule i in the charged (neutral) state and geometry. ΔU corresponds to the adiabatic ionization potential (or electron affinity) of molecule i, as shown in Figure 4. For one-component systems and negligible conformational changes ΔEint = 0, while it is significant for donor–acceptor systems.

In Alq3, significant conformational changes (see section ) lead to a Gaussian distribution of ΔEint with a small variance of σint = 0.01 eV. The internal energy disorder is small compared to the electrostatic (including polarization) energetic disorder and hence does not affect the charge carrier mobility.

Spatial Correlations of Energetic Disorder

Long-range, e.g., electrostatic and polarization, interactions often result in spatially correlated disorder,(81) which affects the onset of the mobility-field (Poole–Frenkel) dependence.[30,82,83]

To quantify the degree of correlation, one can calculate the spatial correlation function of E and E at a distance rwhere ⟨E⟩ is the average site energy. C(r) is 0 if E and E are uncorrelated and 1 if they are fully correlated. For a system of randomly oriented point dipoles, the correlation function decays as 1/r at large distances.(84)

For systems with spatial correlations, variations in site energy differences, ΔE, of pairs of molecules from the neighbor list are smaller than variations in site energies, E, of all individual molecules. Since only neighbor list pairs affect transport, the distribution of ΔE rather than that of individual site energies, E, should be used to characterize energetic disorder.

For Alq3, the spatial correlation function of the electrostatic contribution to site energies, which is calculated for 512 molecules, is shown in Figure 7. It qualitatively reveals strong correlations due to the large dipole moment of the meridional isomer of Alq3 of approximately 4 Debye. Quantitatively, this result is not converged with respect to the system size, and bigger systems will exhibit even longer-ranged correlations. The inset of Figure 7 shows that distributions of ΔE for all and neighbor-list-only pairs are clearly different. Note that respective distributions of internal site energies are identical, indicating that this type of disorder is spatially uncorrelated.

Figure 7

Electrostatic site energy correlation function (eq 11) calculated for pairs from the neighbor list without (blue) and with polarization effects from the Thole model (red) as well as using a parametrized distance-dependent ε(r) according to eq 8 (dashed line). Inset: Gaussian fits to electrostatic site energy distributions for all pairs (σ = 0.30 eV) and for pairs from the neighbor list (σ = 0.21 eV).

Solving the Master Equation

Having determined the list of conjugated segments (hopping sites) and charge transfer rates between them, the next task is to solve the master equation, which describes the time evolution of the systemwhere Pα is the probability of the system to be in a state α at time t and Ωαβ is the transition rate from state α to state β. A state α is specified by a set of site occupations, {α}, where α = 1 (0) for an occupied (unoccupied) site i, and the matrix Ω̂ can be constructed from rates ω.

In particular, for a system with only one charge carrier, each state is uniquely characterized by the index i of the site the carrier occupies. In other words, only states of type i ≡ {0, ..., 0, α = 1, 0, ..., 0} are possible, and the corresponding probabilities P and the transition rates Ω are identical to site occupation probabilities p and the transfer rates ω, respectively. Equation 12 can then be rewritten asand can be solved using linear algebra. While being efficient for stationary, low charge carrier density cases (one charge carrier per simulation box), this approach can become unstable for systems with high energetic disorder, where rates vary by several orders of magnitude.

In more general cases, such as multiple charge carriers, expressing the state picture (eq 12) in terms of site occupations is required because of an extremely large total number of states. For multiple charge carriers, the master equation can still be rewritten in terms of occupation probabilities (see the Supporting Information) by assuming only site-blocking charge–charge interactions and by using a mean-field approximation.(85) The analogue of eq 13 becomes, however, nonlinear and requires special solvers. If, in addition, several different types of carriers, such as holes, electrons, and excitons, are present in the system and their creation/annihilation processes take place, it is practically impossible to link state and site occupation probabilities and the corresponding rates.

Instead, the solution of eq 12 can be obtained by using kinetic Monte Carlo (KMC) methods. KMC explicitly simulates the dynamics of charge carriers by constructing a Markov chain in state space and can find both stationary and transient solutions of the master equation. The main advantage of KMC is that only states with a direct link to the current state need to be considered at each step. Since these can be constructed solely from current site occupations, extensions to multiple charge carriers (without the mean-field approximation), site-occupation dependent rates (needed for the explicit treatment of Coulomb interactions), and different types of interacting particles and processes are straightforward.

To optimize memory usage and efficiency, a combination of the variable step size method(86) and the first reaction method is implemented as described in the Supporting Information.

Extrapolation to Nondispersive Mobilities

Predictions of charge-carrier mobilities in partially disordered semiconductors rely on charge transport simulations in systems which are only several nanometers thick. As a result, simulated charge transport might be dispersive for materials with large energetic disorder,[87,88] and simulated mobilities are system-size-dependent. In time-of-flight (TOF) experiments, however, a typical sample thickness is in the micrometer range, and transport is often nondispersive. In order to link the simulation and experiment, one needs to extract the nondispersive mobility from simulations of small systems, where charge transport is dispersive at room temperature.

Such extrapolation is possible if the temperature dependence of the nondispersive mobility is known in a wide temperature range. For example, one can use analytical results derived for one-dimensional models.[82,89,90] The mobility-temperature dependence can then be parametrized by simulating charge transport at elevated temperatures, for which transport is nondispersive even for small system sizes. This dependence can then be used to extrapolate to the nondispersive mobility at room temperature.(32)

For Alq3, the charge carrier mobility of a periodic system of 512 molecules was shown to be more than 3 orders of magnitude higher than the nondispersive mobility of an infinitely large system.(32) Furthermore, it was shown that the transition between the dispersive and nondispersive transport has a logarithmic dependence on the number of hopping sites N. Hence, a brute-force increase of the system size cannot resolve the problem for compounds with large energetic disorder σ, since N increases exponentially with σ2.

Macroscopic Observables

Spatial distributions of charge and current densities can provide better insight into the microscopic mechanisms of charge transport. If O is an observable which has the value Oα in a state α, its ensemble average at time t is a sum over all states weighted by the probability Pα to be in a state α at time tIf O does not explicitly depend on time, the time evolution of ⟨O⟩ can be calculated asIf averages are obtained from KMC trajectories, Pα = sα/s, where sα is the number of Markov chains ending in the state α after time t, and s is the total number of chains.

Alternatively, one can calculate time averages by analyzing a single Markov chain. If the total occupation time of the state α is τα, thenwhere τ = ∑α τα is the total time used for time averaging.

For ergodic systems and sufficient sampling times, ensemble and time averages should give identical results. In many cases, the averaging procedure reflects a specific experimental technique. For example, an ensemble average over several KMC trajectories with different starting conditions corresponds to averaging over injected charge carriers in a time-of-flight experiment. In what follows, we focus on the single charge carrier (low concentration of charges) case.

Charge Density

For a specific type of particles, the microscopic charge density of a site i is proportional to the occupation probability of the site, pwhere, for an irregular lattice, the effective volume V can be obtained from a Voronoi tessellation of space. For reasonably uniform lattices (uniform site densities), this volume is almost independent of the site, and a constant volume per cite, V = V/N, can be assumed. In the macroscopic limit, the charge density can be calculated using a smoothing kernel function, i.e., a distance-weighted average over multiple sites. Site occupations p can be obtained from eq 14 or eq 16 by using the occupation of site i in state α as an observable.

If the system is in thermodynamic equilibrium, that is without sources or sinks and without circular currents (and therefore no net flux), a condition known as detailed balance holdsIt can be used to test whether the system is ergodic or not by correlating log p and the site energy E. Indeed, if λ = λ, the ratios of forward and backward rates are determined solely by the energetic disorder, ω/ω = exp(−ΔE/kBT) (see eq 1).

Current

If the position of the charge, , is an observable, the time evolution of its average ⟨⟩ is the total current in the systemSymmetrizing this expression, we obtainwhere = – . Symmetrization ensures equal flux splitting between neighboring sites and the absence of local average fluxes in equilibrium. It allows one to define a local current through site i asA large value of the local current indicates that the site contributes considerably to the total current. A collection of such sites thus represents most favorable charge pathways.

Figure 8 illustrates site currents for amorphous Alq3 for a system of 512 molecules. The distribution of currents is very inhomogeneous, and some pathways are sampled more frequently than the others, which is a direct consequence of a rough and correlated energy landscape.(91)

Figure 8

Isosurface of the current density for amorphous Alq3 (512 molecules, external field F = 107V/m, distance-dependent dielectric constant, DFT-based transfer integrals). Currents have filamentary structure due to large correlated energetic disorder.(91) The red streamtraces depict interpolated charge pathways.

Mobility and Diffusion Constant

For a single particle, e.g., a charge or an exciton, a zero-field mobility can be determined by studying particle diffusion in the absence of external fields. Using the particle displacement squared, Δ2, as an observable, we obtainHere, is the coordinate of the site i; Dγδ is the diffusion tensor, γ,δ = x, y, z; and d = 3 is the system dimension. Using the Einstein relationone can, in principle, obtain the zero-field mobility tensor μγδ. Equation 22, however, does not take into account the use of periodic boundary conditions when simulating charge dynamics. In this case, the simulated occupation probabilities can be compared to the solution of the Smoluchowski equation with periodic boundary conditions (see the Supporting Information for details).

Alternatively, one can directly analyze time-evolution of the KMC trajectory and obtain the diffusion tensor from a linear fit to the mean square displacement, Δ = 2dDγδt.

The charge carrier mobility tensor, μ̂, for any value of the external field can be determined either from the average charge velocity defined in eq 19or directly from the KMC trajectory. In the latter case, the velocity is calculated from the unwrapped (if periodic boundary conditions are used) charge displacement vector divided by the total simulation time. Projecting this velocity on the direction of the field yields the charge carrier mobility in this particular direction. In order to improve statistics, mobilities can be averaged over several KMC trajectories and MD snapshots.

For Alq3, the field dependence of the mobility (Poole–Frenkel plot) is shown in Figure 9 for a system of 4096 molecules. To illustrate the role of disorder and correlations, we also show the field dependence for a system without energetic disorder (top panel) and without correlated energetic disorder (randomly shuffled site energies). Energetic disorder reduces the value of mobilty by 6 orders of magnitude. The Poole–Frenkel behavior for small fields can only be observed if correlated disorder is taken into account. Note that, for a system with such large energetic disorder, the absolute values of (nondispersive) mobility are systematically overestimated due to significant finite size effects (see section and ref (32)). The experimentally measured value of the hole mobility at small fields lies between 10–9 and 10–8cm2 V–1 s–1.(92)

Figure 9

Poole–Frenkel plots for a system of 4096 molecules. Mobilities were calculated by averaging over ten, 0.1-s-long (10–5 s for no disorder), KMC runs and six different field directions. Transfer integrals were calculated using DFT; energetic disorder is based on the distance-dependent dielectric constant fitted to the site energy distribution of the Thole model.

Implementation

The toolkit is implemented using modular concepts introduced earlier in the versatile object-oriented toolkit for coarse-graining applications (VOTCA).(42) The VOTCA structures are adapted to reading atomistic trajectories, mapping them onto conjugated segments and rigid fragments, and substituting (if needed) rigid fragments with the optimized copies.

The molecular orbital overlap module calculates electronic coupling elements between conjugated segments from the corresponding molecular orbitals. It relies on the semiempirical INDO Hamiltonian and molecular orbitals in the format provided by the GAUSSIAN package. An alternative, density-functional-based approach, has interfaces to the GAUSSIAN and TURBOMOLE packages. An interface to the TINKER package is provided for calculations of electrostatic and polarization contributions to energetic disorder.

The kinetic Monte Carlo module reads in the neighbor list, site coordinates, and hopping rates and performs charge dynamics simulations using either periodic boundary conditions or charge sources and sinks.

The toolkit is written as a combination of modular C++ code and scripts. The data transfer between programs is implemented via a state file or database, which is also used to restart simulations. Analysis functions and most of the calculation routines are encapsulated by using the observer pattern,(93) which allows the implementation of new functions as individual modules.

Summary

To summarize, we have presented a toolkit for developing and testing methods for charge transport simulations in disordered organic semiconductors. The core of the toolkit is based on a reader and a postprocessor of atomistic trajectories, a fast molecular orbital overlap calculations library, and a kinetic Monte Carlo code. To illustrate its functionality, we have studied charge transport in amorphous tris-(8-hydroxyquinoline)aluminum, a typical organic semiconductor. The source code of the toolkit is available under the Apache license (www.votca.org).