Flexible Surface Hopping Approach to Model the Crossover from

May 17, 2013 - crossover from hopping to band-like transport regimes. This approach is ... polaron model can reproduce the band-to-hopping crossover...
0 downloads 0 Views 421KB Size
Letter pubs.acs.org/JPCL

Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-like Transport in Organic Crystals Linjun Wang* and David Beljonne Laboratory for Chemistry of Novel Materials, University of Mons, Place du Parc 20, B-7000 Mons, Belgium S Supporting Information *

ABSTRACT: Two distinct pictures are usually evoked when modeling charge transport in organic crystals, that is, band and hopping models, the signature of which is conveyed by a characteristic temperature dependence of mobility. Here, we present a novel flexible surface hopping approach compliant with general Hamiltonians that is able to grasp the crossover from hopping to band-like transport regimes. This approach is applied to solve a one-dimensional mixed quantum−classical model and to calculate the temperature dependence of charge mobility along with the degree of charge spatial localization. It is found that the roles of both local and nonlocal electron−phonon couplings strongly depend on the intrinsic charge localization strength.

SECTION: Energy Conversion and Storage; Energy and Charge Transport

A

restricted to small systems. It is widely accepted that the temperature dependence of mobility in conjugated organic crystals is dominated by low-frequency nuclear vibrations.6 It thus seems reasonable to adopt a scheme where a quantum mechanical treatment is reserved for the electrons while the nuclei are treated classically. However, at present, there is very limited cross-talking between the communities applying mixed quantum−classical dynamics (MQCD) and charge-transport simulations. MQCD is often carried out using mean-field (MF) approaches or by means of surface hopping (SH) algorithms,9 which differ in the description of the classical equations of motion. The nuclear dynamics at the MF level is governed by the gradient of the energy expectation value from the timedependent electron wave function, namely, a single average potential energy surface (PES). In contrast, the SH approach can incorporate nonadiabatic transitions between different adiabatic PESs and thus is superior to MF when nuclear motion on different adiabatic surfaces may lead to divergent trajectories. In recent years, there have been some attempts in the transport community to implement MF theory coupled to semiclassical model Hamiltonians to study charge dynamics in one-dimensional molecular stacks or two-dimensional planes.10,11 In parallel, MQCD groups have implemented SH algorithms to model charge-transfer processes in molecular dimers.12 There are, however, a number of caveats. On one hand, the MF approximation is obviously not appropriate in the limit of

long with the rapid development of organic electronics, charge transport in organic crystals has been the focus of many investigations in the past few decades.1 As pointed out by Pope and Swenberg 20 years ago, there exists a great difficulty in formulating a general analytical description of charge transport in organic crystals because many relevant factors, namely, intermolecular electronic couplings, phonon energies, thermal energies, and electron−phonon couplings (EPCs), are of similar magnitudes.2 It is only in limiting cases that transport models can be derived from perturbation theory, (i) the band model for large electronic couplings, where the EPCs are regarded as a perturbation, the charge is delocalized over the whole crystal, and the mobility decreases with temperature and (ii) the hopping model for large EPCs, where the electronic couplings are treated as a perturbation, the charge is completely localized on individual molecules, and the mobility increases with temperature.3 In intermediate regimes, the most widely used approach is the polaron model, which considers the electron and its associated deformation on the surrounding nuclei together as a quasi-particle (through a polaron transformation) moving in a polaron band or by polaron hopping.4,5 Previous ab initio studies have shown that the polaron model can reproduce the band-to-hopping crossover phenomenon in naphthalene crystals.6,7 However, the polaron model is also a perturbation theory, and the widespread used thermal average approximation for the polaron coupling is not appropriate at high temperatures.8 Moreover, the polaron model is often derived from specific model Hamiltonians and can hardly be extended to general cases; thus, the universality of the polaron model is questionable. Numerical simulation of charge transport is also a challenge. Full quantum dynamics are computationally expensive and © XXXX American Chemical Society

Received: April 24, 2013 Accepted: May 17, 2013

1888

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894

The Journal of Physical Chemistry Letters

Letter

localized charges because the mean wave function of many localized states tends to be delocalized; thus, SH should be used instead in such cases. In addition, the use of a SH approach offers the extra benefit that the charge localization length can be extracted from the simulations, which provides useful insight into the transport mechanism. On the other hand, injecting charge-transfer rates between molecular dimers in a chargetransport model would inevitably yield a hopping picture and thus cannot characterize properly band-like transport. To our knowledge, there have been no extensive studies so far aiming at implementing SH formalisms for charge transport in large supramolecular systems due to the difficulty in dealing with the high density of adiabatic PESs. Besides, electron couplings are local interactions, and thus, most crossings happen between localized adiabatic states that are far away in space with no contribution to the electron dynamics; hence, accurately dealing with these massive crossings is not trivial. In this Letter, we propose a flexible surface hopping (FSH) technique, which treats only a small part of the system in a SH manner and does it in a flexible way in time. In our approach, the problematic crossing problem is avoided because all adiabatic states are spatially close and the computational cost is largely reduced because we only diagonalize a very small Hamiltonian matrix to get all important PESs. Furthermore, we propose a flexible time step technique, which ensures the smoothness of all time-dependent adiabatic states, enabling us to describe SH accurately with relative large time steps. Our FSH algorithm based on a minimum subsystem implementing SH and a maximum time interval for SH dynamics can recover all transport regimes described by existing theories and can be easily extended to more complex Hamiltonians. We consider a one-dimensional model system containing N molecular sites with equal spacing between all nearest neighbors, L, and periodic boundary conditions. Each molecule k has one associated orbital |k⟩ (e.g., the highest occupied molecular orbital for hole transport and the lowest unoccupied molecular orbital for electron transport) and two (classical) intramolecular harmonic vibrational degrees of freedom xk,1 and xk,2 (the corresponding velocities, force constants, and effective masses are vk,1(2), K1(2), and m1(2), respectively). At equilibrium geometry, the electronic coupling between nearest neighbors is −τ, and the onsite energies of all molecules are the same and are therefore set to zero for simplicity. The onsite energy of molecule k is linearly modulated by xk,1 with the local EPC constant α, while the electronic coupling between molecule k and its neighbor k + 1 is linearly modulated by xk+1,2 − xk,2 with the nonlocal EPC constant β. Thus, our Hamiltonian reads H = He + Hn He =

Hn =

∑ k

k

(5)

Here, Vk,1(2) ′ = K1(2)xk,1(2) + dEa/dxk,1(2), with a representing the active adiabatic surface, γ is the friction coefficient, ξ1(2) is a Markovian Gaussian random force with standard deviation (2γm1(2)kBT/Δt)1/2, and Δt is the time interval. We can solve eq 4 through the standard fourth-order Runge−Kutta (RK4) algorithm.13 The stochastic differential equation in eq 5 can be also solved at the RK4 level following the work by Hershkovitz (see the Supporting Information (SI)).14 In addition, through the Hellmann−Feynman theorem,15 we obtain αpki pkj dijk ,1 = Ej − Ei (6) dijk ,2 =

β[pki (pk − 1, j − pk + 1, j ) + pkj (pk − 1, i − pk + 1, i )] Ej − Ei

(7)

dEi = αpki2 dxk ,1

(8)

dEi = 2βpki (pk − 1, i − pk + 1, i ) dxk ,2

(9)

d2Ei dxk2,1 d2Ei dxk2,2

= 2α ∑ djik ,1pki pkj j≠i

(10)

= 2β ∑ djik ,2[pki (pk − 1, j − pk + 1, j ) j≠i

+ pkj (pk − 1, i − pk + 1, i )]

(11)

For a given temperature T, the initial nuclear coordinates and velocities uniformly follow Boltzmann (Gaussian) distributions (the variances are kBT/K1(2) and kBT/m1(2) for xk,1(2) and vk,1(2), respectively).10 The charge is initially placed on the central molecule of the stack, and the corresponding nuclear coordinate xN/2,1 is displaced by −α/K1 in the charged geometry. In principle, the evolution of the electronic wave functions, nuclear coordinates, and velocities can be carried out according to eqs 4 and 5. However, we notice that organic crystals often involve relatively weak intermolecular couplings as well as large thermal disorder, so that the charge is generally delocalized over only a few molecules. Therefore, we have developed here a theoretical approach where only this subsystem is described using SH dynamics while the rest of the system is simply treated classically. In addition, we provide flexibility for the SH subsystem along the dynamics, through adding and/or removing neighboring molecules to/from the subsystem at any time step. For example, in order to check if a molecular orbital (|j⟩) should be added, we regard the active state of the SH subsystem (|Φa⟩ = ∑i pia|i⟩) and this test state together as a two-level system. Their onsite energies and coupling can be calculated through

(2)



(4)

k

m1(2)xk̈ ,1(2) = −Vk′ ,1(2) − γm1(2)xk̇ ,1(2) + ξ1(2)

k

1 (m1vk2,1 + K1xk2,1) + 2

j≠i

where = ⟨Φi|dΦj/dxk,1(2)⟩ are nonadiabatic couplings. The nuclear equation of motion is modeled by the Langevin equation

(1)

× (|k⟩⟨k + 1| + |k + 1⟩⟨k|)

∑ cj ∑ (xk̇ ,1·dijk ,1 + xk̇ ,2·dijk ,2)

dk,1(2) ij

∑ αxk ,1|k⟩⟨k| + ∑ [−τ + β(xk+ 1,2 − xk ,2)] k

1 ciEi − iℏ

ci̇ =

1 (m2vk2,2 + K 2xk2,2) 2 (3)

Diagonalizing He, we can get the adiabatic states in terms of the original diabatic orbitals, {Φi = ∑k pki|k⟩}, and the corresponding energies, {Ei}. We define the electron wave function as a linear expansion of these adiabatic states, Φ = ∑i ciΦi. According to the Schrödinger equation, we obtain 1889

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894

The Journal of Physical Chemistry Letters

Letter

Figure 1. Temperature dependence of the (A) charge mobility and (B) inverse participation ratio (IPR) for model stacks. τ varies as indicated with K1 = 14500 amu/ps2, m1 = 250 amu, α = 3500 cm−1/Å, and β = 0 cm−1/Å. The Marcus result and a power law T−1 dependence are shown as guides to the eyes.

⟨Φa |He|Φa ⟩ = α ∑ pi2a xi ,1

gij = Δt ·

i

+

∑ 2[−τ + β(xi+ 1,2 − xi ,2)]pi+ 1,a pia i

⟨j|He|j⟩ = αxj ,1

(15)

If gij < 0, it is reset to 0. We generate a uniform random number, ζ, and make a surface hop if ζ < gij. If Ei + ∑k(m1v2k,1 + m1v2k,1)/2 > Ej, a velocity adjustment is made to conserve the total energy

(12) (13)

A ⎡ vk′ ,1(2) = vk ,1(2) + dijk ,1(2) ·⎢ −1 + B ⎣

⟨Φa |He|j⟩ = αpja xj ,1 + [−τ + β(xj + 1,2 − xj ,2)]pj + 1,a + [−τ + β(xj ,2 − xj − 1,2)]pj − 1,a

2Re[cic*j ∑k (vk ,1dijk ,1 + vk ,2dijk ,2)] ci*ci

1 + 2(Ei − Ej)

B ⎤ ⎥ A2 ⎦ (16)

(14)

k,2 k,1 k,1 where A = ∑k(m1vk,1dk,1 ij + m2vk,2dij ) and B = ∑k(m1dij ·dij + k,2 2 m2dk,2 ·d ); the hop is rejected when E < E − A /2B. Suppose ij ij i j that the hop takes place, the trajectory then evolves on the PES Ej. On the basis of a recent work by Bittner et al.,17 decoherence is achieved by collapsing the electronic wave function to the new PES after every successful and unsuccessful hop. For each realization m, the time evolution of the active state Φ(m) a (t) is restored and used to calculate the mean-squared displacement (MSD) after a total number of M realizations

When the ratio between the absolute coupling, abs{⟨Φa|He|j⟩}, and the absolute energy difference, abs{⟨Φa|He|Φa⟩ − ⟨j|He|j⟩}, is larger than a critical fraction Rc, we add |j⟩ into the new SH subsystem. The same procedure is followed when removing a molecular orbital from the SH subsystem. Meanwhile, the active state is changed to the adiabatic state of the new subsystem that mostly overlaps with the old active state. Obviously, this approach converges to the traditional SH when Rc tends to 0, while alleviating the problem of trivial crossings at finite Rc. In addition, due to the fact that adiabatic states are strongly mixed only in the nonadiabatic regimes, we can use a flexible time interval for the dynamics to increase computational efficiency. Namely, a proper Δt is adjusted at each time step to ensure that the minimum overlap between the adiabatic state at time t and the corresponding state at the next time step t + Δt, min {|⟨Φi(t)|Φi(t + Δt)⟩|}, cannot be smaller than a

MSD(t ) =

1 M

M

∑ [⟨Φ(am)(t )|r 2|Φ(am)(t )⟩ m=1

− ⟨Φ(am)(t )|r |Φ(am)(t )⟩2 ]

(17)

where the elements can be computed with ⟨k|r2|l⟩ = δklk2L2 and ⟨k|r|l⟩ = δklkL. When the MSD shows a linear evolution in time, one enters equilibrium diffusion, and the diffusion coefficient reads

i

critical value Oc. Thus, smaller time steps are automatically generated when the dynamics approaches nonadiabatic regimes to prevent the divergence of eqs 6 and 7. Following Tully’s standard “fewest switches” algorithm,16 the switching probabilities from the active surface i to another surface j are

D=

1 dMSD lim 2 t →∞ dt

(18)

Finally, the charge mobility is obtained by means of the Einstein relation 1890

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894

The Journal of Physical Chemistry Letters

Letter

Figure 2. (A,C) Time evolution of the mean IPR from 10 000 realizations and (B,D) distribution of the IPR at the last snapshot of the calculations. (A) and (B) are calculated for τ = 50 at 300 K, while (C) and (D) are obtained for τ = 800 at 300 K.

μ=

e D kBT

different values for τ, that is, 50, 100, 200, 400, and 800 cm−1, to cover all regimes from weak to strong intermolecular couplings. We carefully checked the convergence of the charge mobility and polaron radius with respect to Rc for different τ (see the SI) and set the corresponding Rc values to 1/15, 1/40, 1/70, 1/130, and 1/150, respectively. From Figure 1A, it is clear that our approach successfully reproduces the expected smooth transition between hopping and band-like regimes as a function of increasing electronic coupling. In detail, mobility is thermally activated for small τ (50 and 100 cm−1), while it monotonously decreases with temperature for large τ (400 and 800 cm−1). In the intermediate regime at around τ = 200 cm−1, the mobility is quite insensitive to temperature. To the best of our knowledge, this is the first MQCD study that can cover both hopping and band-like transport regimes. Note that such a finding was also observed by Fratini and Ciuchi when using a nonperturbative dynamical MF theory to study charge transport.19 On the basis of a similar set of parameters as those used here, they also observed a crossover from hopping to band-like transport in excellent agreement with our results.20 The corresponding temperature dependence of charge localization strength is shown in Figure 1B. It is found that the charge size is strongly correlated with the carrier mobility, namely, the charge gets larger (smaller) when the mobility increases (decreases). Thermally activated hopping happens when the charge size is close to unity, agreeing with the general view on hopping transport. Note that the mean IPR generally converges very rapidly with time at around 300 fs (see Figures

(19)

The charge localization length can be evaluated by the inverse participation ratio (IPR) IPR =

1 M

M

∑ m=1

1 ∑k ⟨k|Φ(am)(t )⟩4

(20)

In all of the numerical calculations below, we adopt molecular stacks containing over 50 molecular sites to avoid the charge interacting with its image. The intermolecular distance is set to be L = 5 Å, consistent with typical values in molecular crystals. We set Oc = 0.9999 to ensure reliable SH dynamics. We have tested different friction coefficients and chose γ = 1 ps−1 to avoid underdamping or overdamping motions. Note that the choice of γ does not influence any conclusion from this study. The maximum time interval is set to be 0.1 fs for the sake of computational stability, and the total simulation time is 2 ps. More than 10000 realizations are performed to get a smooth time evolution of the MSD. We preliminarily investigate the Holstein model, namely, the case with only local EPCs. In light of our recent study on charge transport in molecular crystals,11 we adopt the following set of parameters: K1 = 14500 amu/ps2, m1 = 250 amu, and α = 3500 cm−1/Å. The corresponding vibrational frequency (ω = (K1/m1)1/2) and reorganization energy (λ = α2/K1) are calculated to be about 40 and 1000 cm−1, which are typical for semiconducting molecular materials.10,18 We choose 1891

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894

The Journal of Physical Chemistry Letters

Letter

Figure 3. Temperature dependence of the (A) charge mobility and (B) inverse participation ratio for model stacks. τ and β vary as indicated with K1 = 14500 amu/ps2, m1 = 250 amu, K2 = 14500 amu/ps2, m2 = 250 amu, and α = 3500 cm−1/Å.

Figure 4. Time evolution of (A,C) the mean charge population on each molecular site and (B,D) the mean vibrational coordinate (in Å) of each molecule. The parameters are K1 = 14500 amu/ps2, m1 = 250 amu, α = 3500 cm−1/Å, β = 0 cm−1/Å, and T = 300 K, with τ = 50 cm−1 for (A) and (B) and τ = 800 cm−1 for (C) and (D). The snapshot time interval is 2 fs.

over 10 000 realizations, features a wide distribution especially at large τ, showing the coexistence of localized and delocalized electron states.21 We stress, however, that charge transport

2A and 2C); thus, the reported IPRs are obtained through an average over the time range between 500 fs and 2 ps. As shown in Figures 2B and 2D, the charge carrier size at 2 ps, averaged 1892

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894

The Journal of Physical Chemistry Letters

Letter

FSH approach counts on an equal footing of both intermolecular couplings and EPCs. Fourth, the polaron states obtained through polaron transformation are invariable, while our active charge states are generated automatically and are flexible during the dynamics. Fifth, our time evolution of the charge population gives clear transport channels in real space, which cannot be easily assessed from the polaron model. Finally, our approach can be easily applied to more complex situations, for example, including static disorders, acoustic phonons, asymmetric vibrations, and special topological conformations, which cannot be easily described within the polaron model. Note also that MF calculations typically overestimate the delocalization strength of charge carriers and fail to properly describe the thermal activated charge transport in the hopping regime.24 There are, as well, limitations in the model. (1) The most stringent one relies on the fact that, by design, the FSH approach works better closer to the hopping regime. With the increase of the amount of charge delocalization, the FSH approach is expected to become less efficient and accurate and should fail to describe phenomena associated with pure band transport such as the quantum Hall effect. (2) The present band hopping crossover is different from the crossover from coherent motion at low temperature to classical activated motion at high temperature,4−7 which has been widely studied in polaron theories and is mainly contributed to the quantum effect of nuclear vibrations.25 Because the nuclei are treated classically here, the temperature-dependent mobility at (very) low temperatures cannot be captured by the present model. However, this can be settled by including the vibrational states together with the electronic states (in a multiple particle basis set)26 as the quantum part of the FSH simulation. (3) In this work, we have used a simple dephasing strategy, where the wave function is recollapsed to an appropriate PES after every hop. Though this is a reasonable approximation close to the hopping regime, it would be of interest to study in more details how the dephasing time affects the temperature-dependent charge carrier mobility in organic materials. We further note that various dephasing methods have been suggested in the literature,27 and our FSH can potentially cope with these approaches. To sum up, we have proposed a flexible surface hopping technique, which considers only a fraction of the system for surface hopping dynamics and employs a flexible time interval to bypass the issue of multiple (trivial) crossings at high density of adiabatic states while speeding up numerical calculations. We have shown that it is a convenient and reliable approach to study charge transport in complex organic materials. The model is able to describe the crossover from a thermally activated hopping transport at low electronic coupling to a power law band-like temperature dependence at high couplings. Our approach provides a clear rationale of the relationship between charge mobility and carrier size; in particular, it shows that the hopping-to-band crossover already occurs for charge carriers that spread over only a few sites/molecules. We have also shown how thermal dynamic disorder can either decrease or enhance charge transport depending on the intrinsic localization strength of charge carriers.

follows a band-like temperature dependence even when the charge extends over only two or three molecules. This means that a band-like transport does not necessarily translate into a delocalization of the excess charge over a large amount of molecules as usually implied in band models; it rather indicates that the charge becomes more localized with increasing temperature. Next, we study the role of Peierls terms, that is, the nonlocal EPCs, on charge transport (see Figure 3). For the sake of convenience, we choose the same force constant and effective mass for the corresponding vibrations, namely, K2 = 14500 amu/ps2 and m2 = 250 amu, and consider two different nonlocal EPCs, β = 500 and 1000 cm−1/Å. Two sets of calculations were carried out; one adds nonlocal EPCs to the hopping transport at play with τ = 50 cm−1, while the other one does it in the band-like transport scenario corresponding to τ = 800 cm−1, as discussed above. It is found that the nonlocal EPCs have exactly opposite effects in these two cases; the mobility is increased when accounting for nonlocal EPCs in the hopping regime but reduced in the band-like case. We can understand this observation and the Holstein results consistently based on the intrinsic charge localization strength in the absence of thermal disorder. At zero temperature, the competition between the local EPCs that tend to localize the charge and the intermolecular couplings that favor charge delocalization result in charge states with different localization lengths. For localized charges, increasing the temperature (the local EPC effect is thereby enhanced) or adding nonlocal EPCs induces fluctuations to the localized charge state and increases the possibility for the charge to jump over the barrier; thus, the charge mobility increases. In contrast, if the charge is delocalized, any thermal effect from local or nonlocal EPCs tends to destroy the translational symmetry of the electron Hamiltonian, resulting in smaller delocalization reflected in a lower mobility. Different views have been reported in the literature concerning the role of nonlocal EPCs on charge transport. It was pointed out by Munn and Silbey that incorporating nonlocal EPCs yields higher hopping contributions and increases the charge mobility.22 In contrast, Troisi recently pointed out that nonlocal EPCs localize “dynamically” the charge carrier, resulting in a decreased mobility.10 The results reported here thus unify these two pictures, which appear as the two limiting cases obtained for small and large electronic couplings above. Note that similar findings have been obtained recently using a nonperturbative quantum master equation approach based on the Liouville space hierarchical equations of motion method to study the role of nonlocal EPCs on charge transport in molecular crystals.23 In addition, we compare the time evolution of charge population with that of nuclear coordinates in both hopping and band-like regimes and find that the nuclear vibrations are always closely correlated with the electronic motion (see Figure 4). From this point of view, our FSH algorithm is qualitatively similar to the polaron model. However, we emphasize that our approach is superior over the (small) polaron theory in many aspects. First, the polaron model is based on perturbation theory, while our FSH model is nonperturbative and can thus be applied to a broader parameter range. Second, the polaron model breaks down at high temperatures, where the fluctuations of the polaron couplings are larger than the average couplings,8 which is not an issue in our approach. Third, the temperature dependence of mobility is solely determined by the EPCs in the polaron model,4 while the



ASSOCIATED CONTENT

S Supporting Information *

RK4 approach to solve the stochastic differential equation and convergence of the charge carrier mobility and size with respect 1893

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894

The Journal of Physical Chemistry Letters

Letter

(18) Wang, L. J.; Nan, G. J.; Yang, X. D.; Peng, Q.; Li, Q. K.; Shuai, Z. Computational Methods for Design of Organic Materials with High Charge Mobility. Chem. Soc. Rev. 2010, 39, 423−434. (19) Fratini, S.; Ciuchi, S. Dynamical Mean-Field Theory of Transport of Small Polarons. Phys. Rev. Lett. 2003, 91, 256403. (20) If we apply a phonon frequency of 40 cm−1 to Figure 3 of ref 19, we obtain an electron−phonon coupling constant of 3081 cm−1/Å and a crossover from hopping to band-like transport for τ ≈ 200 cm−1. (21) Fratini, S.; Ciuchi, S. Bandlike Motion and Mobility Saturation in Organic Molecular Semiconductors. Phys. Rev. Lett. 2009, 103, 266601. (22) Munn, R. W.; Silbey, R. Theory of Electronic Transport in Molecular Crystals. III. Diffusion Coefficient Incorporating Nonlocal Linear Electron−Phonon Coupling. J. Chem. Phys. 1985, 83, 1854. (23) Wang, D.; Chen, L. P.; Zheng, R. H.; Wang, L. J.; Shi, Q. Communications: A Nonperturbative Quantum Master Equation Approach to Charge Transport in Organic Molecular Crystals. J. Chem. Phys. 2010, 132, 081101. (24) Wang, L. J.; Beljonne, D. Charge Transport in Organic Semiconductors: Assessment of the Mean Field Theory. J. Chem. Phys., submitted. (25) Nan, G. J.; Yang, X. D.; Wang, L. J.; Shuai, Z. Nuclear Tunneling Effects of Charge Transport in Rubrene, Tetracene, and Pentacene. Phys. Rev. B 2009, 79, 115203. (26) Philpott, M. R. Theory of the Coupling of Electronic and Vibrational Excitations in Molecular Crystals and Helical Polymers. J. Chem. Phys. 1971, 55, 2039−2054. (27) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-Field Dynamics with Stochastic Decoherence (MF-SD): A New Algorithm for Nonadiabatic Mixed Quantum/Classical Molecular-Dynamics Simulations with Nuclear-Induced Decoherence. J. Chem. Phys. 2005, 123, 234106.

to Rc. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank Prof. Sergei Tretiak for helpful suggestions. This work was supported by the Interuniversity Attraction Pole program of the Belgian Federal Science Policy Office (PAI 6/ 27) and FNRS-FRFC. D.B. is a FNRS Research Director.



REFERENCES

(1) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.; Silbey, R.; Brédas, J.-L. Charge Transport in Organic Semiconductors. Chem. Rev. 2007, 107, 926−952. (2) Pope, M.; Swenberg, C. E. Electronic Processes in Organic Crystals; Oxford University Press: New York, 1982. (3) Shuai, Z.; Wang, L. J.; Li, Q. K. Evaluation of Charge Mobility in Organic Materials: From Localized to Delocalized Descriptions at a First-Principles Level. Adv. Mater. 2011, 23, 1145−1153. (4) Hannewald, K.; Stojanović, V. M.; Schellekens, J. M. T.; Bobbert, P. A.; Kresse, G.; Hafner, J. Theory of Polaron Bandwidth Narrowing in Organic Molecular Crystals. Phys. Rev. B 2004, 69, 075211. (5) Cheng, Y.-C.; Silbey, R. J. A Unified Theory for Charge-Carrier Transport in Organic Crystals. J. Chem. Phys. 2008, 128, 114713. (6) Wang, L. J.; Peng, Q.; Li, Q. K.; Shuai, Z. Roles of Inter- and Intramolecular Vibrations and Band-Hopping Crossover in the Charge Transport in Naphthalene Crystal. J. Chem. Phys. 2007, 127, 044506. (7) Wang, L. J.; Li, Q. K.; Shuai, Z. Effects of Pressure and Temperature on the Carrier Transport in Organic Crystal: A FirstPrinciples Study. J. Chem. Phys. 2008, 128, 194706. (8) Troisi, A. Quantum Dynamic Localization in the Holstein Hamiltonian at Finite Temperatures. Phys. Rev. B 2010, 82, 245202. (9) Tully, J. C. Mixed Quantum−Classical Dynamics: Mean-Field and Surface-Hopping. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; World Scientific: Singapore, 1998. (10) Troisi, A.; Orlandi, G. Charge-Transport Regime of Crystalline Organic Semiconductors: Diffusion Limited by Thermal Off-Diagonal Electronic Disorder. Phys. Rev. Lett. 2006, 96, 086601. (11) Wang, L. J.; Beljonne, D.; Chen, L. P.; Shi, Q. Mixed Quantum− Classical Simulations of Charge Transport in Organic Materials: Numerical Benchmark of the Su−Schrieffer−Heeger Model. J. Chem. Phys. 2011, 134, 244116. (12) Landry, B. R.; Subotnik, J. E. Communication: Standard Surface Hopping Predicts Incorrect Scaling for Marcus’ Golden-Rule Rate: The Decoherence Problem Cannot be Ignored. J. Chem. Phys. 2011, 135, 191101. (13) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes; Cambridge University Press: Cambridge, U.K., 1992. (14) Hershkovitz, E. A Fourth-Order Numerical Integrator for Stochastic Langevin Equations. J. Chem. Phys. 1998, 108, 9253. (15) Chernyak, V.; Mukamel, S. Density-Matrix Representation of Nonadiabatic Couplings in Time-Dependent Density Functional (TDDFT) Theories. J. Chem. Phys. 2000, 112, 3572. (16) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061. (17) Bittner, E. R.; Rossky, P. J. Quantum Decoherence in Mixed Quantum−Classical Systems: Nonadiabatic Processes. J. Chem. Phys. 1995, 103, 8130. 1894

dx.doi.org/10.1021/jz400871j | J. Phys. Chem. Lett. 2013, 4, 1888−1894