Detailed analysis of charge transport in amorphous organic thin layer by multiscale simulation without any adjustable parameters

Hopping-type charge transport in an amorphous thin layer composed of organic molecules is simulated by the combined use of molecular dynamics, quantum chemical, and Monte Carlo calculations. By explicitly considering the molecular structure and the disordered intermolecular packing, we reasonably reproduce the experimental hole and electron mobilities and their applied electric field dependence (Poole–Frenkel behaviour) without using any adjustable parameters. We find that the distribution of the density-of-states originating from the amorphous nature has a significant impact on both the mobilities and Poole–Frenkel behaviour. Detailed analysis is also provided to reveal the molecular-level origin of the charge transport, including the origin of Poole–Frenkel behaviour.

tential was used as a model for the van der Waals interaction with tail correction 4 , which is a correction for thermodynamic quantities that takes into account the homogeneous long range van der Waals inter-N N Figure S1. Structural formula of CBP.

1/11
action. Particle-particle-particle-mesh method 5 was used to calculate Coulombic interaction between atoms. To determine the atomic charges, we performed DFT calculation on an optimized structure of isolated CBP molecule and then used the Merz-Singh-Kollman scheme 6 , which is a method for calculating atomic charges to fit the electrostatic potential on molecular surfaces.
As an initial structure of the MD simulation, the DFT-optimized molecules were randomly placed in a large cubic cell. To mimic the vapor deposition process, the MD simulation was performed in a NVT ensemble under 573 K for 10 ps followed by a simulation in a NPT ensemble under 298 K and 1.0 × 10 −4 Pa for 1.0 ns using a Nosé-Hoover thermostat and barostat 7,8 . During the simulation in the NPT ensemble, the simulation cell shrank and finally reached a constant volume. We obtained the cubic cell of the condensed structure with a side length of 14.55 nm and the density of the generated structure was calculated to be 1.04 g cm −3 . Finally, to eliminate the deviation from the stable structure originating from molecular vibration, geometry optimization was performed.

Reorganization energy (λ )
Structure optimizations and the calculations of λ were performed by DFT calculations [Becke threeparameter Lee-Yang-Parr (B3LYP) functional 9 and the 6-31G(d) basis set]. λ is described as λ = λ 1 + λ 2 . λ 1 and λ 2 are described as: where E c n , E c c , E n c , and E n n are the molecular energy of the charged state in the neutral optimized geometry, the charged state in the charged optimized geometry, the neutral state in the charged optimized geometry, and the neutral state in the neutral optimized geometry, respectively 10 . In most cases, λ is calculated for an isolated molecule. However, interactions with neighbouring molecules in the aggregate structure are considered to affect λ . Here, we introduced a quantum mechanics / molecular mechanics (QM/MM) approach 11 to quantify such environmental influences. The λ including the influences is denoted by λ aggr .
For 45 randomly sampled molecules (out of 4000 molecules in the aggregate structure), we calculated λ aggr by the QM/MM method. DFT and the Dreiding force field were used for the QM and MM regions, respectively. The neighbouring molecules within 30Å were included and their geometries were frozen during the calculation. Partial charges of the atoms in the MM region were incorporated into the QM Hamiltonian. The partial charges were set to be same values as those used in MD calculation. The calculated average value of λ aggr was 0.105 ± 0.009 eV for holes and 0.373 ± 0.027 eV for electrons (95% confidence interval), respectively. We also calculated the reorganization energy for the optimized isolated CBP molecule, λ isolated . λ isolated (0.129 eV and 0.516 eV for hole and electron, respectively) is larger than λ aggr , owing to the absence of the steric effect from neighbouring molecules.

Electronic coupling (H i j )
H i j was calculated using the extended Hückel method as described in previous work 10,12,13 . Calculation was performed for all pairs with a centre-to-centre distance within 25Å (167,993 pairs). The cutoff distance (25Å) is validated by H i j versus the intermolecular distance plots (Fig. S2), which show that the molecular pairs with a centre-to-centre distance longer than 25Å have a negligibly small H i j (<10 −4 meV). Figure

Energetic disorder
We calculated the energetic disorder by considering distributed electrostatic interaction with and polarization effect of neighbouring molecules. E 0/+/− i,neighbour , the sum of the two terms, is then: Here, the superscripts 0, +, and − represent neutral, positively charged, and negatively charged states of is an electrostatic interaction energy between the atomic partial charges of molecule i and neighbouring molecules, where q i,k and q j,k ′ are an atomic partial charge of atom k in molecule i and that of atom k ′ in molecule j, respectively, and r r r i j,kk ′ is a vector connecting atom k in molecule i to atom k ′ in molecule j.

5/11
Atomic partial charges of charged and neutral molecules were obtained by DFT calculation and the Merz-Singh-Kollman scheme (same method as used for the MD simulation). E 0/+/− i,charge-charge was calculated taking into account neighbouring molecules with a centre-to-centre distance within 30Å. E 0/+/− i,charge-dipole is a polarization effect of neighbouring molecules, which effectively moderates the interaction between atomic partial charges. This effect can be regarded as the interaction between atomic partial charges in molecule i and induced dipole moments on atoms in neighbouring molecules. Hence, E 0/+/− i,charge-dipole can be approximately written as: where µ µ µ j,k ′ is an induced dipole moment on atom k ′ in neighbouring molecule j. µ µ µ j,k ′ can be described where I isolated and A isolated are the ionization potential and electron affinity of an isolated CBP molecule, respectively. Therefore, the difference of the Gibbs free energy associated with the hole/electron transfer from molecule i to j, ∆G +/− i j , can be expressed as: where q is the charge of the carrier, F F F is an externally applied electric field (|F F F| = F) and x x x i j is a vector connecting the centre of molecule i to the centre of molecule j (displacement of the carrier). The distri-

6/11
bution of the energy level obtained in this work was Gaussian (Fig. S5)  Solid lines show the Gaussian probability density function, which is zero-centred and whose standard deviation is set to be the sample standard deviation of the calculated values.

Kinetic Monte Carlo simulation
Using rate constants for the charge transfer from molecule i to j, k +/− i j , calculated with the Marcus equation (the superscript +/− denotes that the quantities are for hole/electron transport): Hole, x-axis Hole, y-axis Hole, z-axis Electron, x-axis Electron, y-axis Electron, z-axis 100 molecules 1000 molecules 8000 molecules a b c d Figure S7. Calculated µ for hole and electron transport with various size of system. µ was calculated for three orthogonal axes with (a) 100, (b) 1000, (c) 4000, and (d) 8000 CBP systems.