Challenges for the Accurate Simulation of Anisotropic Charge

Article pubs.acs.org/JPCC

Challenges for the Accurate Simulation of Anisotropic Charge Mobilities through Organic Molecular Crystals: The β Phase of merTris(8-hydroxyquinolinato)aluminum(III) (Alq3) Crystal Shiwei Yin,*,† Lanlan Li,† Yongmei Yang,† and Jeffrey R. Reimers*,‡ †

Key Laboratory for Macromolecular Science of Shaanxi Province, School of Chemistry & Chemical Engineering, Shaanxi Normal University, Xi'an City, People's Republic of China, 710062 ‡ School of Chemistry F11, The University of Sydney, NSW 2006 Australia S Supporting Information *

ABSTRACT: Quantitative agreement has been found between observed and calculated charge mobilities through organic conductors, despite the use of many assumptions in the calculations, including: the relative strength of the intermolecular electronic coupling to the reorganization energy driving charge localization, the treatment of site variability in the material, the involvement of tunneling processes during charge hopping between sites, the use of weak-coupling-based perturbation theory to determine hopping rates, the residence times for charges on sites, the effect of the large field strengths used in experimental studies, the general appropriateness of simple one-dimensional diffusion modeling approaches, and the involvement of molecular excited states of the ions. We investigate the impact of these assumptions, concluding that all may be very significant. In some cases, methodological options are considered, and optimum procedures are determined, showing that (i) the use of Koopmans' theorem to estimate intermolecular couplings in solids is problematic and (ii) the correct expression for the residence lifetime of a charge on a crystal site. These conclusions are drawn from simulations of anisotropic charge mobilities through the β phase of mer-tris(8-hydroxyquinolinato)aluminum(III) (Alq3) crystal, a material commonly used in OLED applications. Calculations are compared that determine mobilities at finite applied field from drift velocities through either semianalytical solutions of the master equation or else kinetic Monte Carlo simulations, as well as those that determine mobilities from multidimensional diffusion coefficients at zero field by Monte Carlo and those that analytically solve simplified onedimensional diffusion models. For crystalline Alq3 itself, the calculations predict electron mobilities that are 4−6 orders of magnitude larger than those predicted by similar methods for amorphous Alq3, in agreement with experimental findings. This work vindicates recent theories describing the poor mobilities of the amorphous material, forming a complete basic picture for Alq3 conductivity.

1. INTRODUCTION Because of the potential application in transistors, photodiodes, solar cells, and chemical sensors,1−6 π-conjugated organic semiconducting materials have attracted much attention from both academia and industry. The performance of these devices strongly depends on their ability to transport charge within the organic material, quantified usually via charge-carrier mobilities. The influence of molecular packing on mobilities has been a subject of great interest for many years. Various models have been proposed linking microscopic molecular properties to charge mobility, these being based typically on either strong coupling or weak coupling approximations. The strong-coupling limit arises when the intermolecular electronic coupling [electron transfer (ET) integral] driving charge transport is far larger than the electron-vibration coupling (reorganization energy) driving charge localization. This allows the charge to delocalize over many molecules, with © 2012 American Chemical Society

the material developing a band of densely spaced energy levels within which the mobile charges are located; this effect can be quantified using say a Holstein−Peierls model.7−13 Charge mobilities are then expressed in terms of typically either the thermally averaged velocity−velocity tensors or the effective mass tensor of the charge. The band structures of the appropriate frontier orbital, either the highest occupied molecular orbital (HOMO) for hole transport or the lowest unoccupied molecular orbital (LUMO) for electron transport, assuming either a constant relaxation time or a constant collision-free distance.8,14−16 In the strong-coupling limit, mobilities usually increase at low temperature as the associated ordering increases delocalization lengths. However, Troisi and Received: April 18, 2012 Revised: June 13, 2012 Published: June 14, 2012 14826

dx.doi.org/10.1021/jp303724r | J. Phys. Chem. C 2012, 116, 14826−14836

The Journal of Physical Chemistry C

Article

nonadiabatic limit, but is this approximation valid?58−61 Approximate analytical solutions to the diffusion equation provide significant chemical insight,29 but do the evoked approximations reflect physical reality?27 Two different forms for the hopping rate are commonly used in diffusion-equation approaches,24,27−29,62,63 but which form is most appropriate? Also, two generic types of methods are most commonly applied to determine the intermolecular electronic coupling strengths,19,21,30 but does the computationally attractive much simpler scheme deliver reliable results? Is the commonly used assumption of homogeneity throughout a molecular crystal valid?27 Does the large electric-field strength typically used in experimental studies modulate the mobilities? Only the properties of the ground states of molecular ions are included in modeling scenarios, but could excited states, known to play critical roles in many photochemical charge-transport processes,64 also be important? Here, we investigate these questions, demarking challenges for the accurate simulation of anisotropic charge mobilities through molecular crystals. Charge transport through the β phase of mer-tris(8-hydroxyquinolinato)aluminum(III) (Alq3) crystals65 is taken as a case study for these investigations. This material is very important as both an electron carrier and a hole carrier in organic light-emitting diodes (OLEDs).65,66 The mobility of electrons in crystalline Alq3 is very high (0.23 ± 0.03 cm2 V−1 s−1) but has been measured only in its direction of greatest mobility.67 Amorphous materials have been studied in much greater detail both experimentally68−71 and computationally,22,72−77 but the mobilities are reduced by 4−6 orders of magnitude. Calculations have revealed basic properties of Alq3 materials;16,30,78−81 while those on crystalline Alq3 indicate that the mobility should be strongly anisotropic,16,30 no quantitative predictions have been made. Calculations of the amorphous phase indicate that a combination of large disorder combined with nonoptimal intermolecular alignments dominates the observed low mobilities,7,23 attaining semiquantitative agreement with observed data.77 Here, we unite these studies by applying similar computational methods for the calculation of crystalline anisotropic mobilities. This system serves as a testbed for the understanding of how the many approximations used in mobility calculations manifest in practical molecular conductors.

Orlandi have shown that this is not necessarily the case if structural variations or thermal motions strongly modulate the electronic coupling.17,18 When the ET integral is far less than the reorganization energy, the charge carriers localize onto a single molecule in the weak-coupling limit.19 Charge transport therefore proceeds via an incoherent hopping process that requires thermal activation,20 with a large number of such processes required for macroscopic charge transport. The rate constants for individual reactions in this process are usually modeled19,21−36 using Marcus−Hush ET theory37 in either its traditional hightemperature classical38−41 or extended low-temperature semiclassical forms.42−45 In the classical form, mobilities increase with increasing temperature, but in the semiclassical form, tunneling contributes to the reaction rate as well as thermally activated processes, and it is hence possible for mobilities to decrease with increasing temperature because of thermal inhibition of the tunneling.45 On the basis of the deduced hopping rates, the charge mobility can be evaluated by setting up a macroscopic network of reactions in the presence of an external electric field and hence determining the drift velocity from the solution of the associated master equation.27,46−48 Recently, this solution has become available via efficient and accurate semianalytical means,27,46−48 but traditionally, it has been obtained fully numerically using kinetic Monte Carlo (MC) simulations.22,27,46−50 These techniques are very powerful and can deal with important features such as crystal imperfections and thermal fluctuations. Nevertheless, master-equation approaches are always complex, and it is also desirable to determine mobilities in the weakcoupling limit using simpler and more intuitive methods that offer clear physical insight.29 In the limit of zero-applied electric field strength, the master equation collapses into a simple diffusion equation as a result of the Einstein−Smoluchowski relationship51,52 e μ= D kBT (1) where μ is the mobility, D is the diffusion coefficient, e is the magnitude of the carrier charge, and T is the temperature, provided that all of the molecules in the crystal are symmetrically related to each other, and there are no imperfections in the crystal.27,53 Diffusion equations can be solved for example by MC simulation45 or by using simplified mean-field methods if the sites are homogeneous.28,54 By ignoring interactions between diffusion pathways, simple onedimensional diffusion equations can be generated, leading to analytical expressions for the mobility in terms of fundamental molecular properties29 that are very useful in aiding synthetic design. Many aspects of the first-principles modeling of anisotropic charge mobilities are yet to be fully clarified, however. While approaches are available for the strong-coupling and weakcoupling limits and there is much debate as to the limit that is most appropriate, do typical experimental conditions satisfy either of these limits? If not, should more general methods55−57 developed for other fields such as exciton transport during photosynthesis be applied instead? The weak-coupling approximation is usually solved assuming the high-temperature limit of Marcus−Hush theory, but is the more general lowtemperature formulation more appropriate?45 Also, the couplings are assumed to be so small that Golden Rule perturbation theory may be used to calculate rates in the

2. THEORY AND METHODS a. Charge Mobilities from Solutions to the Master Equation. The kinetics of charge transport through a material containing many possible residence sites for the charges is given by the master equation47,49,50,82 dpi dt

⎡ ⎤ = −∑ ⎣Wij pi(1 − pj ) − Wji pj(1 − pi )⎦ j

(2)

where Wij is the hopping rate from site i to site j in the crystal, pi is the probability of a charge residing on site i, and 1 − pi reflects the Coulomb penalty factor preventing two or more carriers simultaneously occupying the same site. If an electric field F is applied to the crystal, the charges will drift accordingly, and the mobility can be determined from the drift velocity v as the linear response of the motion to the perturbation: v μ= (3) |F | 14827



Article

Charge mobilities may be obtained either by solving directly for the steady-state solution dpi dt

D=

(10)

where ⟨|(t) − l(0)| ⟩ is the mean-squared displacement of the charge carriers at time t from their initial position. The linear response to an applied field determined from solution of the master equation and the field-free variance of the motion are in general related by the Fluctuation−Dissipation Theorem,86,87 expressed here through the Einstein−Smoluchowski relation, eq 1. MC simulations, different in nature from the kinetic-MC schemes used in section 2a to solve the master equation, are typically used to solve the diffusion equation numerically,45 and we employ the method of Shuai et al.13,45 Also, approximate analytical solutions to the diffusion equation can be generated by projecting each diffusional displacement onto the desired transport direction, collapsing the original multidimensional diffusion process onto an effective one-dimensional diffusion model. This process requires the specification of the residence lifetime for a charge on each site in the crystal, and two methods are currently in common use for this purpose. One method is to use the residence lifetime τF (eq 5) that arises naturally from the semianalytical solution27,48−50,84 of the master equation. The subscript “F” is used to indicate that this lifetime embodies the full set of reactive processes available to each site. However, when considering diffusion from site i to site j, another method is to include only the individual rate process leading to this transition, giving a lifetime of28,29,31,45 2

=0

(4)

of the master equation or else through use of kinetic MC methods to solve the dynamics of this equation given some particular initial conditions. Given a full set of hopping rates Wij, the steady-state charge densities pi can be obtained from the master equation using an efficient iterative procedure.48,82,83 Once these are determined, the solution of the master equation becomes analytical,27,48−50,84 indicating an average residence time of the charge on site i of

τF =

1 ∑j Wij

(5)

and a drift velocity of v=

⟨|l(t ) − l(0)|2 ⟩ 1 lim 2n t →∞ t

∑ Wij pi(1 − pj )R ij ,F ij

(6)

where the hopping distance between all pairs of sites projected onto the direction of the applied electric field vector is

R ij , F =

R ij·F (7)

|F |

τI = 1/Wik

and Rij is the vector from atom i to atom j. The traditional fully numerical route to solving the master equation is through kinetic MC simulation of charge dynamics. In this approach, a random number is used to choose the hopping pathway for a charge, and the drift velocity is estimated from the distance walked in the conduction direction after a very large number of random-walk steps.49,50,85 We implement this procedure following previous MC simulations49 by arbitrarily setting one site i as the starting site for the molecular charge. The probability of hopping of this charge to its near neighbor is given by49

Pik =

In this work, we compare both of these approaches. Use of the full set of processes leads to the common expression for the diffusion coefficient27 DF = =

1 2

∑ pi ∑ R ik2 ,FWik

1 2

∑ R ik2 ,FWik

i

k

(12)

k

for equivalent lattice sites. However, use of the individual process lifetime leads instead to28,29,31,45

Wik ∑j Wij

(8)

DI =

where the sum is over all possible one-step hoppings from site i. In each step of the MC procedure, random numbers are used to choose the destination site k based on all available Pik probabilities. The time δt taken to make this transition is also chosen at random from an exponential distribution about the residence lifetime τF: δt = −τF ln X =

(11)

−ln X ∑j Wij

1 2

1 = 2

∑ pi ∑ R ik2 ,FWikPik i

k

∑ R ik2 ,FWikPik

(13)

k

for equivalent sites. Stehr et al.27 have argued that such expressions are fundamentally flawed when the rate constants Wik are inequivalent as they then embody oscillating chemical reactions that do not lead to actual charge motion. However, Wen et al.29 have utilized eq 13 to derive analytical expressions for anisotropic charge mobilities in general molecular conductors that are of significant use for the design of new materials. We consider whether eq 12 is a more appropriate starting point for such analytical expressions, as this approach would flow naturally from the known solution27,48−50,84 to the master equation, eq 2. Because of the significant approximations inherently associated with the introduction of an effective onedimensional diffusion model, we consider this effect by simply replacing τF by τI during kinetic MC solution of the master

(9)

where X is a random number (0 < X ≤ 1). After a large number of random-walk steps (typically >107), the total distance hopped in the direction of the field d and total elapsed time t are summed, and the drift velocity is evaluated as v = d/t. We repeat each MC simulation 10 times and report the average mobility obtained from eq 3. b. Charge Mobilities from Solution of the Diffusion Equation. Even in the absence of an applied electric field, the charge carriers will thermally diffuse through the conductor with a diffusion constant defined by the variance of the motion 14828



Article

dielectric constants, and so the external contribution is almost always neglected in practice.89−91 Indeed, neglecting this term has been argued to deliver enhanced results for molecular conductors because of a fortuitous cancellation of errors that includes neglect of tunneling contributions to the rate constants.90 We utilize this feature in subsequent rate constant and mobility calculations. Here, we evaluate the internal relaxation using traditional means as21

equation itself. The expression for the random time required to hop from site i to site j, eq 9, then becomes δt = −τI ln X =

−ln X Wij

(14)

c. General System Parameters. All simulations are performed using periodic boundary conditions to simulate an infinite molecular crystal, using a supercell of size 10 × 10 × 10.84 We also use a charge density of 0.001 charges per Alq3 molecule. This is sufficiently low such that the factor (1 − pj) in the master equation can be ignored, as is common practice, although here we fully include it. d. Evaluation of the Microscopic Charge-Transfer Rate Constants. The rate of nonadiabatic intermolecular charge processes can be derived from Fermi's golden rule and the Franck−Condon approximation as40 Wij =

2π |Hij|2 DWFCF ℏ

* * λeint / h = (Ee / h − Ee / h) + (E − E)

where the subscript e/h refers to either molecular anions (for electron transport) or molecular cations (for hole transport), E is the energy of neutral state at its own gas-phase optimized geometry, E* is the energy of neutral state at the gas-phase optimized geometry of the charged species, E*e/h is the energy of the charged state at its own gas-phase optimized geometry, and Ee/h is the energy of the charged state at the gas-phase optimized geometry of the neutral species. Gas-phase geometries are used as this procedure is empirically known to give optimum results in conjunction with other evoked approximations.90 These energies are calculated, as is common practice, by the B3LYP density functional method92 combined with the 6-31G* basis set93 using Gaussian.94 The expected accuracy of the traditional Marcus−Hush equation, eqs 17 and 18, arising from its use of a classical treatment of nuclear notion, is gauged by evaluating the contribution of each molecular vibration to the total internal reorganization energy. This is done by curvilinear projection of the changes in nuclear geometry induced by charge transfer onto the normal modes of vibration, again evaluated using Gaussian at the B3LYP/6-31G* level. Tunneling corrections may be directly evaluated from this data, if required.42−45,64 f. Evaluation of the Intermolecular Couplings Driving Charge Transfer. The couplings Hij need to be evaluated for every pair of molecules in the crystal between which one-step charge transfer is feasible. In the primitive unit cell of the β phase of the Alq3 crystal,65 two molecules related by inversion symmetry reside. We define the molecule whose fractional coordinate of the Al atom is (0.2229, −0.1942, 0.2404) as A(0 0 0) and the other molecule as B(0 0 0), its Al fractional coordinates being (−0.2229, 0.1942, −0.2404). In this notation, A and B refer to the two molecules in the primitive unit cell, while (0 0 0) is the displacement index of the interacting molecule along the a, b, and c crystallographic vectors of the unit cell. The charge-transfer couplings are calculated by a dimer direct-coupling method using the one-electron approximation.95 A variety of similar methods are also available for this purpose.19,22,23,27,29,30,45,72,96 Alternate computationally more expensive multielectron methods based say on the interpretation of calculated physical properties such as transition moments are also in common use, for example, the generalized Mulliken−Hush method of Cave and Newton.97−99 However, much cruder methods that evoke Koopmans' theorem (KT) to interpret splittings between paired molecular orbitals are also in common use,30 and herein, we focus on the significant advantages offered by direct-coupling schemes, schemes requiring similar computational effort to KT. The direct-coupling scheme used herein evaluates the couplings Hij between two localized diabatic orbitals, either monomer HOMO (for hole carriers) or LUMO (for electron carriers) orbitals. These are determined using a fragment-orbital

(15)

where Hij is the effective charge-transfer integral describing the interaction between molecule-localized initial (i) and final ( f) electronic states within the crystal, and DWFCF denotes the density-weighted Franck−Condon factor. The nonadiabatic limit remains valid only if the couplings are sufficiently small that the probability of reacting per molecular collision is low. Also, in the high temperature limit in which the frequencies ω of all vibrations coupled to the charge transfer satisfy

ℏω ≪1 kBT

(16)

so that tunneling can be ignored, this leads to the classical hightemperature Marcus−Hush expression38−41 Wij =

1/2 (ΔEij + λ)2 Hij2 ⎛ π ⎞ ⎜ ⎟ exp − 4λkBT ℏ ⎝ λkBT ⎠

(17)

where ΔEij = Ej − Ei is the free-energy change when the charge is transferred from site i to site j and λ is the reorganization energy released subsequent to each charge transfer process. The reorganization energy is also a measure of the strength of vibronic coupling connecting the electronic and vibrational motions.88 All of these parameters may be evaluated using quantum chemistry calculations. When an electric field is applied, the site energies are modulated by the Stark effect so that the hopping rates become 1/2 Hij2 ⎛ π ⎞ (ΔEij − qR ij·F + λ)2 Wij = ⎜ ⎟ exp − 4λkBT ℏ ⎝ λkBT ⎠

(19)

(18)

where q is the carrier charge. This equation assumes that the electronic coupling Hij is independent of the field strength.85 e. Evaluation of the Reorganization Energies. The reorganization energy is usually partitioned into contributions from the intramolecular vibrations39 (the internal reorganization energy) and the intermolecular phonons38 (the external reorganization energy). The former term is related to the relaxation of molecular geometries during the charge-transfer process, while the latter is related to the relaxation of the polarization of the surrounding medium. While the external contribution dominates most charge-transfer processes in solution and in biological environments, the situation is reversed for most molecular solids because of typically low 14829



Article

Table 1. Contributions λi (in cm−1) > 5 cm−1 from Specific Vibrational Modes of Frequency ωi (in cm−1) to the Internal Reorganization Energy Calculated B3LYP/6-31G* for Alq3 Monomer Neutral and Ionic Species Pertaining to HoleTransfer (HT) and ET Processes (for Full Details, See the Supporting Information)

approach. First, separate density functional theory (DFT) calculations are performed for each isolated molecule to determine the orbitals ϕL/H and ϕL/H on molecules i and j, i j respectively, where the superscript L/H refers to either the HOMO or LUMO orbital, as appropriate. The Kohn−Sham matrix FKS for the molecular dimer is then calculated and transformed into the fragment-orbital basis using TijL/H = ⟨ϕi L/H|FKS|ϕjL/H⟩ =

L/H KS ∑ ciL/H , μ c j , ν Fμν μν

(20)

L/H L/H where ci,μ and cj,v are the fragment molecular-orbital coefficients for atomic basis functions μ or molecule i and ν on molecule j. As the fragment atomic-orbital bases are nonorthogonal, the effective coupling is expressed as100

HijL/H =

TijL/H − SijL/H(TiiL/H − T jjL/H)/2 1 − (SijL/H)2

(21)

where SijL/H =

HT neutral

HT cation

ωi

λi

ωi

λi

ωi

ET neutral λi

ωi

ET anion λi

409 472 553 559 661 663

12.2 7.3 21.6 9.8 10.5 5.8

363 406 421 465 530 562 776

5.7 5.9 7.1 5.4 7.4 9.5 7.0

26 44 409 422 472 553 766 1263 1439 1512 1638

5.7 5.1 5.4 10.8 5.7 6.8 13.3 6.1 6.5 11.5 5.7

363 406 421 465 530 562 776

5.7 5.9 7.1 5.4 7.4 9.5 7.0

L/H ∑ ciL/H , μ c j , ν Sμν μν

reorganization energy, with the associated errors identified as canceling out those arising from use of a classical treatment of the high-frequency intramolecular modes.45 For amorphous Alq3, the intermolecular reorganization energy has been recently estimated at 0.11−0.18 or 0.25 eV, depending upon the computational method.77 It would thus appear that a fortuitous cancellation of errors would indeed occur for Alq3, with the low-frequency intermolecular motions “making up” for the quantization-induced lost portion of the intermolecular reorganization energy. b. Calculation and Significance of Intermolecular Electronic Couplings. The B3LYP/6-31G* calculated intermolecular couplings for the most probable one-step hops are listed in Table 2, where they are compared with values calculated by Lin et al.30 at the Hartree−Fock (HF)/6-31G* level. Lin et al. used two different types of calculations, a simple one based on KT that can be applied only when the charges localize on just one of the three organic ligands in an Alq3 complex, and a more complex direct-coupling calculation (slightly different to ours96) otherwise. As Alq3 is synthesized in its mer isomer,65 the molecule has no symmetry, and the three 8-hydroxyquinolinato ligands each have different bond lengths and angles. The charge therefore may either localize onto the ligand that affords the lowestenergy ion or else, if the interligand coupling is larger than the stability differences for the ion on each ligand, then the ion could delocalize over two or three of the ligands. Simplistic approaches such as KT rely on the orbitals taking the same form in a pair of split energy levels and so can only be applied if the ionic charge is strongly localized on one ligand, whereas direct-coupling schemes are of general applicability; herein, we only apply a direct-coupling scheme. Also, the HF method represents a low-level of theory that overestimates charge localization in situations like this and is rarely applied for quantitative analysis of charge-transfer processes. It is thus expected that our new results would be systematically more accurate and less susceptible to uncontrolled errors. The results presented in Table 2 show that both our results and those of Lin et al.30 depict the same qualitative scenario: the ET couplings are highly anisotropic, attaining values an order of magnitude larger than those for HT. However, the details of the calculated couplings show significant differences

(22)

is the associated transformation of the atomic-orbital overlap matrix Sμν. All calculations are performed at the B3LYP/631G* level.

3. RESULTS AND DISCUSSION a. Calculation and Significance of the Reorganization Energies. The internal reorganization energies λ are evaluated at the B3LYP/6-31G* level and partitioned into contributions λi from each vibrational mode of a neutral Alq3 molecule, and it is appropriate ion using curvilinear coordinates via the DUSHIN program;101 these contributions are related to the Huang−Rhys factors Si for each mode by λi = ℏωiSi

(23)

where ωi are the vibration frequencies. Full results are provided in the Supporting Information, including optimized geometries, energies, normal modes, projected displacements, and mode reorganization energies. The calculated total internal reorganization energies for electrons and holes in Alq3 are 0.269 and 0.234 eV, respectively, results that are quite similar to those from previous calculations.22,30 As compared to reorganization energies for conjugated organic molecules extracted from experimental data, estimates at this and similar levels of theory are usually accurate to within 10%.64,102 The vibrational motions with contributions λi > 5 meV are listed in Table 1 for both hole transfer (HT) and ET. The largest single contribution is just 0.022 eV, and it is clear that a very large number of modes are coupled to both types of charge-transfer process. Most significantly, for both ET and HT, less than 15% of the reorganization energy is associated with modes of frequency less than kBT/ℏ = 206 cm−1. Hence, use of the Marcus−Hush expression, eq 17, derived from classical statistical mechanics, is not valid at room temperature. This is expected to be a general phenomenon for organic molecular conductivity, the effects of which have only recently been quantified;45 it is, however, a well-known effect for photochemical charge separation and recombination processes in solution.64 During mobility calculations for molecular conductors, it is also customary to neglect the effects of the intermolecular 14830



Article

Table 2. Calculated Al−Al Distances R, Electronic Couplings Hij, and Charge-Transfer Rates Wij between Alq3 Molecule i at Crystal Lattice Location A(000) and Molecule j at Various Other Locations as Evaluated Using B3LYP/6-31G* by DC Methods (This Work) and Using HF/6-31G* via Both DC and KT Methods by Lin et al.30 Hij for hole transport (meV) j

R (Å)

B3LYP DC

HF DC

A(±1 0 0) A(0 ±1 0) B(0 0 0) B(1 0 0) B(0 −1 0) B(0 0 0 1) B(1 −1 0) B(1 0 1) B(1 −1 1)

8.44 10.25 8.96 10.07 7.91 8.10 9.11 7.71 11.27

2.43 1.37 0.26 0.005 3.89 7.66 5.89 0.82 15.6

2.46 0.40

HT KT

B3LYP DC Wij (s−1)

Hij for electron transport (meV) B3LYP DC

HF DC

3.09 1.27 188 136 0.78 10.71 6.83 4.3 0.0

0.28 1.18 289 194

0.34 0.041 4.94 3.32 3.13 12.1 15.9

HT KT

hole 2.13 6.77 2.44 5.20 2.20 5.46 1.25 2.79 8.81

0.30 2.08 8.37 3.05 0.014

× × × × × × × × ×

1010 109 108 107 1011 1010 1011 109 1011

electron 2.29 3.86 8.49 4.40 2.75 1.46 1.12 4.43 0

× × × × × × × ×

1010 109 1013 1013 1011 109 1011 1010

that would have profound effects for calculated anisotropic mobilities, and so, we prefer to use the fully self-consistent B3LYP/6-31G* values in subsequent calculations. In a more general context, we note that the use of KT to estimate couplings is highly unreliable, with, for example, the HF/6-31G* values of Lin et al.30 obtained using their directcoupling scheme being considerably more similar to our B3LYP/6-31G* results than are their KT values; see Table 2. While Lin et al. noted some clear examples to which KT was inapplicable and hence used direct-coupling instead, our results highlight other undetected examples of this same phenomenon. For example, the largest charge-transfer coupling calculated for holes arises for channel A(0 0 0)−B(1 −1 1) of 15.6 meV, which is very close to Lin et al.'s KT value of 15.9 meV, whereas our value of 0.82 meV for the channel A(0 0 0)−B(1 0 1) is very much smaller than their KT value of 12.1 meV. Figure 1 Figure 2. Occupied frontier orbitals of the Alq3 dimer A(0 0 0)−B(1 0 1).

Table 3. Relative Excited-State Energies (in eV) for Alq3 Anions and Cations Evaluated Using State-Averaged CASSCF/6-31G* state

cation

anion

1 2

0.27 0.44

0.23 0.47

calculations were performed using MOLPRO103 with the 631G* basis set and involve active spaces including either one electron or one hole in three orbitals (the appropriate ligand HOMO or LUMO orbitals). The first excited states are predicted to be at energies of 0.23−0.27 eV above the ground state. Thermal activation allows excited states to be accessed to facilitate charge-transfer processes, and from the results just discussed for the A(0 0 0)−B(1 0 1) HT channel, the electronic coupling can sometimes be orders of magnitude larger for excited-state processes than for the associated ground-state ones. Such effects can easily mask the effects causes by the need for thermal activation, allowing excited-state conductivity to dominate the mobility through unfavored directions in an organic molecular crystal. As ions often have low-energy excited states, this process is well-known in photochemical charge-recombination reactions,64 but Alq3 would seem to provide a near worst-case example with its three weakly coupled inequivalent 8-hydroxyquinolinato ligands. In particular, this situation would appear to be

Figure 1. HOMO and HOMO-1 orbitals of Alq3 monomer.

shows the HF/6-31G* orbitals for the HOMO and HOMO-1 orbitals of the Alq3 monomer, orbitals identified as localizing on different 8-hydroxyquinolinato ligands. However, Figure 2 shows the occupied frontier orbitals for the A(0 0 0)−B(1 0 1) dimer, and these clearly appear as pairs of mixed monomer orbitals from the previous figure. Hence, the KT method is inappropriate for the evaluation of the intermolecular coupling for this channel. Unfortunately, the problems exposed by the comparison of Figures 1 and 2 are not fully resolved by the simple adaption of either direct-coupling schemes or more advanced methods such as generalized Mulliken−Hush theory.97−99 The observation of such orbital mixing indicates the presence of low-lying excited states of the Alq3 ions. In Table 3, this effect is quantified using complete-active-space self-consistent field (CASSCF) calculations for the relative energies of the ground-states and the first two excited states for both hole and electron transport. These 14831



Article

Table 4. Calculated ET and HT Mobilities for Conduction through Alq3 Crystals in the ab Plane at Angle θ to the a-Axis, Determined by Solving for Either the Drift Velocity at Finite Field Strength F or the Diffusion Constant at Zero Field Using Either Analytical Techniques or Numerical MC Simulations from drift velocity (cm−2 V−1 s−1) HT

from diffusion constant (cm−2 V−1 s−1)

ET

ET

anal.

anal.

MC τF

MC τI

anal.

MC

anal. τF

anal. τI

F = 105 V cm−1

F = 105 V cm−1

F = 105 V cm−1

F = 105 V cm−1

F = 105 V cm−1

F=0

F=0

F=0

θ (°)

eq 7

eq 7

eq 9

eq 14

eq 7

eq 11

eq 12

eq 13

0 15 30 45 60 75 90 105 120 135 150 165

0.0374 0.0228 0.0123 0.0088 0.0129 0.0231 0.0368 0.0500 0.0592 0.0618 0.0571 0.0464

4.023 3.748 3.021 2.006 1.034 0.286 0.022 0.290 1.042 2.033 3.031 3.753

4.034 3.755 3.026 2.025 1.035 0.299 0.037 0.301 1.037 2.034 3.029 3.750

0.403 0.381 0.302 0.202 0.103 0.030 0.004 0.030 0.105 0.204 0.299 0.381

4.250 3.692 2.558 1.420 0.606 0.164 0.030 0.180 0.672 1.572 2.780 3.868

4.030 3.666 3.068 2.004 0.964 0.276 0.000 0.264 1.034 1.996 2.956 3.610

4.023 4.286 5.141 6.363 7.617 8.572 8.975 8.711 7.853 6.634 5.377 4.422

1.816 1.541 1.755 2.400 3.304 4.224 4.916 5.191 4.978 4.332 3.427 2.507

particularly significant for the electron mobility through amorphous Alq3, as these excited-state energies are less than the site-energy perturbations and the large coupling strengths manifest in Table 2 are known to be inhibited by poor geometrical alignment in amorphous Alq3,22,77 therefore making the most efficient reaction pathways directly available only for excited-state processes. Table 2 also highlights a second general problem for the calculation of mobilities in molecular crystals that has not been properly explored: for the most significant channels, the calculated couplings are 136 and 188 meV at the B3LYP/631G* level and correspondingly 194 and 289 meV evaluated HF/6-31G*, large values as compared to the calculated internal reorganization energy (for electron transport) of 269 meV. This large ratio voids the applicability of the weak-coupling approximation but is not be sufficiently large for say the chargedelocalized limit of the Holstein−Peierls model7−13 to apply. Accurate determination of charge-transfer rates in this regime thus requires a full-quantum treatment of the electronic and nuclear motions, motions that are anticipated to be occurring on similar time scales. While a variety of approaches of this type have been developed for exciton transport using say Feynman path-integral formulations,55−57 these methods have not yet been applied to organic molecular conductivity. c. Evaluation and Significance of Charge-Transfer Rate Constants. The field-free charge-transfer rate constants evaluated using eq 17 based on the calculated electronic couplings from Table 2 and use of only the total intramolecular reorganization energies are also given in Table 2. The largest rate constants are predicted to be 8 × 1013 s−1, in excess of the estimated “collision rates” associated with charge transfer. Use of nonadiabatic Golden Rule perturbation schemes such as embodied in eqs 15, 17, and 18 are therefore inappropriate;61 this is a common problem for calculations of mobilities in organic conductors.60 As a result, the largest calculated rate constants will significantly overestimate actual transfer rates and hence mobilities. More general schemes that tread adiabatic and nonadiabatic reactions equivalently are available.58−61 d. Simulation and Significance of Anisotropic Electron Mobilities. Table 4 and Figure 3 show electron mobilities

Figure 3. Anisotropic electron mobility calculated in the ab plane for the β phase of Alq3. The bold black line obtained using the semianalytical solution to the master equation at F = 105 V/cm agrees very well with the red points evaluated from full diffusion equation. The green solid squares and blue solid squares come from solution of approximate one-dimensional diffusion equations.

calculated for conduction through the ab plane of the β phase of Alq3 crystals at angle θ to the a-axis based on rate constants Wij from Table 2, possibly modified to account for finite field strengths. The mobilities are evaluated by up to six different methods that evaluate either the drift velocity at finite field or the diffusion coefficient at zero field. Both formally exact semianalytical and MC approaches are used to calculate the drift velocity, while formally exact MC and approximate analytical methods are used for the diffusion coefficient. The results show that the mobilities determined from the semianalytical solution to the steady-state master equation (eq 6) are in excellent agreement with those obtained from the numerical kinetic MC drift-velocity simulations and the MC diffusion-coefficient simulations (eq 11), provided that internally consistent hopping times τF (eq 9) are used in the kinetic MC simulations, as has been previously noted.27 However, both one-dimensional analytical approximations to 14832



Article

terms of structural inhomogeneities.22,77 In Table 4 is given the anisotropic electron mobilities calculated for the crystal without structural disorder at field strengths of both 105 and 106 V cm−1; more results are given in the Supporting Information. Significant effects are predicted, most noticeably a large reduction of the mobility for θ = 60°. For the direction of highest mobility, θ = 0°, the mobility is predicted to increase slightly from 4.02 cm2 V−1 s−1 at F = 105 V cm−1 to 4.25 cm2 V−1 s−1 at 106 V cm−1. Experimentally,67 the mobility in this direction is known to increase by a similar amount but in much greater proportion, from 0.23 ± 0.03 cm2 V−1 s−1 at F = 0 to 0.32 ± 0.06 cm2 V−1 s−1 at F = 106 V cm−1. It is possible that inhomogeneities even within the crystalline phase contribute to both the observed lower mobility and the enhanced electricfield susceptibility. The physical origin of the calculated large electronic mobility anisotropy, as well as the failure of one-dimensional diffusionequation models to adequately describe this effect, is interpreted using a crystal-structure representation shown in Figure 4. In this figure, the various molecular layers crossing the ab plane are marked, as well as the nearest-neighbor intermolecular couplings from Table 2.

the diffusion coefficient, obtained using either the hopping lifetimes τF that arise naturally from the steady-state master equation solution (eq 12) or the commonly used alternative τI (eq 13), deliver very poor approximations to both the magnitude and the anisotropy of the mobility. Neither of these approaches provide a satisfactory method; hence, it is difficult to select an optimum one. However, the generic failure of analytical one-dimensional approximations to the diffusion equation has been attributed to overcounting of chemical 2 processes because of the appearance of the sign-less Rik,F terms 27 in eqs 12 and 13. Hence, these equations are generally expected to only overestimate mobilities; use of τF overestimates or nearly overestimates the calculated mobilities at all angles, whereas use of τI always leads to lower mobilities than does τF, mobilities that can be much less than the exactly obtained ones. We further investigate the optimum residence time by using τI in place of τF during evaluation of the drift velocity by MC. The results shown in Table 4 indicate a 10-fold reduction in the calculated mobilities resulting from use of the approximate residence lifetime. Indeed, if all sites in the crystal are equivalent, then it is in general possible to show that

νF =N νI

(24)

where νF and νI are the drift velocities obtained from the kinetic MC simulations using eqs 5 and 11, respectively, and N is the total number of reaction processes available to each site. In the present example, N = 10 such processes are described in Table 2, based on a truncation, which only includes processes with Hij > 0.01 meV. Including weaker longer-range direct interactions therefore makes the evaluation of νI numerically unstable, with νI → 0 in the limit that all direct hoppings are included. This is a highly undesired feature. The largest value of the electron mobility calculated for transport with the ab crystal plane is 4.25 cm2 V−1 s−1 for conduction in the a direction. This is an order of magnitude larger than the observed67 maximum anisotropic mobility of 0.23 ± 0.03 cm2 V−1 s−1. Given the likely impact of key approximations such as the neglect of the intermolecular reorganization energy, the use of a classical description of the intramolecular vibrations, the use of the Golden Rule expression for large couplings, and the neglect of excited states and nonpairwise contributions to the intermolecular couplings, agreement to within an order of magnitude is encouraging. This is particularly so given that the observed 4−6 order of magnitude reduction of the electron mobility in amorphous Alq3 is reproduced by similar calculations performed taking full account of the irregular structure of the amorphous phase.22,77 Also shown in Table 4 are mobilities calculated for hole transport in crystalline Alq3. As is expected based on the calculated couplings alone,30 the calculated hole mobilities are 2 orders of magnitude less than those for electron conduction. While this parallels results observed for amorphous Alq3, our calculated mobilities are many orders of magnitude larger for the β phase, and so, analogies are difficult to draw, with details of the anisotropy identified as crucial aspects for the amorphous material.22,77 The hole mobility of the crystalline phase is yet to be measured. Application of electric fields is observed to profoundly influence the mobility of both crystalline67 and amorphous68−71 Alq3. For the amorphous material, this effect is interpreted in

Figure 4. Intermolecular electronic couplings mapped on the ab plane for the β phase of Alq3 crystal. The two molecules in the primate unit cell Alq3 are represented by solid circles (A) and open ellipses (B).

When the electric field is applied along the a vector, the electrons will move from layer a1 to a2 and then to a3. From Figure 4, we see that the fastest ET rate from layer a1 to a2 is through channel A(0 0 0)−B(0 0 0) whose coupling is 188 meV, while the fastest rate from layer a2 to a3 is through channel A(0 0 0)−B(1 0 0) whose coupling is 136 meV. Because the largest electronic couplings between the a1−a2 and the a2−a3 layers are quite large and similar, an electron can readily move from layer a1 to layer a2 then to layer a3 and continuously transfer through the material. When the electric field is applied along the b direction, the electron will transfer from layer b1 to b2 and to b3. The largest coupling from layer b1 to b2 is for channel B(−1 1 0)−A(0 0 0), just 6.8 meV, while those for b2 to b3 are again the strongly coupled channels A(0 0 0)−B(0 0 0) and A(0 0 0)−B(1 0 0). As a result, the net drift velocity will be controlled by the b1 to b2 process and is much 14833



Article

smaller than that along the a direction. However, the onedimensional diffusion equation focuses on the fast “localized diffusion”27 between layers b2 and b3 that does not contribute to the net process, hence severely overestimation the diffusion coefficient in the b direction.

recommended. Simple one-electron based direct coupling (DC) schemes should always be used in preference to this approach, although in general more advanced methods such as generalized Mulliken−Hush theory should be applied, if technically feasible. Finally, we show that the use of one-dimensional diffusion equations leads to very poor results for crystals with strong conduction anisotropies such as Alq3. The combination of tightly coupled and loosely coupled connections leads to prolific localized-diffusion events27 that do not contribute to net charge transport. We also investigate the common practice28,29,31,45 of using an empirical expression for the residence time of a charge on a site in one-dimensional diffusion models. It is shown that they can appear to give improved results through cancellation of errors with the localized-diffusion error. However, when applied to model kinetics without use of one-dimensional diffusion approximations, this approach leads to unstable results, which, in the limit of allowing all sites to directly couple with each other, blocks charge transport completely. Such approximations are thus not advised.

4. CONCLUSIONS Our calculations for the electron mobility of the β crystalline phase of Alq3 predict values 4−6 orders of magnitude larger than those from related calculations22,72−77 for amorphous Alq3, in good qualitative agreement with experimental observations.67−71 These results thus provide the missing link between calculations suggesting that electron mobility in crystalline Alq3 is dominated by some special high-conductivity paths16,30 and calculations suggesting that anisotropy dominates the mobilities of the amorphous phase.22,77 A comprehensive picture describing the electron and hole mobilities through Alq3 materials is therefore established. Many challenges are identified for the full quantitative prediction of mobilities in organic conductors, however. The ratios of the calculated intermolecular couplings to reorganization energies for the most conductive channels are too small for the strong-coupling limit to apply, meaning that nuclear relaxations compete with charge transport on the same time scale. Accurate calculations thus require the coupled solution of both nuclear and electronic motions using say Feynman pathintegration techniques55−57 rather than the application of Marcus−Hush theory.38−41 However, the calculated couplings are so large that the rate constants calculated using nonadiabatic perturbation theory (eqs 15, 17, and 18) exceed likely collision rates, demanding computational methods also appropriate for the adiabatic limit.58−61 Furthermore, partitioning of the intramolecular relaxation energy into contributions from each vibrational mode indicates that high-frequency vibrations dominate charge-transfer relaxation processes, voiding the conditions of applicability of standard Marcus−Hush equation, thus demanding the use of quantum-mechanical treatments for the vibrational motions.42−45 As is customary, we neglect the intermolecular relaxation energy in our treatment, but recent calculations77 indicate that this contribution could well exceed the intramolecular one. While the combination of these two approximations is known to often lead to useful results via the cancellation of two large errors,90 both aspects need to be properly treated in any accurate calculation. Calculations of mobilities in organic conductors usually proceed by ignoring excited states of the ions transporting the charges. This can be a poor approximation in general as organic conductors often have closely spaced energy levels. In particular, we find this to be a very poor approximation for Alq3, as this molecule has its charges stored in one of three possible near-degenerate ligand-localized states per molecule. These effects may be particularly important for electron transport through amorphous Alq3 as the local geometries in the amorphous material block the most highly conductive intermolecular pathways for the ground state, making them instead available only to excited-state processes. The energies of the excited states are also compatible to the energy differences induced by site variations in the amorphous material, making them readily accessible for charge transport. While use of KT is commonplace for the evaluation of electronic couplings between molecules in molecular crystals, this method is found to give results of variable quality and is not

■

ASSOCIATED CONTENT

S Supporting Information *

Optimized coordinates and normal modes of Alq3 in its neutral, cationic, and anionic forms, as well as the projections of the geometry changes on the normal modes and hence mode reorganization energies. Numerical analyses of the MC solution to the diffusion equation, details of the electric field dependence of the electron mobility, and the full form of ref 94. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (S.Y.) or Jeffrey.Reimers@ sydney.edu.au (J.R.R.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank National Science Foundation of China (Grants 20833004 and 21173138) and the Australian Research Council (DP110102932) for financial support.

■

ABBREVIATIONS Alq3, mer-tris(8-hydroxyquinolinato)aluminum(III); DFT, density functional theory; HF, Hartree−Fock; CASSCF, completeactive-space self-consistent field; KT, Koopmans' theorem; DC, direct coupling; ET, electron transfer; HT, hole transfer; HOMO, highest occupied molecular orbital; LUMO, lowest unoccupied molecular orbital

■

REFERENCES

(1) Burroughes, J. H.; Bradley, D. D. C.; Brown, A. R.; Marks, R. N.; Mackay, K.; Friend, R. H.; Burns, P. L.; Holmes, A. B. Nature 1990, 347, 539−541. (2) An, H.; Chen, B.; Hou, J.; Shen, J.; Liu, S. J. Phys. D: Appl. Phys. 1998, 31, 1144−1148. (3) Baldo, M. A.; Lamansky, S.; Burrows, P. E.; Thompson, M. E.; Forrest, S. R. Appl. Phys. Lett. 1999, 75, 4−6. (4) Cölle, M.; Brütting, W. Phys. Status Solidi A 2004, 201, 1095− 1115.

14834



Article

(5) Zaumseil, J.; Sirringhaus, H. Chem. Rev. 2007, 107, 1296−1323. (6) Gunes, S.; Neugebauer, H.; Sariciftci, N. S. Chem. Rev. 2007, 107, 1324−1338. (7) Hannewald, K.; Bobbert, P. A. Appl. Phys. Lett. 2004, 85, 1535− 1537. (8) Hannewald, K.; Bobbert, P. A. Phys. Rev. B 2004, 69, 752121− 7521212. (9) Ortmann, F.; Bechstedt, F.; Hannewald, K. J. Phys.: Condens. Matter 2010, 22, 465802. (10) Ortmann, F.; Bechstedt, F.; Hannewald, K. New J. Phys. 2010, 12, 023011. (11) Wang, L. J.; Peng, Q.; Li, Q. K.; Shuai, Z. J. Chem. Phys. 2007, 127, 044506. (12) Wang, L. J.; Li, Q. K.; Shuai, Z. J. Chem. Phys. 2008, 128, 194706. (13) Wang, L.; Nan, G.; Yang, X.; Peng, Q.; Li, Q.; Shuai, Z. Chem. Soc. Rev. 2010, 39, 423−434. (14) Cheng, Y. C.; Silbey, R. J.; Filho, D. A. d. S.; Calbert, J. P.; Cornil, J.; Bredas, J. L. J. Chem. Phys. 2003, 118, 3764−3774. (15) Troisi, A.; Orlandi, G. J. Phys. Chem. B 2005, 109, 1849−1856. (16) Yang, Y.; Geng, H.; Yin, S.; Shuai, Z.; Peng, J. J. Phys. Chem. B 2006, 110, 3180−3184. (17) Troisi, A.; Orlandi, G. Phys. Rev. Lett. 2006, 96, 086601− 086604. (18) Troisi, A. Adv. Mater. 2007, 19, 2000−2004. (19) Brédas, J. L.; Norton, J. E.; Cornil, J.; Coropceanu, V. Acc. Chem. Res. 2009, 42, 1691−1699. (20) Schein, L. B.; Duke, C. B.; McGhie, A. R. Phys. Rev. Lett. 1978, 40, 197−200. (21) Brédas, J. L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chem. Rev. 2004, 104, 4971−5003. (22) Kwiatkowski, J. J.; Nelson, J.; Li, H.; Bredas, J. L.; Wenzel, W.; Lennartz, C. Phys. Chem. Chem. Phys. 2008, 10, 1852−1858. (23) Brédas, J. L.; Calbert, J. P.; Da Silva Filho, D. A.; Cornil, J. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 5804−5809. (24) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.; Silbey, R.; Brédas, J. L. Chem. Rev. 2007, 107, 926−952. (25) Berlin, Y. A.; Hutchison, G. R.; Rempala, P.; Ratner, M. A.; Michl, J. J. Phys. Chem. A 2003, 107, 3970−3980. (26) Hutchison, G. R.; Ratner, M. A.; Marks, T. J. J. Am. Chem. Soc. 2005, 127, 2339−2350. (27) Stehr, V.; Pfister, J.; Fink, R. F.; Engels, B.; Deibel, C. Phys. Rev. B 2011, 83, 155208. (28) Deng, W.-Q.; Goddard, W. A., III J. Phys. Chem. B 2004, 108, 8614−8621. (29) Wen, S.-H.; Li, A.; Song, J.; Deng, W.-Q.; Han, K.-L.; Goddard, W. A. J. Phys. Chem. B 2009, 113, 8813−8819. (30) Lin, B. C.; Cheng, C. P.; You, Z.-Q.; Hsu, C.-P. J. Am. Chem. Soc. 2005, 127, 66−67. (31) Song, Y.; Di, C.; Yang, X.; Li, S.; Xu, W.; Liu, Y.; Yang, L.; Shuai, Z.; Zhang, D.; Zhu, D. J. Am. Chem. Soc. 2006, 128, 15940−15941. (32) Yang, X.; Li, Q.; Shuai, Z. Nanotechnology 2007, 18, 424029. (33) Yin, S.; Yi, Y.; Li, Q.; Yu, G.; Liu, Y.; Shuai, Z. J. Phys. Chem. A 2006, 110, 7138−7143. (34) Sancho-García, J. C.; Pérez-Jiménez, A. J. J. Chem. Phys. 2008, 129, 024103. (35) Sancho-García, J. C.; Pérez-Jiménez, A. J. Phys. Chem. Chem. Phys. 2009, 11, 2741−2746. (36) Sancho-García, J. C.; Pérez-Jiménez, A. J.; Olivier, Y.; Cornil, J. Phys. Chem. Chem. Phys. 2010, 12, 9381−9388. (37) Marcus, R. A. Rev. Mod. Phys. 1993, 65, 599−610. (38) Marcus, R. A. J. Chem. Phys. 1956, 24, 966−978. (39) Hush, N. S. J. Chem. Phys. 1958, 28, 962−972. (40) Kubo, R.; Toyozawa, Y. Prog. Theo. Phys. 1955, 13, 160−182. (41) Levich, V. G.; Dogonadze, R. R. Dokl. Acad. Nauk. SSSR Ser. Fiz. Khim. 1959, 124, 123−126. (42) Jortner, J. J. Chem. Phys. 1976, 64, 4860−4867. (43) Warshel, A. J. Phys. Chem. 1982, 86, 2218.

(44) Lin, S. H.; Chang, C. H.; Liang, K. K.; Chang, R.; Shiu, Y. J.; Zhang, J. M.; Yang, T.-S.; Hayashi, M.; Hsu, F. C. Adv. Chem. Phys. 2002, 121, 1−88. (45) Nan, G.; Yang, X.; Wang, L.; Shuai, Z.; Zhao, Y. Phys. Rev. B 2009, 79, 115203. (46) Yu, Z. G.; Smith, D. L.; Saxena, A.; Martin, R. L.; Bishop, A. R. Phys. Rev. Lett. 2000, 84, 721. (47) Pasveer, W. F.; Cottaar, J.; Tanase, C.; Coehoorn, R.; Bobbert, P. A.; Blom, P. W. M.; de Leeuw, D. M.; Michels, M. A. J. Phys. Rev. Lett. 2005, 94, 206601. (48) Yin, S.; Lv, Y. Org. Electron. 2008, 9, 852−858. (49) Bässler, H. Phys. Status Solidi B 1993, 175, 15−56. (50) Chatten, A.; Tuladhar, S.; Choulis, S.; Bradley, D.; Nelson, J. J. Mater. Sci. 2005, 40, 1393−1398. (51) Einstein, A. Ann. Phys. (Weinheim, Ger.) 1905, 17, 549−560. (52) Smoluchowski, M. Bull. Int. l'Acad. Sci. Cracovie 1906, 202. (53) Landsberg, P. T. Eur. J. Phys. 1981, 2, 213−219. (54) Datta, A.; Mohakud, S.; Pati, S. K. J. Chem. Phys. 2007, 126, 144710. (55) Huo, P.; Coker, D. F. J. Chem. Phys. 2012, 136, 115102. (56) Huo, P.; Coker, D. F. J. Phys. Chem. Lett. 2011, 2, 825−833. (57) Bonella, S.; Montemayor, D.; Coker, D. F. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6715−6719. (58) Hush, N. S. J. Electroanal. Chem. 1999, 460, 5. (59) Zhao, Y.; Liang, W. Chem. Soc. Rev. 2012, 41, 1075−1087. (60) Shuai, Z.; Wang, L.; Li, Q. Adv. Mater. 2011, 23, 1145−1153. (61) Reimers, J. R.; Hush, N. S. Chem. Phys. 1989, 134, 323. (62) Liu, Y. H.; Xie, Y.; Lu, Z. Y. Chem. Phys. 2010, 367, 160−166. (63) Bisquert, J. Phys. Chem. Chem. Phys. 2008, 10, 3175−3194. (64) Lee, S.-H.; Larsen, A. G.; Ohkubo, K.; Cai, Z.-L.; Reimers, J. R.; Fukuzumi, S.; Crossley, M. J. Chem. Sci. 2012, 3, 257−269. (65) Brinkmann, M.; Gadret, G.; Muccini, M.; Taliani, C.; Masciocchi, N.; Sironi, A. J. Am. Chem. Soc. 2000, 122, 5147−5157. (66) Tang, C. W.; VanSlyke, S. A. Appl. Phys. Lett. 1987, 51, 913− 915. (67) Drew, A. J.; Pratt, F. L.; Hoppler, J.; Schulz, L.; Malik-Kumar, V.; Morley, N. A.; Desai, P.; Shakya, P.; Kreouzis, T.; Gillin, W. P.; Kim, K. W.; Dubroka, A.; Scheuermann, R. Phys. Rev. Lett. 2008, 100, 116601. (68) Naka, S.; Okada, H.; Onnagawa, H.; Yamaguchi, Y.; Tsutsui, T. Synth. Met. 2000, 111, 331−333. (69) Malliaras, G. G.; Shen, Y.; Dunlap, D. H.; Murata, H.; Kafafi, Z. H. Appl. Phys. Lett. 2001, 79, 2582−2584. (70) Murata, H.; Malliaras, G. G.; Uchida, M.; Shen, Y.; Kafafi, Z. H. Chem. Phys. Lett. 2001, 339, 161−166. (71) Fong, H. H.; So, S. K. J. Appl. Phys. 2006, 100, 094502/ 094501−094502/094505. (72) Li, H.; Bredas, J.-L.; Lennartz, C. J. Chem. Phys. 2007, 126, 164704. (73) Nagata, Y.; Lennartz, C. J. Chem. Phys. 2008, 129, 034709. (74) Lukyanov, A.; Lennartz, C.; Andrienko, D. Phys. Status Solidi A 2009, 206, 2737−2742. (75) Difley, S.; Wang, L.-P.; Yeganeh, S.; Yost, S. R.; Voorhis, T. V. Acc. Chem. Res. 2010, 43, 995−1004. (76) Nagata, Y. ChemPhysChem 2010, 11, 474−479. (77) Fuchs, A.; Steinbrecher, T.; Mommer, M. S.; Nagata, Y.; Elstner, M.; Lennartz, C. Phys. Chem. Chem. Phys. 2012, 14, 4259−4270. (78) Johansson, N.; Osada, T.; Stafström, S.; Salaneck, W. R.; Parente, V.; Dos Santos, D. A.; Crispin, X.; Brédas, J. L. J. Chem. Phys. 1999, 111, 2157−2163. (79) Halls, M. D.; Schlegel, H. B. Chem. Mater. 2001, 13, 2632− 2640. (80) Kushto, G. P.; Iizumi, Y.; Kido, J.; Kafafi, Z. H. J. Phys. Chem. A 2000, 104, 3670−3680. (81) Martin, R. L.; Kress, J. D.; Campbell, I. H.; Smith, D. L. Phys. Rev. B 2000, 61, 15804−15811. (82) Yu, Z. G.; Smith, D. L.; Saxena, A.; Martin, R. L.; Bishop, A. R. Phys. Rev. B 2001, 63, 085202. 14835



Article

(83) Pasveer, W. F.; Cottaar, J.; Bobbert, P. A.; Michels, M. A. J. Synth. Met. 2005, 152, 157−160. (84) Yin, S.; Yang, Y.; Lv, Y. Synth. Met. 2010, 160, 1241−1246. (85) Olivier, Y.; Lemaur, V.; Bredas, J. L.; Cornil, J. J. Phys. Chem. A 2006, 110, 6356−6364. (86) Nyquist, H. Phys. Rev. 1928, 32, 110−113. (87) Weber, J. Phys. Rev. 1956, 101, 1620−1626. (88) Reimers, J. R.; Hush, N. S. Chem. Phys. 2004, 299, 79. (89) Brovchenko, I. V. Chem. Phys. Lett. 1997, 278, 355−359. (90) Norton, J. E.; Brédas, J. L. J. Am. Chem. Soc. 2008, 130, 12377− 12384. (91) McMahon, D. P.; Troisi, A. J. Phys. Chem. Lett. 2010, 1, 941− 946. (92) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (93) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257−2261. (94) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (95) Kubar, T.; Woiczikowski, P. B.; Cuniberti, G.; Elstner, M. J. Phys. Chem. B 2008, 112, 7937−7947. (96) Zhang, L. Y.; Friesner, R. A.; Murphy, R. B. J. Chem. Phys. 1997, 107, 450−459. (97) Hush, N. S. Prog. Inorg. Chem. 1967, 8, 391−444. (98) Reimers, J. R.; Hush, N. S. J. Phys. Chem. 1991, 95, 9773. (99) Cave, R. J.; Newton, M. D. Chem. Phys. Lett. 1996, 249, 15−19. (100) Löwdin, P.-O. J. Chem. Phys. 1950, 18, 365−375. (101) Reimers, J. R. J. Chem. Phys. 2001, 115, 9103−9109. (102) Reimers, J. R.; Hush, N. S. J. Chem. Phys. 2003, 119, 3262− 3277. (103) Knowles, P. J.; Lindh, R.; Manby, F. R.; Schütz, M.; et al. MOLPRO, version 2010.1, A Package of Ab Initio Programs; University of Birmingham: Birmingham, 2010.

14836


Challenges for the Accurate Simulation of Anisotropic Charge

Recommend Documents