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 . 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 . in parallel ,
increased device efficiencies could be achieved , e.g. , by optimizing light in- and out - coupling and introducing tandem concepts .
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 , 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 .
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 .
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 .
alternatively , charge transfer theories can be used to evaluate rates from quantum chemical calculations . 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 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 , ij = ijint + ijout is the reorganization energy , which is a sum of intra- and intermolecular ( outer - sphere ) contributions , eij is the site - energy difference , or driving force , and jij 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 ) .
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 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 ) .
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 ) .
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
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 . 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 ,
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 .
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 .
in alq3 , a metal - coordinated compound , a charge carrier is delocalized over all three ligands .
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 . 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 .
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 .
the electronic transfer integral jij 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 , jij 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 , it is essential to calculate jij 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 .
this molecular orbital overlap ( moo ) method has been applied successfully to study charge transport , for instance , in discotic liquid crystals , polymers,(36 ) or partially disordered organic crystals .
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 . 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 ,
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 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 .
the reorganization energy ij 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 , ijint , which is due to a reorganization of nuclear coordinates of the two molecules forming the charge transfer complex , and intermolecular ( outersphere ) , ijout , which is due to the relaxation of the environment . in what follows , we discuss how these contributions can be calculated . 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 , ijint .
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 , uinc 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 , ujcn 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 ,
thus , ijint is a property of the charge transfer complex and not of a single molecule . potential energy surfaces of ( a ) donor and ( b ) acceptor in charged and neutral states . after the change of the charge state
if all vibrational modes are treated classically , the total internal reorganization energy and the internal energy difference of the electron transfer reaction are ijint = icn + jnc and eijint = ui  uj , respectively . in alq3
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 ijint 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 ijint 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 . 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 .
assuming that charge transfer is much slower than electronic polarization but much faster than nuclear rearrangement of the environment , ijout can be calculated from the electric displacement fields created by the charge transfer complex:(67)where 0 is the permittivity of free space , di , f(r ) 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 , v 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 v outside of the two spheres of radii ri and rj centered on molecules i and j leads to the classical marcus expression for the outersphere reorganization energy : where rij 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 can not provide quantitative estimates , since charge distributions are rarely spherically symmetric .
the difference of the displacement fields at the position of an atom bk outside the charge transfer complex ( molecule k  i , j ) can be expressed aswhere qain ( qaic ) 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 , ijout / 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 resulting distribution of ijout / cp is shown in figure 5 , together with a fit to eq 5 .
the fit yields ralq3 = 0.57 nm and predicts negative ijout for separations smaller than this radius , which is unphysical .
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 ijout 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 ijout < 0.08 ev , which are small compared to ijint .
similar results have been reported for other organic semiconductors and different methods for computing ijout . a charge transfer reaction between molecules
i and j is driven by the site energy difference , eij = ei  ej . since the transfer rate , ij , depends exponentially on eij ( 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 .
the contribution to the total site energy difference due to an external electric field f is given by eijext = q(frij ) , where q = e is the charge and rij = ri  rj 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 < 10 v / m , this term is always smaller than 0.1 ev .
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 , eijel can be computed from the site energies(31)where raibk = |rai  rbk| is the distance between atoms ai and bk 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 eijel is constructed on the basis of the neighbor list as described in section .
distributions of site energy differences without ( blue ) and with ( red ) polarization effects for pairs from the neighbor list .
the dashed line corresponds to a parametrized distance - dependent (r ) according to eq 8 .
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 = 25 for typical organic semiconductors uniformly scale all site energies but are not capable of describing polarization effects on a microscopic level .
the contribution to eiel 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 (raibk ) 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.(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 can not 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 , fai(0 ) , created by all atomic partial charges ( s = 1 , nearest image convention ) , the induced dipole moments , ai(0 ) , are computed . during this first step ,
. induced dipole moments are then iteratively refined as ai(k+1 ) = fai(k)ai + ( 1)ai(k ) , where ai 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 .
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 ,
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  , 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  .
this agrees with dft calculations ( b3lyp functional and 6 - 311g(d , p ) basis set ) , yielding  = 55.2 .(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 , which is significantly smaller than the inverse screening length of water , 3 nm .
the contribution to the site energy difference due to different internal energies ( see figure 4 ) can be written aswhere uicc(nn ) is the total energy of molecule i in the charged ( neutral ) state and geometry .
ui 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
eijint = 0 , while it is significant for donor  acceptor systems . in alq3 , significant conformational changes ( see section )
lead to a gaussian distribution of eijint with a small variance of  = 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 .
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 . to quantify the degree of correlation
, one can calculate the spatial correlation function of ei and ej at a distance rijwhere e is the average site energy .
c(rij ) is 0 if ei and ej 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 , eij , of pairs of molecules from the neighbor list are smaller than variations in site energies , ei , of all individual molecules . since only neighbor list pairs affect transport , the distribution of eij rather than that of individual site energies , ei , 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 eij 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
. 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 ) .
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 , { i } , where i = 1 ( 0 ) for an occupied ( unoccupied ) site i , and the matrix  can be constructed from rates ij . 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 , i = 1 , 0 , ... , 0 } are possible , and the corresponding probabilities pi and the transition rates ij are identical to site occupation probabilities pi and the transfer rates ij , respectively .
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 .
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 , 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 .
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 can not resolve the problem for compounds with large energetic disorder  , since n increases exponentially with .
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 .
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 ,
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 . for a specific type of particles , the microscopic charge density of a site
i is proportional to the occupation probability of the site , piwhere , for an irregular lattice , the effective volume vi 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 , vi = 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 pi 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 pi and the site energy ei . indeed ,
if ij = ji , the ratios of forward and backward rates are determined solely by the energetic disorder , ji/ij = exp(eij / kbt ) ( see eq 1 ) . if the position of the charge , r , is an observable , the time evolution of its average r is the total current in the systemsymmetrizing this expression , we obtainwhere rij = ri  rj .
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 .
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 ) isosurface of the current density for amorphous alq3 ( 512 molecules , external field fz = 10v / 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 . 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 , ri , as an observable , we obtainhere
, ri 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 , ri,ri, = 2ddt .
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 f yields the charge carrier mobility in this particular direction . in order to improve statistics ,
, 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 ) .
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 and 10 cm v s.(92 ) poole 
mobilities were calculated by averaging over ten , 0.1-s - long ( 10 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 .
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 .
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 ) .
 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
Methods
Macroscopic Observables
Implementation
Summary