Computational Study of Excess Electron Mobility in High-Pressure

Apr 18, 2016 - In recent years, excess electron transfer in organic liquids has attracted increasing interest owing to the emerging class of liquid or...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCC

Computational Study of Excess Electron Mobility in High-Pressure Liquid Benzene Masahiro Sato,† Akiko Kumada,*,† Kunihiko Hidaka,† Toshiyuki Hirano,‡ and Fumitoshi Sato‡ †

Department of Electrical Engineering and Information Systems, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan ‡ Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro, Tokyo 153-8505, Japan S Supporting Information *

ABSTRACT: In recent years, excess electron transfer in organic liquids has attracted increasing interest owing to the emerging class of liquid organic semiconductors. In this study, to achieve a comprehensive understanding of electron conduction in liquids, we investigate hopping electron conduction in liquids from an atomistic viewpoint. High-pressure liquid benzene is chosen as a simple model system. Hopping electron mobility is computed using a combination of molecular dynamics simulations, quantum chemical calculations, and kinetic Monte Carlo methods. The computed electron mobility is in good agreement with the experimental values. Because the amplitude of the intermolecular vibration observed in liquids is larger compared to that in solids, the effect of dynamic disorder on electron mobility is investigated. The time scale of the change in electronic couplings due to the rotation of molecules is comparable to that of the electron residence time at each benzene molecule at the absence of change in the arrangement of a benzene dimer. Thus, the effect of dynamic disorder is evaluated by a Monte Carlo-based analysis. It is shown that the effect of dynamic disorder on electron mobility is small even though electronic couplings between molecules fluctuate by more than an order of magnitude.



INTRODUCTION Organic electronics have attracted considerable attention in recent years.1 In particular, liquid organic semiconductors (LOSs) have emerged as a novel class of organic semiconductors.2 Because charge transfer is one of the most fundamental properties for designing LOS devices, many studies have investigated the charge carrier mobility in organic liquids.3−5 At the same time, dielectric liquids have been used for insulating electrical equipment for more than 100 years.6 It is known that the electrical conduction of excess electrons in liquids triggers the electrical breakdown of the liquid. Therefore, electron mobilities in various liquids were measured several decades ago.7−10 Excess electron mobility in nonpolar liquids is classified as follows according to the magnitude of mobility: (I) quasi-free mobility11−13 (μ ≥ 10 cm2/(V·s)), (II) two-state mobility,14−16 and (III) hopping mobility9,10 (μ < 10 cm2/(V·s)). When equilibrium of quasi-free and trapped states of excess electrons M + e− ⇌ M −

described by quasi-free electron mobility and hopping electron mobility, quasi-free electron transfer and hopping electron transfer are of primary interest for understanding electron conduction in liquids. Quasi-free electron mobility is well described by deformation potential theory, which led to the successful treatment of electron mobility in Xe, Ar, Kr, and methane.11 However, hopping transfer of carriers in liquids, that is, carrier transfer via localized sites, has not been studied adequately from an atomistic viewpoint. Revealing the nature of electron hopping transfer enables us to quantitatively describe electron mobilities in liquids in a unified manner. Furthermore, it is thought that charge transport in liquid organic semiconductors is dominated by intermolecular hopping. Therefore, we evaluated hopping electron mobility in high-pressure liquid benzene as a simple model system for organic liquids. It is known that hexafluorobenzene (C6F6) exhibits hopping electron transfer, which is reasonable since electrons are bound to the molecules and form negative ions owing to the large electron affinity between 1 and 2 eV.17 In contrast, electron affinity of benzene ranges from −1.12 ± 0.03 to −1.15 ± 0.03 eV,18 and it is unfavorable to form anions in the gas phase.9 Nevertheless, it is well-known that benzene molecular anion is

(1)

is achieved, electron mobility is described by two-state mobility μtwo = μtrap P + μqf (1 − P), where M and M− represent neutral molecules and anions, respectively; P denotes the probability of localization of electrons; and μtrap and μqf represent the mobility in localized and extended states, namely, quasi-free and hopping mobility, respectively.12 Because two-state mobility is © XXXX American Chemical Society

Received: February 16, 2016 Revised: March 29, 2016

A

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

Article

The Journal of Physical Chemistry C

to the dispersion force.29,30 In addition, because ab initio postHartree−Fock calculations are computationally demanding, at most a few tens of molecules can be handled. Moreover, we are interested in the time evolution of the system when computing the mobility in liquids. The molecular structure of liquid benzene has been studied extensively by Fu et al.,35 and it has been shown that the Optimized Potentials for Liquid Simulations All-Atom (OPLS-AA)33,34 force field reliably describes the configuration of benzene molecules and various thermodynamic properties. Therefore, the structure of liquid benzene was computed using classical MD simulations. All classical MD simulations were performed using the GROMACS program package.31,32 In addition to the OPLS-AA force field, we have adopted the Generalized AMBER Force Field (GAFF)36 for the purpose of comparisons. We have followed the computational procedures used by Fu et al.35 Table 1 shows a brief overview of the equilibrium MD simulation. See ref 35 and the Supporting Information for more details.

stable in the condensed phase owing to the solvent-induced stabilization,19 and some portion of electrons is expected to attach to molecules.9 In fact, it is shown that there is a threshold size that can form stable benzene anion cluster: 53 molecules, depending on the formation method of the cluster anions.19 With regard to high pressure-liquid benzene, since reaction 1 shows large negative volume change, the equilibrium shifts toward M− as the pressure increases. Electron mobility decreases with increasing pressure and levels off at a pressure of around 2000 bar when the equilibrium completely shifts to the M− side.9,10 Therefore, it is thought that electron mobility in liquid benzene under high pressure is described by hopping electron transfer. M− + M → M + M−

(2)

This hopping electron mobility in the saturation region, i.e., lower limit of two-state electron mobility with respect to pressure, is simply referred to as “hopping electron mobility” in this article. We have investigated excess electron hopping mobility in liquid benzene using the Marcus theory, which is extensively used to examine the charge transfer process in solid organic semiconductors,20−22 biosystems,23−25 a catalytic environment,26 and so on. We discuss the quasi-free to hopping transition from an atomistic viewpoint by investigating excess electron localization characteristics in terms of electronic coupling and reorganization energy. From a structural viewpoint, solids and liquids have similar density. However, considering the time evolution of the structure, it is expected that the effect of dynamic disorder plays a more important role in liquids than in solids. The effect of dynamic disorder on hopping electron mobility is studied using Monte Carlo (MC)-based analysis.

Table 1. Conditions for MD Simulation number of molecules boundary condition time step MD duration temperature pressure Coulomb interaction restraints

500 periodic (3D) 1 fs 10 ns (NPT product run) Nosé−Hoover thermostat37,38 Parrinello-Rahman isotropic barostat39 Particle-mesh Ewald (PME) method40 P-LINCS algorithm41,42

Computation of Charge Transfer Rate. Pairs of benzene molecules (dimers) were taken from the liquid benzene structure obtained by the MD simulation described above. The electron transfer rates between these pairs were computed using the Marcus theory. For practical purposes, the hopping rate between molecules for which the distance between centroids is larger than the cutoff radius rc was set to 0, thereby reducing the computational complexity from O(Nmol2) to O(Nmol). Here, Nmol is the number of molecules in a cell. The determination of the appropriate cutoff radius is discussed below. All QM computations were performed using Gaussian 09, revision D.01.43 As shown in Figure 1, the charge transfer between molecules “D” and “A,” which is considered a transition between diabatic states, is understood on the basis of Fermi’s golden rule.44 Assuming the Condon approximation, for the high-temperature limit, the formula is reduced to the well-known Marcus formula:45



COMPUTATIONAL PROCEDURE The computational approach is based on those used to simulate charge carrier dynamics in crystalline and amorphous solid semiconductors.21,27,28 In this study, we combine molecular dynamics (MD) simulations, quantum chemical (QM) calculations, and kinetic Monte Carlo (KMC) simulations to investigate excess electron transport in liquid benzene. The configuration of liquid benzene at various temperatures and pressures is generated by MD simulations. Given the configuration, the electron hopping rates between pairs of benzene molecules are computed by QM calculations using the Marcus theory. The effect of dynamic disorder is studied by computing the time evolution of hopping rates between molecules. Because intermolecular vibrations are slow owing to weak intermolecular forces, quantum effects are negligible.27 Thus, the effect of dynamic disorder is treated classically by computing the hopping rates between pairs of benzene molecules taken from successive snapshots of classical MD simulations. Because the effect of dynamic disorder is moderate (vide inf ra), we use a snapshot of MD simulations for the structure of benzene molecules to compute the mobility. Hopping electron mobility is evaluated by KMC simulations from the configuration of the benzene molecules and the hopping rates between them. Generating Liquid Benzene Structure. The microscopic configuration is of great importance in computing the hopping rates between molecules. An amorphous structure can be obtained either by MD or MC methods based on QM or classical calculations. However, with regard to QM methods, the density functional theory has some drawbacks with respect

2π ⟨|Hif |2 ⟩(FCWD) ℏ ≡ ⟨φi|Ĥ |φf ⟩

kif = Hif

(FCWD) ≡

⎛ (λ + ΔG 0)2 ⎞ 1 exp⎜ − ⎟ 4λkT ⎠ ⎝ 4πλkT

(3)

where kif is the electron transfer rate constant, φi and φf are the diabatic initial and final states, respectively, ΔG0 is the free energy difference between the initial and the final states, λ is the reorganization energy (energy difference owing to relaxation of geometry), and FCWD is the density of states weighted B

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

Article

The Journal of Physical Chemistry C

Figure 1. Schematic diagram of electron hopping between molecules “D” and “A.” The solid and broken lines show the relation between adiabatic states and diabatic states for the reactant and product, respectively. The figure on the right shows the schematic diagram of the four-point method. The value ΔG* ≡ (ΔG0 + λ)2/4λ corresponds to the activation energy. Within the framework of four-point method, ΔG0 = 0.

Franck−Condon factor. It should be noted that ⟨|Hif |2⟩ is the mean square (ensemble-averaged) coupling, which is used instead of |Hif |2 in conformationally flexible systems.46,47 Computation of Electronic Coupling. Kubas et al.48,49 examined the electronic couplings of positively and negatively charged π-conjugated organic dimers with high-level ab initio calculations such as the multireference configuration interaction MRCI+Q level of theory to assess the performance of DFTbased approaches. For the fragment-orbital (FO) method,50,51 they have evaluated the effect of the electron number for constructing the Hamiltonian. It is concluded that although the correct number of electrons is 2N + 1 for electron transfer and 2N − 1 for hole transfer, the method using 2N electrons shows better agreement (by a few percentage points) with the benchmarks. The deviation from the benchmark is at most a few tens of percentage points.48,49 The integer N is the number of electrons in a neutral monomer. Thus, we have adopted the FO method and fragment charge difference (FCD) method52 using 2N electrons. The electronic coupling was computed at PW91PW91/6-31G*53 level of theory, which well describes the electronic couplings.56 As is well-known, benzene has doubly degenerate π* LUMO. Thus, the charge localized states cannot be estimated by a twostate approach,54,55 that is, the electronic coupling of the benzene dimer for electron transfer cannot be computed from two adiabatic states. For such cases, the electronic coupling is obtained as a root-mean-square of the electronic couplings for all degenerate states. The effective electronic coupling Heff IF , which is substituted into the Marcus formula, is obtained by further considering the multiplicity of the final state as49 eff |HIF |=

N2 N1N2

N1

N2

i

f

the LUMO + 1s of each molecule, LUMO−LUMO, LUMO− LUMO + 1, LUMO + 1−LUMO, and LUMO + 1−LUMO + FO FO FO 1, respectively denoted as HFO LL , HLL1, HL1L, and HL1L1, were computed as follows where 1

HIFFO

=

Jif − 2 (ei − ef )Sif 1 − Sif 2

(5)

ei = ⟨φi|Ĥ |φi⟩ Jif = ⟨φi|Ĥ |φf ⟩

(6)

Here, Sif is the overlap between orbitals i and f.50,51 The orbitals i and f belong to molecules “D” and “A,” respectively. With regard to the FCD method, a multistate method that includes several states for a single charge localized site is applied. Rust et al. discussed a multistate approach for the general Mulliken-Hush method,57 and this approach has been shown to be valid for the FCD method.58 We use this method to handle degenerate states by adopting four states that are mainly formed by the interaction between LUMOs and LUMO + 1s of each benzene molecule. Because the benzene monomer has a stiff molecular structure, the quasi-degenerate LUMO and LUMO + 1 (π* states) are clearly distinct from LUMO + 2 (Rydberg-like state); the energy difference between LUMO + 2 and LUMO + 1 was more than 2 eV for all benzene molecules. Because the electronic coupling between benzene molecules is at most around 0.2 eV (see below), LUMO, LUMO + 1, LUMO + 2, and LUMO + 3 are clearly distinct from LUMO + 4 and LUMO + 5 for the same reason. Nevertheless, LUMO and LUMO + 1 of benzene molecules are not completely degenerate owing to slight structural deformation. This enabled us to simply include LUMO, LUMO + 1, LUMO + 2, and LUMO + 3 in the multistate treatment. With a unitary transformation U1 that diagonalizes the donor−acceptor charge difference matrix

∑ ∑ Hif 2 (4)

|Heff IF |

Note that is the effective electronic coupling between molecules m and n, and |Hif | is the electronic coupling between diabatic states i and f, and N1 and N2 are the degrees of degeneracy of the diabatic states for molecules I and F, respectively. With regard to the FO method, because the LUMO and LUMO + 1 of a single molecule (monomer) is degenerate, the electronic coupling between the LUMOs and C

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

Article

The Journal of Physical Chemistry C ⎛ ΔQ 11 ⎜ ⎜ ΔQ 21 ΔQ = ⎜ ⎜ ΔQ 31 ⎜ ⎜ ΔQ ⎝ 41

ΔQ 12 ΔQ 13 ΔQ 14 ⎞ ⎟ ΔQ 22 ΔQ 23 ΔQ 24 ⎟ ⎟ ΔQ 32 ΔQ 33 ΔQ 34 ⎟ ⎟ ΔQ 42 ΔQ 43 ΔQ 44 ⎟⎠

ΔG 0 = ΔGD0 + ΔGA0 = [ED(Q D) − E D−(Q D−)] + [E A−(Q A−) − EA (Q A)] λ = λD + λA = [ED(Q D−) − ED(Q D)] + [E A−(Q A) − E A−(Q A−)]

(7)

(13)

where ED(QD) and ED (QD) are the ground-state energy of the neutral and negatively charged state, respectively, of the “D” molecule at its optimized geometry in the neutral state. Similarly, ED(QD−) and ED−(QD−) are the energy of the neutral and negatively charged state, respectively, at its optimal geometry in the negatively charged state.59,60 Therefore, within the framework of four-point method, ΔG0 = 0. Estimation of Hopping Electron Mobility. KMC simulations were conducted using the coordinates of 500 benzene molecules computed with MD simulations and hopping rates computed by QM calculations. The KMC procedure is based on Bortz’s method.61 Hereafter, kmn denotes ) the transition rate of electrons from site m to n. Let Ω(m m ′ ≡ ∑mn ′ kmn and Ω(0) = 0 be the cumulative transition probability, m where m′ = 1, 2, ..., Nm are the hopping sites accessible from site m. A brief overview of Bortz’s KMC procedure is given below: (I) generate random number u1, u2 ∈ (0,1); (II) if

where ΔQ ij =



∫r∈D ρij (r) dr − ∫r∈A ρij (r) dr

(8)

the corresponding Hamiltonian is accordingly transformed as ⎛ E1 ⎞ ⎜ ⎟ E ⎜ ⎟ 2 U1T⎜ ⎟U1 = E3 ⎜ ⎟ ⎜ ⎟ E4 ⎠ ⎝

⎛ H11 ⎜ ⎜ H21 ⎜ ⎜ H31 ⎜ ⎝ H41

H12 H13 H14 ⎞ ⎟ H22 H23 H24 ⎟ ⎟ H32 H33 H34 ⎟ ⎟ H42 H43 H44 ⎠

(9)

where T denotes the transpose of the matrix, and Ei are the eigenvalues of the adiabatic states. The quantity ρij is a one-particle density (i = j) or a transition density between |i⟩ and |j⟩ (i ≠ j).58 The diagonalized donor− acceptor charge difference matrix UT1 ΔQU1 has values close to 1 and −1, which corresponds to the charge-transfer (CT) states D− + A (CT1) and D + A−(CT2), respectively. Thus, the CT states are divided into two subspaces to block-diagonalize each subspace to remove the ambiguity due to near-degeneracy;58 the right-hand side of eq 9 is transformed into

Ω(mm ′− 1)/Ω(mNm) < u1 ≤ Ω(mm ′)/Ω(mNm)

holds, we chose path m′ as the next hopping site and increment the time by Δt = −ln(u2)/ΩNmm.61 The time t and hopping site m is updated until t reaches a certain time. The trajectory of electron hopping transfer as a function of each time step r(t), which is obtained through this procedure, is used to compute the diffusion coefficient D from the Einstein relation:

⎛ HCT1 HCT1,CT2 ⎞ T ⎜ ⎟U Ublock ⎜H HCT2 ⎟⎠ block ⎝ CT1,CT2 ⎛ E H̃ CT1,CT2 ⎞ CT1 ⎟ = ⎜⎜ ⎟ ̃ H E CT2 ⎠ ⎝ CT1,CT2

D = lim

t →∞

(10)

(15)



RESULTS AND DISCUSSION Structure of Liquid Benzene. The electron mobility in liquid benzene is measured under pressures of 1000−2500 bar and temperatures of 300−400 K. Thus, we have parametrized the pressure and temperature dependence of the density of liquid benzene using a Tait-like equation

(11)

The matrix U is a unitary transformation that diagonalizes the CT subspace. Thus, the diabatic states are estimated from the adiabatic states as ⎛ φ1 ⎞ ⎛ ψ1 ⎞ ⎜ ⎟ ⎜ ⎟ ψ ⎜ φ2 ⎟ T⎜ 2⎟ ⎜ φ ⎟ = (U1Ublock ) ⎜ ψ ⎟ ⎜⎜ 3 ⎟⎟ ⎜⎜ 3 ⎟⎟ ⎝ φ4 ⎠ ⎝ ψ4 ⎠

⟨|r(t ) − r(0)|2 ⟩ 2dt

where the brackets ⟨ ⟩ denote the ensemble average, and d is the dimension of the system. The electron mobility μ is computed using the Einstein relation μ = eD/kBT, where e is the elementary charge. The three-dimensional periodic boundary condition used for MD simulations was applied for KMC simulations. Random numbers were generated using Matlab 2014b.62

where E is a diagonal matrix, HCT1 is a subspace of the CT1 state, and

⎛UCT1 0 ⎞ ⎟ Ublock = ⎜ ⎝ 0 UCT2 ⎠

(14)

−1 ⎛ P + B (T ) ⎞ ρ(P , T ) = A(T )⎜1 − C ln ⎟ B (T ) ⎠ ⎝

(16)

to compare the experimental values and those computed from MD simulations. A(T), B(T), and C are fitting parameters. Equation 16 is slightly modified from the Tait equation used by Gibson and Kincaid63 to cover the range up to 3000 bar and 470 K. The details are described in the Supporting Information. Table 2 shows the relation between experimental values,63−65 including those interpolated using eq 16, and the computed values of the density of liquid benzene. As shown in Table 2, the OPLS-AA values and GAFF values shows good agreement

(12)

Computation of Franck−Condon Weighted Density of States. The reorganization energies and free energy differences were calculated using the four-point method.59,60 As shown in Figure 1, these values can be calculated directly from the adiabatic potential energy surface. The total reorganization energy and free energy difference are written as D

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

Article

The Journal of Physical Chemistry C

and temperature range, which is in agreement with the comparable density for all conditions. It should be noted that for 1 atm and 298 K, RDF computed by OPLS-AA and GAFF shows close agreement with the experimental values.35 Electron Hopping Rates between Pairs of Benzene Molecules. Electronic Coupling between Benzene Molecules. Figure 3 shows the adiabatic and diabatic states for the benzene dimer. The diabatic states are computed from the adiabatic states using eq 12. The dimer is taken from a snapshot of the MD simulation at 1500 bar and 355 K. The charge localized states are successfully estimated by including four states. As shown in Figure 4, the electronic couplings computed by the FO method show excellent agreement with those computed by the FCD method. This strongly indicates that the weak coupling approximation holds well for this system. As is wellknown, electronic coupling is strongly correlated with the overlap between interacting orbitals.66 The relation between effective electronic coupling |Heff IF | and effective overlap

Table 2. Computed Densities of Liquid Benzene at Various Pressures and Temperatures (Unit: 103 kg/m3)a P (bar)b

T (K)

exptlc

OPLS-AAd

GAFFd

1 1200 1500 1500 2000 2500

300 315 335 355 375 395

0.872 0.932 0.932 0.919 0.929 0.939

0.862 0.932 0.928 0.914 0.923 0.925

0.848 0.924 0.923 0.903 0.919 0.924

a

Differences between computed and experimental values are 2% and 3% for OPLS-AA and GAFF force fields, respectively. bIn order to facilitate comparison with experimental data, bar is used instead of SI unit. cEstimated from experimental values. See text for details. d Average over 10000 steps: 10 ps.

with the experimental values over a wide range of pressures and temperatures. The structures obtained by both OPLS-AA and GAFF force fields are used in the following sections. Figure 2

eff |SIF |

=

N2 N1N2

N1

N2

i

f

∑ ∑ Sif 2 (17)

is shown in the inset of Figure 4. As seen in the histogram in the right, the maximum value of electronic coupling is around 100−200 meV, and the mode is around 50 meV. Figure 5 shows the relation between the electronic coupling and the distance between benzene molecules. As expected from the insensitivity of RDF to temperature, temperature has negligible effect on the electronic coupling distribution (inset of Figure 5a). The electronic coupling decays roughly exponentially with the distance between the centroid of benzene molecules d. However, it is notable that even for constant distance, the electronic couplings vary more than 10-fold owing to the angular degrees of freedom. As shown in Figure 5b, this is attributable to two reasons. First, although the distance between the centroids of molecules “D” and “A” is constrained, the distance of closest approach between the atoms of molecule “D” and those of molecule “A” dmin can change. Because the spatial correlation of the benzene dimer decays with distance, the variation of dmin and therefore that of the electronic coupling becomes larger when d increases. Second, the overlap integral of the interacting orbital is strongly dependent on the phase, nodal plane of benzene monomer orbitals, and geometrical arrangement of the dimer. As shown in Figure

Figure 2. Radial distribution function for the distance between centroid of benzene molecules under various temperatures and pressures.

shows the corresponding radial distribution function (RDF). There is little difference between the OPLS-AA and the GAFF values. RDF does not show a significant change in this pressure

Figure 3. Adiabatic states (LUMO, LUMO + 1, LUMO + 2, and LUMO + 3) of benzene dimer and diabatic states estimated by FCD method and block diagonalization with four states. Isovalue of the contour is 0.02. MOs were visualized using Gauss View program.67 E

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

Article

The Journal of Physical Chemistry C

Figure 4. Relation between electronic couplings |Heff IF | computed by the FCD and FO methods. The histogram on the right shows the frequency distribution of the electronic coupling values. Note that for |Heff IF | < 20 meV the histogram is affected by the cutoff radius: rc = 0.8 nm; see also Figure eff 5. The inset shows the correlation between the electronic couplings computed by the FO method |Heff IF | and the effective overlap |SIF |.

concluded that with an augmented basis, the 2E2u state is no longer a ground state. In addition, they found that a satisfactory description of the molecular structure of the benzene anion can be obtained using the CCSD/cc-pVTZ level.71 Thus, to explore electron hopping transfer via localized states, we computed the reorganization energy for the 2E2u state at the CCSD/ccpVTZ72,73 level in benzene solvent. We used the integral equation formalism polarizable continuum model (IEFPCM)74,75 to account for solvation, thereby including the outer-sphere reorganization term. The computed reorganization energy was λ = 468 meV, which is comparable to the value calculated at MP2/6-311G*.68 The extent of delocalization of electrons can be assessed by comparing the values of electronic coupling and reorganization energy: (I) complete localization (|Hif | = 0), (II) weakly coupled and localized (2|Hif | ≪ λ), and (III) strongly coupled and delocalized (2|Hif | > λ);76 class (II) corresponds to the carrier hopping regime. The reorganization energy of the Rydberg-like state is several orders of magnitude smaller than the electronic couplings for the 2E2u states, which are around 100 meV. This will presumably become more prominent when compared with the electronic couplings of the Rydberg-like states, although we could not determine the electronic coupling and reorganization energy of the Rydberg-like state. This is because the additional diffuse function results in probing continuum-like states71 (see also Supporting Information). However, the reorganization energy for the 2E2u state is larger compared to the electronic couplings for the 2E2u states. The reorganization energies of oligoacenes are comparable (e.g., | Hif | = 110 and λ = 190 meV in the case of naphthalene56) or even smaller than the electronic couplings (e.g., |Hif | = 130 and λ = 99 meV in the case of pentacene56). However, it is assumed

5b, the electronic couplings vary by an order of magnitude even for comparable d and dmin values owing to the phase alignment of interacting orbitals. The inset of Figure 5b shows the typical arrangement of the benzene dimer with large (dimer no. 1) and small (dimer no. 2) electronic coupling. Reorganization Energy of Benzene Molecule. The electronic coupling is computed by the DFT method because a detailed comparison between the DFT method and high-level ab initio calculations has already been conducted by Kubas et al.48,49 For the computation of the reorganization energy, because the electronic states of the benzene anion are sensitive to the theory levels, we used a theory level that has been analyzed in detail (see below). The reorganization energies of benzene and naphthalene were computed by Klimkans ̅ and Larsson68 and Datta et al.69 The energies for naphthalene were found to be 185 and 234 meV at MP2/6-31G*70 and B3LYP/ 6-31G++(d,p) levels, respectively, which is in reasonable agreement, whereas those for benzene were found to be 400 and 3.9 meV at MP2/6-311G* and B3LYP/6-31G++(d,p) levels, respectively, which differ by 2 orders of magnitude. This is because the electronic state that Datta et al. obtained is a Rydberg-like state, whereas that obtained by Klimkans ̅ is the 2 E2u state in which excess electrons are in π*-type orbitals such as 2B1u and 2Au (see Supporting Information). In ref.,71 the Rydberg-like state is referred to as a discretized free electron state coupled to the states of the neutral molecule. Because excess electrons in the Rydberg-like state are weakly bound to the molecule, the extraction of an electron has little effect on the geometry of the molecule. Therefore, the reorganization energy computed for the Rydberg-like state is drastically underestimated compared to that for the 2E2u state. Bazante et al. investigated the nature of the benzene anion in detail and F

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

Article

The Journal of Physical Chemistry C

at a certain site is longer than or comparable to the time-scale of the fluctuation of electronic coupling, the framework using the thermally averaged coupling ⟨|HIF|2⟩ to account for the thermal fluctuation of electronic couplings is no longer effective. Figure 6 shows the distribution of the residence time of electrons for a single snapshot. The “static” residence time at

Figure 6. Residence time distribution of excess electrons at each molecule. The residence time is computed for the OPLS-AA structure by the FCD method. N

molecule m is defined as Δtstm = 1/Ωmm. This is equivalent to the time step of KMC simulations except for the −ln(u2) factor, which is set to unity. It should be noted that the average of −ln(ui) is unity. Within the framework of Marcus theory, Δtst is proportional to the inverse of the FCWD term. Because the FCWD term increases with temperature, Δtst and the dispersion of the static residence time probability decrease as the temperature increases. The time-scale of the static residence time is of the order of 10−100 fs. To evaluate the time-scale of molecular interchange and fluctuation of electronic couplings, we computed the time evolution of the distance between centroids and the transfer rates between benzene molecules for a series of snapshots. We focused on a single molecule m and molecules adjacent to it: all molecules whose distance of closest approach to molecule m during 10 ps is within 0.8 nm. The trajectories of benzene molecules for 10 ps, which is sufficiently longer than the static residence time, were extracted from an additional NPT MD simulation. Therefore, conformational changes that occur on longer time scales can be considered static.81 Snapshots were taken every 100 fs. As shown in Figure 7a, frequent molecular interchange occurs between nearest neighbors, and therefore, the corresponding electronic couplings fluctuate by more than an order of magnitude (Figure 7b). As seen in Figure 7b, the electronic coupling fluctuates by an order of magnitude along with changes in distance between molecules (d) on a time scale of 1 ps, which is of the order of the typical intermolecular vibrational frequency of benzene82 and other nonpolar molecules.27 In addition to this time scale, it is seen that the electronic coupling changes several-fold within a few hundred femtoseconds, which is comparable to the static residence time of electrons. This fluctuation occurs without the occurrence of change in d, presumably owing to changes in the roll and twist angles between benzene molecules, as discussed above. Because

Figure 5. Electronic couplings computed by the FCD method are plotted versus (a) distance between centroid of benzene molecules and (b) distance of closest approach between carbon atoms in each molecule. The electronic couplings are computed for the liquid benzene structure computed using the OPLS-AA force field. The cutoff radius is 0.9 nm. The inset in (a) shows the frequency distribution of the electronic coupling for dimers with a distance of 0.54−0.56 nm. The inset in (b) shows the LUMO of two dimers; the monomer LUMOs of dimer no. 1 are in phase with each other.

from both the experimental results77 and the computational analysis27,56 that hole transfer in oligoacenes shows hopping conduction. Because the reorganization energy and electronic coupling of benzene for electron transfer are comparable to those of oligoacenes for hole transfer, the electron transfer in benzene via the 2E2u state occurs in the hopping regime. Therefore, it is reasonable to think that the transition from quasi-free electron conduction to hopping electron conduction is caused by the destabilization of the Rydberg-like state owing to the increase in pressure. Effect of Dynamic Disorder. We evaluated the effect of nonlocal electron−phonon coupling (dynamic disorder): effect of variation of electronic couplings due to modulations of relative coordination of molecules.79,80 In solids, the pairs of nearest neighbors are fixed; therefore, the contribution of dynamic disorder to electron mobility is directly evaluated by computing the time evolution of the electronic couplings between certain pairs of molecules.27 However, this approach is not valid for liquids owing to molecular interchanges; specifically, the time course of change in neighboring molecules must be taken into account. If the residence time of an electron G

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

Article

The Journal of Physical Chemistry C

Figure 8. Breakdown of Ω̂ at (i), (ii), and (iii) in Figure 7. The bracket denotes the ensemble average, that is, Ω̂ = ∑Nsites k/⟨Ω⟩. Hopping rates whose contribution is less than 1% are not shown in this figure. Number of sites whose contribution is above 1% roughly coincides with the first coordination number of liquid benzene, namely, 12 molecules.

Figure 8, the large value of Ω̂ at (iii) is mostly due to the large hopping rate between site 1. Similarly, Ω̂ at (i) is small because the hopping rate between site 1 is small. Although d between site 1 at (iii) is smaller than that at (i), the hopping rate at (iii) is smaller than (i) by an order of magnitude. This confirms the fact that changes in the molecular orbital alignment are responsible for the fluctuation of the electronic coupling on a time scale of 100 fs. Nevertheless, it is clearly seen that the fluctuation of ΩN(t) is far more moderate compared to the change in the square of electronic coupling. This is counterintuitive to some extent, however. As shown in Figure 7b, the maximum value of the electronic coupling is rather stable compared to the individual electronic couplings; the fluctuation ranges from a few to a few tens of percentage points. Therefore, at least one readily accessible path exists, and this path contributes to the insensitivity of ΩN(t) to the thermal motion of molecules. The results show that the consideration of the time evolution of the hopping rates has little effect on the electron residence time. Further, to include the effect of change in hopping distance for the corresponding hopping rates, we introduce a quantity

Figure 7. Thermal fluctuation of (a) distance between centroids of benzene molecules d, (b) electronic couplings, and (c) D̂ value defined in the text. Snapshots of the geometry are captured every 100 fs from the NPT simulation with the OPLS-AA force field at 395 K and 2500 bar. Molecules whose distance of closest approach during 10 ps (dmin) is smaller than 0.5 nm are indicated by solid lines. The colors of the solid line in (a) correspond to those in (b).

the electronic coupling has a large effect on the hopping rates (k ∝ |H|2) and the electron hopping and fluctuation of electronic couplings occur on comparable time scales, it may be necessary to apply a time-dependent algorithm to evaluate the electron residence time with reasonable accuracy. According to the dynamical MC algorithms proposed by Prados et al.,83 the problem that the transition rates are timedependent is solved by simply modifying the time step used in the KMC simulations from Δt = −ln(u2)/ΩNmm to Δt that holds −ln(u 2) =

∫0

Δt

dt ′ΩmNmt ′

Dm(t ) ≡

(∑n dmn(t )kmn(t )/ΩmNm(t ))2 1/ΩmNm(t )

(19)

whose unit is identical to that of the diffusion constant, which is proportional to the mobility. The denominator represents the characteristic time scale, that is, time between hopping computed for a certain snapshot, and the numerator represents the square of the characteristic length, which is defined as an expected value of the hopping distance. Figure 7c shows the time evolution of D̂ ≡ D/⟨D⟩. The brackets ⟨ ⟩ denote the average of snapshots. As shown in Figure 7c, the time evolution of D̂ resembles that of ΩN. This is not surprising because the fluctuation of the characteristic hopping distance is moderate compared to that of electronic couplings. The standard deviation of D̂ (σD) is 0.30. Thus, the electron mobilities are evaluated reasonably, at least with respect to the order of magnitude, even when dynamic disorder is neglected.

(18)

Thus, to evaluate the effect of time-dependent hopping rates on N the residential time, we computed Ωmm(t), as shown in Figure Nm 7c. It should be noted that Ωm (t) in the figure is normalized by the mean value for comparison. Because we focus on the time evolution of the hopping rates and hopping distances between a single molecule m and molecules adjacent to it, the suffix m is abbreviated. The breakdown of Ω̂ at (i), (ii), and (iii), represented in Figure 7, is shown in Figure 8. As seen from H

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

Article

The Journal of Physical Chemistry C Hopping Electron Mobility. We carried out KMC simulations for snapshots taken from MD simulations to compute the hopping electron mobility. Figure 9 shows the

mobility increases owing to the increase in the number of accessible hopping sites, while rc is smaller than 7.5 Å. The electron mobility is adequately described when rc is larger than rmin = 7.5 Å. As shown in Figure 2, RDF takes its first minimum c at rm1 = 7.5 Å, which roughly coincides with rmin c . This is reasonable because the integral to the first minimum of the RDF over the radius is the first coordination number, that is, taking the radius that gives the first minimum of RDF as rc resembles the nearest neighbor interaction approach. Figure 11 shows the experimental9 and computed hopping electron mobility in liquid benzene under various temperatures

Figure 9. Mean square displacement and its xyz component versus time. The value is averaged over 5000 trajectories. The liquid benzene structure and hopping rates are computed at 1200 bar and 315 K. The inset shows the typical trajectory of excess electrons. The simulation time for the graphical representation of the electron trajectory is 500 ps.

relation between the mean square displacement and time t that are computed from the trajectory of electrons. We ignored the effect of dynamic disorder, as it was shown that it has little effect on the mobility. Because the time between hopping events was roughly of the order of 10−100 fs, the single KMC simulation time was set to 10 ps so that a sufficient number of hopping events could be carried out. We confirmed that extending the simulation time did not affect the computed mobilities. The isotropy of the simulated mean square displacement indicates that the simulation cell size is sufficiently large. To investigate the validity of the cutoff radius rc, the hopping electron mobility was computed by varying rc. The statistical error of the simulated mobility is evaluated by the method proposed by Nan et al.78 Figure 10 shows that the electron

Figure 11. Experimental9 and computed values of hopping electron mobility in high-pressure liquid benzene. Error bars represent the statistical error. It should be noted that the statistical error for KMC simulations is evaluated by the method proposed by Nan et al. and that the error bar for snapshot dependence shows the standard deviation for mobilities computed from 10 snapshots. The effect of the force fields used in MD simulations (see legend) is negligible.

and pressures. The electron mobilities computed from snapshots of the liquid benzene structure predicted using the OPLS-AA and GAFF force fields show little difference. The differences are comparable to the statistical error of the KMC simulations, which is at most 10% of the mobility. The mobility roughly coincides with those estimated from μest = ∑NmmolDm/ Nmol and the Einstein relation, which is around 0.3 cm2/(V s) for 355 K. The snapshot dependence of electron mobility was evaluated by computing the mobility using 10 snapshots extracted from MD simulation runs at 355 K and 1500 bar. The mean value over 10 snapshots μ̅s and corresponding standard deviation σμs were 0.176 and 0.005 cm2/(V s), respectively. The snapshot dependence is comparable to the statistical error. Thus, although D̂ shows variations owing to local structural fluctuations (Figure 7), the simulation cell size is sufficiently large to capture macroscopic electron transfer properties such as the bulk mobility. The electron mobility computed by the KMC simulation is in reasonable agreement with the experimental values regardless of the temperature and pressure. Nevertheless, there is a 30−40% deviation between experimental and computed values. The main sources of error are probably the use of (1) the high-temperature approximation, (2) approximation employed by four-point method, and (3) approximations for predicting electronic couplings. Thus, in order to increase the accuracy of the computation, one may (1) go beyond Marcus theory and treat vibrational modes quantum mechanically,27 (2) go beyond four-point method by evaluating

Figure 10. Computed electron mobility as a function of cutoff radius rc. The broken line shows the interval average from rc = 0.75 to 0.9 nm. I

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

Article

The Journal of Physical Chemistry C the effect of variation of the local electric field,28 and (3) compute electronic couplings at a more accurate level of theory.

fellowship. We would like to thank Editage (www.editage.jp) for English language editing.





CONCLUSIONS We have investigated excess electron conduction in highpressure liquid benzene as a simple model system of liquid organic semiconductor and organic dielectric liquids. The electronic coupling and reorganization energy for electron transfer via the 2E2u state of the benzene molecule were comparable to those for hole transfer in oligoacenes, which supports the experimental assumption that hopping electron transfer occurs in liquid benzene. The electron hopping mobility was evaluated by employing a multiscale approach with a combination of MD simulations, QM calculation, and KMC simulations. Although the model adopted in this study is simplified to a certain extent, the computed electron mobility was in good agreement with the experimental values over a wide range of temperatures and pressures. To explore the charge transfer property specific to liquids, namely, the effect of large intermolecular vibration, we investigated the effect of the dynamic disorder on electron mobility. The effect of the dynamic disorder is studied by treating the time evolution of accessible hopping sites, hopping distances, and corresponding hopping rates. The transfer integral fluctuates by an order of magnitude along with changes in distance between molecules on the time scale of 1 ps, which is of the order of the typical intermolecular vibrational frequency of nonpolar molecules both in liquid and in amorphous and crystalline solid phases. In addition to this time scale, it is seen that the electronic coupling changes several times within a few hundred femtoseconds, which is comparable to the static residence time of electrons. It was shown that these fluctuations were due to the rotational motion of benzene molecules. Because electron hopping and fluctuation of electronic couplings occur on comparable time scales, we evaluated the effect of dynamic disorder on the electron mobility by a MC-based analysis. We demonstrated that the effect of the thermal motion of molecules is moderate even though frequent molecular interchange occurs between nearest neighbors and that the electronic coupling fluctuates by more than an order of magnitude.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b01581. Force field parameters for MD simulations, regression curve of density of experimental liquid benzene, basis set dependence of electronic states, and reorganization energy of benzene anion (PDF)



REFERENCES

(1) Mei, J.; Diao, Y.; Appleton, A. L.; Fang, L.; Bao, Z. Integrated Materials Design of Organic Semiconductors for Field-Effect Transistors. J. Am. Chem. Soc. 2013, 135, 6724−6746. (2) Babu, S. S. Paradigms shift when solvent-less fluids come into play. Phys. Chem. Chem. Phys. 2015, 17, 3950−3953. (3) Plint, T. G.; Kamino, B. A.; Bender, T. P. Charge Carrier Mobility of Siliconized Liquid Triarylamine Organic Semiconductors by Timeof-Flight Spectroscopy. J. Phys. Chem. C 2015, 119, 1676−1682. (4) Kamino, B. A.; Bender, T. P.; Klenkler, R. A. Hole Mobility of a Liquid Organic Semiconductor. J. Phys. Chem. Lett. 2012, 3, 1002− 1006. (5) Ribierre, J.-C.; Zhao, L.; Inoue, M.; Schwartz, P.-O.; Kim, J.-H.; Yoshida, K.; Sandanayaka, A. S.; Nakanotani, H.; Mager, L.; Mery, S.; Adachi, C. Low threshold amplified spontaneous emission and ambipolar charge transport in non-volatile liquid fluorene derivatives. Chem. Commun. 2016, 52, 3103−3106. (6) Schmidt, W. F. Electronic Conduction Processes in Dielectric Liquids. IEEE Trans. Electr. Insul. 1984, EI-19, 389−418. (7) Baxendale, J. H.; Bell, C.; Wardman, P. Pulse Radiolysis Observations of Solvated Electrons in Liquid Hydrocarbons. Chem. Phys. Lett. 1971, 12, 347−348. (8) Munoz, R. C.; Holroyd, R. A. The Effect of Temperature and Pressure on Excess Electron Mobility in n-Hexane, 2,2,4-Trimethylpentane, and Tetramethylsilane. J. Chem. Phys. 1986, 84, 5810−5815. (9) Itoh, K.; Holroyd, R. Effect of Pressure on the Electron Mobility in Liquid Benzene and Toluene. J. Phys. Chem. 1990, 94, 8850−8854. (10) Itoh, K.; Ninomiya, S.; Nishikawa, M.; Holroyd, R. Charge Transport at High Pressures in Toluene-n-Pentane Mixtures. J. Phys. Chem. 1994, 98, 1974−1977. (11) Basak, S.; Cohen, M. H. Deformation-Potential Theory for the Mobility of Excess Electrons in Liquid Argon. Phys. Rev. B: Condens. Matter Mater. Phys. 1979, 20, 3404−3414. (12) Holroyd, R. A.; Schmidt, W. F. Transport of Electrons in Nonpolar Fluids. Annu. Rev. Phys. Chem. 1989, 40, 439−468. (13) Chandler, D.; Leung, K. Excess Electrons in Liquids: Geometrical Perspectives. Annu. Rev. Phys. Chem. 1994, 45, 557−591. (14) Huang, S. S.-S.; Freeman, G. R. Electron Transport in Gaseous, Critical, and Liquid Benzene and Toluene. J. Chem. Phys. 1980, 72, 2849−2855. (15) Knoesel, E.; Bonn, M.; Shean, J.; Wang, F.; Heinz, T. F. Conductivity of Solvated Electrons in Hexane Investigated with Terahertz Time-Domain spectroscopy. J. Chem. Phys. 2004, 121, 394− 404. (16) Shkrob, I. A.; Sauer, M. C., Jr. Photo-Stimulated Electron Detrapping and the Two-State Model for Electron Transport in Nonpolar Liquids. J. Chem. Phys. 2005, 122, 134503. (17) Nylkos, L.; van den Ende, C. A. M.; Warman, J. M.; Hummel, A. High Mobility Excess Electrons in the Electron-Attaching Liquid Hexafluorobenzene. J. Phys. Chem. 1980, 84, 1154−1155. (18) Hajgato, B.; Deleuze, M. S.; Tozer, D. J.; De Proft, F. A Benchmark Theoretical Study of the Electron Affinities of Benzene and Linear Acenes. J. Chem. Phys. 2008, 129, 84308. (19) Mitsui, M.; Nakajima, A. Formation of Large Molecular Cluster Anions and Elucidation of Their Electronic Structures. Bull. Chem. Soc. Jpn. 2007, 80, 1058−1074. (20) Friederich, P.; Symalla, F.; Meded, V.; Neuman, T.; Wenzel, W. Ab initio Treatment of Disorder Effects in Amorphous Organic Materials: Toward Parameter Free Materials Simulation. J. Chem. Theory Comput. 2014, 10, 3720−3725. (21) Yavuz, I.; Martin, B. N.; Park, J.; Houk, K. N. Theoretical Study of the Molecular Ordering, Paracrystallinity, And Charge Mobilities of Oligomers in Different Crystalline Phases. J. Am. Chem. Soc. 2015, 137, 2856−2866. (22) Noriega, R.; Salleo, A.; Spakowitz, A. J. Chain Conformations Dictate Multiscale Charge Transport Phenomena in Disordered

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +81 3 5841 6759. Fax: +81 3 5841 6725. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was partly supported by JSPS Grants-in-Aid for Scientific Research Grant Number 15J05161 and IEEE DEIS J

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

Article

The Journal of Physical Chemistry C Semiconducting Polymers. Proc. Natl. Acad. Sci. U. S. A. 2013, 8, 16315. (23) Heck, A.; Woiczikowski, P. B.; Kubař, T.; Giese, B.; Elstner, M.; Steinbrecher, T. B. Charge Transfer in Model Peptides: Obtaining Marcus Parameters from Molecular Simulation. J. Phys. Chem. B 2012, 116, 2284−2293. (24) Heck, A.; Woiczikowski, P. B.; Kubavr, T.; Welke, K.; Niehaus, T.; Giese, B.; Skourtis, S.; Elstner, M.; Steinbrecher, T. B. Fragment Orbital Based Description of Charge Transfer in Peptides Including Backbone Orbitals. J. Phys. Chem. B 2014, 118, 4261−4272. (25) Genereux, J. C.; Barton, J. K. Mechanisms for DNA Charge Transport. Chem. Rev. 2010, 110, 1642−1662. (26) Yoon, Y.; Wang, Y.-G.; Rousseau, R.; Glezakou, V.-A. Impact of Nonadiabatic Charge Transfer on the Rate of Redox Chemistry of Carbon Oxides on Rutile TiO2(110) Surface. ACS Catal. 2015, 5, 1764−1771. (27) Wang, L.; Li, Q.; Shuai, Z.; Chen, L.; Shi, Q. Multiscale Study of Charge Mobility of Organic Semiconductor with Dynamic Disorders. Phys. Chem. Chem. Phys. 2010, 12, 3309−3314. (28) Rühle, V.; Lukyanov, A.; May, F.; Schrader, M.; Vehoff, T.; Kirkpatrick, J.; Baumier, B.; Andrienko, D. Microscopic Simulations of Charge Transport in Disordered Organic Semiconductors. J. Chem. Theory Comput. 2011, 7, 3335−3345. (29) Becke, A. D. Perspective: Fifty Years of Density-Functional Theory in Chemical Physics. J. Chem. Phys. 2014, 140, 18A301. (30) Klimes, J.; Michaelides, A. Perspective: Advances and Challenges in Treating Van Der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901. (31) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible and free. J. Comput. Chem. 2005, 26, 1701−1718. (32) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (33) Jorgensen, W. L.; Tirado-Reves, J. The OPLS Force Field for Proteins. Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (34) Jorgensen, W. L.; Severance, D. L. Aromatic-Aromatic Interactions: Free Energy Profiles for the Benzene Dimer in Water, Chloroform, and Liquid Benzene. J. Am. Chem. Soc. 1990, 112, 4768− 4774. (35) Fu, C.-F.; Tian, S. X. A Comparative Study for Molecular Dynamics Simulations of Liquid Benzene. J. Chem. Theory Comput. 2011, 7, 2240−2252. (36) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general AMBER force field. J. Comput. Chem. 2004, 25, 1157−1174. (37) Nosé, S. Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (38) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (39) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (40) Essman, U.; Perela, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (41) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (42) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (43) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al.Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2010.

(44) Dirac, P. A. M. The Quantum Theory of the Emission and Absorption of Radiation. Proc. R. Soc. London, Ser. A 1927, 114, 243− 265. (45) Marcus, R. A. Electron Transfer Reactions in Chemistry. Theory and Experiment. Rev. Mod. Phys. 1993, 65, 599−610. (46) Prytkova, T. R.; Krnikov, I. V.; Beratan, D. N. Coupling Coherence Distinguishes Structure Sensitivity in Protein Electron Transfer. Science 2007, 315, 622−625. (47) Fornari, R. P.; Aragó, J.; Troisi, A. A Very General Rate Expression for Charge Hopping in Semiconducting Polymers. J. Chem. Phys. 2015, 142, 184105. (48) Kubas, A.; Fruzsina, G.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic Couplings for Molecular Charge Transfer: Benchmarking CDFT, FODFT, and FODFTB Against High-Level ab initioCalculations. J. Chem. Phys. 2014, 140, 104105. (49) Kubas, A.; Gajdos, 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. II. Phys. Chem. Chem. Phys. 2015, 17, 14342−54. (50) Valeev, E. F.; Coropceanu, V.; Filho, D. A.; da, S.; 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. (51) Wolter, M.; Elstner, M.; Kubar̆, T. Charge transport in Desolvated DNA. J. Chem. Phys. 2013, 139, 125102. (52) Voityuk, A. A.; Rösch, N. Fragment Charge Difference Method for Estimating Donor-Acceptor Electronic Coupling: Application to DNA π-stacks. J. Chem. Phys. 2002, 117, 5607−5616. (53) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 6671−6687. (54) Cave, R. J.; Newton, M. D. Calculation of Electronic Coupling Matrix Elements for Ground and Excited State Electron Transfer Reactions: Comparison of the Generalized Mulliken?Hush and Block Diagonalization Methods. J. Chem. Phys. 1997, 106, 9213−9226. (55) Newton, M. D. Quantum Chemical Probes of Electron-Transfer Kinetics: the Nature of Donor-Acceptor Interactions. Chem. Rev. 1991, 91, 767−792. (56) Deng, W.-Q.; Goddard, A. G., III Predictions of Hole Mobilities in Oligoacene Organic Semiconductors from Quantum Mechanical Calculations. J. Phys. Chem. B 2004, 108, 8614−8621. (57) Rust, M.; Lappe, J.; Cave, R. J. Multistate Effects in Calculations of the Electronic Coupling Element for Electron Transfer Using the Generalized Mulliken-Hush Method. J. Phys. Chem. A 2002, 106, 3930−3940. (58) Yang, C.-H.; Hsu, C.-P. A Multi-State Fragment Charge Difference Approach for Diabatic States in Electron Transfer: Extension and Automation. J. Chem. Phys. 2013, 139, 154104. (59) 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. (60) Wu, Q.; Voorhis, T. V. Direct Calculation of Electron Transfer Parameters through Constrained Density Functional Theory. J. Phys. Chem. A 2006, 110, 9212−9218. (61) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. A New Algorithm for Monte Carlo Simulation of Ising Spin Systems. J. Comput. Phys. 1975, 17, 10−18. (62) MATLAB 2014b; The MathWorks Inc.: Natick, MA, 2014. (63) Gibson, R. E.; Kincaid, J. F. The Influence of Temperature and Pressure on the Volume and Refractive Index of Benzene. J. Am. Chem. Soc. 1938, 60, 511−518. (64) Deul, R.; Rosenzweig, S.; Franck, E. U. The Dielectric Constant and Density of Benzene to 400 °C and 3000 bar. Ber. Bunsenges. Phys. Chem. 1991, 95, 515−519. (65) dos Santos, F. J. V.; de Castro, C. A. N. Viscosity of Toluene and Benzene Under High Pressure. Int. J. Thermophys. 1997, 18, 367−378. K

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

Article

The Journal of Physical Chemistry C (66) Troisi, A.; Orlandi, G. Hole Migration in DNA: a Theoretical Analysis of the Role of Structural Fluctuations. J. Phys. Chem. B 2002, 106, 2093−2101. (67) Dennington, R.; Keith, T.; Millam, J. GaussView, Ver. 5; Semichem Inc.: Shawnee Mission, KS, 2009. (68) Klimkns, A.; Larsson, S. Reorganization Energies in Benzene, Naphthalene, and Anthracene. Chem. Phys. 1994, 189, 25−31. (69) Datta, A.; Mohakud, S.; Pati, S. K. Electron and Hole Mobilities in Polymorphs of Benzene and Naphthalene: Role of Intermolecular Interactions. J. Chem. Phys. 2007, 126, 144710. (70) Head-Gordon, M.; Pople, J. A.; Frisch, M. J. MP2 Energy Evaluation by Direct Methods. Chem. Phys. Lett. 1988, 153, 503−506. (71) Bazante, A. P.; Davidson, E. R.; Bartlett, R. J. The Benzene Radical Anion: A Computationally Demanding Prototype for Aromatic Anions. J. Chem. Phys. 2015, 142, 204304. (72) Bartlett, R. J.; Purvis, G. D. Many-Body Perturbation Theory, Coupled-Pair Many-Electron Theory, and the Importance of Quadruple Excitations for the Correlation Problem. Int. J. Quantum Chem. 1978, 14, 561−581. (73) Dunning, T. H., Jr. Gaussian Basis Sets For Use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (74) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic Interaction of a Solute with a Continuum. A Direct Utilization of ab initio Molecular Potentials for the Prevision of Solvent Effects. Chem. Phys. 1981, 55, 117−129. (75) Scalmani, G.; Frisch, M. J. Continuous surface charge polarizable continuum models of solvation. I. General formalism. J. Chem. Phys. 2010, 132, 114110. (76) Chang, J.-F.; Sirringhaus, H.; Giles, M.; Heeney, M.; McCulloch, I. Relative Importance of Polaron Activation and Disorder on Charge Transport in High-Mobility Conjugated Polymer Field-Effect Transistors. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 205204. (77) Schein, L. B.; McGhi, A. R. Band-Hopping Mobility Transition in Naphthalene and Deuterated Naphthalene. Phys. Rev. B: Condens. Matter Mater. Phys. 1979, 20, 1631−1639. (78) Nan, G.; Yang, X.; Wang, L.; Shuai, Z.; Zhao, Y. Nuclear Tunneling Effects of Charge Transport in Rubrene, Tetracene, and Pentacene. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 115203. (79) Coropceanu, V.; Cornil, J.; Filho, D. A.; da, S.; Olivier, Y.; Silbey, R.; Bredas, J.-L. Charge Transport in Organic Semiconductors. Chem. Rev. 2007, 107, 926−952. (80) Li, Y.; Coropceanu, V.; Bredas, J.-L. Nonlocal Electron-Phonon Coupling in Organic Semiconductor Crystals: The Role of Acoustic Lattice Vibrations. J. Chem. Phys. 2013, 138, 204713. (81) Fornari, R. P.; Troisi, A. Theory of Charge Hopping Along a Disordered Polymer Chain. Phys. Chem. Chem. Phys. 2014, 16, 9997− 10007. (82) van der Avoird, A.; Podeszwa, R.; Szalewiez, K.; Leforestier, C.; van Harrevelt, R.; Bunker, P. R.; Schnell, M.; von Helden, G.; Meijer, G. Vibration-Rotation-Tunneling States of the Benzene Dimer: an ab initio Study. Phys. Chem. Chem. Phys. 2010, 12, 8219−8240. (83) Prados, A.; Brey, J. J.; Sanchez-Rey, B. A Dynamical Monte Carlo Algorithm for Master Equations with Time-Dependent Transition Rates. J. Stat. Phys. 1997, 89, 709−734.

L

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