Evaluation of Electronic Coupling in Solids from Ab Initio Periodic

Article pubs.acs.org/JPCC

Evaluation of Electronic Coupling in Solids from Ab Initio Periodic Boundary Condition Calculations: The Case of Pentacene Crystal and Bilayer Graphene Alessandro Biancardi and Marco Caricato* Department of Chemistry, University of Kansas, 1251 Wescoe Hall Drive, Lawrence, Kansas 66045, United States S Supporting Information *

ABSTRACT: Understanding the modulation of the electronic coupling in molecular crystals and two-dimensional materials is crucial from a fundamental point of view as well as for the development of organic electronics. In this work, we present a first-principles quantum-mechanical method for the calculation of the electronic coupling (or transfer integrals) between fragments or layers, using density functional theory with periodic boundary conditions (DFT-PBC) within the Γ-point approximation. This method is applied to two periodic systems: crystalline pentacene and a bilayer graphene film. For the former system, we find that the inclusion of the solid environment affects the interfragment electronic couplings, with changes of the order of 10%. However, we confirm the qualitative trends obtained with the “isolated molecular dimer” model. For the graphene film, we show how the interlayer coupling changes with the relative position of two π-stacked layers. Interestingly, we find that particle−particle coupling is large even for configurations that are not perfectly stacked.

1. INTRODUCTION An increased understanding of both electronic coupling and intermolecular interactions in organic molecular crystals and two-dimensional materials is crucial for materials1,2 and biomolecular science.3 For instance, organic-based materials have received much attention because of their unusual properties, e.g., flexibility, lightness, and tunability, which make them competitive with traditional inorganic materials for electronic devices. In particular, the reliable evaluation of charge transport rates is a key component in this research field.4−9 The electronic coupling, J, also known as transfer integral, is a fundamental term in Marcus formula for the hopping charge transfer rate:10,11 k=

|J |2 ℏ

⎛ (ΔE − λ)2 ⎞ π exp⎜ − ⎟ 4λkBT ⎠ λkBT ⎝

materials. Indeed, computed electronic couplings can be semiquantitatively correlated with measured carrier mobilities.15 In this respect, Brédas and co-workers4,18,19 have pioneered the development of techniques for the ab initio quantum-mechanical calculation of electronic coupling between isolated pairs of fragments (which we will call “isolated molecular dimer” model, or IMD model). However, as mentioned above, it is desirable to consider the effect of the extended solid in the evaluation of the coupling. Therefore, in this work, we go beyond the IMD model4,18−21 to compute the interfragment/layer electronic coupling in periodic solids. We derive the coupling expression for periodic systems within the Γ-point approximation, using density functional theory with periodic boundary conditions (DFTPBC) to determine the electronic structure of the system. Our approach directly includes polarization effects from the extended solid, and it allows the treatment of extended systems that cannot be divided into molecular fragments (e.g., a graphene layer). Additionally, our code can handle a general definition of (finite or infinite size) fragments, so that the coupling can be computed between multiple fragment pairs at once. We use this method to study the interfragment/layer electronic coupling of a single layer of crystalline pentacene along its ab-plane, and π-stacked bilayer graphene. Pentacene has attracted considerable interest due to the high charge

(1)

where ΔE is the difference between the site energies, and λ is the reorganization energy. Although the coupling is not the only factor governing the transfer process, it is often sufficient for a qualitative description of such a process. In spite of extensive efforts, however, a general microscopic understanding of the hole/particle transfer modulation is still unsatisfactory.12−14 A factor that is difficult to determine experimentally is the modulation of the crystal packing on the transfer.15,16 An important example is bilayer graphene, where the strength of the interlayer electronic coupling is still under debate.17 Computational modeling can help in achieving a better microscopic understanding of the interplay between structure and transfer modulation as well as in guiding the design of new © XXXX American Chemical Society

Received: June 14, 2016 Revised: July 21, 2016

A

DOI: 10.1021/acs.jpcc.6b06011 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C mobility exhibited by its crystals,22 and it is now a prototypical organic conducting system in both experimental and theoretical studies. We use this system to investigate the effect of crystalline packing on the coupling, which can affect the preferential pathway of charge transport.23,24 In this case, we find that intermolecular interactions in the crystal are weak enough that the IMD approach provides a qualitatively correct picture,25 although it underestimates |J|2 by about 10%. Bilayer graphene is another prototypical system for twodimensional few-layer materials that show unusual properties26 with promising applications in organic electronics.27−32 A graphene layer consists of sp2 carbon atoms32 whose electronic properties originate from the structure of the energy levels near the Fermi energy.33−35 These properties can be effectively modulated both by size quantization effects (e.g., by cutting, bending, and folding) and by layer stacking. In bilayer graphene, different stacking orders may occur: the most stable is the AB (Bernal) stacking, although the less common AAstacking has also been reported36 (including in graphite37). In this work, we evaluate the change in interlayer coupling by translating one of the layers along the lattice vectors. As expected, the hole coupling is largest when the two layers are perfectly stacked (AA), while the coupling decreases rapidly when the geometry is moved away from this optimal orientation. However, the particle coupling is surprisingly large not only in the AA geometry but also along the a1⃗ − a2⃗ translation direction from the AA geometry (where a1⃗ and a2⃗ are the lattice vectors). The paper is organized as follows. Section 2 describes the theory for the calculation of the coupling. The results for pentacene and graphene are described in section 3. Finally, section 4 summarizes our findings and reports some concluding remarks.

KS matrix of the supersystem projected on the basis of the fragments:19,21

FijAB = ⟨ϕi A |F |̂ ϕjB⟩

(2)

where F̂ is the Fock operator, is the ith orbital of fragment A, and ϕBj is the jth orbital of fragment B. Within the single band approximation, i.e., assuming only one relevant orbital per fragment, the hole/particle transfer takes place within the highest-occupied/lowest-unoccupied orbitals, respectively. The expression in eq 2 can be evaluated by using a projective technique, i.e., by using the resolution of the identity, 1̂ = ∑l | ϕl⟩ ⟨ϕl|, twice: ϕAi

FijAB =

∑ ⟨ϕi A |ϕl⟩⟨ϕl|F |̂ ϕm⟩⟨ϕm|ϕjB⟩ = ∑ γilAεlγljB lm

l

(3)

where ϕl and ϕm are the supersystem canonical orbitals. The Fock matrix is diagonal in that basis, thus εl are the supersystem orbital energies, and γBlj is the projector onto the jth orbital of fragment B. For isolated molecules, it is convenient to expand the orbitals over atomic centered Gaussian functions. The same expansion will be also used in the PBC development, although this treatment can be equally applied to plane wave basis sets as well. In the atomic orbital basis {χα}, the projector takes the form γmjB = ⟨ϕm|ϕjB⟩ =

∑ CαmSαβCβBj (4)

αβ

where Sαβ = ⟨χα|χβ⟩ is the overlap element between the α and β atomic orbitals, and Cαm/CBβj are the expansion coefficients of the orbitals of the supersystem/fragment, respectively. The fragment orbitals used for the projection of the supersystem orbitals (eq 2) are not orthonormal in general, and their overlap can be written as

2. THEORY The electronic coupling is a fundamental factor governing charge transfer processes that occur with a hopping mechanism, and the transfer rate can be evaluated using Marcus formula, see eq 1.10,11 An efficient and rigorous way to obtain the coupling from electronic structure calculations is through the Fock/ Kohn−Sham (KS) matrix reconstruction method in the basis of the fragment orbitals.4,19,21,38,39 In this approach, the supersystem where the transfer occurs is divided into a donor and an acceptor fragment, and the electronic coupling is obtained as the off-diagonal element of the Fock/KS matrix projected on the space of the orthogonal fragment orbitals. Thus, a meanfield approximation is introduced to describe the electron− electron interaction in the supersystem and in the fragments. More importantly, the supersystem and, thus, the fragments are isolated and do not feel any polarization effect from the solid state environment. In this work, we extend this approach in two simultaneous directions: (i) to include the effect of the environment through periodic boundary conditions, and (ii) by considering an arbitrary number of fragments so that the two-body coupling can be computed simultaneously for all pairs of fragments. In order to describe these developments, we review the IMD approach in section 2.1, while the extension to PBC is presented in section 2.2. 2.1. Electronic Couplings between Molecular Orbitals. Once two fragments A and B are defined out of the supersystem, the coupling between two orbitals of the fragments is related to the off-diagonal elements of the Fock/

SijAB = ⟨ϕi A |ϕjB⟩ =

∑ ⟨ϕi A |ϕl⟩⟨ϕl|ϕjB⟩ = ∑ γilAγljB l

l

(5)

Therefore, in order to maintain the local character of the monomer orbitals,19 we apply a Löwdin orthogonalization:40 F̃

AB

= (SAB)−1/2 F(SAB)−1/2

(6)

Finally, the electronic couplings are given by the off-diagonal elements of the transformed Fock matrix: AB

JijAB ≡ Fij̃ , i ∈ A, j ∈ B

(7)

2.2. Electronic Couplings between Γ-Point Crystal Orbitals. In the case of a periodic system, we can obtain the couplings within the unit cell by rewriting eq 2 as a sum over each lattice vector Rn: FijAB =

1 M

M

∑ ∫ ψi A*(r + R n)Fψ̂ jB(r + R n) dr n

(8)

where M is the number of cells considered within the summation over Rn. The fragment orbitals can be expanded in k space as ψ jA(r + R n) =

1 N

N

∑ ϕjA,k (r) e−ik·R k

n

(9)

obtaining B


Article

The Journal of Physical Chemistry C FijAB

1 = MN

M

∑∑e n

are periodic systems where packing41,42 and periodicity43 modulate the hole/particle couplings. All DFT and DFT-PBC calculations were performed with a development version of the GAUSSIAN suite of programs44−46 using the HSE0647−52 functional. HSE06 is a range-separated hybrid functional with exact exchange in the short range, and it was developed specifically for PBC calculations as a good compromise between accuracy and computational effort (as the long-range decay of the exchange considerably increases the computational cost of the PBC procedure). Since it includes exact exchange, this functional does not require any ad hoc parameter for the band gap. For the basis set, we use atomic centered Gaussian functions. In particular, we use the 6-31G basis for pentacene so that we can compare directly with the recent work of Li et al.,25 while for graphene we use 6-31G(d), which contains an extra set of polarization d functions. For pentacene, we also computed the couplings with the 6-31G(d,p) basis set, which includes a set of polarization functions on each atom. These results show that the effect of the polarization functions is rather small, as it leads to a systematic decrease of |J|2 by at most 4 × 10−4 eV2 compared to the 6-31G results. Therefore, we only report the data with the larger basis set in the Supporting Information, see Tables S1 and S2. This choice of basis sets is convenient since we based the implementation of our method on the crystal orbital energies and coefficients obtained with the GAUSSIAN program. However, as mentioned above, the implementation can be extended to plane wave DFT in the future. 3.1. Pentacene Crystal. Pentacene has become a standard system to study charge transport in organic molecular crystals due to its high charge-carrier mobility (>1 cm2 V−118) and simple molecular structure. We focus on transport within the ab plane since |J|2 is 1 order of magnitude smaller in the direction perpendicular to the ab plane, as previously noted by Troisi13 and confirmed by our own calculations. Our goal is 2-fold: (i) quantify the effect of the extended solid on the coupling between two pentacene units, and (ii) investigate whether the coupling should involve only pairs of monomer units as donor and acceptor moieties, or if these moieties should include extended layers of pentacene units. In either case, our method will provide information on the effects of the crystalline structure on the coupling. The pentacene layer is built using the experimental crystal structures obtained by Haas et al.53 at different temperatures. The unit cell contains two pentacene molecules, whose periodic repetition forms the well-known herringbone pattern.54 We replicate the unit cell building a 3 × 3 × 1 supercell, containing a single layer of crystalline pentacene along its ab plane, as shown in Figure 1. This size of supercell ensures convergence in the value of the electronic coupling within the Γ-point approximation, and for the energy calculation we use a uniform

N i(k − k ′)·R n

∫

̂ B (r) ϕi ,Ak*(r)F ϕ j,k′

dr

k, k ′

(10)

Within the Γ-point approximation, eq 10 reduces to FijAB ,Γ =

∫ ϕi ,AΓ*(r)Fϕ̂ jB,Γ(r) dr = ⟨ϕi ,AΓ|F|̂ ϕjB,Γ⟩

ϕAi,Γ

(11)

ϕBj,Γ

where and are the Γ-point canonical fragment orbitals, which are real functions. This expression is similar to the IMD model in eq 2. Thus, we can apply the same procedure as in section 2.1 to obtain the couplings. The final expression for the couplings is AB

JijAB ≡ Fij̃ , Γ , i ∈ A, j ∈ B ,Γ

(12)

where AB

FΓ̃ = (S ΓAB)−1/2 FΓ (S ΓAB)−1/2

(13)

SAB Γ

The overlap matrix and the corresponding projectors and orbital coefficients have expressions similar to those for molecules, see eqs 4 and 5, except that now the orbitals are the Γ-point crystal orbitals. So far, we have worked with two fragments, so that the orbital coefficient matrix of the fragments can be written as the direct sum of the orbital coefficient matrices of the individual fragments: C frag Γ

⎛C A 0 ⎞ Γ ⎟ = ⎜⎜ B⎟ ⎝ 0 CΓ ⎠

(14)

and the total projector matrix (compare with eq 4) is written as γΓ = C†Γ SC frag Γ

(15)

However, there is no reason why the number of fragments must be limited to two. We can define an arbitrary number N of fragments so that the matrix of the fragment coefficients Cfrag Γ can be extended to

C frag Γ

⎛C A 0 ⋯ 0 ⎞ ⎜ Γ ⎟ ⎜ 0 CB ⎟ Γ ⎟ =⎜ ⎜⋮ ⋱ ⋮ ⎟ ⎜⎜ ⎟⎟ ⋯ C ΓN ⎠ ⎝ 0

(16)

In this way, the coupling between every pair of fragments can be computed directly using eqs 6 and 13 in the case of the IMD and PBC approaches, respectively. It is important to note that the Γ-point approximation is only applied in the evaluation of the coupling, see eq 11, but the calculation of the crystal orbitals and their energy is performed in the full k space both for the supersystem and for the fragments. Hence, the polarization effect of the extended solid on the coupling is introduced through the polarization of the Γpoint orbitals. Our implementation also allows a mixed procedure, where the supersystem calculation is performed with PBCs, while the fragments can be treated as isolated molecules (as in the pentacene case) or as periodic systems (as in the graphene case).

3. RESULTS AND DISCUSSION In this section, we consider the interfragment/layer couplings within crystalline pentacene and a bilayer graphene film. Both

Figure 1. Left: Top view of the ab-layer pentacene supercell and translation vectors, a⃗ and b⃗. Right: Transfer pathways t1, t2, t3, and t4 between pentacene units, as defined in ref 25. C


Article

The Journal of Physical Chemistry C grid of 6 × 6 irreducible k-points. According to the energies of the highest occupied crystal orbital (HOCO) and lowest unoccupied crystal orbital (LUCO), the gaps are 1.25 eV, 1.33 eV, and 1.41 eV for the crystal structures obtained at 120 K, 293 K, and 414 K, respectively. These values are qualitatively in agreement with the previous result of 1.1 eV reported by Endres et al.,55 using plane-wave DFT (Perdew−Burke− Ernzerhof functional56) and the experimental crystal structures obtained by Campbell and co-workers,57 and smaller than the experimental one of 1.82 eV reported for a thin pentacene film at room temperature.58 This underestimation has been related to the approximate nature of current functionals.55 We start our analysis from the t1−t4 couplings for isolated pentacene dimers, defined in Figure 1, for a direct comparison with the results reported by Li et al.25 in their investigation of the temperature effect on the charge transport in pentacene crystals. Note that the temperature effect is only represented by the different experimental crystal structures, and there is no direct temperature effect in the coupling calculations. Table 1 Table 1. Comparison of |J|2 (eV2 × 1000) along the Transfer Pathways t1, t2, t3, and t4, See Figure 1, Using Isolated Pentacene Pairs hole t1

t2

B3LYP25 HSE06

2.1 1.9

1.8 1.7

B3LYP25 HSE06

1.6 1.5

1.3 1.2

B3LYP25 HSE06

1.4 1.3

1.0 0.9

particle t3

t4

T = 120 K 4.5 10.3 4.4 10.1 T = 293 K 3.0 7.8 2.9 7.7 T = 414 K 2.1 6.5 2.0 6.3

t1

t2

t3

t4

3.0 2.9

2.8 2.7

9.4 9.4

8.9 8.9

2.6 2.5

2.4 2.3

8.0 7.9

7.5 7.5

2.3 2.3

2.1 2.0

6.5 6.5

6.3 6.2

reports |J|2 computed with the IMD model and the B3LYP functional in ref 25, and with the HSE06 functional. This data validates our choice of functional as the values of the coupling are nearly identical to those reported in the literature: the strongest channels for both hole and particle transport are t3 and t4, probably due to a more favorable overlap of the π electron cloud. The t1 and t2 couplings are larger for the electron transport than for the hole transport, likely due to the larger spatial extent of the unoccupied orbitals, but their overall values remain smaller than the t3 and t4 couplings. The magnitude of the coupling decreases at higher temperature, and the overall coupling is similar for hole and particle transport.25,59 In order to assess the polarization effect of the crystal packing, we computed the electronic coupling between all possible pairs of pentacene units within the supercell in Figure 1 using the method described in section 2. We evaluate the coupling by using crystal orbital energies and coefficients obtained with the GAUSSIAN program, by performing a PBC calculation for the supercell and two individual calculations of the pentacene units shown in Figure 1. The latter is an important simplification that can be used in this case since more computationally intensive PBC calculations on the fragments would not affect the value of the coupling. The pairwise coupling for the 18 fragments at 120 K is reported in Figure 2. The figure includes all possible t1−t4 couplings between nearby units as well as the couplings between more

Figure 2. |J|2 (eV2 × 1000) for hole and particle transport between the 18 pentacene units contained in the supercell by using the experimental crystal structure obtained at 120 K. The numeration of the fragment is also reported.

distant units, although the latter are very weak (e.g., the 1−5, 3−4, 1−7, 1−8, and 2−7 couplings). The color coding shows that the value of the coupling is consistently the same for geometrically equivalent pairs, e.g., the 1−4 and 8−11 couplings, or the 1−3 and 7−9 couplings. These trends are D


Article

The Journal of Physical Chemistry C

the 9 possible combinations of orbital couplings when a fragment contains three (replicated) units, see Figure 3, and when each fragment is composed of one pentacene unit, see Figure 2. The table reports the mean |J|2 for nearest-neighbor and second nearest-neighbor rows of pentacenes for transport along both a⃗ and b⃗ directions. Not surprisingly, the single- and multiple-unit pentacene fragments are equivalent, consistently with the weak interaction picture discussed above. Nonetheless, from a methodological point of view, our method allows us to define infinite fragments in which the mutual polarization between the pentacene orbitals (used in eq 12) is introduced through the PBC procedure. 3.2. Bilayer Graphene. We build the graphene bilayer based on experimental graphite distances:60−63 1.42 Å for the C−C distance and 3.36 Å for the interlayer distance. We consider a rectangular and a rhombic supercell of approximately 200 carbon atoms each. However, since the results are equivalent, here we discuss the rectangular one and we report the results for the rhombic supercell in Figures S3 and S4. The rectangular supercell is defined by the lattice vectors (in Å) X⃗ = (14.76,0) and Y⃗ = (0,17.04), and a uniform grid of 8 × 6 irreducible k-points. The supercell is shown in Figure 4 together with the single layer crystal orbitals involved in the evaluation of the coupling. Since the two highest occupied COs and the two lowest unoccupied COs are degenerate in the single layer, we will consider the mean |J|2 for each set, i.e., the average of four values for the hole transfer as well as for the particle transfer. The various orientations are obtained as translations from the perfectly stacked AA configuration along the symmetry vectors of graphene, a⃗1 and a2⃗ , which connect the centers of two adjacent hexagons. The center distance is 2.46 Å, and the scan includes 8 × 8 grid points. The translation vectors and three relevant configurations are shown in Figure 5: the AA configuration (maximum coupling), the AB configuration (the most energetically stable), and an overbond configuration (OB, corresponding to the center of the surfaces in Figure 6). Figure 6 reports the change in the coupling with the relative position of the two layers at the same interlayer distance, for both hole and particle transfer. In terms of hole transfer, see the plot in Figure 6i, the coupling is largest at the AA configuration, and it decreases when moving in any direction. This trend can be rationalized by recognizing that the AA configuration provides the best overlap between the fragments and the bilayer orbitals, see eq 12, and moving away from this configuration reduces the overlap and thus the coupling strength. This behavior is in agreement with other computational studies and opposite to what was observed experimentally.64,65 Such discrepancy has been attributed to a distortion of the AA configuration in real systems, either via relative rotation of the layers or via folding (with related increase in interlayer distance), since AA is energetically unfavored.64,65 The coupling in the particle transfer, however, shows a different trend, see the plot in Figure 6ii. While the AA and AB configurations still correspond to large and small coupling, respectively, the OB configuration reveals a large value of the coupling that is comparable to that of AA. In fact, the coupling remains large for translations in the a1⃗ −a2⃗ direction from the AA configuration. This is somewhat surprising, and it requires a more detailed analysis of the contributions to the coupling in eq 12. First of all, we can limit the analysis to only eight bilayer orbitals, which provide the largest contributions to the summation over l in eqs 3 and 12 for all three configurations

consistent also at the other two temperatures, 293 and 414 K, whose results are reported in Figures S1 and S2. The actual values of the four unique couplings, which we can still label t1− t4, are reported in Table 2 at the three temperatures. The data Table 2. Unique |J|2 Values (eV2 × 1000) from PBC Calculations between Pentacene Pairs of the Supercell by Using the Experimental Crystal Structures at 120, 293, and 414 K (See Figure 2) hole 120 K 293 K 414 K

particle

t1

t2

t3

t4

t1

t2

t3

t4

1.7 1.1 0.8

1.4 1.3 1.2

4.8 3.2 2.2

10.9 8.3 6.8

2.9 2.3 2.0

2.7 2.5 2.2

10.2 8.6 6.9

9.8 8.2 6.8

in Figure 2 and Table 2 indicates that the effect of the extended solid on the coupling is weak, but non-negligible. In fact, a comparison between the values in Tables 1 and 2, i.e., between the isolated dimers and the solid, shows differences of about 10% for the larger values of |J|2 at 120 K: tIMD = 8.9 × 10−3 eV2 4 PBC IMD −3 2 −3 vs t4 = 9.8 × 10 eV , or t3 = 9.4 × 10 eV2 vs tPBC = 10.2 3 × 10−3 eV2 for the particle transport. Changes for the hole transport are smaller because occupied orbitals are more compact and thus are less influenced by long-range intermolecular interactions. Smaller values of the coupling show small changes, and the polarization effect decreases at higher temperature, as a result of a complex interplay between changes in intermolecular distances and angles.53 The second aspect that we investigate is the definition of the donor and acceptor fragments. Indeed, one may imagine that the hole/particle is delocalized across multiple units. Therefore, we test this possibility by defining the fragments as infinite 1D chains of pentacene units as shown in Figure 3. Contrary to the case discussed above, here PBC calculations are required also for the donor and acceptor fragments. Figure 3a defines the fragments along the b⃗ lattice vector so that the transport is investigated along the a⃗ direction. In Figure 3b, the definitions of the fragments and transport direction are inverted. In Figure 3, the numeration refers to the ij fragment orbitals for which we compute the coupling, see eq 12, rather than to the fragments themselves; for instance, one fragment now includes units 1, 2, and 3, and each of them contributes one orbital (either HOCO or LUCO) so that for this fragment we consider three orbitals, which in first approximation can be considered linear combination of the orbitals of the isolated units. There is no coupling between orbitals within each fragment since the full KS-PBC equations are solved, and the corresponding blocks in Figure 3 are empty. Figure 3a shows that the coupling is strong only between adjacent fragments (1−3 with 4−6, or 4−6 with 7−9), and it goes to zero for fragments that are further apart. However, while the coupling is completely symmetrical for the particle transport, the hole transport shows stronger coupling between the 1−3 and 4−6 fragments than between the 4−6 and 7−9 fragments. Along the b⃗ transport direction, there is non-negligible coupling also between second nearest-neighbor fragments, e.g., 1−3 with 7−9. This is likely due to the fact that these fragments are spatially closer in this configuration that in case a. As in the previous case, and consistently with the individual-fragments option in Figure 2, the coupling for the particle transport is stronger than for the hole transport. In order to compare the single-unit and multiple-unit fragment definition, we report in Table 3 the mean |J|2 for E


Article


Figure 3. |J|2 (eV2 × 1000) for hole and particle transport between infinite linear fragments by using the experimental crystal structure obtained at 120 K. Each fragment in the supercell contributes three orbitals, and the numeration refers to the orbitals in each fragment. (a) The fragments are defined along the b⃗ direction. (b) The fragments are defined along the a⃗ direction.

energy. According to the same notation, we will refer to the occupied and unoccupied COs of the fragments in Figure 4 as (Hf1,Hf2) and (Lf1,Lf2), respectively. For the purposes of this discussion, we can focus on the fragment-bilayer COs with the largest overlap, although the values for all possible combinations are also reported in Tables S3−S6. For the ⟨Lf1|F̂|Lf1⟩

AA, AB, and OB. These are the highest four occupied COs and the four lowest unoccupied COs, and their energy diagram is shown in Figure 7, while their isodensity surfaces are reported in Figure S5. In order to simplify the notation, we will call the bilayer orbitals (H1b,Hb2,H3b,Hb4) for the occupied set, and (Lb1,Lb2,Lb3,Lb4) for the unoccupied set, in order of increasing F


Article

The Journal of Physical Chemistry C Table 3. Comparison of the Mean |J|2 (eV2 × 1000) between Nearest-Neighbor (Ji,i+1) and Second Nearest-Neighbor Rows/Columns (Ji,i+2) for the Single-Unit (1) and MultipleUnit (3) Fragment Definition along the Transport Direction a⃗ and b⃗ (see Figures 2 and 3) hole |Ji,i+1| , |Ji,i+1|2, |Ji,i+2|2, |Ji,i+2|2, 2

(1) (3) (1) (3)

particle

a⃗

b⃗

a⃗

b⃗

5.2 5.2 0.0 0.0

5.2 5.2 0.6 0.5

6.6 6.7 0.0 0.0

6.6 6.7 1.0 0.9

Figure 7. Energy diagram of the eight bilayer graphene orbitals that provide significant contributions to the coupling in Figure 6.

coupling in AA (FAA L1L1), the largest m contributions come from f b b f b the ⟨L1|H1⟩εH1⟨H1|L1⟩ and ⟨Lf1|Lb1⟩εLb1⟨Lb1|Lf1⟩ terms: 2.15 and −1.71 eV, respectively. For OB (FOB L1L1), the largest values are f b b f f b b from the ⟨L1|H1⟩εH1⟨H1|L1⟩ and ⟨L1|L4⟩εLb4⟨Lb4|Lf1⟩ terms: −2.14 and 1.57 eV, respectively. These values are similar because both the orbital energies and overlap are similar, see Figure 7 and Tables S3−S6. An equivalent behavior is observed for the FAA L2 L2 and FOB couplings, while the mixed couplings F and F L2L2 L1L2 L2 L1 are zero for both configurations. For AB, all m terms with these eight orbitals are close in magnitude but opposite in sign so they tend to cancel out. These trends should now be compared to those for the hole coupling. For the ⟨Hf1|F̂ |Hf1⟩ coupling in f AA (FAA H1H1), the largest m contributions come from the ⟨H1| Hb3⟩εHb3⟨Hb3|Hf1⟩ and ⟨Hf1|Lb3⟩εLb3⟨Lb3|Hf1⟩ terms: 2.42 and −1.82 eV, respectively. For OB (FOB H1H1), the largest values are from the f b b f b ⟨H1|H4⟩εH4⟨H4|H1⟩ and ⟨Hf1|Lb1⟩εLb1⟨Lb1|Hf1⟩ terms: 2.07 and −2.01 eV, respectively. For AA, there is a large difference in energy between the L and H sets so that the magnitude of the m terms is different and there is less cancellation. On the other hand, in OB, the splitting of the orbital energies due to the reduction in symmetry implies a smaller energy difference between the L and H orbitals that mostly contribute to the hole coupling (Hb4 and Lb1, see Figure 7). Hence, these contributions are closer in magnitude and cancel out more. The same

Figure 4. Left: Top view of the rectangular supercell for the AB bilayer graphene with translation vectors X⃗ and Y⃗ . Right: Isodensity surfaces of the single layer orbitals involved in the hole/particle transfer process.

Figure 5. (i) Lattice vectors a⃗1 and a⃗2 used to create the rhombic 8 × 8 displacement grid. Three relevant configurations are also shown: (ii) AA, (iii) AB, and (iv) overbond (OB, where the center of each hexagonal ring is on top of the center of a C−C bond of the other layer).

Figure 6. Mean |J|2 for hole (i) and particle (ii) transfer as a function of the displacement between the graphene layers in the rectangular supercell. The dashed lines used in defining the rhombic 8 × 8 displacement grid are also shown at the bottom. G



■

OB considerations apply to FAA H2H2 and FH2H2, while the mixed

Article

ASSOCIATED CONTENT

S Supporting Information *

couplings FH1H2 and FH2H1 are again zero. Interestingly, even in the hole transfer (as it is in the particle transfer), the value of the overlap for the AA and OB configurations is equivalent, see Tables S3−S6, and the difference between the two is only determined by the energy difference between the COs. For the AB configuration, comparable m terms of opposite sign tend to cancel each other out, providing an overall small magnitude of the coupling also in the hole transfer.

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b06011. Coupling data and figures, isodensity surfaces, and raw data for terms in eq 12 (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 785-864-6509.

4. CONCLUSIONS

Notes

The authors declare no competing financial interest.

We propose a method based on first-principles electronic structure calculations to evaluate the electronic coupling between a donor and an acceptor moiety in solids. This method is an extension to solid systems through periodic boundary conditions within the Γ-point approximation of the Fock/Kohn−Sham (KS) matrix reconstruction method applicable to isolated molecules.19,21 The donor and acceptor moieties can be represented by either finite or infinite fragments, and an arbitrary number of fragments can be defined so that all possible donor−acceptor pair couplings are obtained in one step. We applied this method to two prototypical systems for charge transport: crystalline pentacene and bilayer graphene. For the former system, we show that the isolated molecular dimer approach provides a qualitatively correct description of the coupling for all possible pentacene pairs. However, the polarization induced by the extended solid on the pentacene electron density influences the value of the coupling, with changes of the order of 10% for the largest values of |J|2. For the graphene bilayer, we investigate the dependence of the interlayer coupling on the relative orientation of the layers by performing a rigid translation scan from the perfectly overlapped AA configuration. Our calculations show that, for the hole transfer, the AA configuration is the most favorable, suggesting that the opposite observation in experiments may be due to a deformation of the sample.64,65 For the particle transfer, we find that translations along the a1⃗ −a2⃗ direction from the AA configuration also provide large values of the coupling. A detailed analysis of the various contributions to the coupling reveals that the bilayer orbital energies and the overlap with the single layer COs are both equally important. Indeed, it is the splitting of the orbital energies in the OB configuration (due to symmetry reduction compared to AA) that produces a small value of the coupling in the hole transfer, and not the orbital overlap (which is similar in AA and OB). The value of the current approach, based on the Γ-point approximation, is it simplicity, computational efficiency, and convergence to the k space integration limit for large enough supercells. This approach can provide important qualitative understanding of the electronic coupling in extended systems. In fact, polarization effects of the extended solid on the orbitals used in the coupling evaluation are included since the COs are obtained from a full PBC calculation. Accurate results for smaller supercells will be obtained by lifting the Γ-point approximation, and computing the coupling in the full k space. This development is in progress and will be presented in a future publication.

■

ACKNOWLEDGMENTS The authors would like to thank Prof. Wai-Lun Chan for useful discussions. The authors gratefully acknowledge support from the National Science Foundation under Grant No. OIA1539105, and under Award No. EPS-0903806 and matching support from the State of Kansas through the Kansas Board of Regents. Startup funds from the University of Kansas are also gratefully acknowledged.

■

REFERENCES

(1) Anthony, J. E. Organic electronics: addressing challenges. Nat. Mater. 2014, 13, 773−775. (2) Lanzani, G. Materials for bioelectronics: Organic electronics meets biology. Nat. Mater. 2014, 13, 775−776. (3) Blumberger, J. Recent advances in the theory and molecular simulation of biological electron transfer reactions. Chem. Rev. 2015, 115, 11191−11238. (4) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.; Silbey, R.; Bredas, J.-L. Charge Transport in Organic Semiconductors. Chem. Rev. 2007, 107, 926−952. (5) Troisi, A. Charge transport in high mobility molecular semiconductors: classical models and new theories. Chem. Soc. Rev. 2011, 40, 2347−2358. (6) Winkler, J. R.; Gray, H. B. Electron flow through metalloproteins. Chem. Rev. 2014, 114, 3369−3380. (7) Genereux, J. C.; Barton, J. K. Mechanisms for DNA charge transport. Chem. Rev. 2010, 110, 1642−1662. (8) Mikkelsen, K. V.; Ratner, M. A. Electron tunneling in solid-state electron-transfer reactions. Chem. Rev. 1987, 87, 113−153. (9) Fox, M. A. Introduction-electron transfer: a critical link between subdisciplines in chemistry. Chem. Rev. 1992, 92, 365−368. (10) Marcus, R. Chemical and electrochemical electron-transfer theory. Annu. Rev. Phys. Chem. 1964, 15, 155−196. (11) Marcus, R. A. Electron transfer reactions in chemistry. Theory and experiment. Rev. Mod. Phys. 1993, 65, 599. (12) Karl, N. Charge carrier transport in organic semiconductors. Synth. Met. 2003, 133−134, 649−657. (13) Troisi, A.; Orlandi, G. Dynamics of the intermolecular transfer integral in crystalline organic semiconductors. J. Phys. Chem. A 2006, 110, 4065−4070. (14) Troisi, A. Charge transport in high mobility molecular semiconductors: classical models and new theories. Chem. Soc. Rev. 2011, 40, 2347−2358. (15) Troisi, A. Theories of the charge transport mechanism in ordered organic semiconductors. Adv. Polym. Sci. 2009, 223, 259−300. (16) Yavuz, I.; Martin, B. N.; Park, J.; Houk, K. Theoretical Study of the Molecular Ordering, Paracrystallinity, And Charge Mobilities of Oligomers in Different Crystalline Phases. J. Am. Chem. Soc. 2015, 137, 2856−2866. (17) Gonzalez, J.; Santos, H.; Pacheco, M.; Chico, L.; Brey, L. Electronic transport through bilayer graphene flakes. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 195406. H


Article

The Journal of Physical Chemistry C (18) Coropceanu, V.; Li, H.; Winget, P.; Zhu, L.; Brédas, J.-L. Electronic-Structure Theory of Organic Semiconductors: ChargeTransport Parameters and Metal/Organic Interfaces. Annu. Rev. Mater. Res. 2013, 43, 63−87. (19) Valeev, E. F.; Coropceanu, V.; da Silva Filho, D. A.; Salman, S.; Bredas, J.-L. Effect of electronic polarization on charge-transport parameters in molecular organic semiconductors. J. Am. Chem. Soc. 2006, 128, 9882−9886. (20) Kirkpatrick, J. An approximate method for calculating transfer integrals based on the ZINDO Hamiltonian. Int. J. Quantum Chem. 2008, 108, 51−56. (21) Baumeier, B.; Kirkpatrick, J.; Andrienko, D. Density-functional based determination of intermolecular charge transfer properties for large-scale morphologies. Phys. Chem. Chem. Phys. 2010, 12, 11103− 11113. (22) Cornil, J.; Beljonne, D.; Calbert, J.-P.; Bredas, J.-L. Interchain interactions in organic π-conjugated materials: impact on electronic structure, optical response, and charge transport. Adv. Mater. 2001, 13, 1053−1067. (23) Hutchison, G. R.; Ratner, M. A.; Marks, T. J. Intermolecular charge transfer between heterocyclic oligomers. Effects of heteroatom and molecular packing on hopping transport in organic semiconductors. J. Am. Chem. Soc. 2005, 127, 16866−16881. (24) Haddon, R.; Chi, X.; Itkis, M.; Anthony, J.; Eaton, D.; Siegrist, T.; Mattheus, C.; Palstra, T. Band electronic structure of one-and twodimensional pentacene molecular crystals. J. Phys. Chem. B 2002, 106, 8288−8292. (25) Li, Y.; Coropceanu, V.; Bredas, J.-L. Thermal narrowing of the electronic bandwidths in organic molecular semiconductors: impact of the crystal thermal expansion. J. Phys. Chem. Lett. 2012, 3, 3325−3329. (26) Abergel, D.; Chakraborty, T. Electron correlations in bilayer graphene. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 161409. (27) Bottari, G.; Suanzes, J. A.; Trukhina, O.; Torres, T. Phthalocyanine- carbon nanostructure materials assembled through supramolecular interactions. J. Phys. Chem. Lett. 2011, 2, 905−913. (28) Umeyama, T.; Imahori, H. Photofunctional hybrid nanocarbon materials. J. Phys. Chem. C 2013, 117, 3195−3209. (29) Malig, J.; Jux, N.; Guldi, D. M. Toward Multifunctional Wet Chemically Functionalized Graphene-Integration of Oligomeric, Molecular, and Particulate Building Blocks that Reveal Photoactivity and Redox Activity. Acc. Chem. Res. 2013, 46, 53−64. (30) Matte, H. S. S. R.; Subrahmanyam, K. S.; Rao, K. V.; George, S. J.; Rao, C. N. R. Quenching of fluorescence of aromatic molecules by graphene due to electron transfer. Chem. Phys. Lett. 2011, 506, 260− 264. (31) Dresselhaus, M. What’s Next for Low-Dimensional Materials? Mater. Res. Lett. 2014, 2, 1−9. (32) Meunier, V.; Souza Filho, A. G.; Barros, E. B.; Dresselhaus, M. S. Physical properties of low-dimensional sp2-based carbon nanostructures. Rev. Mod. Phys. 2016, 88, 025005. (33) Wallace, P. R. The band theory of graphite. Phys. Rev. 1947, 71, 622. (34) Novoselov, K. S.; Jiang, Z.; Zhang, Y.; Morozov, S.; Stormer, H.; Zeitler, U.; Maan, J.; Boebinger, G.; Kim, P.; Geim, A. Roomtemperature quantum Hall effect in graphene. Science 2007, 315, 1379−1379. (35) Aoki, M.; Amawashi, H. Dependence of band structures on stacking and field in layered graphene. Solid State Commun. 2007, 142, 123−127. (36) Liu, Z.; Suenaga, K.; Harris, P. J.; Iijima, S. Open and closed edges of graphene layers. Phys. Rev. Lett. 2009, 102, 015501. (37) Lee, J.-K.; Lee, S.-C.; Ahn, J.-P.; Kim, S.-C.; Wilson, J.; John, P. The growth of AA graphite on (111) diamond. J. Chem. Phys. 2008, 129, 234709−234709. (38) Newton, M. D. Quantum chemical probes of electron-transfer kinetics: the nature of donor-acceptor interactions. Chem. Rev. 1991, 91, 767−792. (39) Senthilkumar, K.; Grozema, F.; Bickelhaupt, F.; Siebbeles, L. Charge transport in columnar stacked triphenylenes: Effects of

conformational fluctuations on charge transfer integrals and site energies. J. Chem. Phys. 2003, 119, 9809−9817. (40) Löwdin, P.-O. On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J. Chem. Phys. 1950, 18, 365−375. (41) Giri, G.; Verploegen, E.; Mannsfeld, S. C.; Atahan-Evrenk, S.; Kim, D. H.; Lee, S. Y.; Becerril, H. A.; Aspuru-Guzik, A.; Toney, M. F.; Bao, Z. Tuning charge transport in solution-sheared organic semiconductors using lattice strain. Nature 2011, 480, 504−508. (42) Arjona-Esteban, A.; Krumrain, J.; Liess, A.; Stolte, M.; Huang, L.; Schmidt, D.; Stepanenko, V.; Gsänger, M.; Hertel, D.; Meerholz, K.; et al. Influence of Solid-State Packing of Dipolar Merocyanine Dyes on Transistor and Solar Cell Performances. J. Am. Chem. Soc. 2015, 137, 13524−13534. (43) Dvorak, M.; Oswald, W.; Wu, Z. Bandgap opening by patterning graphene. Sci. Rep. 2013, 3, 2289. (44) 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. et al. Gaussian 09, Revision D.37; Gaussian Inc.: Wallingford, CT, 2009. (45) Izmaylov, A. F.; Scuseria, G. E.; Frisch, M. J. Efficient evaluation of short-range Hartree-Fock exchange in large molecules and periodic systems. J. Chem. Phys. 2006, 125, 104103. (46) Kudin, K. N.; Scuseria, G. E. Linear-scaling density-functional theory with Gaussian orbitals and periodic boundary conditions: Efficient evaluation of energy and forces via the fast multipole method. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 61, 16440. (47) Heyd, J.; Scuseria, G. E. Efficient hybrid density functional calculations in solids: assessment of the Heyd-Scuseria-Ernzerhof screened Coulomb hybrid functional. J. Chem. Phys. 2004, 121, 1187− 1192. (48) Heyd, J.; Scuseria, G. E. Assessment and validation of a screened Coulomb hybrid density functional. J. Chem. Phys. 2004, 120, 7274− 7280. (49) Heyd, J.; Peralta, J. E.; Scuseria, G. E.; Martin, R. L. Energy band gaps and lattice parameters evaluated with the Heyd-ScuseriaErnzerhof screened hybrid functional. J. Chem. Phys. 2005, 123, 174101. (50) Izmaylov, A. F.; Scuseria, G. E.; Frisch, M. J. Efficient evaluation of short-range Hartree-Fock exchange in large molecules and periodic systems. J. Chem. Phys. 2006, 125, 104103. (51) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 2006, 125, 224106− 224106. (52) Henderson, T. M.; Izmaylov, A. F.; Scalmani, G.; Scuseria, G. E. Can short-range hybrids describe long-range-dependent properties? J. Chem. Phys. 2009, 131, 044108. (53) Haas, S.; Batlogg, B.; Besnard, C.; Schiltz, M.; Kloc, C.; Siegrist, T. Large uniaxial negative thermal expansion in pentacene due to steric hindrance. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 205203. (54) Moulton, B.; Zaworotko, M. J. From molecules to crystal engineering: supramolecular isomerism and polymorphism in network solids. Chem. Rev. 2001, 101, 1629−1658. (55) Endres, R.; Fong, C.; Yang, L.; Witte, G.; Wöll, C. Structural and electronic properties of pentacene molecule and molecular pentacene solid. Comput. Mater. Sci. 2004, 29, 362−370. (56) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (57) Campbell, R.; Robertson, J. M.; Trotter, J. The crystal and molecular structure of pentacene. Acta Crystallogr. 1961, 14, 705−711. (58) Park, S.; Kim, S.; Kim, J.; Whang, C.; Im, S. Optical and luminescence characteristics of thermally evaporated pentacene films on Si. Appl. Phys. Lett. 2002, 80, 2872−2874. (59) Cornil, J.; Bredas, J.-L.; Zaumseil, J.; Sirringhaus, H. Ambipolar transport in organic conjugated materials. Adv. Mater. 2007, 19, 1791− 1799. I


Article

The Journal of Physical Chemistry C (60) Neto, A. H. C.; Guinea, F.; Peres, N.; Novoselov, K. S.; Geim, A. K. The electronic properties of graphene. Rev. Mod. Phys. 2009, 81, 109. (61) Crespo, E.; Luque, F.; Barrenechea, J.; Rodas, M. Influence of grinding on graphite crystallinity from experimental and natural data: implications for graphite thermometry and sample preparation. Mineral. Mag. 2006, 70, 697−707. (62) Tsukamoto, J.; Matsumura, K.; Takahashi, T.; Sakoda, K. Structure and conductivity of graphite fibres prepared by pyrolysis of cyanoacetylene. Synth. Met. 1986, 13, 255−264. (63) Birowska, M.; Milowska, K.; Majewski, J. Van der Waals density functionals for graphene layers and graphite. Acta Phys. Pol., A 2011, 120, 845−848. (64) Berashevich, J.; Chakraborty, T. Interlayer repulsion and decoupling effects in stacked turbostratic graphene flakes. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 033403. (65) Rakhmanov, A. L.; Rozhkov, A. V.; Sboychakov, A. O.; Nori, F. Instabilities of the AA-Stacked Graphene Bilayer. Phys. Rev. Lett. 2012, 109, 206801.

J


Evaluation of Electronic Coupling in Solids from Ab Initio Periodic

Recommend Documents