Multiscale Simulations on Charge Transport in Covalent Organic

Subscriber access provided by UNIV OF NEWCASTLE

Article

Multi-Scale Simulations on Charge Transport in Covalent Organic Frameworks: Including Dynamics of Transfer Integrals from FMO-DFTB/LCMO Hirotaka Kitoh-Nishioka, Kai Welke, Yoshio Nishimoto, Dmitri G. Fedorov, and Stephan Irle J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b05779 • Publication Date (Web): 19 Jul 2017 Downloaded from http://pubs.acs.org on July 21, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Multi-Scale Simulations on Charge Transport in Covalent Organic Frameworks: Including Dynamics of Transfer Integrals from FMO-DFTB/LCMO Hirotaka Kitoh-Nishioka,∗,†,‡ Kai Welke,† Yoshio Nishimoto,¶ Dmitri G. Fedorov,§ and Stephan Irle∗,†,∥,⊥ Department of Chemistry, Graduate School of Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan, Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan, Fukui Institute for Fundamental Chemistry, Kyoto University, 34-4 Takano Nishihiraki-cho, Sakyo-ku, Kyoto 606-8103, Japan, Research Center for Computational Design of Advanced Functional Materials (CD-FMat), National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan, and WPI-Research Initiative-Institute of Transformative Bio-Molecules, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan E-mail: [email protected]; [email protected] Phone: +81 29 853 6284. Fax: +81 29 853 6406

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Covalent organic frameworks (COFs) are potential candidates for applications in optoelectronic devices and solar cells due to their ability to transport charge through their aromatic molecular units. The highly ordered π-conjugated TP-COF, consisting of pyrene and triphenylene functional units alternately linked in a mesoporous hexagonal skeleton, is known as the first semiconducting COF. In this paper, we investigate the transport of holes as charge carriers in TP-COF through the π-stacked pyrene units with a multi-scale technique, which combines classical molecular dynamics simulations, quantum chemical calculations, and carrier dynamics simulations. To efficiently estimate the charge transfer integrals from quantum chemical calculations, we developed the FMO-DFTB/LCMO approach by combining the fragment molecular orbital (FMO), density-functional tight-binding (DFTB), and linear-combination of fragment molecular orbitals (LCMO) methods. We observed that the thermal motions of TP-COF cause substantial fluctuations of the transfer integrals. To evaluate the charge carrier diffusion, we performed Ehrenfest dynamics and kinetic Monte Carlo simulations, including the fluctuations of the transfer integrals. Using both simulation approaches, we obtained high carrier mobilities of ca. 2 cm2 V−1 s−1 . We found that the characteristics of charge transport in COFs are similar to oligoacene crystals, suggesting a common mechanism associated with “band-like” transport.

∗

To whom correspondence should be addressed Nagoya University ‡ Tsukuba University ¶ Fukui Institute § AIST ∥ WPI-ITBM ⊥ Current address: Chemical Sciences and Engineering Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 U.S.A. †


Page 2 of 54

Page 3 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION Covalent organic frameworks (COFs) are novel crystalline porous polymers. Since the first work by Yaghi and co-workers in 2005, 1 COFs have attracted much attention to their notable structural features and functionalities, 2–6 including automatic self-integration, low mass density, high thermal stability, gas adsorption and storage, heterogeneous catalysis, semi- and photoconductivity, etc. Two-dimensional (2D) COFs are candidates for the application in optoelectronic devices, transistors, and solar cells because their covalent framework is arranged in 2D sheets with atomic thickness. The sheets can stack on top of each other in an (almost) eclipsed fashion with inter-layer distances of 3.4–4.0 ˚ A. Therefore, 2D COFs allow for a long-range π-orbital overlap in the stacking direction, leading to a path for exciton and charge carrier transport with high mobilities. 2,6,7 Jiang and co-workers 8 reported the first synthesis of a semiconducting COF called TP-COF via the co-condensation reaction of 2,3,6,7,10,11-hexahydroxytriphenylene (HHTP) and pyrenediboronic acid (PDBA), which adopts a belt shape and leads to alternately linked pyrene and triphenylene units in a mesoporous hexagonal skeleton, as shown in Figure 1(a). The electrical conductivity of TP-COF measured across a 10 µm gap between platinum electrodes confirmed (1) a linear current-voltage (I-V) curve, (2) an increase in current upon doping with electron-accepting iodine, and (3) a capability for a repetitive on-off current switching at room temperature. These observations 8 show that TP-COF is a p-type semiconductor, where the charge carrier “hole” probably transports via the π-stacked pyrene units owing to their lower ionization potential 9 compared to the triphenylene units. Other semiconducting COFs have been studied as well, including pyrene COF (PPyCOF), 10 nickel-phthalocyanine based COFs (NiPc-COF 11 and 2D-NiPc-BTDA COF 12 ), other metallo-Pc based COFs (CoPc-, CuPc- and ZnPc-COFs), 13 metallo-porphyrin based COFs (MP-COFs), 14 metal-free porphyrin based COFs (COFs-366 and -66 15 ), and others. 16 The carrier mobilities of several semiconducting COFs were measured by laser flashphotolysis time-resolved microwave conductivity (FP-TRMC). For example, NiPC-COF 11 is 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

a p-type semiconductor with hole mobilities of 1.3 cm2 V−1 s−1 . On the other hand, 2D-NiPcBTDA COF with a different linkage unit 12 is a n-type semiconductor with electron mobilities of 0.6 cm2 V−1 s−1 . COF-66 and COF-366 are p-type semiconductors with hole mobilities as high as 8.1 and 3.9 cm2 V−1 s−1 , respectively. 15 These mobilities are significantly higher than those of discotic liquid crystal semiconductors consisting of phthalocyanine, triphenylene, or pyrene derivatives (10−3 –10−1 cm2 V−1 s−1 ) 17–20 and are comparable to those of amorphous silicon thin-film transistors (0.5–1 cm2 V−1 s−1 ). 21 Therefore, for future application in optoelectronic devices, it is important to understand the microscopic origins of these high mobilities of 2D COFs. Recently, Patwardhan et al. 22 developed a method to calculate the relative stability of localized and delocalized charges in isolated π-stacked triphenylene dimers. They varied twist angles and intermolecular distances to investigate the hole-transport ability of tryphenylenebased COF (T-COF). 1 Consequently, they observed large transfer integrals (> 0.4 eV) and significant charge delocalization within the dimer, when the triphenylene units were eclipsed at intermolecular distances of 3.4–3.6 ˚ A and zero twist angle. Such large transfer integrals significantly exceed those calculated for oligoacene single-molecule crystals (< 0.1 eV, see Table 2 in Ref. 23.) The charge transport in oligoacene crystal appears to be “band-like”, 24–27 with decreasing carrier mobility upon increasing temperature. In addition, some oligoacene crystals, such as pentacene and rubrene, exhibit charge mobilities of a few tens cm2 V−1 s−1 at room temperature. 28–30 Form these results, Patwardhan et al. 22 expected that T-COF can exhibit “band-like” charge transport with a mobility of the order of 10 cm2 V−1 s−1 at room temperature. However, based on a static point of view, their study neglects the electronphonon (e-p) interaction that regulates the charge transport dynamics in conjunction with the amplitude of the transfer integral. 23 Multi-scale simulations 31–40 are prominent methods for understanding the CT dynamics in 2D COFs from a microscopic perspective beyond the static picture. These methods usually combine quantum chemical (QC) calculations of the CT parameters on structures


Page 4 of 54

Page 5 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


obtained from classical molecular dynamics (MD) simulations with subsequent carrier dynamics simulations based on these parameters. The carrier dynamics simulations are often performed using either the Ehrenfest dynamics 31–36,41–46 or kinetic Monte Carlo (KMC) approach. 34,35,37–40 These techniques have been widely used for studying of the CT in organic semiconductors including discotic liquid crystals 33,36,39 and organic solid crystals. 31,32,36–38,40 In particular, some of the studies based on the Ehrenfest dynamics 31,32 or KMC simulations 37,38 clarified the impact of dynamic disorder arising from thermal fluctuations of the transfer integral (associated with non-local e-p interactions) on the CT and reproduced the order of the carrier mobility at room temperature and the experimentally observed “bandlike” behavior. To efficiently obtain the transfer integrals for the multi-scale simulations, we propose a new approach called FMO-DFTB/LCMO by combining the fragment molecular orbital (FMO) 47–51 and self-consistent-charge density-functional tight-binding (SCC-DFTB) methods, 52 using DFTB parameters 53,54 tuned for the CT in organic molecules. This new approach makes use of the FMO-DFTB method, 55,56 followed by transfer integral calculations 57–60 based on the linear-combination of fragment MOs (FMO-LCMO 61,62 ). We apply this approach to the CT dynamics in TP-COF. The transfer integrals of TPCOF calculated from FMO-DFTB/LCMO fluctuates considerably over time with the amplitude of fluctuations being similar to the average value of the transfer integral. The power spectrum of the transfer integral shows a notable peak at low-frequency modes around 25 cm−1 . These results resemble those calculated for oligoacene crystals, 32,37,63 rather than discotic liquid crystals. 33 Both the Ehrenfest dynamics and KMC simulations take into account the dynamics of the transfer integrals and predict high carrier mobilities of ca. 2 cm−1 at room temperature, which partially supports the prediction by Patwardhan et al. 22



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2. METHODS 2.1 Molecular Dynamics Simulations The initial structure for MD simulations is a 3 × 3 × 4 supercell with theoretically derived atomic positions and lattice parameters for the inclined zigzag (Iz ) stacked TP-COF 64 (Figure 1(b)). Based on SCC-DFTB optimizations, 64 the Iz stacking is energetically favored compared to the eclipsed (AA) stacking and its simulated X-ray diffraction (XRD) pattern agrees with the experiment. 8 We performed MD simulations of a 3 × 3 × 4 supercell with periodic boundary conditions using TINKER 7.1. 65 The supercell consists of 108 pyrene, 72 triphenylene, and 216 boronic acid units, with the total of 4 968 atoms. Standard MM3 parameters 66 were used in addition to parameters for boronic acids in COFs. 67,68 More details are provided in the Supporting Information. For a minimal model for TP-COF (one pyrene unit linked to one triphenylene unit via one boronic acid, see Figure S1(a) in the Supporting Information), geometry optimizations and normal modes using these molecular mechanics (MM) force field parameters can reproduce those obtained by B3LYP 69–71 with cc-pVDZ. (Figures S1(b) and (c) in the Supporting Information). The MD simulations were performed under NVT conditions using the Bussi–Parrinello thermostat 72 and the RESPA integrator 73 with the volume fixed to the value provided in Ref. 64 at all temperatures. The integration time step was set to 2.0 fs using the RATTLE 74 algorithm. Particle Mesh Ewald (PME) was used for the long-range electrostatic interactions with a real-space cutoff of 12 ˚ A. The overall system was heated to the target temperatures 200 K and 300 K at a rate of 50 K per 10 ps and equilibrated during the following 500 ps. Structures were collected every 60 fs over the following 492 ps, yielding 213 = 8 192 snapshots at each temperature. These snapshots were used for subsequent electronic structure calculations for determining the dynamics of the transfer integral.


Page 6 of 54

Page 7 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2 FMO-DFTB/LCMO Approach to the Transfer Integral Calculation For the efficient estimation of transfer integrals in large organic semiconductor systems, we propose a novel method, called FMO-DFTB/LCMO, by making use of the following three approaches; (1) FMO-LCMO for the treatment of charge transport phenomena, 59 (2) FMODFTB to be able to treat large systems, 55,56 and (3) the unconfined parameter set of DFTB and the proposed uniform scaling factor to the transfer integral. 53,54 First, we briefly review the FMO-LCMO method, 57–62 whose details have been described in Refs. 61 and 62. The current work is based on the FMO2 47,48 version with fragment dimer exchange corrections, although the extension to trimer corrections (FMO3 49 ) is straightforward. 60,62 In FMO the total molecule is first divided into N fragments. The electronic structure of each fragment is then solved self-consistently under the electrostatic potential (ESP) from all other fragments. (MO)

The pth orbital of fragment I and the corresponding MO energy are denoted by φIp and ϵIp

respectively, which are optimized self-consistently under the ESP from the other fragments. (MO)

The subscript Ip of ϵIp

refers to φIp . Next, the electronic structure of each fragment dimer is

solved under the ESP from all other fragments determined above. The resultant ath orbital (MO)

of fragment dimer IJ and the corresponding MO energy are denoted by φIJ a and ϵIJa , respectively. FMO-LCMO provides intra- and inter-fragment Hamiltonian matrix elements as (total)

HIp,Iq = ⟨φIp |hI |φIq ⟩ +

∑ {⟨

⟩ ⟨ ⟩} φIp |hIJ |φIq − φIp |hI |φIq ,

(1)

J̸=I (total)

HIp,Jq = ⟨φIp |hIJ |φJq ⟩ (I ̸= J).


(2)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 54

Here, “effective” one-electron Hamiltonians, such as Fock and Kohn-Sham operators, of fragment monomers I and dimers IJ are written as

hI =

∑

(MO)

|φIp ⟩ϵIp

⟨φIp |,

(3)

p∈I

hIJ =

∑

(MO)

IJ |φIJ a ⟩ϵIJa ⟨φa |.

(4)

a∈IJ

Since the dimer Hamiltonian matrices in eqs 1 and 2 are projected to the monomer orbital space, the total Hamiltonian matrix H total is represented in the fragment monomer MO basis. When bond-detached atoms (BDA) are present, one should remove from eqs 1, 2, 3, and 4 the φIp and φIJ a with huge MO energies that are an effect of the hybrid orbital projection (HOP) operators 61,62 as described later. The MOs and corresponding energies for the total system can be obtained by solving a generalized eigenvalue problem of H total with the overlap matrix among monomer fragment molecular orbitals, SIp,Jq = ⟨φIp |φJq ⟩. 57,58,60–62 The FMO method usually employs the electrostatic dimer (ES-DIM) approximation 48 that avoids self-consistent field calculations of far separated dimers. When the ES-DIM ap{⟨ ⟩ ⟨ ⟩} proximation is used for dimer IK, one should remove the term φIp |hIK |φIq − φIp |hI |φIq (total) Kq

from the summation in eq 1 and set HIp,

to zero. 62

Ref. 59 demonstrated that the diagonal and non-diagonal matrix elements of H (total) can be used for the site energies ϵIp and transfer integrals in CT systems. To consider the nonorthogonality between the fragment monomer MOs in those cases, one should correct the transfer integrals by using the Löwdin’s symmetric transformation as follows: 59 ( TIp,

Jq

=γ (

(total) HIp, Jq

− SIp,

(total) Jq

− SIp,

= γ HIp,

( Jq

(total) HIp, Ip

+

(total) HJq, Jq

)

) 2 /2 /(1 − SIp,

) 2 Jq (ϵIp + ϵJq ) /2 /(1 − SIp,

Jq ),

Jq )

(5)

where γ is a uniform scaling factor, which is always set to unity for FMO-HF and FMODFT calculations. Hereafter the “corrected transfer integral” TIp,


Jq

is just called “transfer

Page 9 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


integral”. The site energies (eq 1) and transfer integrals, (eq 5) naturally include the orbital relaxation from the self-consistent ESP and the pairwise exchange interactions. Note that (MO)

the site energy ϵIp is different from the fragment monomer MO energy ϵIp (MO)

reflects only the polarization effect. 59 (Nevertheless, the use of ϵIp

, the latter

instead of ϵIp in eq 5

may be useful to study CT. 75 ) This study employs the SCC-DFTB method for the calculation of the fragments’ electronic structures within the framework of FMO, which is called FMO-DFTB. 55,56 The standard DFTB parameters are derived using a slightly compressed atomic basis set, which is suitable for many molecular and solid-state applications. However, using the standard DFTB parameters significantly underestimates the transfer integrals between the fragment MOs at greater distances. 76 Reducing the confinement in the parameterization procedure for DFTB can approximate transfer integrals between distant fragment MOs obtained from DFT calculations with double-ζ basis sets. 76 In this study, we used an approach similar to the recent HAB11 53 and HAB7− 54 benchmark studies as follows: (1) standard confined DFTB parameters are used for the fragment monomer calculations; (2) the same parameters are also used for the intra-layer fragment-dimer calculations; (3) for inter-layer fragment-dimer DFTB calculations, much less confined parameters are used for the X–Y atomic-pair interactions, where X belongs to one fragment and Y belongs to another. On the other hand, the standard confined parameters are used for the X–Y atomic-pair interactions where both X and Y belong to the same fragment; (4) a uniform scaling factor, γ in eq 5, is set to 1.54 for the calculation of hole transfer integrals by following Ref. 53. The unconfined parameters, called “8-∞”, are identical to those in Refs. 53 and 54. Parameters related to boron were prepared following the same procedure as reported in Ref. 53. Figure 2(a) shows the TP-COF model on which we performed the FMO-DFTB calculations. The coordinates are taken from the central portion of each MD snapshot, which is marked by a dashed circle in Figure 1(b). The dangling bonds produced by the extraction of the TP-COF model from the MD snapshots are capped by hydrogen atoms, as denoted



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

by bold ‘H’ letters in Figure 2(b). The TP-COF model consists of 4 pyrene, 8 triphenylene, and 8 boronic acid units, totaling 344 atoms. In the FMO-DFTB calculations, the TP-COF model is divided into 12 fragments, as shown in Figure 2(a). For fragmenting intra-layer moieties, we detached B–C bonds at the sp2 boron atoms, thereby treating the boron atoms as BDAs. In this study, we simply used the DFTB atomic orbitals of boron for the HOP operators. (2s, 2px , 2py atomic orbitals of boron in Figure 2(b). We confirmed that the FMO-DFTB and FMO-DFTB/LCMO calculations with these HOPs can approximate the total energy and canonical MOs from the full SCC-DFTB calculation, as described below. As the confined DFTB parameters, we used the mio-0-1 set 52 for X–Y element-pair interactions (X, Y = C, O, H), the rscc-materials set 77 for B–X (X = C, O) interactions, and the borg-0-1 set 78 for B–X (X = B, H) interactions. The effects of the mixing of such DFTB parameters on the transfer integral calculations are discussed below. In this study, we calculate the transfer integrals between the highest occupied fragment MOs of adjacent pyrene units of the TP-COF model with eq 5, which are denoted by T1,2 , T2,3 , and T3,4 . Hereafter, the four indices in eq 5 are always simplified as two indices by dropping orbital’s indices because we only treat the highest occupied MO each pyrene molecule/unit for the hole transfer integral calculations. Note that in eq 4 we take into account all dimer MOs except those with huge MO energies arising from the HOPs. The detail of such purifications related to BDAs in the FMO-LCMO method is given in Refs. 61,62. The ES-DIM approximation 48 was used with the threshold value of 2.0, which negligibly affects the results as shown below. All FMO calculations were performed using the GAMESS program, 79 which was locally modified to allow for FMO-DFTB calculations with mixed confined and unconfined parameters.


Page 10 of 54

Page 11 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3 Reorganization Energy In addition to the time-dependent transfer integral obtained from the mixed MD and FMODFTB/LCMO approaches, the reorganization energy λ is a key parameter for the carrier dynamics simulations performed in this study. The reorganization energy consists of intramolecular and outer-sphere (or external) contributions. The outer-sphere reorganization energy plays the dominant role in solution-phase and biological charge transfer reactions. In contrast, the intra-molecular reorganization energy becomes dominant in most molecular solids because of their low dielectric constants. This study, therefore, focuses only on the intra-molecular reorganization energy, similar to CT studies on molecular crystals of oligoacenes and their derivatives. 32,34,35,37,38,40 (The impact of the external reorganization energy on CTs in oligoacene crystals has been investigated by McMahon and Troisi. 80 ) Since TP-COF is regarded as an infinite 2D-sheet polymer, the four model systems in Figure 2(c) are used for the estimation of λ of the pyrene moiety as charge carrier in TP-COF. The intramolecular reorganization energy for hole transfer can be obtained from the adiabatic potential energy surfaces of the neutral and cationic charged model molecules as 23,81 ∗ λ = (E+1 − E+1 ) + (E0∗ − E0 ),

(6)

where E0 and E+1 are the energy of the neutral state and the energy of the cation molecular state, respectively. E0∗ is the energy of the neutral molecule at the optimal cation geometry ∗ and E+1 is the energy of the cation state at the optimal geometry of the neutral molecule. ∗ The optimized geometries of the neutral molecule are used for the estimation of E0 and E+1 ,

while the optimized geometries of the cation are used for the estimation of E+1 and E0∗ . All geometry optimizations were performed under D2h symmetry restriction with the Gaussian 09 Revision D.01 suite of programs. 82



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 54

2.4 Carrier Dynamics Simulations To simulate the carrier transport through the columnar stacked pyrene moieties of TP-COF, we performed Ehrenfest dynamics and kinetic Monte Carlo simulations with the charge transport parameters obtained from the QC calculations. Both simulations calculate the diffusion coefficient D of the carrier charge. The carrier mobility µ can be calculated by means of the Einstein-Smoluchowski relation as 32,37,38,42,43

µ=

e D. kB T

(7)

2.4.1 Ehrenfest Dynamics Simulations We briefly review the Ehrenfest dynamics (ED) 32,42,43 used in this study. We consider a one-dimensional N ED molecular stack with a periodic boundary, in which each molecule j is associated with one electronic state |j⟩. In this study, |j⟩ represents the highest occupied molecular orbital, HOMO, of each pyrene moiety of TP-COF in the stacking direction. The Ehrenfest dynamics simulation employs the Su–Schrieffer–Heeger (SSH)-type Hamiltonian; 32,35,42,43,83 H(t) =

∑ j

αxj,1 |j⟩⟨j| +

∑

[−τ + β(xj+1,2 − xj,2 )](|j⟩⟨j + 1| + |j + 1⟩⟨j|)

j

∑ 1 ∑ 1 1 1 2 2 + ( m1 vj,1 + m1 (ω1 xj,1 )2 ) + ( m2 vj,2 + m2 (ω2 xj,2 )2 ), 2 2 2 2 j j

(8)

where we consider the following two kinds of harmonic nuclear vibrations (indexed with k = 1, 2) on each molecule j with displacement xj,k , velocity vj,k , effective mass mk , and frequency ωk : (1) the vibration indexed with k = 1 is associated with the Holstein local electron-phonon (e-p) coupling and (2) the vibration indexed with k = 2 is associated with the Peierls nonlocal e-p coupling. In eq 8, the vibration xj,1 linearly modulates the site energy on |j⟩ with the local e-p coupling constant α, leading to the term αxj,1 |j⟩⟨j|. On the other hand, (xj+1,2 − xj,2 ) linearly modulates the transfer integral between the molecules j and j + 1 with 12 ACS Paragon Plus Environment

Page 13 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the non-local e-p coupling constant β, leading to the term β(xj+1,2 −xj,2 )(|j⟩⟨j+1|+|j+1⟩⟨j|); −τ represents the transfer integral at the equilibrium geometry xj,2 = xj+1,2 = 0. Here, ∑ we express the wave function describing the charge carrier as ψ(t) = j ψj |j⟩. The first ˙ = H(t)ψ(t)/i¯ derivative of ψ(t), the time-dependent Schrödinger equation, ψ(t) h, is written as ψ˙ j = {αxj,1 ψj − τ (ψj+1 + ψj−1 )

(9)

+ β[(xj,2 − xj−1,2 )ψj−1 + (xj+1,2 − xj,2 )ψj+1 ]}/i¯h. The motions of the nuclear vibrations xj,k are obtained by solving a set of Newton’s equation mk v˙ j,k = −∂⟨ψ(t)|H(t)|ψ(t)⟩/∂xj,k , which are explicitly written as: v˙ j,1 = (−m1 ω12 xj,1 − αψj∗ ψj )/m1 ,

(10)

∗ ∗ v˙ j,2 = [−m2 ω22 xj,2 − β(ψj∗ ψj+1 + ψj+1 ψj − ψj−1 ψj − ψj∗ ψj−1 )]/m2 .

(11)

and

The mean squared displacement (MSD) is written as

MSD(t) =

M 1 ∑ (m) [⟨ψ (t)|r2 |ψ (m) (t)⟩ − ⟨ψ (m) (t)|r|ψ (m) (t)⟩2 ], M m=1

(12)

where M is the number of the realizations. The matrix elements can be computed from the relations ⟨j|r2 |k⟩ = δj,k j 2 L2 and ⟨j|r|k⟩ = δj,k jL, where L is set to the equilibrium distance between the closest neighbor molecules. In this study, L = 3.6 ˚ A was taken from the DFTB optimized structure. 64 In normal diffusive regime, there is a linear relationship between MSD and the simulation time. Therefore, the diffusion coefficient in the one-dimensional case can be evaluated through the time derivative of the MSD as

D=

1 dMSD(t) lim . 2 t→∞ dt 13

ACS Paragon Plus Environment

(13)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Here, |j⟩ are represented by the highest occupied fragment MOs of the pyrene units of the TP-COF model as illustrated in Figure 2(a). The parameters associated with the transfer integrals, τ , β, and ω2 are, therefore, estimated from the FMO-DFTB/LCMO calculations performed on the TP-COF model. On the other hand, the Holstein-type e-p coupling con√ stant α is associated with the reorganization energy λ through the equation α = ω1 m1 λ. 32 By following Ref. 32, we assume that the C-C stretching mode dominates the reorganization energy and set m1 and ω1 to 6 amu and 1400 cm−1 respectively. The single effective mode can approximate the CT in organic semiconductors and electron transfer reactions in organic and biological molecules. 84 The system size N ED is set to 600 sites with periodic boundary condition (|N ED +1⟩ = |1⟩ and |0⟩ = |N ED ⟩ ). 32,42,43 The initial nuclear positions {xj,k (t = 0)} and velocities {vj,k (t = 0)} for each realization m at temperature T were randomly chosen to satisfy the Boltzmann distribution for a set of non-interacting harmonic oscillators. 32,42,43 We take the center site |N ED /2⟩ = |300⟩ as the initial wavefunction ψ (m) (t = 0) for each realization m. 42,43 The numerical integrations of eqs 9, 10, and 11 are performed using the fourth-order Runge–Kutta (RK4) algorithm 85 with a time interval of 0.025 fs. We carried out M = 1000 realizations to obtain the MSD. The simulation settings are summarized in Table 1.

2.4.2 Kinetic Monte Carlo Simulations In addition to the Ehrenfest dynamics simulations based on the SSH Hamiltonian, we performed kinetic Monte Carlo (KMC) simulations to examine the charge carrier diffusion. Assuming the incoherent hopping regime for the carrier diffusion, KMC enables us to solve a set of kinetic master equations for the time-evolution of the population of the carrier charge ∑ at each molecule, ∂Pj (t)/∂t = i (ki,j (t)Pi (t) − kj,i (t)Pj (t)), where Pj is the occupation number of the carrier charge at molecule j and kj,i is the charge transfer rate from j to its neighboring molecule i. 43 Here, we again consider a one-dimensional molecular stack, in which each molecule j is associated with one electronic state |j⟩. Hence, this KMC approach


Page 14 of 54

Page 15 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


is actually an one-dimensional random-walk simulation with forward and backward hopping rates, kj,j+1 and kj−1,j , respectively, from site j. 37,38,43 The rate kj,j+1 (t) is expressed by Fermi’s golden rule: 84,86

kj,j+1 (t) =

2π KMC 2 |T (t)| FC(ϵj,j+1 (t)), h ¯ j,j+1

(14)

KMC where FC is the Franck–Condon factor, ϵj,j+1 is the site energy difference, and Tj,j+1 (t) is

the time-dependent transfer integral between sites j and j +1 in the KMC simulations. As in the SSH Hamiltonian used in the Ehrenfest dynamics simulations, |j⟩ represents the HOMO of each pyrene moiety of TP-COF in the π-stacking direction. To take into account both local and non-local e-p couplings in the KMC simulations, we follow the works by Shuai and co-workers. 37,38 Although their studies 37,38 used multiple quantum modes with the CT rate formula by Nan et al., 87 our study employs a single effective quantum mode ⟨ω⟩ as the local e-p interaction one with the FC factor derived by Jortner 84,86

FC(ϵ) =

√ n ¯ + 1 p/2 1 exp[−S(2¯ n + 1)]Ip (2S n ¯ (¯ n + 1))[ ] . h ¯ ⟨ω⟩ n ¯

(15)

Here, p = ϵ/¯h⟨ω⟩ is the effective-frequency normalized energy gap, and S = λ/¯ h⟨ω⟩ is the effective-frequency normalized reorganization energy (or Huang–Rhys factor). n ¯ = 1/(exp[⟨ω⟩/kB T ] − 1) is the averaged quantum number of the effective mode at temperature T . Ip (x) is the modified Bessel function of the first kind of the pth order. In the high-temperature limit, kB T ≫ h ¯ ⟨ω⟩, eq 15 can be approximated as the well-known semiclassical Marcus formula 88 as follows: 84,86 [ ] 1 (−ϵ − λ)2 FC(ϵ) = √ exp − , 4λkB T 4πλkB T

(16)

We again regard |j⟩ as the highest occupied fragment MOs of the pyrene units of the TP-COF model as illustrated in Figure 2(a). We take into account the dynamic disorder



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 54

of the transfer integral in the KMC simulations by following Refs. 37,38. The time-course transfer integral Tj,j+1 (t) calculated with FMO-DFTB/LCMO for the MD snapshots with the time interval ∆t is written as

Tj,j+1 (∆t × h) = ⟨Tj,j+1 ⟩ + Vj,j+1 (∆t × h) (h = 1, 2, · · · , N snap ),

(17)

where N snap = 213 represents the total number of MD snapshots. A discrete Fourier transformation of Vj,j+1 (t) provides

Vkj,j+1

=

1 N snap

snap N ∑

h=1

[

] k×h exp 2πi snap Vj,j+1 (∆t × h). N

(18)

KMC Using eq 18, we can express the transfer integral Tj,j+1 (t) for arbitrary times t and arbitrary

pairs j and j + 1 in the KMC simulations as 37,38 (

N snap /2

KMC Tj,j+1 (t)

∑

2πk t + φj,j+1 cos = ⟨Tj,j+1 ⟩ + 2 0,k snap N ∆t k=0 ( ) N snap /2 ∑ 2πk t j,j+1 j,j+1 ImVk sin +2 + φ0,k , snap ∆t N k=0

)

ReVkj,j+1

(19)

where the phase factor φj,j+1 is chosen randomly to make the transfer integral of each 0,k (j, j + 1) site-pair fluctuate independently. At each realization for the KMC simulation, the phase factor for each Fourier component of each (j, j + 1) site-pair is randomly set as rωk ∆tN snap , where r is a random number uniformly distributed in [0, 1]. 37,38 In the KMC simulation, 37,38 we first create a one dimensional system consisting of N KMC = 100 sites with periodic boundary condition and choose the site j = N KMC /2 = 50 as the starting center of charge at each realization. The hopping rates for the carrier charge at site j to the forward site j + 1 and to the backward site j − 1 are calculated using eq 14. The forward hopping probability is kj,j+1 /(kj,j+1 + kj,j−1 ) and the backward hopping probability is kj,j−1 /(kj,j+1 + kj,j−1 ). By using a random number r uniformly distributed in [0, 1], we 16 ACS Paragon Plus Environment

Page 17 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


determine which hopping is accepted. In the case of r ≤ kj,j+1 /(kj,j+1 + kj,j−1 ), the forward hopping is chosen, otherwise the backward hopping is chosen. After the determination of the next site of the charge, the time t of the KMC simulation is incremented by the residence time of the charge on site j, 1/(kj,j+1 + kj,j−1 ). Then, the transfer integrals are updated at the new time (eq 19). The KMC simulation at each realization terminates when the time exceeds the alloted total time. The diffusion constant can be calculated from the obtained MSD(t). using eq 13. As in the Ehrenfest dynamics simulations, for the single effective mode in eq 15, we assume that the C–C stretching mode dominantly contributes to the reorganization energy and set ⟨ω⟩ to 1400 cm−1 . 32,84 We carried out M = 5000 realizations to obtain the MSD. The simulation settings are summarized in Table 1.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3. RESULTS and DISCUSSION 3.1 Benchmarks of DFTB, FMO-DFTB, and FMO-DFTB/LCMO Calculations First, we confirmed that DFTB with the mio-0-1 52 parameter set provides adequate shapes of the HOMOs of isolated pyrene and triphenylene molecules, compared with B3LYP/6-31G(d), as shown in Figure S2 in the Supporting Information. For the benchmark, the geometry of each molecule was optimized using B3LYP/cc-pVTZ. The HOMO energies of pyrene and triphenylene calculated with DFTB are 5.62 and 6.02 eV respectively, corresponding to the site energy of each isolated molecule for the “hole” transfer. The ionization potentials (IPs) of pyrene and triphenylene as determined by experiment are 7.41 and 7.88 eV respectively. 9 Therefore the HOMO energy differences of pyrene and triphenylene agree well between calculation and experiment, although the calculated HOMO energies are significantly lower than the experimental IPs. Table S1 in the Supporting Information additionally lists the HOMO energies calculated with B3LYP/6-31G(d) and LC-BLYP/6-31G(d). Next, we checked the applicability of the atomic orbitals of boron to the HOPs in the FMO-DFTB calculations. For the benchmark, we adopted a model molecule consisting of one pyrene, two triphenylene, and two boronic acid units (Figure 2(b)), where two BDAs are involved in the FMO calculations. The model was first optimized using full DFTB (without fragmentation). Subsequently, the total energies of FMO-DFTB (−116.1446387 hartree) and DFTB (−116.1445937 hartree) were found to be very similar. Thus, the fragmentation results in a negligible error of 0.045 milli-hartree. Finally, we calculated the MO energy spectrum of the same model molecule by diagonalizing the FMO-DFTB/LCMO Hamiltonian matrix,. Here, in eqs 1, 2, 3, and 4, we consider all monomer and dimer MOs except those with anomalous MO energies arising from HOPs. In Figure S3 in the Supporting Information, it can be seen that the obtained occupied MO energy spectrum is consistent with full DFTB. Although the fragmented bonds in this study 18 ACS Paragon Plus Environment

Page 18 of 54

Page 19 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are closer than those in typical FMO calculations, 50 FMO-LCMO reasonably recovers the MOs of the entire molecule in a wide energy range using the electronic structures of fragments. 57–62 The HOMO of the entire molecule is localized on the pyrene fragment, while HOMO-1 and HOMO-2 of the molecule are delocalized over the two triphenylene fragments. As described in Section 2.2, this study focuses on the transfer integrals (eq 5) between the highest occupied fragment MOs of adjacent pyrene fragments of the TP-COF model (Figure 2(a)). Since FMO-DFTB/LCMO reasonably approximates the occupied MOs of the TP-COF model, we can expect that the calculated transfer integral adequately includes the effects of the triphenylene fragments covalently linked via the boronic acids.

3.2 Assessment of Hole Transfer Integrals of the Pyrene Homodimer Here we show the performance of the “8-∞” DFTB parameters for the estimation of hole transfer integrals of the pyrene homodimer because they are not included in the HAB11 and HAB7− databases. 53,54 The reference data was obtained using eq 5 with LC-BLYP 89 because it can reproduce transfer integrals from high-level ab initio calculations like MRCI+Q and NEVPT2. 59 For the LC-BLYP calculations, we used the cc-pVTZ basis set for carbon atoms and cc-pVDZ for hydrogens. The range-separation parameter for was set to 0.33. For comparison, we also performed FMO-DFT/LCMO calculations with B3LYP and PBE 90 with the 6-31G(d) basis sets. First, the structure of one neutral pyrene molecule was optimized using B3LYP/cc-pVDZ under D2h symmetry. Then we made a cofacial pyrene dimer with D2h symmetry by superposing one optimized monomer on top of another with a displacement along the z-axis perpendicular to the planar molecular planes, as shown in the inset of Figure 3(a). The FMO calculations treat one pyrene monomer as one fragment. Figure 3(a) shows the resultant hole transfer integrals Ti,j for the cofacial pyrene dimer as a function of intermolecular distance, d. Figure 3(b) shows the Ti,j in cofacial pyrene dimers when one pyrene monomer is translated 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

along its long axis, as shown in the inset of Figure 3(b). Here we designate the FMODFTB/LCMO calculation with a uniform scaling factor γ (eq 5) of 1.54 as sDFTB. In the calculations labelled DFTB on the other hand, γ = 1. The resultant Ti,j from sDFTB (red closed circles in Figure 3) are in reasonable agreement with the reference values (black plus symbols) compared to B3LYP and PBE with 6-31G(d). Therefore, we can expect that the application of the “8-∞” DFTB parameters with γ = 1.540 to the FMO-DFTB/LCMO approach can provide good-quality transfer integrals associated with the hole transport through the pyrene units of TP-COF. For the rest of this work we always use sDFTB.

3.3 Transfer Integral Fluctuations in TP-COF To reproduce the thermal fluctuations of TP-COF, we performed classical MD simulations for the 3 × 3 × 4 supercell (Figure 1(b)) with periodic boundary conditions. The dynamics of the transfer integrals were investigated at target temperatures 200 K and 300 K, although we mainly focus on the results at T = 300 K. We plot the total energy and the coordinate rootmean-square deviation (RMSD) of the atomic coordinates as a function of time in Figure S4 in the Supporting Information. To address the dynamics of the transfer integral in TP-COF, we performed FMODFTB/LCMO calculations on the TP-COF model as shown in Figure 2(a). The coordinates were extracted from MD snapshots (dotted circle in Figure 1(b)). We especially focus on T2,3 , which is the transfer integral between the HOMOs of fragments 2 and 3 calculated with eq 5, because it represents the transfer integral responsible for the hole transport through the columnar-stacked pyrene units in TP-COF. ( See Figure 2(a). ) Figures 4(a) and (b) show T2,3 as a function of time and the normalized probability distributions of T2,3 , respectively. The average values ⟨Tj,j+1 ⟩ and standard deviations σ{Tj,j+1 } of the transfer integrals calculated with eq 5 for MD snapshots at the two target temperatures are listed in Table 2. σ{Tj,j+1 } is larger than ⟨Tj,j+1 ⟩, reflecting large thermal fluctuations of the transfer integrals. Figure 4(b) shows that T2,3 at 300 K seems to follow a normal distribution (green 20 ACS Paragon Plus Environment

Page 20 of 54

Page 21 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


solid line). On the other hand, T2,3 at 200 K somewhat deviates from a normal distribution (red solid line), which is probably due to insufficient sampling. The distributions of T1,2 and T3,4 resemble that of T2,3 (Figure S5 in the Supporting Information), which supports, that the fluctuation of Ti,j between adjacent π-stacked pyrene units in TP-COF indeed follow a normal probability distribution. Now, we show the correlation between the fluctuations of Ti,j using the mutual correlation coefficient between the two transfer integrals

rab,cd =

⟨(Ta,b − ⟨Ta,b ⟩) (Tc,d − ⟨Tc,d ⟩)⟩ . σ{Ta,b }σ{Tc,d }

(20)

When Ta,b and Tc,d fluctuate synchronously, rab,cd approaches 1 or −1. On the other hand, when Ta,b and Tc,d fluctuate independently, rab,cd ≈ 0. Using Ti,j at T = 300 K, we obtained correlation coefficients r12,23 of 0.077, r23,34 of 0.082, and r12,34 of 0.025. The negligible correlation between the transfer integral fluctuations have been observed in oligoacene crystals as well. 63 The lack of the correlation supports the assumption of random initial phases (eq 19) used for the KMC simulation. Next we calculated the power spectrum of T2,3 at 300 K by discrete Fourier transformation of the auto-correlation function A(t) = ⟨V2,3 (t)V2,3 (0)⟩. (See eq 17 for the definition of the thermal deviation amplitude Vi,j (t).) We employed a time resolution ∆t of 60 fs and the data of A(t) from −∆t × (211 − 1) to ∆t × 211 with the application of a Welch window function. 85 As shown in Figure 5(a), only low-frequency modes under 65 cm−1 contribute to the modulation of the transfer integral with strong peaks around 25 cm−1 . We also performed a discrete Fourier transformation of V2,3 (t) without using a window function with the same time resolution of 60 fs, using all N snap = 213 data points. Figure 5(b) shows the resulting real (ReV ) and imaginary (ImV ) parts, which were used for reproducing KMC (t) in the KMC simulation (eq 19). Naturally, the time-dependent transfer integral Tj,j+1

ReV and ImV exhibit features similar to the power spectrum shown in Figure 5(a).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The properties of the transfer integral dynamics of TP-COF, such as a normal probabilitylike distribution and a distinctive peak in the power spectrum around 25 cm−1 , resemble features observed in oligoacene crystals. 32,37,38,63 ¿From these properties, Troisi and coworkers 31,32 assumed that the modes contributing to the transfer integral dynamics (namely, Peierls’ modes) behave as simple harmonic oscillators, as represented by xj,2 in the SSH Hamiltonian (eq 8). Therefore, we use a simple harmonic oscillator with a frequency ω2 of 25 cm−1 for the Ehrenfest dynamics simulations as described in Section 2.4.1. In contrast, the power spectrum of the transfer integral fluctuations obtained for discotic liquid crystals shows broad and structureless bands in the 0 – 200 cm−1 region. 33 In that case, the simple harmonic oscillator in eq 8 should be replaced by multiple Langevin oscillators reproducing its power spectrum. 33,91 As explained in Section 2.2, we combined four different kinds of Slater–Koster (SK) parameters in the FMO-DFTB calculations, which are the unconfined “8–∞” and standard confined “mio-0-1”, “rscc-materials”, and “borg-0-1” parameters. As shown in Section 3.2, the combination of mio-0-1 with “8–∞” can provide reasonable hole transfer integrals in the pyrene homodimer. However, the mixing of different DFTB parameter sets is still a source of concern. To check the effects of this combination on the transfer integrals, we used only the “matsci-0-3” parameter set 77 as the confined parameters. Figure S6 in the Supporting Information shows the comparison between T2,3 calculated with “matsci-0-3” only and T2,3 calculated with the mixed SK parameters. We can see that there are negligible differences between them. This study employs the ES-DIM approximation in FMO calculations. The effect of the approximation on the transfer integrals is introduced via the term (ϵIp + ϵJq ) in eq 5. Figure S7 in the Supporting Information shows the correlation between T2,3 calculated with and without ES-DIM approximation. We can see that the both calculations produced almost identical transfer integrals.


Page 22 of 54

Page 23 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.4 Reorganization Energy Using eq 6 for the four molecules I to IV illustrated in Figure 2(c), we estimated the intramolecular reorganization energy λ of the CT in TP-COF. Geometry optimizations were performed with restricted and unrestricted B3LYP/cc-pVDZ methods. For I to IV, the values for λ were 0.1505, 0.2148, 0.09987, and 0.1512 eV respectively. ¿From these results, we set λ = 0.15 eV for the carrier dynamics simulations in TP-COF with the Ehrenfest dynamics and KMC approaches.

3.5 Carrier Dynamics from Ehrenfest Dynamics Simulations In this subsection, we analyze the carrier dynamics in TP-COF obtained by Ehrenfest dynamics simulations. Table 3 lists the representative parameters used for the SSH Hamiltonian (eq 8). The transfer integral τ was set to the average value from the FMO-DFTB/LCMO calculations for the MD snapshots: τ = ⟨T2,3 ⟩ = 38.07 meV or 307 cm−1 (Table 2). The mass of Peierls’ mode, m2 , was set to the mass of an isolated pyrene molecule, 200 amu. The frequency of Peierls’ mode, ω2 , was set to 25 cm−1 , at which the strong peaks of the power spectrum of T2,3 were observed, (Figure 5(a)). Assuming that only the single mode xj,2 yields a “normal probability”-like distribution of T2,3 shown in Figure 4(b), the non-local √ e-p coupling β can be calculated from σ{T2,3 }/ 2kB T /m2 (ω2 )2 . 32 The determination of the parameters related to Holstein’s mode (m1 , ω1 , and α) was explained in Section 2.4.1; the √ local e-p coupling α was calculated from ω1 m1 λ, where λ was set to 0.15 eV as described above. These values are also summarized in Table 3. Figure 6(a) shows the MSD(t)/L2 at T = 300 K. The MSD(t) calculated with m2 = 2000 amu and m2 = 200 amu are almost identical, indicating that the results are independent of m2 . Keeping m1 and λ constant, we varied ω1 , as shown in Figure 6(a). As a result, the reduction in ω1 between 800 cm−1 and 1 400 cm−1 slightly decreased the MSD(t). Therefore, we regard µ calculated with ω1 = 1400 cm−1 as representative, similar to previous studies. 32 MSD(t)/L2 can be well fitted using a straight line a × t + b. We confirmed that the number 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of reliazation and ED sites (namely, M = 1000 and N ED = 600) listed in Table 1 is sufficient to obtain a converged MSD(t)/L2 , as shown in Figure S8(a) in the Supporting Information. The diffusion coefficient D can be calculated via eq 13 using the slope a as L2 ×a/2. Using eq 7 with the calculated D, we estimate a carrier mobility µ of 1.93 cm2 V−1 s−1 . The calculated µ for TP-COF is close to the lower limit of 1.3 cm2 V−1 s−1 estimated by experiments in NiPc COF. 11 The Ehrenfest dynamics approach treats the nuclear vibrations classically. Since the frequency of ω2 is estimated as 25 cm−1 , we can expect that Peierls’ mode works classically in the carrier dynamics in TP-COF. On the other hand, we set the frequency of the Holstein’s mode, ω1 , to as high as 1400 cm−1 . In such a high frequency regime, the quantum effect of the nuclear vibration on the carrier dynamics might become apparent, which is discussed below. It is also known that Ehrenfest dynamics fails to reproduce the Boltzmann distribution of quantum states. 92 To avoid this problem, Kranz and Elstner used Boltzmann-corrected Ehrenfest dynamics to study singlet exciton diffusion in anthracene crystal. 93 Another prominent technique for mixed quantum-classical dynamics simulations of charge/exciton transport is the fewest-switches surface-hopping with adaptions to handle trivial crossing. 43,93,94 Applying these techniques to CT simulations in TP-COF may be advantageous in future extensions of this work.

3.6 Carrier Dynamics from KMC Simulations We next show the carrier dynamics in TP-COF observed by using KMC simulations. The rate constant for each incoherent hopping step (eq 14) was calculated using either a a quantummechanically treated (eq 15), or a classical-limit Franck–Condon factor (eq 16). Hereafter we call the former QM CT and the latter Marcus’ CT. We used the following values for the parameters in eqs 15 and 16; the site energy difference ϵj,j+1 was set to 0, T was set to 300 K, λ was set to 0.15 eV, and ⟨ω⟩ was set to 1 400 cm−1 , as in the Ehrenfest dynamics 24 ACS Paragon Plus Environment

Page 24 of 54

Page 25 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulations. Moreover, we accounted for the dynamics of the transfer integral by using eq 19. Similar to the Ehrenfest dynamics simulations, the values of the parameters used for eq 19 were determined from the statistical parameters of T2,3 from the FMO-DFTB/LCMO calculations; ⟨Tj,j+1 ⟩ was set to 38.07 meV. The ReVkj,j+1 and ImVkj,j+1 used are plotted in Figure 5(b). The values of MSD(t)/L2 involving the QM CT hopping are plotted in Figure 6(b). We confirmed that the number of realization and KMC sites (namely, M = 100 and N KMC = 5000) listed in Table 1 is sufficient to obtain a converged MSD(t)/L2 , as shown in Figure S8(b) in the Supporting Information. We fitted MSD(t)/L2 to a simple straight line a × t + b, which estimates the diffusion coefficient D at L2 × a/2 via eq 13. As a result, eq 7 with the obtained D predicts a high carrier mobility of µ = 2.067 cm2 V−1 s−1 . The values of MSD(t)/L2 involving Marcus’ CT hopping are also plotted in Figure 6(b). which results in the µ of 1.025 cm2 V−1 s−1 . We can see that the quantum effects of the nuclear high-frequency mode (nuclear tunneling effects) double the predicted carrier mobility in TP-COF. This study incorporates the dynamics of the transfer integrals from FMO-DFTB/LCMO calculations for the MD snapshots into the KMC simulations via eq 19. To clarify the impact of such the dynamic disorder of the transfer integral on the carrier diffusion, we also performed the KMC KMC simulations involving the QM CT by replacing Tj,j+1 (t) with the static ⟨T2,3 ⟩ and

plotted the resultant MSD(t)/L2 in Figure 6(b) for comparison. Consequently, we can see that the dynamic disorder of the transfer integral increases the carrier mobility, similar to previous theoretical studies. 38 The carrier mobilities are summarized in Table 4.

3.7 Comparison of Ehrenfest Dynamics Results with KMC Results For comparison, the MSD(t)/L2 from the Ehrenfest dynamics simulations (red line in Figure 6(a)) is shown as a black line in Figure 6(b). Interestingly, µ = 2.067 cm2 V−1 s−1 from the KMC simulations involving QM CT is close to µ = 1.93 cm2 V−1 s−1 from the Ehrenfest dynamics simulations. Note that the assumptions, on which the Ehrenfest dynamics approach 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

with the SSH Hamiltonian (eq 8) and the KMC approach with the QM CT hopping (eq 15) are based, differ from each other. The former allows coherent motions of the carrier (the wavepacket of the carrier can delocalize over several pyrene units), but ignores the nuclear tunneling effect. On the other hand, the KMC approach takes into account nuclear tunneling via eq 15, but ignores the coherent motions of the carrier because the incoherent hopping model is used (the carrier always localizes on a single molecule for each hopping step). As a result, the dynamic disorder of the transfer integral affects the carrier mobilities in Ehrenfest dynamics simulations in an opposite way compared to KMC simulations. In Ehrenfest dynamics simulations, the disorder (namely, increasing β) enhances localization of the charge carrier, leading to the reduction of the carrier mobility. 27,31,32,42,43 In KMC on the other hand, the carrier mobility at room temperature is increased, as shown in Figure 6(b) and presented in the previous theoretical studies. 38 To provide clarification about mechanisms of CT in COFs, both ED and KMC approaches are applied in this study. Spectroscopic experiments, such as electron spin resonance and charge modulated spectroscopy, could be very useful as they may be able to provide details of the localization of charge carriers, 27,38 which could give evidence for a verification which approach is more appropriate for further studies on semiconductive 2D-COFs. As mentioned above, the characteristics of the obtained CT parameters of TP-COF are similar to those of oligoacene crystals. As a result, the Ehrenfest dynamics simulations with the CT parameters listed in Table 3 indicate “band-like” transport behavior in TP-COF, as shown in Figure S9 in the Supporting Information. On the other hand, KMC simulations with the data of the time-dependent transfer integral from the FMO-DFTB/LCMO calculations perfomed on the MD trajectories at 300 K and with the nuclear tunneling effect yield a nearly temperature-independent carrier mobility in TP-COF because the FC factor, estimated by eq 15 with ⟨ω⟩ = 1400 cm−1 and λ = 0.15 eV, is almost independent of the temperature. Note that the conventional “hopping” regime 23,27 in which the carrier mobility is larger for higher temperature, arises mainly from the temperature dependence of the semi-


Page 26 of 54

Page 27 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


classical Marcus formula (see eq 16), which corresponds to the high-temperature limit of eq 15. The fact that eq 15 yields a temperature-independent FC factor supports the band-like behavior of CT in TP-COF. For an accurate examination of the temperature dependence of carrier mobilites within the framework of the KMC approach, the time-dependent transfer integral at each target temperature should be obtained using the MD trajectory and multiple quantum modes should be taken into consideration, as demonstrated by Shuai and co-workers. 37,38,87 Even though KMC simulations are based on the hopping model, with multiple quantum modes they can reproduce the band-like behavior depending on the CT parameters. 37,38,87 Our results from the two types of simulations support the band-like transport behavior in TP-COF rather than the conventional hopping mechanism. In a recent theoretical study by Zhao and co-workers, 95 a time-dependent wavepacket diffusion method was developed that essentially incorporates the carrier coherence and tunneling effects simultaneously, and they argued that the band-like behavior arises from the nuclear tunneling effect rather than the carrier coherence motion. However, this interesting issue on the origin of the band-like behavior is left for a future study.

3.8 Advantages of FMO-DFTB/LCMO Here we highlight some advantages of using FMO-DFTB/LCMO for studying the CT in COFs and other organic semiconductors. We performed QC calculations on a TP-COF model (Figure 2(a)) to obtain the thermal motion of the hole transfer integral between adjacent πstacked pyrene units. Since the model system is small, FMO-DFTB/LCMO has little advantage over full DFTB in computational cost. Regarding CT studies using full DFT or DFTB calculations, the fragment orbital (FO) approach is usually employed. 36,45,46,53,54,59,76,96 In these conventional FO approaches, the MOs of isolated fragments are regarded as basis states for the estimation of the site energies and transfer integrals. In the cases of the oligoacene crystals, the fragments contributing to the CT are spatially well-defined, probably guaranteeing the choice of the basis states in the FO approaches. On the other hand, 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in the case of COFs consisting of organic building blocks linked via two or three covalent bonds, it is difficult to define reasonable basis states for describing the CT. When the conventional FO approach is used for TP-COF and the basis states are taken from the MOs of isolated pyrene units with hydrogen-capped dangling bonds, the effect of the two adjacent triphenylene units and boronic acid linkers are neglected. If one pyrene unit linked to two adjacent triphenylene units via boronic acid linkers on the same sheet is treated as one fragment (molecule III in Figure 2(c)), other artificial boundary effects arise from the four adjacent boronic acid linkers and pyrene units. By making use of the HOP scheme, FMODFTB/LCMO can circumvent the difficulty in the definition of the basis states as described in Section 2.2. The FMO φIp appear naturally as basis states. As shown in Figure S3 in the Supporting Information, the MO energy spectrum obtained from the diagonalization of the FMO-DFTB/LCMO Hamiltonian matrix supports the validity of the basis states in our approach. DFTB using standard confined parameters is not suitable for the calculation of transfer integrals, as we described above. Since the TP-COF model used in this study is rather small, FMO-DFT calculations performed on several thousands of MD snapshots are still feasible. So far, FMO–LCMO at the Hartree-Fock/DFT level of theory have provided reasonable transfer integrals when BDAs were not involved or when a minimal atomic basis set was used. 57,59 However, the application of FMO-DFT/LCMO with double-ζ basis sets to TPCOF faces the problem 60,62 of linear dependent basis set functions associated with BDAs on boron atoms. Although the three–body correction (FMO3) remarkably reduces the error associated with BDAs, 60 it increases the computational cost. Therefore, in addition to its low computational cost, FMO-DFTB/LCMO has an advantage in the definition of the basis states describing the CT in semiconduncting COFs, compared to the FO approach based on full DFTB and FMO-LCMO based on DFT with double-ζ basis sets. The FMO-DFTB/LCMO approach provides both site energies ϵIp (eq 1) and transfer integrals Ti,j (eq 5). Because of the crystal’s uniformity and the use of a single effective


Page 28 of 54

Page 29 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mode for Holstein’s local e-p interactions, this study neglects the thermal motion of the site energies, in line with previous theoretical works. 31–35,37,38 In their studies on CT in liquid crystals and DNA, Elstner and co-workers 36,45,46 incorporated the dynamics of the site energies into the carrier dynamics simulations, where both the fast local e-p coupling modes and the slow relaxation modes are included. In addition to soft matter, systems associated with spatially inhomogeneous electronic properties, such as bulk heterojunction solar cells, require the estimation of the site energy from QC calculations for describing their CT properties. The FMO-LCMO approach at DFT level of theory provided accurate site energies (eq 1) in DNA including the orbital relaxation from environmental self-consistent ESP and pairwise exchange interactions. 59 Corresponding accuracy tests of the site energies from FMO-DFTB/LCMO calculations with “8-∞” DFTB parameters are left for future work. The secondary objective of this study is to probe the FMO-DFTB/LCMO method itself towards future applications to larger CT systems. ¿From this study, we expect that FMODFTB/LCMO works well for other COFs sharing common building blocks with TP-COF. However, since the accuracy tests reported in this paper are quite limited, we need further examinations to show if the FMO-DFTB/LCMO is generally applicable to describe CT in various systems.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4. CONCLUDING REMARKS In this study, we investigated the “hole” carrier transport through π-stacked pyrene units in TP-COF using a multi-scale approach combining classical MD simulations, QC calculations, and carrier dynamics simulations. To efficiently estimate the transfer integral from the QC calculations, we have developed the FMO-DFTB/LCMO method with the “8-∞” DFTB parameter set and demonstrated that the method can provide accurate hole transfer integrals for the π-stacked pyrene dimer (Figure 3). For π-stacked pyrene units in realistic TP-COF structures, we found that the calculated transfer integrals fluctuate considerably over time (Figure 4(a)). The standard deviation of the transfer integral, σ{Ti,j }, is larger than the corresponding average value, ⟨Ti,j ⟩, (Table 2), and Ti,j is normally distributed, as shown in Figure 4(b). The power spectrum of Ti,j (Figure 5) shows that only low-frequency modes under 65 cm−1 contribute to the fluctuations with a strong peak around 25 cm−1 . To incorporate these Ti,j fluctuations into carrier dynamics simulations, we used the Ehrenfest dynamics and KMC approaches with both Holstein’s local and Peierls’ non-local e-p couplings estimated from the QC calculations. Interestingly, both approaches provided similar carrier mobilities of ca. 2 cm2 V−1 s−1 without using adjustable CT parameters, which to the best of our knowledge is the first theoretical prediction in COFs. Although Ehrenfest dynamics and KMC simulations produced similar carrier mobilities in TP-COF, there is a large difference in them regarding the treatment of coherent motion of the carriers, leading to the opposite dependence of the resultant carrier mobility on the dynamic disorder of the transfer integral: the dynamic disorder reduces the carrier mobility in the former, but increases the mobility in the latter (Figure 6(b)). Future spectroscopic experiments on the charge localization may be able to discern whether the Ehrenfest dynamics or KMC approach is better for simulating CT in COFs and designing future 2D-COFs with a desirable semiconductivity. CT properties, including Ti,j fluctuations, and high carrier mobilities in TP-COF, are calculated in a “first-principle” manner in this study; they are similar to the properties of 30 ACS Paragon Plus Environment

Page 30 of 54

Page 31 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


oligoacene crystals rather than liquid crystals. Consequently, the carrier dynamics simulations with the CT properties in COFs support the band-like transport behavior (Figure S9 in the Supporting Information) rather than the conventional hopping mechanism, observed in oligoacene crystals. This study indicates a necessity of experimental studies on the temperature dependence of the carrier mobility in COFs to obtain further physical insight into their CT mechanism. Although measurements of the carrier mobilities have been performed on a limited number of semiconducting COFs, some of them reported carrier mobilities of ca. 1 cm2 V−1 s−1 and above. 11,12,15 The calculated carrier mobility of TP-COF of 2 cm2 V−1 s−1 is consistent with these experiments, suggesting that this study on TP-COF elucidates general features of CT in semiconducting COFs. To perform multi-scale simulations with FMO-DFTB/LCMO on semiconducting COFs, for which measured carrier mobilities are available, new classical MM force fields need to designed in order to perform MD simulations. Alternatively, integrating FMO-DFTB 56,97 into Ehrenfest MD and KMC can provide a consistent QM approach to study charge transfer in materials.

Acknowledgement S.I. acknowledges partial support by a CREST grant from JST. H.K.-N. also acknowledges support from Collaborative Research Program for Young Scientists of ACCMS and IIMC, Kyoto University. K.W. acknowledges support from the JSPS.

Supporting Information Available Details of the force-field preparation, additional structural, vibrational, and one-electron properties, dependence of transfer integrals, and results of additional carrier dynamics simulations. This material is available free of charge via the Internet at http://pubs.acs.org/.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Côté, A. P.; Benin, A. I.; Ockwig, N. W.; O’Keeffe, M.; Matzger, A. J.; Yaghi, O. M. Porous, Crystalline, Covalent Organic Frameworks. Science 2005, 310, 1166–1170. (2) Feng, X.; Ding, X.; Jiang, D. Covalent Organic Frameworks. Chem. Soc. Rev. 2012, 41, 6010–6022. (3) El-Kaderi, H. M.; Hunt, J. R.; Mendoza-Cortés, J. L.; Côté, A. P.; Taylor, R. E.; O’Keeffe, M.; Yaghi, O. M. Designed Synthesis of 3D Covalent Organic Frameworks. Science 2007, 316, 268–272. (4) Ding, S.-Y.; Wang, W. Covalent Organic Frameworks (COFs): From Design to Applications. Chem. Soc. Rev. 2013, 42, 548–568. (5) Xiang, Z.; Cao, D. Porous Covalent-Organic Materials: Synthesis, Clean Energy Application and Design. J. Mater. Chem. A 2013, 1, 2691–2718. (6) Dogru, M.; Bein, T. On the Road towards Electroactive Covalent Organic Frameworks. Chem. Commun. 2014, 50, 5531–5546. (7) Colson, J. W.; Dichtel, W. R. Rationally Synthesized Two-Dimensional Polymers. Nat. Chem. 2013, 5, 453–465. (8) Wan, S.; Guo, J.; Kim, J.; Ihee, H.; Jiang, D. A Belt-Shaped, Blue Luminescent, and Semiconducting Covalent Organic Framework. Angew. Chem. Int. Ed. 2008, 47, 8826– 8830. (9) Montalti, M.; Credi, A.; Prodi, L.; Gandolfi, M. T. Handbook of Photochemistry. Thrid Edition; CRC Press, 2006; Chapter 7, pp 493–527. (10) Wan, S.; Guo, J.; Kim, J.; Ihee, H.; Jiang, D. A Photoconductive Covalent Organic Framework: Self-Condensed Arene Cubes Composed of Eclipsed 2D Polypyrene Sheets for Photocurrent Generation. Angew. Chem. Int. Ed. 2009, 48, 5439–5442. 32 ACS Paragon Plus Environment

Page 32 of 54

Page 33 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(11) Ding, X.; Guo, J.; Feng, X.; Honsho, Y.; Guo, J.; Seki, S.; Maitarad, P.; Saeki, A.; Nagase, S.; Jiang, D. Synthesis of Metallophthalocyanine Covalent Organic Frameworks That Exhibit High Carrier Mobility and Photoconductivity. Angew. Chem. Int. Ed. 2011, 50, 1289–1293. (12) Ding, X.; Chen, L.; Honsho, Y.; Feng, X.; Saengsawang, O.; Guo, J.; Saeki, A.; Seki, S.; Irle, S.; Nagase, S. et al. An n-Channel Two-Dimensional Covalent Organic Framework. J. Am. Chem. Soc. 2011, 133, 14510–14513. (13) Ding, X.; Feng, X.; Saeki, A.; Seki, S.; Nagai, A.; Jiang, D. Conducting Metallophthalocyanine 2D Covalent Organic Frameworks: The Role of Central Metals in Controlling π-Electronic Functions. Chem. Commun. 2012, 48, 8952–8954. (14) Feng, X.; Liu, L.; Honsho, Y.; Saeki, A.; Seki, S.; Irle, S.; Dong, Y.; Nagai, A.; Jiang, D. High-Rate Charge-Carrier Transport in Porphyrin Covalent Organic Frameworks: Switching from Hole to Electron to Ambipolar Conduction. Angew. Chem. Int. Ed. 2012, 51, 2618–2622. (15) Wan, S.; Gándara, F.; Asano, A.; Furukawa, H.; Saeki, A.; Dey, S. K.; Liao, L.; Ambrogio, M. W.; Botros, Y. Y.; Duan, X. et al. Covalent Organic Frameworks with High Charge Carrier Mobility. Chem. Mater. 2011, 23, 4094–4097. (16) Chen, Y.; Cui, H.; Zhang, J.; Zhao, K.; Ding, D.; Guo, J.; Li, L.; Tian, Z.; Tang, Z. Surface Growth of Highly Oriented Covalent Organic Framework Thin Film with Enhanced Photoresponse Speed. RSC Adv. 2015, 5, 92573–92576. (17) Warman, J. M.; Van De Craats, A. M. Charge Mobility in Discotic Materials Studied by Pr-Trmc. Mol. Cryst. Liq. Cryst. 2003, 396, 41–72. (18) Senthilkumar, K.; Grozema, F. C.; Bickelhaupt, F. M.; Siebbeles, L. D. A. Charge Transport in Columnar Stacked Triphenylenes: Effects of Conformational Fluctuations on Charge Transfer Integrals and Site Energies. J. Chem. Phys. 2003, 119, 9809–9817. 33 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(19) Fujikake, H.; Murashige, T.; Sugibayashi, M.; Ohta, K. Time-of-Flight Analysis of Charge Mobility in a Cu-Phthalocyanine-Based Discotic Liquid Crystal Semiconductor. Appl. Phys. Lett. 2004, 85, 3474–3476. (20) Sienkowska, M. J.; Monobe, H.; Kaszynski, P.; Shimizu, Y. Photoconductivity of Liquid Crystalline Derivatives of Pyrene and Carbazole. J. Mater. Chem. 2007, 17, 1392–1398. (21) Gleskova, H.; Hsu, P.; Xi, Z.; Sturm, J.; Suo, Z.; Wagner, S. Field-Effect Mobility of Amorphous Silicon Thin-Film Transistors under Strain. J. Non-Cryst. Solids 2004, 338, 732 – 735. (22) Patwardhan, S.; Kocherzhenko, A. A.; Grozema, F. C.; Siebbeles, L. D. A. Delocalization and Mobility of Charge Carriers in Covalent Organic Frameworks. J. Phys. Chem. C 2011, 115, 11768–11772. (23) 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. (24) Warta, W.; Karl, N. Hot Holes in Naphthalene: High, Electric-Field-Dependent Mobilities. Phys. Rev. B 1985, 32, 1172–1182. (25) Cheng, Y. C.; Silbey, R. J.; da Silva Filho, D. A.; Calbert, J. P.; Cornil, J.; Brédas, J. L. Three-Dimensional Band Structure and Bandlike Mobility in Oligoacene Single Crystals: A Theoretical Investigation. J. Chem. Phys. 2003, 118, 3764–3774. (26) Hannewald, K.; Bobbert, P. A. Ab Initio Theory of Charge-Carrier Conduction in Ultrapure Organic Crystals. Appl. Phys. Lett. 2004, 85, 1535–1537. (27) Troisi, A. Charge Transport in High Mobility Molecular Semiconductors: Classical Models and New Theories. Chem. Soc. Rev. 2011, 40, 2347–2358. (28) Jurchescu, O. D.; Baas, J.; Palstra, T. T. M. Electronic Transport Properties of Pentacene Single Crystals upon Exposure to Air. Appl. Phys. Lett. 2005, 87, 052102. 34 ACS Paragon Plus Environment

Page 34 of 54

Page 35 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(29) Podzorov, V.; Menard, E.; Rogers, J. A.; Gershenson, M. E. Hall Effect in the Accumulation Layers on the Surface of Organic Semiconductors. Phys. Rev. Lett. 2005, 95, 226601. (30) Takeya, J.; Yamagishi, M.; Tominari, Y.; Hirahara, R.; Nakazawa, Y.; Nishikawa, T.; Kawase, T.; Shimoda, T.; Ogawa, S. Very High-Mobility Organic Single-Crystal Transistors with In-Crystal Conduction Channels. Appl. Phys. Lett. 2007, 90, 102120. (31) 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. (32) Troisi, A. Prediction of the Absolute Charge Mobility of Molecular Semiconductors: The Case of Rubrene. Adv. Mater. 2007, 19, 2000–2004. (33) Troisi, A.; Cheung, D. L.; Andrienko, D. Charge Transport in Semiconductors with Multiscale Conformational Dynamics. Phys. Rev. Lett. 2009, 102, 116602. (34) Vehoff, T.; Baumeier, B.; Troisi, A.; Andrienko, D. Charge Transport in Organic Crystals: Role of Disorder and Topological Connectivity. J. Am. Chem. Soc. 2010, 132, 11702–11708. (35) Vehoff, T.; Chung, Y. S.; Johnston, K.; Troisi, A.; Yoon, D. Y.; Andrienko, D. Charge Transport in Self-Assembled Semiconducting Organic Layers: Role of Dynamic and Static Disorder. J. Phys. Chem. C 2010, 114, 10592–10597. (36) Heck, A.; Kranz, J. J.; Kubaˇr, T.; Elstner, M. Multi-Scale Approach to Non-Adiabatic Charge Transport in High-Mobility Organic Semiconductors. J. Chem. Theory Comput. 2015, 11, 5068–5082. (37) Wang, L.; Li, Q.; Shuai, Z.; Chen, L.; Shi, Q. Multiscale Study of Charge Mobility of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Organic Semiconductor with Dynamic Disorders. Phys. Chem. Chem. Phys. 2010, 12, 3309–3314. (38) Geng, H.; Peng, Q.; Wang, L.; Li, H.; Liao, Y.; Ma, Z.; Shuai, Z. Toward Quantitative Prediction of Charge Mobility in Organic Semiconductors: Tunneling Enabled Hopping Model. Adv. Mater. 2012, 24, 3568–3572. (39) Kirkpatrick, J.; Marcon, V.; Nelson, J.; Kremer, K.; Andrienko, D. Charge Mobility of Discotic Mesophases: A Multiscale Quantum and Classical Study. Phys. Rev. Lett. 2007, 98, 227402. (40) Yin, S.; Li, L.; Yang, Y.; Reimers, J. R. Challenges for the Accurate Simulation of Anisotropic Charge Mobilities through Organic Molecular Crystals: The β Phase of mer-Tris(8-hydroxyquinolinato)aluminum(III) (Alq3) Crystal. J. Phys. Chem. C 2012, 116, 14826–14836. (41) Hultell, M.; Stafström, S. Polaron Dynamics in Highly Ordered Molecular Crystals. Chem. Phys. Lett. 2006, 428, 446 – 450. (42) Wang, L.; Beljonne, D.; Chen, L.; Shi, Q. Mixed Quantum-Classical Simulations of Charge Transport in Organic Materials: Numerical Benchmark of the Su-SchriefferHeeger Model. J. Chem. Phys. 2011, 134, 244116. (43) Wang, L.; Beljonne, D. Charge Transport in Organic Semiconductors: Assessment of the Mean Field Theory in the Hopping Regime. J. Chem. Phys. 2013, 139, 064316. (44) Ciuchi, S.; Fratini, S.; Mayou, D. Transient Localization in Crystalline Organic Semiconductors. Phys. Rev. B 2011, 83, 081202. (45) Kubaˇr, T.; Elstner, M. Coarse-Grained Time-Dependent Density Functional Simulation of Charge Transfer in Complex Systems: Application to Hole Transfer in DNA. J. Phys. Chem. B 2010, 114, 11221–11240. 36 ACS Paragon Plus Environment

Page 36 of 54

Page 37 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(46) Kubaˇr, T.; Elstner, M. Efficient Algorithms for the Simulation of Non-Adiabatic Electron Transfer in Complex Molecular Systems: Application to DNA. Phys. Chem. Chem. Phys. 2013, 15, 5794–5813. (47) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: An Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701 – 706. (48) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Use of Approximate Electrostatic Potential. Chem. Phys. Lett. 2002, 351, 475 – 480. (49) Fedorov, D. G.; Kitaura, K. The Importance of Three-Body Terms in the Fragment Molecular Orbital Method. J. Chem. Phys. 2004, 120, 6832–6840. (50) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2012, 14, 7562–7577. (51) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. ElectronCorrelated Fragment-Molecular-Orbital Calculations for Biomolecular and Nano Systems. Phys. Chem. Chem. Phys. 2014, 16, 10310–10344. (52) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260–7268. (53) Kubas, A.; Hoffmann, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic Couplings for Molecular Charge Transfer: Benchmarking CDFT, FODFT, and FODFTB against High-Level Ab Initio Calculations. J. Chem. Phys. 2014, 140, 104105. (54) Kubas, A.; Gajdos, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic Couplings for Molecular Charge Transfer: Benchmarking CDFT, FODFT and 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FODFTB against High-Level Ab Initio Calculations. II. Phys. Chem. Chem. Phys. 2015, 17, 14342–14354. (55) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2014, 10, 4801–4812. (56) Nishimoto, Y.; Nakata, H.; Fedorov, D. G.; Irle, S. Large-Scale Quantum-Mechanical Molecular Dynamics Simulations Using Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. J. Phys. Chem. Lett. 2015, 6, 5034– 5039. (57) Nishioka, H.; Ando, K. Electronic Coupling Calculation and Pathway Analysis of Electron Transfer Reaction Using Ab Initio Fragment-Based Method. I. FMO-LCMO Approach. J. Chem. Phys. 2011, 134, 204109. (58) Kitoh-Nishioka, H.; Ando, K. Fragment Molecular Orbital Study on Electron Tunneling Mechanisms in Bacterial Photosynthetic Reaction Center. J. Phys. Chem. B 2012, 116, 12933–12945. (59) Kitoh-Nishioka, H.; Ando, K. Charge-Transfer Matrix Elements by FMO-LCMO Approach: Hole Transfer in DNA with Parameter Tuned Range-Separated DFT. Chem. Phys. Lett. 2015, 621, 96 – 101. (60) Kitoh-Nishioka, H.; Ando, K. FMO3-LCMO Study of Electron Transfer Coupling Matrix Element and Pathway: Application to Hole Transfer between Two Tryptophans through Cis- and Trans-Polyproline-Linker Systems. J. Chem. Phys. 2016, 145, 114103. (61) Tsuneyuki, S.; Kobori, T.; Akagi, K.; Sodeyama, K.; Terakura, K.; Fukuyama, H. Molecular Orbital Calculation of Biomolecules with Fragment Molecular Orbitals. Chem. Phys. Lett. 2009, 476, 104 – 108.


Page 38 of 54

Page 39 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(62) Kobori, T.; Sodeyama, K.; Otsuka, T.; Tateyama, Y.; Tsuneyuki, S. Trimer Effects in Fragment Molecular Orbital-Linear Combination of Molecular Orbitals Calculation of One-Electron Orbitals for Biomolecules. J. Chem. Phys. 2013, 139, 094113. (63) Troisi, A.; Orlandi, G. Dynamics of the Intermolecular Transfer Integral in Crystalline Organic Semiconductors. J. Phys. Chem. A 2006, 110, 4065–4070. (64) Lukose, B.; Kuc, A.; Frenzel, J.; Heine, T. On the Reticular Construction Concept of Covalent Organic Frameworks. Beilstein J. Nanotechnol. 2010, 1, 60–70. (65) Ponder, J. W. Tinker - Software Tools for Molecular Design, version 7.1. http:// dasher.wustle.edu/tinker/. (66) Allinger, N. L.; Li, F.; Yan, L.; Tai, J. C. Molecular Mechanics (MM3) Calculations on Conjugated Hydrocarbons. J. Comput. Chem. 1990, 11, 868–895. (67) Schmid, R.; Tafipolsky, M. An Accurate Force Field Model for the Strain Energy Analysis of the Covalent Organic Framework COF-102. J. Am. Chem. Soc. 2008, 130, 12600–12601. (68) Amirjalayer, S.; Snurr, R. Q.; Schmid, R. Prediction of Structure and Properties of Boron-Based Covalent Organic Frameworks by a First-Principles Derived Force Field. J. Phys. Chem. C 2012, 116, 4921–4929. (69) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (70) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627. (71) Hertwig, R. H.; Koch, W. On the Parameterization of the Local Correlation Functional. What is Becke-3-LYP? Chem. Phys. Lett. 1997, 268, 345 – 351. 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(72) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (73) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. J. Chem. Phys. 1992, 97, 1990–2001. (74) Andersen, H. C. Rattle: A ”Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24 – 34. (75) Fujita, T.; Haketa, Y.; Maeda, H.; Yamamoto, T. Relating Stacking Structures and Charge Transport in Crystal Polymorphs of the Pyrrole-Based π-Conjugated Molecule. Org. Electron. 2017, 49, 53 – 63. (76) Kubaˇr, T.; Woiczikowski, P. B.; Cuniberti, G.; Elstner, M. Efficient Calculation of Charge-Transfer Matrix Elements for Hole Transfer in DNA. J. Phys. Chem. B 2008, 112, 7937–7947. (77) Frenzel, J.; Oliveira, A. F.; Jardillier, N.; Heine, T.; Seifert, G. Semi-Relativistic, SelfConsistent Charge Slater-Koster Tables for Density-Functional Based Tight-Binding (DFTB) for Materials Sciece Simulations; 2004-2009. (78) Grundkötter-Stock, B.; Bezugly, V.; Kunstmann, J.; Cuniberti, G.; Frauenheim, T.; Niehaus, T. A. SCC-DFTB Parametrization for Boron and Boranes. J. Chem. Theory Comput. 2012, 8, 1153–1163. (79) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347–1363. (80) McMahon, D. P.; Troisi, A. Evaluation of the External Reorganization Energy of Polyacenes. J. Phys. Chem. Lett. 2010, 1, 941–946.


Page 40 of 54

Page 41 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(81) Nelsen, S. F.; Blackstock, S. C.; Kim, Y. Estimation of Inner Shell Marcus Terms for Amino Nitrogen Compounds by Molecular Orbital Calculations. J. Am. Chem. Soc. 1987, 109, 677–682. (82) 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 C01. Gaussian, Inc.: Wallingford, CT, 2009. (83) Su, W. P.; Schrieffer, J. R.; Heeger, A. J. Solitons in Polyacetylene. Phys. Rev. Lett. 1979, 42, 1698–1701. (84) Jortner, J. Temperature Dependent Activation Energy for Electron Transfer between Biological Molecules. J. Chem. Phys. 1976, 64, 4860–4867. (85) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T.; Mezey, P. G. Numerical Recipes in C, 2nd ed.; Cambridge University Press: Cambridge, United Kingdom, 1992. (86) DeVault, D. Quantum-Mechanical Tunneling in Biological Systems; Cambridge University Press: Cambridge, 1984. (87) Nan, G.; Yang, X.; Wang, L.; Shuai, Z.; Zhao, Y. Nuclear Tunneling Effects of Charge Transport in Rubrene, Tetracene, and Pentacene. Phys. Rev. B 2009, 79, 115203. (88) Marcus, R. A. Electron Transfer Reactions in Chemistry. Theory and Experiment. Rev. Mod. Phys. 1993, 65, 599–610. (89) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A Long-Range-Corrected Time-Dependent Density Functional Theory. J. Chem. Phys. 2004, 120, 8425–8433. (90) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(91) Troisi, A.; Cheung, D. L. Transition from Dynamic to Static Disorder in OneDimensional Organic Semiconductors. J. Chem. Phys. 2009, 131, 014703. (92) Parandekar, P. V.; Tully, J. C. Detailed Balance in Ehrenfest Mixed Quantum-Classical Dynamics. J. Chem. Theory Comput. 2006, 2, 229–235. (93) Kranz, J. J.; Elstner, M. Simulation of Singlet Exciton Diffusion in Bulk Organic Materials. J. Chem. Theory Comput. 2016, 12, 4209–4221. (94) Wang, L.; Beljonne, D. Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-like Transport in Organic Crystals. J. Phys. Chem. Lett. 2013, 4, 1888–1894. (95) Han, L.; Ke, Y.; Zhong, X.; Zhao, Y. Time-Dependent Wavepacket Diffusion Method and Its Applications in Organic Semiconductors. Int. J. Quant. Chem. 2015, 115, 578– 588. (96) Valeev, E. F.; Coropceanu, V.; da Silva Filho, D. A.; Salman, S.; Brédas, J.-L. Effect of Electronic Polarization on Charge-Transport Parameters in Molecular Organic Semiconductors. J. Am. Chem. Soc. 2006, 128, 9882–9886. (97) Nishimoto, Y.; Fedorov, D. G. Three-Body Expansion of the Fragment Molecular Orbital Method Combined with Density-Functional Tight-Binding. J. Comput. Chem. 2017, 38, 406–418.


Page 42 of 54

Page 43 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Simulations settings for the Ehrenfest dynamics (ED) and KMC approaches ED KMC Number of Realization M 1 000 5 000 Number of Sites N ED = 600 N KMC = 100 Time Increment 0.025 fs 1/(kj,j+1 + kj,j−1 )a b ˚ Lattice Length L (A) 3.6 3.6 a See the main body of the text.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Average values and standard deviations of transfer integrals obtained from MD (in meV) (j, j + 1) (2, 3) (1, 2) (3, 4)

T = 300 K ⟨Tj,j+1 ⟩ σ{Tj,j+1 } 38.07 102.0 50.42 96.55 46.86 101.1

T = 200 K ⟨Tj,j+1 ⟩ σ{Tj,j+1 } 76.63 92.19 81.27 90.43 88.91 89.20


Page 44 of 54

Page 45 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Su-Schrieffer-Heeger (SSH) model parameters for the Ehrenfest dynamics approach. Parameter Values meaning −1 τ 307 cm (=38.07 meV) transfer integral m1 6 amu mass of Holstein’s mode ω1 1 400 cm−1 frequency of Holstein’s mode a α 20 543 cm−1 /˚ A Holstein’s local e-p coupling m2 200 amu mass of Peierls’ mode −1 ω2 25 cm frequency of Peierls’ mode βb 2 054 cm−1 /˚ A Peierls’ non-local e-p coupling √ a α√ = ω1 m1 λ with the reorganization energy λ of 0.15 eV b β = σ/ 2kB T /m2 (ω2 )2 with the standard deviation σ of 823 cm−1 .



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 4: Hole mobility through columnar-stacked pyrene moieties in TP-COF obtained form different carrier diffusion simulations at 300 K. Numbers are obtained by fitting a straight line to MSD(t) in Figure 6 Method mobility (cm2 V−1 s−1 ) Monte Carlo simulations QM CT with dynamic disorder 2.067 QM CT with static ⟨T2,3 ⟩ 1.690 Marcus’ CT with dynamic disorder 1.025 Ehrenfest dynamics simulations 1.930


Page 46 of 54

Page 47 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

HO

PDBA

OH

HHTP OH

HO B

B

HO

OH

HO

OH HO

OH

TP-COF BO O

OB O O B O

O B O O BO

O OB

OB O

BO O O B O

O B O O BO

O OB

OB O

BO O O B O

O B O O BO

(b)

O OB

3 × 3 × 4 Super-cell

Figure 1: (a) Chemical Structure of TP-COF and (b) 3 × 3 × 4 supercell used for MD simulations under periodic boundary condition. The central region marked by a dashed open circle is used as the TP-COF model on which we performed the FMO-DFTB calculations



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 54

(a) 8

12

4

7

3

11 T2,3

6

10

2

5

9

1

(b) H

H

H

H H

H

H

H

H

O

H

O B

H

H

H H

H

H

O H

H

H

B

O

H

H

H

H

H

H

H

H

X-axis

Y-axis Z-axis

(c)

I

III O B

O B O

O

II O B O

IV

B O O

O B O

O O

B

O B

O

O

B O

O B O

O O B

Figure 2: (a) TP-COF model extracted from the MD snapshots, its fragmentation, and their numbering for the FMO-DFTB calculations. (b) Top view of the single-layered TPCOF model for the FMO-DFTB calculations. (c) Four molecules for the estimation of the reorganization energy.


Page 49 of 54

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




(a) Transfer Integral |Ti,j | (meV)

700 600

d 500 400 300 200 LC-BLYP sDFTB DFTB 100 B3LYP PBE 0 3.2 3.4

3.6

3.8

4

4.2

Intermolecular Distance, d (Å) (b) 500

Transfer Integral |Ti,j | (meV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 50 of 54

LC-BLYP sDFTB DFTB B3LYP PBE

400

300

200

3.4Å

100

0

0

0.5

1

1.5

2

Lateral Displacement (Å) Figure 3: Transfer integrals from DFTB and DFT calculations for pyrene homodimer cations. (a) Distance dependence for the perfect cofacial dimer and (b) dependence on the lateraldisplacement along the long-axis.


Page 51 of 54

(a) Transfer integral (meV)

500 400 300 200 100 0 -100 -200 -300 0

50

100 150 200 250 300 350 400 450

Time (ps)

(b) 0.09

200K 300K

0.08 0.07

Probability

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.06 0.05 0.04 0.03 0.02 0.01 0

-200

-100

0

100

200

300

400

Transfer integral (meV) Figure 4: (a) Time-dependence of transfer integrals, T2,3 , calculated by the FMODFTB/LCMO method with eq 5 for the MD snapshots at 300 K and (b) normalized probability distributions of T2,3 at 200 K (red open circle) and 300 K (green closed circle). The red and green solid lines represent the fitted Gaussian functions.



Power Spectrum (arb. unit)

(a)

0

20

40

60

Frequency (cm-1)

80

100

(b) Re V (meV)

10 8 6 4 2 0 -2 -4 -6 -8 8 6 4 2 0 -2 -4 -6 -8

Im V (meV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 52 of 54

0

50

100

Frequency

150 (cm-1)

200

Figure 5: (a) Power spectrum density of the time-course transfer integral, T2,3 , calculated by the FMO-DFTB/LCMO method with eq 5 for the MD snapshots at 300 K and (b) Fourier transformation of the thermal deviation amplitude V2,3 at 300 K.


Page 53 of 54

(a) 8000

MSD(t) / L2

7000

ω1=1400 cm-1; α= 20543. cm-1/Å m2=200 m2=2000

6000 5000 4000 3000 2000 1000 0

0

m2=200 ω1=1000 cm-1; α= 14673. cm-1/Å ω1=800 cm-1; α= 11739. cm-1/Å

20

40

60

t (ps)

80

100

(b) 9000

KMC Simulations QM CT QM with static 7000 Marcus’ CT 8000

MSD(t) / L2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6000 5000 4000 3000 2000 Ehrenfest Dynamics Simulations

1000 0

0

20

40

60

80

100

t (ps) Figure 6: (a) Time (t) evolution of the mean squared displacement normalized by the squared equilibrium distance between the π-stacked pyrene units, MSD(t)/L2 , from the mean field approach at 300 K with different m2 and ω1 parameters and (b) Time evaluation of MSD(t)/L2 from Monte Carlo simulations based on the QM CT rate formula, eq 15, with dynamic disKMC (t) (red), the QM CT rate formula with static ⟨T2,3 ⟩ (green), and the Marcus’ order of Tj,j+1 CT rate formula, eq 16, with dynamic disorder (blue), which are compared to the Ehrenfest dynamics result (black; m2 and ω1 are set to 200 amu and 1400 cm−1 , respectively).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry


Page 54 of 54

Multiscale Simulations on Charge Transport in Covalent Organic

Recommend Documents