Influence of Structural Dynamics on Polarization Energies in

McMahon and Troisi also applied this methodology to anthracene, tetracene, ...... Zheng , Saadullah G. Aziz , Veaceslav Coropceanu , and Jean-Luc Bré...
0 downloads 0 Views 1MB Size
Download PDF
20678

J. Phys. Chem. C 2010, 114, 20678–20685

Influence of Structural Dynamics on Polarization Energies in Anthracene Single Crystals† Nicolas G. Martinelli,‡,§,| Julien Ide´,‡,§ Roel S. Sa´nchez-Carrera,|,# Veaceslav Coropceanu,| Jean-Luc Bre´das,‡ Laurent Ducasse,§ Fre´de´ric Castet,§ Je´roˆme Cornil,‡,| and David Beljonne*,‡,| Laboratory for Chemistry of NoVel Materials, UniVersity of Mons, Place du Parc 20, BE-7000 Mons, Belgium, Institut des Sciences Mole´culaires, UMR CNRS 5255, UniVersite´ de Bordeaux, Cours de la Libe´ration 351, FR-33405 Talence, France, and School of Chemistry and Biochemistry and Center for Organic Photonics and Electronics, Georgia Institute of Technology, 901 Atlantic DriVe NW, Atlanta, Georgia 30332-0400, U.S. ReceiVed: June 24, 2010; ReVised Manuscript ReceiVed: September 7, 2010

Normal mode sampling and molecular dynamics simulations are coupled to a valence-bond/Hartree-Fock approach to evaluate the impact of the lattice and molecular vibrations on site energies in anthracene single crystals. The calculations are conducted in the temperature range 0-400 K and show substantial contributions from high-frequency modes, which calls for a quantum-mechanical model even at room temperature. External reorganization energies are also obtained from these modeling studies and found to be much smaller than their internal counterparts. Implications for charge transport in organic single crystals are discussed. 1. Introduction

H)

Over the last twenty years, the field of organic electronics has attracted tremendous interest.1 However, despite the good performance of organic electronic devices, the chargetransport mechanisms are not fully understood at the microscopic level yet. In highly ordered systems (i.e., molecular crystals of high purity at very low temperature) in which the electronic coupling between adjacent molecules is strong, charge transport operates in a band regime where the carrier wave functions are delocalized.2,3 In most practical cases, however, a pure band regime does not apply.4-6 In amorphous materials, charge transport occurs in the hopping regime, in which the charge carriers are localized over a single (or a few) molecule(s) and migrate across the solid through a sequence of hopping events. It is generally considered that a similar mechanism applies in most molecular organic crystals and thin films at room temperature; there, the localized nature of the charge carriers is triggered by both static disorder due to the presence of impurities7 or structural defects8 and dynamic disorder induced by lattice inter- and intramolecular vibrations. Lattice vibrations are responsible for fluctuations in the amplitude of the electronic parameters that govern the chargetransport properties, namely, the site energies of the carriers (i.e., the ionization potential (IP) or electron affinity (EA) of a molecule within the crystal) and the electronic coupling terms (i.e., the transfer integrals) between adjacent molecules, that impact the ability for a charge to hop from a molecule to a neighbor.9 These two parameters define the Hamiltonian for itinerant holes or electrons in the tight-binding approximation: †

Part of the “Mark A. Ratner Festschrift”. University of Mons. Universite´ de Bordeaux. | Georgia Institute of Technology. # Present address: Research and Technology Center North America, Robert Bosch LLC, Cambridge, Massachusetts 02142 ‡ §

∑ a+n anεn(t) + ∑ am+anJmn(t) n

(1)

m*n

Here an+/an are the creation/annihilation operators for an electron (hole) on molecular site n, Jmn is the transfer integral between sites m and n, and εn is the site energy. Both Jmn and εn depend on time t as described below. The impact of the dynamical disorder on the electronic couplings (referred to as nonlocal or Peierls electron-vibration interaction) has been investigated by several groups who followed various computational strategies.9-13 Recent studies on naphthalene single crystals have pointed to significant interactions with both inter- and intramolecular vibrations and have emphasized the role of zero-point fluctuations.14,15 The site energies, being local in nature, are also significantly affected by crystal vibrations. Their fluctuations with time, corresponding to what is known as a local Holstein-type16,17 electron-vibration mechanism, have received much less attention. In the case of holes, the site energy is quantified by the ionization potential of the crystal (IPcry, defined as the difference between the relaxed energy of the positively charged system and the relaxed energy of the neutral system). Similarly, in the case of electrons, the site energy is defined by the electron affinity of the crystal (EAcry, obtained by considering a negatively charged system): 0 IPcry ) E+ cry - Ecry 0 EAcry ) Ecry - Ecry

(2)

The modulations of the on-site energy as a function of the lattice vibrations can be, to a first approximation, divided into an intramolecular component and an intermolecular component. The intermolecular component, commonly referred to as the polarization energy, arises from the interaction of localized electrons (holes) with the induced and permanent multipole moments of surrounding molecules. The polarization energy can be quantified by the difference between the crystal and gas phase ionization potentials/electron affinities, IPcry and IPgas/EAcry and

10.1021/jp105843t  2010 American Chemical Society Published on Web 09/27/2010

Polarization Energies in Anthracene EAgas, respectively. While in most recent studies11,18-22 the emphasis was placed on the internal contribution to the local electron-vibration interaction, which arises from the dependence of the intramolecular component of the site energy on intramolecular vibrations, the external contribution to the local electron-vibration interaction associated with the modulation of polarization energy by lattice vibrations has been less studied. However, a quantitative assessment of this contribution is needed for a reliable theoretical description of charge-carrier mobilities in organic materials23 and is the focus of this work. Recently, several computational approaches have been worked out to compute the polarization energy. For example, Tsiper and Soos have used a mixed quantum/classical model to investigate the polarization energies of anions in anthracene and 3,4,9,10-perylen-tetracarboxydianhydride (PTCDA);24 in their approach, individual molecules are treated at a quantummechanical level and are subjected to classical external fields created by all other molecules; the stabilization energies of the neutral lattice, of singly charged systems, and of systems in which charge-transfer pairs are present were calculated. Verlaak et al. have investigated similar effects in anthracene and pentacene single crystals using a microelectrostatic approach in which molecules are coarse-grained into electrical multipoles;25 they studied the effect of structural defects (grain boundaries) on the electronic state distribution and the energy levels of a hole on molecules at and near the defect. Kwiatkowski et al. have evaluated the effect of vibronic coupling on the charge-transport parameters (transfer integrals and site energies) in crystalline naphthalene by means of density functional theory (DFT) methods, treating the vibrations at both classical and quantum levels.14 The magnitude of the polarization energy for localized charge carriers has further been evaluated recently by Norton and Bre´das in oligoacenes using hybrid quantum-mechanics/molecular mechanics (QM/MM) schemes with a polarizable force field as well as QM methods with a polarizable continuum model;26 the polarization energies for anionic and cationic naphthalene clusters, as well as the IP, EA, and reorganization energies were estimated. A similar methodology has been recently applied by McMahon and Troisi to evaluate the external reorganization energies for oligoacenes.27 Very recently, Nagata has combined polarizable molecular dynamics simulations to a charge response kernel method to evaluate site energy disorder in amorphous Alq3.28 We also recently proposed a quantum-chemical approach based on a fragment orbital (valence bond/Hartree-Fock, VBHF) formalism to evaluate the polarization energies and transfer integrals in anthracene clusters.29 The modulation of the polarization energy by lattice vibrations was first discussed by Gosar and Choi,30 Vilfan,31 and Friedman32 using a model that accounts only for the induced dipole moments and interaction with acoustic phonons. This model provided an estimate of the lattice relaxation energy of about 40 meV in anthracene. A later study by Brovchenko based on a force-field approach resulted in a smaller lattice relaxation energy of about 15 meV.33 A yet smaller relaxation energy of about 2.5 meV was derived by McMahon and Troisi27 using a polarizable force field. For the case of the external relaxation energy for holes in naphthalene, both Norton and Bre´das and McMahon and Troisi computed relaxation energies values of about 8 meV. McMahon and Troisi also applied this methodology to anthracene, tetracene, pentacene, and rubrene. In the present work, we have combined molecular dynamics (MD) simulations with VB/HF calculations to assess the impact of lattice fluctuations on site energies in the anthracene single

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20679 crystal. In addition, we have also investigated the contribution of optical vibrations (phonons) to the external relaxation energy. 2. Theoretical Methodology In the harmonic approximation and when only the linear electron-vibration interactions are considered, the Hamiltonian of the lattice vibrations and the dependence of site energy on normal-mode coordinates are given by15

H)

1 2

∑ pωj(Pj2 + Qj2)

(3)

j

ε ) ε0 +

∑ εjQj + ...

(4)

j

where ωj, Qj, and Pj denote the frequency and the dimensionless coordinates and momentum of vibrational mode j, respectively; ε0 is the site energy at the equilibrium geometry, and the εj terms represent the linear electron-phonon coupling constants. The relaxation energy is calculated on the basis of εj as9

λrel )

ε2

j ∑ 2pω j

(5)

j

The determination of the coupling constants εj further allows us to calculate the temperature dependence of the variance in site energies. The average site energy and its variance are given by15,34

〈ε〉QM ) ε0

σQM2 )

∑ j

εj2 pωj coth 2 2kBT

(6)

(7)

The QM subscript indicates that the vibrations are treated at the quantum-mechanical level. The classical limit, or hightemperature limit, yields34

〈ε〉Cl ) ε0 σCl2 )

k T

B ) 2λrelkBT ∑ εj2 pω j

(8)

(9)

j

The reorganization energy, used in the framework of Marcus theory for hopping transport,35 describes the energetic barrier associated with the nuclear relaxation of the molecules during a charge-transfer process; it can be approximated as twice the relaxation energy given in eq 5.9 The contribution of the molecules involved in the charge-transfer process, i.e., the internal reorganization energy, can be calculated using eq 5 from the intramolecular modes of a single molecule as well as by using a method based on adiabatic potential energy surfaces, where the internal relaxation energy is calculated as36

λint,adiab ) E(1)(M) - E(0)(M) + E(1)(M•+) - E(0)(M•+) (10)

20680

J. Phys. Chem. C, Vol. 114, No. 48, 2010

Martinelli et al.

Figure 1. View of the 2D anthracene cluster used for the VB/HF calculations. At the bottom right, the chemical structure of anthracene is displayed.

Here, E(0)(M) and E(0)(M+•) denote the ground-state energy of the neutral state and of the charged state, respectively; E(1)(M) is the energy of the neutral molecule in the optimized geometry of the ion, and E(1)(M+•) is the energy of the charged state at the optimized geometry of the neutral molecule. The consistency between these approaches has been demonstrated.22 We consider here the modulation of IP and EA by lattice fluctuations in the anthracene single crystal. Classical MD simulations were performed on molecular aggregates, with geometric structures saved at regular time steps along the MD trajectory. Single-point VB/HF calculations of site energies were then carried out on the basis of the MD snapshots. The distributions of site energies derived from these calculations were fitted to a Gaussian function to estimate the standard deviations and the related reorganization energies by means of eq 9. Note that, in these calculations, the excess charge is assumed to be localized on a single molecule. The delocalization of the charge has been assessed by Hultell and Strafstro¨m for the pentacene single crystal37 and was estimated to be roughly 4 molecules at low temperature for a transfer integral of 50 meV, while neglecting the nondiagonal disorder. In the case of anthracene, the electron-phonon coupling is higher (the reorganization energy is 50% larger for anthracene than pentacene),22 which therefore should further reduce the polaron size. A pronounced charge localization at room temperature is also predicted by advanced theories of charge transport in the case of anthracene.38 Thus, localization of the charge on a single molecule, albeit being a limiting case, appears to be a reasonable approximation. Moreover, the approach used to calculate the linear electron-phonon couplings is general and does not assume any charge localization. The local electron-phonon couplings were estimated only for the optical modes using the normal-mode coordinates derived at the Γ-point. In this case, the unit cells were displaced along

each normal mode and single-point VB/HF calculations of site energy were performed at each geometrical configuration; the electron-phonon coupling constants were then computed by means of numerical differentiations. We refer to the displacement of the unit cell along each normal mode as normal-mode sampling (NMS). We note that the coupling constants in general depend on the phonon wave vector; thus, the present results only provide a first insight into the issue.39 The lattice dynamics at the Γ-point were investigated: (i) by means of MM simulations using the TINKER package40 and the MM3 force field;41 and (ii) by the computation of the lattice vibrations by means of quantum-chemical calculations using the VASP code,42 with the Perdew-Burke-Ernzerhof exchange-correlation functional, a plane-wave basis set (300 eV cutoff) and projector augmented wave potentials.43,44 These selfconsistent calculations were carried out with an 8 × 8 × 8 k-point mesh. While the quantum-chemical geometry optimizations were performed by constraining the cell parameters to the experimental values (a ) 8.414 Å, b ) 5.990 Å, c ) 11.095 Å, and β ) 125.29°),45 during the MM geometry optimization the monoclinic unit cell was allowed to relax, resulting in the following parameters: a ) 8.379 Å, b ) 6.001 Å, c ) 11.144 Å, R ) 90°, β ) 124.19°, and γ ) 90°. VB/HF calculations based on the semiempirical Austin Model 1 (AM1) method46 of site energies were then carried out on clusters of 25 molecules (see Figure 1) determined from NMS. This cluster size is chosen as it offers a good compromise between results accuracy and computational time.29 Despite the fact that anharmonic terms of the potential energy surface and the dispersion of the vibrational modes were not considered in this study, the NMS approach has the advantage to provide detailed information on the contribution of each normal mode to the overall fluctuation of site energies. Moreover, the thermal fluctuations can be quantified at any

Polarization Energies in Anthracene

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20681

temperature, in both quantum and classical limits. To avoid any confusion, we will refer to the normal mode sampling in the quantum-mechanical approach as NMSQM and to the normal mode sampling in the classical limit as NMSCl. The MD simulations were conducted on three molecular layers, each layer containing 25 unit cells, using the TINKER package40 and the MM3 force field,41 whose parameters are fitted to structural and vibrational data for hydrocarbon molecular crystals.47 All simulations were carried out in the canonical ensemble (NVT) at 100, 200, 300, and 400 K, using the Andersen thermostat48 and the Beeman algorithm.49 The canonical ensemble was chosen to fix the cell dimensions, as done in the NMS approach. After a thermalization period of 200 ps, 8000 snapshots of the structure were extracted, at a rate of 0.2 fs-1. VB/HF calculations of site energies were then carried out on the structures issued from the dynamics, after scaling the size of the clusters down to 25 molecules as in the NMS approach. Finally, Fourier transforms (FT) of the fluctuating IP and EA were calculated with the R-project package50 to investigate the various vibrational components of the fluctuations. 3. Results and Discussion 3.1. Normal Mode Sampling. We first assessed the validity of the harmonic approximation to describe the potential energy surface and the linear character of the fluctuations in site energies along the normal modes. To do so, we performed single-point calculations of the potential energy at the classical level as a function of the distortions of the molecules in the unit cell along each normal mode. The curve so obtained was then fitted to eq 3, with a momentum equal to zero. The value of the force constant obtained with this fit should result in a ground-state energy of pω/2. These results are collected in Table S1 (see Supporting Information). All the normal modes can be considered harmonic in view of the small error corrections on the fits, and the fitted and calculated energies show also a good agreement. The evolution of IPcry and EAcry has then been calculated for a cluster of 25 molecules along each mode; a linear fit has been performed and the correlation coefficient R2 calculated. Table S2 reports the correlation coefficient for each mode as well as max min max min - IPcry , where IPcry [IPcry ] is the maximum ∆IP ) IPcry [minimum] IP calculated along the distortion for a given normal mode. The same has been done for EA. The modes leading to high ∆IP and ∆EA values display mostly linear electroncoupling constants (R2 ≈ 1). Among the modes leading to large IP fluctuations (i.e., larger than 0.01 eV), three modes at 41, 47, and 235 cm-1 do not show a strict linear behavior (with linear correlation coefficients ranging from 0.54 to 0.83), and their influence has been neglected hereafter. This approximation will be better justified later, as the results show that lowfrequency modes have a much weaker impact than highfrequency modes. In the case of a negative charge, all modes leading to large EA fluctuations are strictly linear and are therefore taken into account. Figure 2a plots the evolution of the standard deviation of the electron-phonon coupling as a function of temperature, from 0 to 400 K, for both NMSQM and NMSCl approaches. In the case of the NMSQM approach at 0 K, the average value is equal to 6.99 [-1.55] eV for IP [EA] and the standard deviation is 0.17 [0.14] eV, respectively. Despite differences between the absolute value of the averages, the standard deviations are of the same order of magnitude. σQM remains almost constant in the temperature range investigated. This is explained by the fact that the main contributions entering summations (7) and (9) at

Figure 2. Evolution of the standard deviation calculated for IP and EA from eqs 7 and 9, using the TINKER eigenvectors for the normal mode sampling. (a) All molecules in the cluster are distorted along the normal modes. (b) Only the charged molecule is distorted; all the surrounding molecules are kept fixed. (c) The charged molecule is kept fixed, only the surrounding molecules are distorted.

300 K are two high-frequency modes at around 1670 cm-1 (C-C stretching), as shown in Table S3. These two modes represent about 70% of the variance calculated at the quantummechanical level for electrons or holes. Other modes at 1280 and 1554 cm-1 (both C-C stretching modes as well) also play a significant role. Since these are vibrational modes characterized by quantum energies much higher than kT at room temperature (∼200 cm-1), the thermal energy is not sufficient to significantly perturb the population distribution at 0 K, so that these modes are in the vibrational ground state. In the case of the NMSCl approach, as expected, σCl vanishes at 0 K and increases with temperature; we also note that σCl reaches only half the value of σQM at 300 K. Again, this overall discrepancy over the temperature range is due to the role of the high-frequency modes at nearly 1670 cm-1. Note that the calculations described above do not allow us to distinguish between the contributions to the overall electronphonon coupling coming from distortions of the charged molecule or from the surrounding neutral molecules. To address this issue, two additional NMS have been performed. In the first case, only the charged molecule can distort while the surrounding molecules are kept frozen (Figure 2b). In the second case, the charged molecule is frozen and the surrounding molecules are distorted (Figure 2c). It appears that the fluctuations of the charged molecule alone accounts for the same contribution to the total variance, as that resulting from the fluctuations of the entire cluster. The detailed analysis of eq 7 shows that the mode contributing the most is the vibration at 1676 cm-1, which provides around 75% of the total variance for both holes and electrons. The fluctuations due to the surrounding molecules lead to a standard deviation of only 0.04 eV. In this case, the quantum-mechanical and classical samplings provide similar variances at room temperature; this is consistent with the fact that these fluctuations are mediated by lowfrequency modes. For holes, modes at 688 (23%), 3074 (13%),

20682

J. Phys. Chem. C, Vol. 114, No. 48, 2010

Martinelli et al. TABLE 1: Comparison of Standard Deviations Calculated from the NMS (Quantum-Mechanical Approach and Classical Approximation) and MD Simulations, at Various Temperatures IP

Figure 3. Evolution of the standard deviation calculated for IP and EA from eqs 7 and 9, using the VASP normal modes. All molecules in the cluster are distorted along the normal modes.

974 (11%), 684 (7%), 704 (6%), and 138 (5%) cm-1 are the most strongly coupled vibrations. In the case of electrons, similar modes at 138 (17%), 688 (17%), 704 (8%), 974 (7%), 1676 (6%), and 3074 (6%) cm-1 contribute the most to eq 7. It is worth stressing that the intramolecular modes of the surrounding molecules (even some C-H vibrations) influence the electronphonon coupling, as evidenced upon freezing the central molecule. Note also that the variance associated to the entire system is not simply the sum of that estimated for the two limiting cases; the VB/HF procedure indeed optimizes the electronic cloud taking into account the environment of each molecule, which depends on the degree of geometric distortion of the molecules. To validate the normal modes calculated with the force field, the normal modes have also been calculated at a quantumchemical level using a DFT approach as implemented in the VASP code. Figure 3 plots the evolution of the standard deviation of the electron-phonon coupling as a function of temperature using the DFT-calculated modes. When compared to Figure 2a, the plots are found to be very similar, which validates the evolution obtained with the force-field approach. A mode-by-mode analysis from the DFT results shows that there are four important modes contributing to the fluctuations of IP and EA: 2 modes at around 1430 cm-1 and 2 modes at 1560 cm-1. They represent 75% of the variance for holes and 67% for electrons. These modes (C-C stretching) are very similar and, after careful analysis, no striking differences were noticed between these modes and the two modes at 1676 cm-1 obtained with the force field (which represents 70% of the variance). The results provided by the two theoretical approaches are thus consistent and justify the use of the force-field approach. 3.2. Molecular Dynamics Simulations. The average values and standard deviations of IP and EA have been estimated at various temperatures with data derived from the MD simulations. A slight change in the average values is observed as a function of temperature, ranging from 6.98 to 6.96 eV for IP and -1.56 to -1.60 eV for EA; these can be considered as negligible when compared to the standard deviation. The latter is reported in Table 1; its evolution with temperature is very similar to the standard deviation obtained with the NMSCl approach. The difference between the MD and NMSCl approaches increases with temperature (as does the nonlinear behavior of the fluctuations) but remains modest, which confirms the validity of the linear electron-phonon coupling approximation. FTs have been performed on the IP and EA traces provided by the MD run to quantify the contribution of each mode to the fluctuations in IP and EA. The FTs at 300 K are plotted in

EA

T

NMSQM

NMScl

MD

NMSQM

NMScl

MD

100 200 300 400

0.169 0.169 0.169 0.170

0.051 0.072 0.088 0.102

0.052 0.075 0.099 0.112

0.140 0.141 0.142 0.144

0.045 0.064 0.078 0.090

0.049 0.071 0.092 0.101

Figure 4a,b and are similar to the results obtained with the NMSCl approach: the modes at around 1670 cm-1 contribute much more significantly than any other mode, and their contributions are related to the fluctuations of the charged molecule. For the sake of comparison, rigid-body calculations have also been performed at 300 K, by extracting structures every 10 fs (Figure 4c,d). An average of 6.99 [-1.56] eV is found for IP [EA], which is similar to the values obtained with a fully relaxed MD. However, the standard deviation is much lower: 0.02 eV for IP and 0.04 eV for EA, instead of 0.10 and 0.09 eV in the fully relaxed cases at 300 K; this again points to the smaller influence of intermolecular modes. 3.3. Calculation of Reorganization Energies. The internal reorganization energy is easier to estimate at the theoretical level than the external reorganization energy since the former requires calculations on a single molecule and the latter on large clusters of molecules. The internal and external reorganization energy parameters have recently been estimated by Norton and Bre´das26 and by McMahon and Troisi27 using polarizable force fields, since this approach allows an excess charge to be localized on a molecule in a molecular cluster. The VB/HF procedure also allows charge localization on molecules, while taking into account explicitly the electronic cloud of all surrounding molecules and their relaxation as an excess charge is brought to the system. The linear couplings calculated in section 3.1 can thus be used to estimate the total reorganization energy via eq 5. If we subtract the internal reorganization energy from this quantity, the external contribution is obtained. Note that we were not able to calculate the total reorganization energy using the method based on adiabatic potential energy surfaces (eq 10) since geometry optimization has not yet been implemented in the VB/HF procedure used here (only molecular orbital relaxation is allowed). The linear electron-phonon couplings calculated with the force field have been used to calculate the total (λtotal) and internal (λint) reorganization energies for holes and electrons (Table 2). A good agreement between the MM and DFTcalculated internal and external reorganization energies is observed (see Table 2). On the basis of eq 9 and assuming that λtotal ) 2λrel,9 λtotal can also be estimated from the variances calculated using the MD approach (Table 1). A total reorganization energy of 346 [299] meV is found for IP [EA], which is slightly higher than the value obtained from the linear electron-phonon coupling constants (303 and 234 meV for IP and EA, respectively). The difference between the values obtained from linear electron-phonon coupling constants and variances obtained from the MD simulations can be explained by the fact that in the first approach, nonlinear coupling constants are not taken into account, which is the case in the MD approach. Thus, eq 9 is not strictly correct in this case since the (limited) nonlinearity is not included on the right-hand side of the equation. The external part of the reorganization energy

Polarization Energies in Anthracene

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20683

Figure 4. Fourier transforms calculated from the MD simulations: (a) and (b) fully relaxed dynamics; (c) and (d) rigid-body constraints. The intensity I is scaled in the range 0 to 1 on each plot.

TABLE 2: Reorganization Energies Calculated from the Linear Couplingsa hole (meV)

electron (meV)

TABLE 3: Reorganization Energies Calculated from the Linear Couplings Using Two Different Partition Methods (Pintra-inter and Pfreq)a hole (meV)

TINKER λtotal λint λext

303 259 44

234 182 52

VASP λtotal λint λext λint,adiab λext,QM/MM27

335 298 37 245 5

264 219 45 187 5

a

λtotal is the total reorganization energy, λint the internal part, and λext the external part (on a cluster of 25 molecules). λint,adiab is the internal part calculated from the adiabatic potential energy curves and λext,QM/MM the external part calculated using a QM/MM approach.

λext ) λtotal - λint is estimated to be 44 and 52 meV for holes and electrons, respectively, which is much lower than the internal part, estimated (using eq 5) to be 259 and 182 meV for holes and electrons. The internal reorganization energies calculated from the adiabatic potential energy surfaces (λint,adiab) evaluated at the AM1 level using the AMPAC 9 package51 amount to 245 and 187 meV for holes and electrons. This further demonstrates the consistency between the two approaches. Note that DFT yields a similar internal reorganization energy for electrons (196 meV at the B3LYP level with the 6-31G** basis set), yet a smaller value for holes (137 meV at the same level).22 The partition of the total reorganization energy λ considered above, i.e., with λint corresponding to the reorganization energy of the two molecules involved in the charge transfer process and λext corresponding to the contribution of the surrounding molecules (this approach will be referred to as Pintra-inter), is ubiquitous in Marcus electron transfer theory but not unique. An alternative partition scheme is based on the frequencies of the vibration normal modes, being of intra- or intermolecular nature. In the latter picture, all contributions from the highenergy vibrational modes (with zero point energy larger or similar to kT) are included in λint, while the low-energy

λint (meV) λext (meV) k (s-1)

electron (meV)

Pintra-inter

Pfreq

Pintra-inter

Pfreq

259 44 1.29 × 1013

299 4 5.14 × 1013

182 52 1.61 × 1013

215 19 3.11 × 1013

a k is the corresponding Marcus-Levich-Jortner rate calculated at room temperature with a transfer integral of 30 meV.

vibrational modes yield λext (we note this partition Pfreq). For electrons, using Pfreq instead of Pintra-inter leads to a decrease of λext from 52 to 19 meV, while λint increases from 182 to 215 meV. For holes, the λextvalue drops from 44 to 4 meV, and the λint increases from 259 to 299 meV (see Table 3). As mentioned earlier, external reorganization energies of (2 × 2.5)) 5 meV for holes27 and of (2 × 15)) 30 meV33 have been estimated in anthracene on the basis of QM/MM calculations and the sole consideration of the intermolecular modes, respectively. These external reorganization energies are consistent with the energies calculated in this work using the Pfreq partition but smaller if Pintra-inter is used instead. In the framework of the Marcus-Levich-Jortner (MLJ) formalism, the transfer rate for a charge jumping from one molecule to another is52 ∞

kif )

1 2π 2 Sn exp(-S) × Vif p n! √4πλextkBT n)0

∑

[

exp -

(λext + npω)2 4λextkBT

]

(11)

with kB the Boltzmann constant, T the temperature, Vif the electronic coupling between the two interacting molecules, S the Huang-Rhys factor defined as S ) λint/pω, pω the energy of an effective mode assisting the charge transfer, and n the vibrational quantum number. We can now assess the influence of the partition scheme on the resulting transfer rates by plugging

20684

J. Phys. Chem. C, Vol. 114, No. 48, 2010

the two set of λint and λext values in eq 11. We find that, if Pfreq is used instead of Pintra-inter, the MLJ transfer rate increases by a factor of 2 for electrons and 4 for holes (when considering room temperature and a typical value of 30 meV for the transfer integral;53 see Table 3). This result shows that substantial variations in the transfer rates (and thus charge-carrier mobilities) result from subtle changes in the partitioning of the total reorganization energy in the MLJ expression and calls for a more general electron transfer theory. 4. Conclusions We have combined the VB/HF methodology to a normalmode sampling approach to estimate the linear coefficients of the electron-phonon couplings in the anthracene single crystal. The evolutions of the standard deviations of these couplings have been calculated between 0 and 400 K within quantummechanical and classical frameworks. The quantum-mechanical approach shows a nearly constant value of the standard deviation, due to dominant contributions from high-frequency modes. The classical approximation underestimates the standard deviation by a factor of 2 at high temperature for the same reason. A good agreement is found between the standard deviations of IP and EA calculated from MD simulations and the NMSCl approach; this confirms that the nonlinear behavior of the electron-phonon coupling is rather limited. However, the underestimation of the standard deviation calculated using the NMSCl approximation in comparison to the value obtained with the NMSQM approach points out that classical samplings do not provide reliable estimates of the impact of dynamical fluctuations on the site energies in these systems, even at room temperature. We have also derived the reorganization energies from the linear coupling constants of IP and EA and compared them to values extracted from adiabatic potential energy surfaces. The internal reorganization energies computed with both approaches are in good agreement, which supports the use of the linear approximation for the calculation of the electron-phonon couplings. External reorganization energies calculated for a cluster of 25 molecules are 2 orders of magnitude smaller than the internal part for holes, and one order smaller for electrons. We have further shown that these energies depend on the partition scheme used to disentangle external and internal contributions, based on either the nature of the vibrational modes, Pintra-inter (with modes involving the molecules that undergo electron transfer reaction contributing to the internal part and modes that involve the surrounding molecules contributing to the external part) or their energies, Pfreq (with highenergy vibrations retained in the internal part and low-energy vibrations in the external part). Our values of the external reorganization energy are consistent with those obtained from other theoretical works in the Pfreq scenario but larger in the Pintra-inter case. The different partition schemes lead to significant changes in the values of the charge transfer rates obtained within the framework of the Marcus-Levich-Jortner theory for hopping. An accurate estimate of the external reorganization energy and its contributions is thus critical to quantitatively describe the charge-transport properties in organic materials. The subtle effects revealed here in the dynamic fluctuations of the electronic properties (namely IP and EA) can hardly be accounted for within the current electron transfer theories and require more elaborate theoretical frameworks.54-56 Acknowledgment. We acknowledge Professors Mark Ratner and Robert Silbey for many fruitful discussions and collabora-

Martinelli et al. tions over the recent years. The work in Mons was partly supported by the European project ONE-P (NMP3-LA-2008212311) and FNRS. The work at Georgia Tech was primarily supported by the National Science Foundation under the MRSEC Program (Award DMR-0819885), as well as partly under the STC Program (Award DMR-0120967). J.C. and D.B. are FNRS Research Fellows; N.M. acknowledges a grant from “Fonds pour la Formation a` la Recherche dans l’Industrie et dans l’Agriculture (FRIA)”; J.I. is grateful to the “Advanced Materials in Aquitaine” program (www.ama-materials.com) for his Ph.D. grant. Supporting Information Available: The fitting details of the harmonic potential energy and the comparison to the expected value are provided for each normal mode in Table S1. A complete list of the linear couplings, maximum energy differences at 300 K, and linear correlation coefficients is provided in Table S2. Finally, the percentage of each mode in the calculation of the variance at 300 K is provided in Table S3. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Hanbook of conducting polymers, 3rd ed.; Skotheim, T. A., Reynolds, J. R., Eds.; CRC Press: Boca Raton, FL, U.S., 2007. (2) Glaeser, R. M.; Berry, R. S. J. Chem. Phys. 1966, 44, 3797. (3) Karl, N. Synth. Met. 2003, 133, 649. (4) Cheung, D. L.; Troisi, A. Phys. Chem. Chem. Phys. 2008, 10, 594. (5) Troisi, A.; Orlandi, G. J. Phys. Chem. A 2006, 110, 4065. (6) Troisi, A.; Orlandi, G. Phys. ReV. Lett. 2006, 96, 086601. (7) Forrest, S. R. Nature 2004, 428, 911. (8) Kwiatkowski, J. J.; Nelson, J.; Li, H.; Bre´das, J. L.; Wenzel, W.; Lennartz, C. Phys. Chem. Chem. Phys. 2008, 10, 185. (9) Coropceanu, V.; Cornil, J.; da Silva, D. A.; Olivier, Y.; Silbey, R. J.; Bre´das, J. L. Chem. ReV. 2007, 107, 2165. (10) Lemaur, V.; da Silva Filho, D. A.; Coropceanu, V.; Lehmann, M.; Geerts, Y.; Piris, J.; Debije, M.; Van de Craats, A.; Senthilkumar, K.; Siebbeles, L. D. A.; Warman, J.; Bre´das, J. L.; Cornil, J. J. Am. Chem. Soc. 2004, 126, 3271. (11) Sa´nchez-Carrera, R. S.; Coropceanu, V.; da Silva Filho, D. A.; Friedlein, R.; Osikowicz, W.; Suess, R. C.; Salaneck, W. R.; Bre´das, J. L. J. Phys. Chem. B 2006, 110, 18904. (12) Wang, L. J.; Peng, Q.; Li, Q. K.; Shuai, Z. J. Chem. Phys. 2007, 127, 044506. (13) Lemaur, V.; Bouzakraoui, S.; Olivier, Y.; Brocorens, P.; Burhin, C.; El Beghdadi, J.; Martin-Hoyas, A.; Jonas, A.; Serban, D. A.; Vlad, A.; Boucher, N.; Leroy, J.; Sferrazza, M.; Melinte, S.; Sergeev, S.; Geerts, Y.; Lazzaroni, R.; Cornil, J.; Nysten, B. J. Phys. Chem. C 2010, 114, 4617. (14) Kwiatkowski, J. J.; Frost, J. M.; Kirkpatrick, J.; Nelson, J. J. Phys. Chem. A 2008, 112, 9113. (15) Coropceanu, V.; Sa´nchez-Carrera, R. S.; Paramonov, P.; Day, G. M.; Bre´das, J. L. Phys. Chem. C 2009, 113, 4679. (16) Holstein, T. Ann. Phys. 1959, 8, 325. (17) Holstein, T. Ann. Phys. 1959, 8, 343. (18) Coropceanu, V.; Kwon, O.; Wex, B.; Kaafarani, B. R.; Gruhn, N. E.; Durivage, J. C.; Neckers, D. C.; Bre´das, J. L. Chem.sEur. J. 2006, 12– 2073. (19) Malagoli, M.; Coropceanu, V.; da Silva, D. A.; Bre´das, J. L. J. Chem. Phys. 2004, 120, 7490. (20) Kato, T.; Yamabe, T. J. Chem. Phys. 2003, 119, 11318. (21) Gruhn, N. E.; da Silva, D. A.; Bill, T. G.; Malagoli, M.; Coropceanu, V.; Kahn, A.; Bre´das, J. L. J. Am. Chem. Soc. 2002, 124, 7918. (22) Coropceanu, V.; Malagoli, M.; da Silva Filho, D. A.; Gruhn, N. E.; Bill, T. G.; Bre´das, J. L. Phys. ReV. Lett. 2002, 89, 275503. (23) Olivier, Y.; Lemaur, V.; Bre´das, J. L.; Cornil, J. J. Phys. Chem. A 2006, 110, 6356. (24) Tsiper, E. V.; Soos, Z. G. Phys. ReV. B 2001, 64, 195124. (25) Verlaak, S.; Heremans, P. Phys. ReV. B 2007, 75, 115127. (26) Norton, J. E.; Bre´das, J. L. J. Am. Chem. Soc. 2008, 130, 12377. (27) McMahon, D. P.; Troisi, A. J. Phys. Chem. Lett. 2010, 1, 941. (28) Nagata, Y. ChemPhysChem 2010, 11, 474. (29) Castet, F.; Aurel, P.; Fritsch, A.; Ducasse, L.; Liotard, D.; Linares, M.; Cornil, J.; Beljonne, D. Phys. ReV. B 2008, 77, 115210. (30) Gosar, P.; Choi, S. Phys. ReV. 1966, 150, 529. (31) Vilfan, I. Phys. Status Solidi B 1973, 59, 351. (32) Friedman, L. Solid State Commun. 1981, 40, 41.

Polarization Energies in Anthracene (33) Brovchenko, I. V. Chem. Phys. Lett. 1997, 278, 355. (34) Kubo, R. Statistical mechanics, an adVanced course with problems and solutions; North Holland Publishing Co.: Amsterdam, New York, 1967. (35) Marcus, R. A. ReV. Mod. Phys. 1993, 65, 599. (36) Bre´das, J. L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chem. ReV. 2004, 104, 4971. (37) Hultell, M.; Stafstro¨m, S. Chem. Phys. Lett. 2006, 428, 446. (38) Ortmann, F.; Bechstedt, F.; Hannewald, K. New J. Phys. 2010, 12, 023011. (39) Hannewald, K.; Stojanovic´, V. M.; Schellekens, J. M. T.; Bobbert, P. A.; Kresse, G.; Hafner, J. Phys. ReV. B 2004, 69, 075211. (40) Ponder. J. W. TINKER: Software Tools for Molecular Design, 4.2 ed.: Washington University School of Medicine: Saint Louis, MO, U.S., 2004. (41) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933. (42) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15. (43) Blochl, P. E. Phys. ReV. B 1994, 50, 17953. (44) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (45) Brock, C. P.; Dunitz, J. D. Acta Crystallogr., Sect. B: Struct. Sci. 1990, 46, 795.

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20685 (46) Dewar, M. J. S.; Zoebisch, E. E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (47) Lii, J. H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 111, 8566. (48) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (49) Beeman, D. J. Comput. Phys. 1976, 20, 130. (50) R. Development Core Team. R: A Language and EnVironment for Statistical Computation; R Foundation for Statistical Computing; Vienna, Austria, 2008. (51) AMPAC 9; Semichem, Inc.: 12456 W 62nd Terrace - Suite D, Shawnee, KS 66216, U.S.; 1992-2008. (52) Jortner, J. J. Chem. Phys. 1976, 64, 4860. (53) Martinelli, N. G.; Olivier, Y.; Athanasopoulos, S.; Ruiz Delgado, M. C.; da Silva Filho, D. A.; Sa´nchez-Carrera, R. S.; Venuti, E.; Della Valle, R. G.; Bre´das, J. L.; Beljonne, D.; Cornil, J. ChemPhysChem 2009, 10, 2265. (54) Nan, G.; Yang, X.; Wang, L.; Shuai, Z.; Zhao, Y. J. Phys. Chem. B 2009, 79, 115203. (55) Ortmann, F.; Bechstedt, F.; Hannewald, K. Phys. ReV. B 2009, 79, 235206. (56) Stafstro¨m, S. Chem. Soc. ReV. 2010, 39, 2484.

JP105843T