Effect of Nuclear Motion on Charge Transport in Fullerenes: A Combined Density Functional Tight Binding—Density Functional Theory Investigation

Introduction

Fullerene-based materials remain the preferred acceptor materials in bulk heterojunction (BHJ) organic solar cells (OSC) (Lee et al., 2008; Park et al., 2009; Cowan et al., 2010; Ganesamoorthy et al., 2017), and fullerenes are efficient electron transport layer materials in planar perovskite solar cells (PSC) (Liang et al., 2015; Nie et al., 2015; Gil-Escrig et al., 2016). Indeed, fullerene-based perovskite solar cells have now achieved power conversion efficiencies (PCE) exceeding 20% which are competitive with PCEs achieved with titania-based cells (Jiang et al., 2017) while being more suitable for large-area production. While addends to the buckyball have to be employed in OSCs to control the BHJ morphology (PCBM is the famous example), bare C60 and C70 [or their mix (Lin et al., 2018)] work well in PSCs (Brédas et al., 2009; Ameri et al., 2013). Fullerenes with addends can also be used in PSCs to control film morphology and band alignment with the perovskite and they also have a significant effect on the charge transport (Pal et al., 2017).

Critical electronic properties of fullerenes in these applications include band alignment with the organic donor or the perovskite, which determines charge separation rate at the interface, and the electron transport rate in bulk and across grain boundaries. The band alignment is controlled by the position of the LUMO (lowest unoccupied molecular orbital), which is centered on the buckyball, even when addends are used (Pal et al., 2017). The LUMO must be lower in energy than the LUMO of the organic donor or the conduction band minimum (CBM) of the perovskite, to provide sufficient driving force for charge separation, where the charge transfer rate can be approximated with the Marcus equation:

where ω_ij is the charge transfer rate between states i and j, λ is the reorganization energy, ΔG_ij is the driving force (difference in Gibbs free energies between the two states), k_B is Boltzmann constant, and T is the temperature. J_ij is the overlap integral between the wavefunctions of states i and j: 〈ψ_i|V_ij|ψ_j〉, where V_ij is the coupling (Coulomb interaction) term. The driving force can be approximated by the differences between corresponding single electron state energies. In general, a larger driving force leads to more efficient charge separation at the price of a lower open circuit voltage V_oc. It is therefore important to optimize the driving force so that it is large enough to provide efficient charge separation but as small as possible to reduce voltage loss. Selection of a particular fullerene (optimal for a given donor) is part of such optimization. The driving force is also important for charge transport in bulk fullerenes, in which case i and j are LUMO (for electron transport) or HOMO (for hole transport) orbitals centered on different molecules. ΔG_ij in this case is non-zero due to reorganization effects as electrons occupy and de-occupy orbitals. The overlap integral J_ij strongly depends on mutual position of molecules.

Ab initio modeling can be used to estimate the LUMO level as well as the reorganization energy; this is relatively easy at the molecular level but is also feasible for solid fullerenes (Ching et al., 1991; Xue et al., 2017). Typically, ab initio models of charge transfer consider equilibrium geometries. Finite temperature nuclear motions are, however, expected to affect the band alignment. Indeed, we have previously shown that room-temperature nuclear dynamics cause changes in LUMO of dyes and polymers by tenths of an eV (Manzhos et al., 2012a,b, 2013; Manzhos, 2013). This can have significant effects on charge separation, especially in interfaces with small driving forces designed for high voltage, as well as on bulk charge transport, via Equation (1). Further, molecular units in organic solids such as fullerenes, which are bound by van der Waals (vdW) forces, can be relatively mobile (Pal et al., 2017). This is expected to have a drastic effect on the overlap integral. Note that as Equation (1) depends (strongly) non-linearly on ΔG_ij and on intermolecular separation, the effects of nuclear motions will not in general average out.

In this work, we study the effects of room-temperature nuclear motions on band structure of C60 and C70 and on electron and hole transport in these fullerenes. We consider transport in bulk, as effects of microstructure are still not feasible to compute ab initio. We perform approximate ab initio (using Density Functional Tight Binding) molecular dynamics and compute expectation values of key quantities entering Equation (1). We show that nuclear motions have a strong effect on both the overlap integral and the exponent in the Marcus equation as well as on the overall charge transport rate, which can become several-fold smaller once the effects of nuclear motions are included.

Methods

Density Functional Theory (DFT) (Hohenberg and Kohn, 1964; Kohn and Sham, 1965) calculations on molecules and dimers for calculating the transfer integral were performed in Gaussian 09 (Frisch et al., 2016) using the 6–31+g(d,p) basis set. The ΔG and reorganization energy were calculated as described in the previous work (Pal et al., 2017). Namely, ΔG estimated from the difference in the redox potentials of the neutral and charged systems: $\Delta G=E_{ox(red)}^{+(-)}-E_{ox(red)}^0$, where $E_{ox(red)}^{+(-)}$ and $E_{ox(red)}^0$ were the oxidation (reduction) potential of ionic and neutral molecules. These potentials were computed as E_ox(red) = E⁺⁽⁻⁾ − E⁰, where E⁺, E⁻, and E⁰ were the energies of the cation, anion, and neutral molecules (Dardenne et al., 2015). The redox potentials were computed with the ωB97XD functional (Chai and Head-Gordon, 2008), as it was reported that a large fraction of exact exchange is desired to accurately compute them (Dardenne et al., 2015; Chen and Manzhos, 2016; Sk et al., 2016). The reorganization energy was calculated using the equation $\lambda ={\lambda }_1+{\lambda }_2=(E_{+(-)}^{*}-E_{+(-)})+(E^{*}-E)$, where $E_{+(-)}^{*}$ is ion energy at the geometry of the neutral molecule, E₊₍₋₎ is ion energy at its optimized geometry, E* is the energy of the neutral molecule at ion's geometry, and E is the energy of the neutral molecule at its optimized geometry. This was done with the B3LYP functional (Becke, 1993).

The transfer integral was computed by DFT/DIPRO (DImer PROjection) approach included in the VOTCA-CTP module (Baumeier et al., 2010; Rühle et al., 2011). The PBE functional (Perdew et al., 1996) was used for transfer integral calculations for computational efficiency (Baumeier et al., 2010). We previously showed (Pal et al., 2017) that this approach results in electron transfer rate for C60 in good agreement with available experimental and computational data reported elsewhere (Günes et al., 2007; Grzegorczyk et al., 2010; Ferguson et al., 2011; Nardes et al., 2012, 2014; Arntsen et al., 2013; Pelzer et al., 2017).

Solid state calculations and ab initio molecular dynamics were performed using the Self-Consistent Charge Density Functional Tight Binding (SCC-DFTB) method (Elstner et al., 1998) in DFTB+ (Aradi et al., 2007). The 3ob-3-1 Slater-Koster parameter set and DFTB3 capability were used (Gaus et al., 2011). Dispersion interactions were modeled with the UFF (universal force field) (Rappé et al., 1992). We previously showed that (dispersion-corrected) DFTB computes molecular and solid-state structures of fullerenes in good agreement with experiment and (dispersion-corrected) DFT, see (Pal et al., 2017, 2018; Lin et al., 2018) for corresponding benchmarks. Specifically, while the DFTB parameterization causes an underestimation of the absolute value of the bandgap, molecular and solid structures as well as relative changes in electronic structure sufficiently reproduce those computed with DFT (Pal et al., 2017, 2018; Lin et al., 2018).

Simulation cells consisted of 8 C60 or C70 molecules. Initial structures were taken from known bcc and hcp structures for C60 and C70, respectively (Quo et al., 1991; Soldatov et al., 2001). The Brillouin zone was sampled at the Gamma point which was sufficient due to the size of the cells (more than 17 Å in each dimension). Molecular dynamics simulations were performed using the Berendsen thermostat/barostat with a time constant of 100 fs at room temperature (300 K). First 100 ps were used for equilibration, and the following 900 ps were used for data analysis. Intermolecular (based on center of mass) distances as well as HOMO and LUMO energy levels were collected during the trajectory.

Results

Effect of Nuclear Motions on Band Structure and Driving Force

The simulation cells of C60 and C70 are shown in Figure 1. The distributions of HOMO and LUMO energy levels for the two fullerenes are shown in Figure 2. The HOMO and LUMO energies at the equilibrium geometry as well as their expectation values over the MD (molecular dynamics) trajectories are listed in Table 1. The expectation values over the trajectories differ from the values at the equilibrium geometry by about 0.1 eV; this is equivalent to an amount of change in ΔG, which was previously associated with changes in electron transfer rate by about a factor of 2 (Koops et al., 2009; Santos et al., 2010). The distributions are broad on the order of 0.1 eV with respect to their expectation value (3σ width, where the standard deviations are about σ = 0.02 eV). The distribution of the exponential of the Marcus equation is shown in Figure 3; it is computed as

where ΔE′ and ΔE″ are deviations from the equilibrium-geometry levels for neighboring molecular units. The values of reorganization energies for electron and hole transport are listed in Table 1.

Figure 1. The simulations cells of C60 (left) and C70 (right).

Figure 2. The distributions of HOMO and LUMO levels in C60 (left) and C70 (right) during molecular dynamics in solid state. Individual molecules are shown next to respective plots.

Table 1. HOMO and LUMO energies at the equilibrium geometry (“Equil.”), their expectation values over MD trajectories (“ < …>”), standard deviations of HOMO and LUMO distributions over MD trajectories (“σ”) as well as reorganization energies λ and driving forces ΔG_eq for electron and hole transport in C60 and C70.

Figure 3. The distributions of the values of the exponent in the Marcus equation (Equation 2) for C60 (left) and C70 (right). The values at the equilibrium geometry and expectation values over molecular dynamics trajectories are indicated by solid and dotted red lines, respectively, with corresponding numbers as inserts. Note the logarithmic scale of the abscissa.

The expectation value of this distribution for C60 and C70 for electron transport (i.e., based on LUMO) is 0.0063 and 0.0068, respectively, compared to values at the equilibrium geometry of 0.0030 and 0.0029, respectively. This change is due to both the different expectation values of LUMO over MD trajectories from their equilibrium values (Table 1) as well as the finite width of the distribution which has a strong effect on E due to the strongly non-linear dependence of E on orbital energies. A similar effect is seen for hole transport (based on HOMO): the expectation value of E over the trajectories for C60 is 2.02 × 10⁻⁴ vs. the equilibrium value of 1.08 × 10⁻⁴ and for C70, 4.85 × 10⁻⁷ vs. the equilibrium value of 1.11 × 10⁻⁷.

Effect of Nuclear Motions on the Overlap Integral and Charge Transfer Rate

Molecular dynamics causes a distribution of inter-molecular separation of C60, C70 units. These distributions (of distances between centers of mass) are shown in Figure 4. The equilibrium geometry C60-C60 and C70-C70 nearest neighbor distances are 9.3 and 11.3 Å, respectively. The distributions have a strong effect on the overlap integral. The overlap integrals as a function of the inter-molecular distance are also shown in Figure 4. The strongly non-linear dependence of the integral on the distance leads to a strong effect of nuclear motions on the overlap integral; for C60, while the value of $J_{ij}^2$ for electron/hole transport at the equilibrium inter-molecular distance is 0.045/0.258 eV, its expectation value over the MD trajectory is 0.0023/2.15 × 10⁻⁶ eV. For C70, these numbers are 0.151/0.002 and 0.0010/0.0012 eV, respectively. Figure 4 was computed by changing the inter-molecular distance starting from nearest C60 and C70 dimers as they are positioned in the crystal at a specific orientation. At selected distances, we computed overlap integrals at other molecular orientations. Not surprisingly, given the strong sensitivity of the overlap integral, different orientations resulted in changes in J_{e, h} of more than an order of magnitude. It is not feasible to fully sample and average the curves in Figure 4 for all possible orientations. Nevertheless, Figure 4 still clearly demonstrates the strong effect of nuclear dynamics on the charge transfer integral and thereby on the rate.

Figure 4. The squared electron transfer integrals (left ordinate axis) for electrons ($J_e^2$, top panels) and holes ($J_h^2$, bottom panels) for C60 (left panels) and C70 (right panels) are shown as black lines vs. intermolecular distance r. The probability density function (pdf ) of r is shown as blue histogram (right ordinate axis).

We also observed that the distributions in HOMO and LUMO energies were uncorrelated with the distribution in inter-molecular distances. This is demonstrated in Figure 5. This was not unexpected in a vdW material like fullerene, where the intramolecular vibrational modes (mostly responsible for HOMO, LUMO changes) have little coupling to inter-molecular motion. The absence of significant correlation permits estimating the expected charge transfer rate with the effect of nuclear motions 〈ω_ij〉 as

where the angular brackets denote expectation values (averages over MD trajectories), e.g., 〈E〉 is the expectation value of the Marcus exponent, Equation (2). The resulting average electron and hole transfer rates for C60 were 6.75 × 10¹¹ and 1.79 × 10⁷ s⁻¹, respectively, and for C70 are 3.32 × 10¹¹ and 2.69 × 10⁷ s⁻¹, respectively. Except for hole transport in C70 (where the equilibrium-geometry charge-transfer integral is very low leading to a very low equilibrium-geometry rate), these rates were significantly lower than the rates computed at the equilibrium geometry, which for C60 were 6.24 × 10¹² and 1.15 × 10¹² s⁻¹, respectively, and for C70 were 2.10 × 10¹³ and 1.16 × 10⁷ s⁻¹, respectively. Rates can thus change by orders of magnitude. Specifically, for hole transport in C70, a significantly higher <E> than its equilibrium value outweighs the dynamic lowering of J_ij, resulting in a computed expectation value for the rate which is higher than its (very low) computed equilibrium value.

Figure 5. Correlations between HOMO, LUMO energies, and intermolecular distances (between centers of mass) during molecular dynamics in C60 and C70. Pearson R² coefficients are also shown.

The rate calculations are very sensitive, due to a very sensitive dependence of the overlap integral on small changes in geometry or orbital distribution as well as due to the exponential dependence on ΔG and the reorganization energy. We therefore believe that such data are reliable mostly on the trend basis.

Conclusions

In this work, we considered effects due to room-temperature atomic and molecular motions on charge transfer rates in solid C60 and C70. We modeled solid C60 and C70 using Density Functional Tight Binding and performed DFTB MD to estimate effects of dynamics on components of the Marcus equation. Room-temperature molecular dynamics cause a spread of frontier orbital energies on the order of 0.1 eV which leads to a distribution of driving forces to charge transfer and an increase by more than a factor of two in the Marcus exponent. We also computed a decrease, in some cases by orders of magnitude, in the expectation value of the overlap integral, leading to an overall decrease in the charge transport rate in most cases.

Considering the high degree of sensitivity of the rates to various approximations (the DIPRO approximation and the level of ab initio treatment within it, the approximate treatment of vdW interactions responsible for the structure of the solid, the classical MD used to model dynamics and the structures used for rate calculations, etc.), the numeric results are likely only semi-quantitative at best. This might have led e.g., to the very low expectation values of hole transfer rates vs. equilibrium geometry rate. Nevertheless, the qualitative trends are likely trustworthy, and they consistently point to a strong effect of nuclear dynamics already at room temperature on the electron and hole transport rate.

The effect of nuclear dynamics should therefore be included when feasible. The calculations are not routine and are multi-stepped, involving solid state MD calculations, dimer generation, and DIPRO model application. We have shown here that effects on the charge transfer integral and on the exponent of the Marcus equation are largely uncorrelated, which significantly simplifies such estimates. Development of approximate schemes, which would allow routine inclusion of dynamics effect on the Marcus rate, is therefore desired, especially those considering orientation distributions of molecules which were not fully considered here. Hopefully this work will help encourage such development.

Author Contributions

SM: conceptualization, project guidance, and text; SA: DFTB calculations and data analysis; AP: charge transport calculations; KY: advisor to SA. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the Ministry of Education of Singapore. KY was supported by MEXT as Priority Issue on Post-K computer (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Prof. Yutaka Matsuo and Dr. Il Jeon of the University of Tokyo for discussions.

References