Chinese Journal of Chemical Physics  2017, Vol. 30 Issue (3): 268-276

#### The article information

Cai-bin Zhao, Zhi-hua Tang, Xiao-hua Guo, Hong-guang Ge, Jian-qi Ma, Wen-liang Wang

Modeling Photovoltaic Performances of BTBPD-PC61BM System via Density Functional Theory Calculations
BTBPD-PC61BM体系光伏性质的密度泛函理论计算
Chinese Journal of Chemical Physics, 2017, 30(3): 268-276

http://dx.doi.org/10.1063/1674-0068/30/cjcp1702016

### Article history

Accepted on: April 7, 2017
Modeling Photovoltaic Performances of BTBPD-PC61BM System via Density Functional Theory Calculations
Cai-bin Zhaoa, Zhi-hua Tanga, Xiao-hua Guoa, Hong-guang Gea, Jian-qi Maa, Wen-liang Wangb
Dated: Received on February 16, 2017; Accepted on April 7, 2017
a. Shaanxi Key Laboratory of Catalysis, School of Chemical and Environmental Science, Shaanxi University of Technology, Hanzhong 723001, China;
b. Key Laboratory for Macromolecular Science of Shaanxi Province, School of Chemistry and Chemical Engineering, Shaanxi Normal University, Xi'an 710062, China
*Author to whom correspondence should be addressed. Cai-bin Zhao, E-mail:zhaocb@snut.edu.cn; Hong-guang Ge, E-mail:gehg@snut.edu.cn, Tel.:+86-916-2641660
Abstract: Designing and fabricating high-performance photovoltaic devices have remained a major challenge in organic solar cell technologies.In this work, the photovoltaic performances of BTBPD-PC61BM system were theoretically investigated by means of density functional theory calculations coupled with the Marcus charge transfer model in order to seek novel photovoltaic systems.Moreover, the hole-transfer properties of BTBPD thin-film were also studied by an amorphous cell with 100 BTBPD molecules.Results revealed that the BTBPDPC61BM system possessed a middle-sized open-circuit voltage of 0.70 V, large short-circuit current density of 16.874 mA/cm2, large fill factor of 0.846, and high power conversion efficiency of 10%.With the Marcus model, the charge-dissociation rate constant was predicted to be as fast as 3.079×1013 s-1 in the BTBPD-PC61BM interface, which was as 3-5 orders of magnitude large as the decay (radiative and non-radiative) rate constant (108-1010 s-1), indicating very high charge-dissociation efficiency (~100%) in the BTBPD-PC61BM system.Furthermore, by the molecular dynamics simulation, the hole mobility for BTBPD thin-film was predicted to be as high as 3.970×10-3 cm2V-1s-1, which can be attributed to its tight packing in solid state.
Key words: BTBPD    PC61BM    Photovoltaic performances    Density functional theory
Ⅰ. INTRODUCTION

In the past 100 years, with the overconsumption for fossil energies (coal, petroleum, and natural gas), environmental pollution problems have received widespread attention, and actively exploring for the clean and renewable energy has being become a hot and focus issue [1, 2]. As one of the most promising long-term solutions for the clean and renewable energy, free-metal photovoltaic technology has attracted intense interests in recent years due to its numerous advantages compared to the commercial inorganic photovoltaic technology, such as low manufacturing cost, flexibility, ease of solvent processing, and large-area capability [3-6]. Previous studies indicated that the high-performance donor materials should meet the following requirements: (ⅰ) narrow optical band gap, (ⅱ) low-lying highest occupied molecular orbital (HOMO) level, and (ⅲ) high hole carrier mobility [7-9]. Unfortunately, electron-donating materials that simultaneously satisfy those three demands are still scarce up to date.

Recently, Cai et al. synthesized a novel small molecule material (BTBPD) with the donor-acceptor-donor (D-A-D) character, and found that the thin-film field effect transistor fabricated with BTBPD had high hole mobility under natural ambient conditions [10]. More interestingly, BTBPD also exhibited the prominent capture to solar radiation, and its strongest absorption peak was found to red-shift to 696 nm in solid state. In a word, all properties suggest that BTBPD should be an excellent electron donor candidate. In this work, taking BTBPD as donor and [6, 6]-phenyl-C $_{61}$ -butyric acid methyl ester (PC $_{61}$ BM) as acceptor, we carried out a systematic theoretical study on the photovoltaic properties for BTBPD-PC $_{61}$ BM system by means of quantum chemistry and molecular dynamics calculations coupled with the incoherent charge hopping model in order to seek novel high-performance photovoltaic systems. In this work, our main objectives are to theoretically explore the applicability of BTBPD as an electron-donating material and estimate the photovoltaic performances of BTBPD-PC $_{61}$ BM system. Theoretical calculations clearly show that BTBPD, as expected, is an excellent electron donor material, and the power conversion efficiency (PCE) of BTBPD-PC $_{61}$ BM system theoretically reaches up to 10%.

Ⅱ. COMPUTATIONAL METHODS

To simplify calculations, the long-alkyl chain (2-ethylhexyl) in BTBPD was replaced with the CH $_3$ , because it has been confirmed that the substituted alkyl in organic compounds has hardly any effect on their electronic structure and optical properties, and merely promotes solubility [11-13]. All stable species were fully optimized without any symmetry constraints by means of density functional theory (DFT) calculations with the B3LYP hybrid functional [14] and the 6-31G(d) basis set, with subsequent frequency calculations to confirm that they were true minima of potential energy surface. Based on the optimized geometries, the UV-Vis spectrum for BTBPD was simulated with the time-dependent density functional theory (TD-DFT) [15, 16] and the B3LYP/6-31G(d) scheme. In order to determine the most reasonable geometry of BTBPD-PC $_{61}$ BM complex, the detailed potential-surface scan was carried out between PC $_{61}$ BM and BTBPD with the CAM-B3LYP-D3(BJ)/6-31G(d) method [17, 18]. As seen in FIG. S1 (supplementary materials), the BTBPD-PC $_{61}$ BM complex was found to be the most stable when the centroids distance of PC $_{61}$ BM and BTBPD is at 8.0 Å, which is in good agreement with the recent studies [19, 20]. Then, in subsequent calculations, the centroids distance of PC $_{61}$ BM and BTBPD was invariably fixed at 8.0 Å. In addition, in this work the influence of molecular orientation was also considered. As is shown in FIG. S2 (supplementary materials), the molecular orientation affects a little on the BTBPD-PC $_{61}$ BM complex. Based on optimized structures for PC $_{61}$ BM, BTBPD, and BTBPD-PC $_{61}$ BM complex, total density of states (TDOS) and partial density of states (PDOS) were visualized with the Multiwfn 3.3.6 software developed by Lu et al. [21, 22]. All quantum chemistry calculations were carried out with the Gaussian 09 software [23].

Ⅲ. RESULTS AND DISCUSSION A. Photovoltaic performances of BTBPD-PC61BM system 1. Electronic properties and open circuit-voltage

BTBPD and PC $_{61}$ BM molecular structures are depicted in FIG. 1. The geometric optimization revealed that BTBPD molecule has a near planar conformation (FIG. S3 in supplementary materials), and the dihedral angle ( $\alpha$ ) between its adjacent units is close to 20 $^{\circ}$ , which indicates its good $\pi$ -conjugated character. FIG. 2 shows the TDOS and PDOS of PC $_{61}$ BM, BTBPD, and BTBPD-PC $_{61}$ BM complex. With the DOS, it is very easy to directly observe the contribution from each substituent to the frontier molecular orbital (FMO). As seen, in PC $_{61}$ BM molecule all density of HOMOs and LUMOs was found to concentrate on the C $_{60}$ sphere in the energy range from -10.0 eV to $\mbox{2.0 eV}$ , and the contribution from the substituent (methyl-4-phenylbutanoate) is very small, meaning the substituent only enhances the C $_{60}$ solubility, and has a little influence on its electronic properties, which agrees well with the previous experimental studies [24, 25]. Moreover, it can be noticed that CH $_3$ contributes very small to the HOMO and LUMO in BTBPD molecule, verifying it is rational to replace 2-ethylhexyl with CH $_3$ in the current work. Interestingly, the benzo[b]thiophene (BT) and bipyrrolylidene-2, 2 $'$ (1H, 1 $'$ H)-dione (BPD) in BTBPD molecule contribute very much to both the HOMO and the LUMO, denoting the HOMO and LUMO of BTBPD delocalize over the molecular skeleton, rather than centralize at a certain molecular fragment, which benefits the rapid charge-transfer between two molecules. As for BTBPD-PC $_{61}$ BM complex, the HOMO and the LUMO exhibit an obvious separation characteristic, and the HOMO completely locates on the BTBPD, while the LUMO mainly centuries on the PC $_{61}$ BM, which suggests the easy formation of BTBPD $^\cdot$ $^+$ -PC $_{61}$ BM $^\cdot$ $^{-}$ charge-separated state. According to the previous study, the open circuit-voltage for organic solar cells (OSCs), $V_{\textrm{oc}}$ , can be estimated with [26]:

 FIG. 1 Molecular structures of BTBPD and PC61BM.
 FIG. 2 Total density of states and partial density of states of (a) PC61BM, (b) BTBPD, and (c) PC61BM-BTBPD complex.
 $\begin{eqnarray}V_{\textrm{oc}} = \frac{1}{e} \left( {\left| {E_{\textrm{HOMO}} \left(\textrm{D}\right)} \right| - \left| {E_{\textrm{LUMO}} \left(\textrm{A}\right)} \right|} \right) - 0.3\end{eqnarray}$ (1)

where $E_{\textrm{HOMO}}$ ( $\textrm{D}$ ) and $E_{\textrm{LUMO}}$ (A) are the HOMO level of donor and the LUMO level of PC $_{61}$ BM, $e$ is the electron charge, and the value of 0.3 is an empirical factor. Then, based on the experiment HOMO (-5.0 eV [10]) for BTBPD and LUMO (-4.0 eV [27, 28]) for PC $_{61}$ BM, the $V_{\textrm{oc}}$ was estimated to be as large as 0.70 V for BTBPD-PC $_{61}$ BM system. More interestingly, the PCE of BTBPD-PC $_{61}$ BM system was predicted to be over 10% (FIG. 3) by means of the Scharber diagram, indicating the BTBPD-PC $_{61}$ BM system is a promising OSC candidate.

 FIG. 3 Predicted PCE for BTBPD-PC61BM cell with the Scharber diagram.
2. Charge binding energy and optical absorption properties

As is well-known, the charge binding energy ( $E_\textrm{b}$ ) is one of the most parameters in photovoltaic devices, which is directly related to the charge separation. Usually $E_\textrm{b}$ is taken as the difference between the transport gap ( $E_\textrm{t}$ ) and the optical band gap ( $E_{\textrm{opt}}$ ). The former is the difference between the adiabatic ionization potential ( $E_{\textrm{AIP}}$ ) and electron affinity ( $E_{\textrm{AEA}}$ ) of donor in the solid state, while the latter is taken as the first-singlet emission energy ( $E_\textrm{m}$ ). Then, the $E_\textrm{b}$ can be calculated as the following expression [29]:

 $\begin{eqnarray}E_\textrm{b} = E_{\textrm{AIP}} - E_{\textrm{AEA}} - E_\textrm{m}\end{eqnarray}$ (2)

As seen in Eq.(2), to calculate $E_\textrm{b}$ , the $E_{\textrm{AIP}}$ and $E_{\textrm{AEA}}$ of donor in solid state firstly should be calculated. Here, the $E_{\textrm{AIP}}$ / $E_{\textrm{AEA}}$ of solid BTBPD was estimated via the scheme reported by Schwenn et al. [30], which has been verified to be an excellent selection that estimates the electronic properties of organic materials in the solid state. Table Ⅰ shows calculated $E_{\textrm{AIP}}$ and $E_{\textrm{AEA}}$ values for BTBPD in gas phase and solid state with two different DFT methods. It can be noted that $E_\textrm{b}$ estimated by two methods is as large as 1.596 and 1.733 eV in gas phase, which is much larger compared with those measured values of 0.2-1.0 eV in many organic materials [31], which can be attributed to the solid stacking effect. Comparing the results in gas phase to the ones in solid state, it can be noticed that in gas phase $E_{\textrm{AIP}}$ is larger, while $E_{\textrm{AEA}}$ is smaller, which is similar to the measured and theoretical results in acenes [32]. According to the calculated $E_{\textrm{AIP}}$ , $E_{\textrm{AEA}}$ , and $E_\textrm{m}$ for the solid BTBPD, the $E_\textrm{b}$ was estimated to be about $\mbox{0.594 eV}$ . The precious study showed that the exciton is unstable when $E_\textrm{b}$ $<$ $k_\textrm{B}T$ , which amounts to 0.025 eV at room temperature [33]. According to the calculated $E_\textrm{b}$ for BTBPD, it can be deduced that the photo-induced exciton in BTBPD is relatively stable, which can be efficiently transported to the BTBPD-PC $_{61}$ BM interface without rapidly being decayed in transit.

Table Ⅰ Calculated $E_{\textrm{AIP}}$ , $E_{\textrm{AEA}}$ , and $E_\textrm{b}$ values in the gas and solid state for BTBPD with two different DFT methods.

As is known to all, the good harvest for solar radiation is essential for efficient dye sensitizers, which determines the short-circuit current density $J_{\textrm{sc}}$ of DSC devices to some extent. To explore reliable DFT methods estimating optical absorption properties for BTBPD, a set of popular DFT methods were tested. As seen in Table Ⅱ, compared with the experimental value, the B3LYP hybrid functional can estimate accurately the excited energy of BTBPD, and the derivation between the theoretical and experimental values is only about 3.0 nm (about 0.008 eV). Moreover, it can be noticed that the strongest absorption in UV-Vis spectrum for the BTBPD molecule can be assigned to the $\pi-\pi^*$ type, and dominated completely by the electron transition of HOMO $\rightarrow$ LUMO ( $\sim$ 100%). In addition, as observed in FIG. 4, the TD-B3LYP/6-31G(d) method can reproduce well the UV-Vis spectrum of BTBPD molecule in solid state, which confirms again the method reliability used in this work. Since the HOMO and LUMO for BTBPD distribute the whole molecule skeleton (FIG. S4 in supplementary materials), the lowest-excited singlet is a typically local excited state, and no obvious charge is transferred in light absorption process.

 FIG. 4 Simulated and experimental absorption spectra for BTBPD in solid state.
Table Ⅱ Calculated excited energies (λmax), molar absorption coefficients ( $\varepsilon$ ), oscillator strengths (f), and main configuration for BTBPD with different DFT methods coupled with the 6-31G(d, p) basis set.
3. Short-circuit current density Jsc, fill factor FF, and PCE η

Short-circuit current density $J_{\textrm{sc}}$ is another key parameter that determines the PCE of OSC devices, which can be expressed as [34, 35]:

 $\begin{eqnarray}J_{\textrm{sc}} = q\int_0^\infty {\eta _{\textrm{EQE}} (λ )} × S(λ )\textrm{d}λ\end{eqnarray}$ (3)

where $S$ ( $λ$ ) is incident photon-to-current conversion efficiency at a fixed wavelength, $e$ is the unit charge, and $\eta _{\textrm{EQE}} (λ )$ is the external quantum efficiency. The $\eta _{\textrm{EQE}} (λ )$ term can be described as the product of $\eta_λ$ (light-harvesting efficiency), $\eta_{\textrm{CT}}$ (charge transfer efficiency), and $\eta_{\textrm{coll}}$ (charge collection efficiency) [36],

 $\begin{eqnarray}\eta _{\textrm{EQE}} = \eta _λ \eta _{\textrm{CT}} \eta _{\textrm{coll}}\end{eqnarray}$ (4)

where $\eta_λ$ can be calculated as $\eta_λ$ =1-10 $^{-f}$ , $f$ is the oscillator strength. Then, to estimate the maximum $J_{\textrm{sc}}$ , we sent the $\eta_{\textrm{CT}}$ =1.0 and the $\eta_{\textrm{coll}}$ =1.0. Our calculation showed that $f$ is about 1.3963 at the lowest-excited singlet state for BTBPD, then yielding the $\eta_λ$ =0.752. FIG. 5 shows that the simulated $\eta_λ$ and $J_{\textrm{sc}}$ with the above-mentioned parameters. As seen, the $J_{\textrm{sc}}$ was estimated to be as high as 16.874 mA/cm $^2$ for the BTBPD-PC $_{61}$ BM system, which can be attributed to its strong spectral response. In addition, it can be noticed that the $\eta_λ$ is as large as 81.7% in visible region. Relatively, BTBPD has a weak capture for ultraviolet radiation ( $\eta_λ$ $\approx$ 57%). For the $FF$ calculation, an approximate scheme can be expressed as [37, 38],

 FIG. 5 Predicted $J_{\textrm{sc}}$ and $\eta_λ$ for BTBPD-PC $_{61}$ BM cell.
 $\begin{eqnarray}FF=\frac{\nu_{\textrm{oc}}-\textrm{ln}(\nu_{\textrm{oc}}+0.72) }{\nu_{\textrm{oc}}+1}\end{eqnarray}$ (5)

where $\nu_{\textrm{oc}}$ is the dimensionless voltage, which can be estimated with the $V_{\textrm{oc}}$ [39, 40],

 $\begin{eqnarray} \nu_{\textrm{oc}}=\frac{qV_{\textrm{oc}}}{nk_\textrm{B}T}\end{eqnarray}$ (6)

where $k_\textrm{B}$ , $T$ and $q$ are Boltzmann constant, temperature (here, we set $T$ =300 K), and elementary charge respectively, $n$ is an ideality factor relating to an ideal ( $n$ =1) or non-ideal ( $n$ $>$ 1) diode [41], organic solar cells typically have ideality factors in the range of 1.5-2.0 due to their inherent disorder [42]. According to the calculated $V_{\textrm{oc}}$ (0.70 V) for the BTBPD-PC $_{61}$ BM system, the $\nu_{\textrm{oc}}$ was estimated to be 27.08 at $n$ =1.0 and 13.54 at $n$ =2, then, the $FF$ for P $_{61}$ BM-BTBPD was predicted to be as high as 0.748 ( $n$ =2.0) and 0.846 ( $n$ =1.0), in excellent agreement with measured values in most OSC devices. According to the previous study, the PCE ( $\eta$ ) of OSC devices can be estimated with the following equation [43, 44]

 $\begin{eqnarray}\eta = \frac{{P_{\max } }}{{P_{\textrm{in}} }} = \frac{{V_{\textrm{oc}} J_{\textrm{sc}} }}{{P_{\textrm{in}} }}FF\end{eqnarray}$ (7)

where $P_{\textrm{max}}$ and $P_{\textrm{in}}$ (=100 mW/cm $^2$ ) are the maximum and incident power respectively. With the calculated $V_{\textrm{oc}}$ , $J_{\textrm{sc}}$ , and $FF$ , the PCE of BTBPD-PC $_{61}$ BM system was predicted to be 8.83% ( $n$ =2.0) and 9.99% ( $n$ =1.0), which is slightly smaller than the value ( $>$ 10%) estimated by the Scharber diagram.

B. Charge dissociation and recombination rates

Generally, the charge transfer in organic photoelectric materials obeys the incoherent charge-hopping mechanism, and the transfer rate constant between donor and acceptor, $k_{\textrm{DA}}$ , can be evaluated via the Marcus model [45, 46],

 $\begin{eqnarray}k_{\textrm{DA}} = \frac{{2\pi }}{h}\sqrt {\frac{\pi }{{λ k_\textrm{B} T}}} \left| {V_{\textrm{DA}} } \right|^2 \exp \left[{ - \frac{{(\Delta G + λ )^2 }}{{4λ k_\textrm{B} T}}} \right]\end{eqnarray}$ (8)

where $λ$ is the total reorganization energy, $V_{\textrm{DA}}$ is the effective charge transfer integral between donor and acceptor, $\Delta G$ is the Gibbs free energy change between the initial and final states, $k_\textrm{B}$ is Boltzmann constant, $h$ is Planck constant, and $T$ is the temperature (here, we set $T$ =300 K).

1. Gibbs free energy change in charge dissociation and recombination

As seen in Eq.(8), the Gibbs free energy change, $\Delta G$ , has a remarkable influence on the $k_{\textrm{DA}}$ . Generally, the $\Delta G$ can be estimated as the energy difference in the final and initial states, accounting for the Coulombic attraction between the two charges in charge-separated state. Thus, for the charge-dissociation, the $\Delta G$ is written as [47],

 $\begin{eqnarray}\Delta G_{\textrm{dis}} \hspace{-0.15cm}&=&\hspace{-0.15cm} E_\textrm{D}^ + (Q^ + ) + E_\textrm{A}^ - (Q^ - ) - E_\textrm{D}^* (Q^* ) - \nonumber \\&&\hspace{-0.15cm}E_\textrm{A}^0 (Q^0 ) + \Delta E_{\textrm{coul}}\end{eqnarray}$ (9)

where $E_\textrm{D}^*(Q^*)$ and $E_\textrm{D}^+(Q^+)$ is the total energy of isolated donor in its equilibrium geometry of lowest singlet-excited and cationic state, $E_\textrm{A}^{-} (Q^{-} )$ and $E_\textrm{A}^0 (Q^0 )$ are the total energies of isolated acceptor in its equilibrium geometries of anionic and ground states, $\Delta$ $E_{\textrm{coul}}$ is the Coulombic attraction between donor and acceptor in charge-separated state, which can be estimated with the flowing equation,

 $\begin{eqnarray}\Delta E_{\textrm{coul}} \hspace{-0.1cm} = \hspace{-0.1cm}\sum\limits_{\textrm{D}^ + } \hspace{-0.1cm}{\sum\limits_{\textrm{A}^ - } {\frac{{q_{\textrm{D}^ + }q_{\textrm{A}^ - } }}{{4\pi \varepsilon _0 \varepsilon _\textrm{s} r_{\textrm{D}^ + \textrm{A}^ - } }}} } \hspace{-0.1cm}-\hspace{-0.1cm}\sum\limits_{\textrm{D}^* } \hspace{-0.1cm}{\sum\limits_\textrm{A} {\frac{{q_{\textrm{D}^*}q_\textrm{A} }}{{4\pi \varepsilon _0 \varepsilon_\textrm{s} r_{\textrm{D}^*\textrm{A}} }}} }\end{eqnarray}$ (10)

where $q_\textrm{D}$ and $q_\textrm{A}$ are the atomic charges on donor and acceptor in their relevant states with a separation distance, $r_{\textrm{DA}}$ . The $\varepsilon_0$ is the vacuum dielectric constant (8.854 $×$ 10 $^{-12}$ F/m), and the $\varepsilon_\textrm{s}$ is the static dielectric constant of medium. Similarly, the Gibbs free energy change ( $\Delta G_{\textrm{rec}}$ ) in charge recombination can also be estimated with the expression similar to Eq.(9) and Eq.(10). Here, the $\varepsilon_s$ is estimated with the Clausius-Mossotti equation [48],

 $\begin{eqnarray}\varepsilon _\textrm{s} = \left( {1 + \frac{{8\pi \overline{\alpha} }}{{3V}}} \right)\left( {1 - \frac{{4\pi \overline{\alpha} }}{{3V}}} \right)^{ - 1}\end{eqnarray}$ (11)

where $V$ is the Connolly molecular volume, $\overline{\alpha}$ is the isotropic component of molecular polarizability, $\bar \alpha = 1/3\sum\limits_{} {_i{\alpha _{ii}}}$ , and the $\alpha_{ii}$ is the diagonal matrix elements of first-order polarizability tensor. Calculations show that the $\varepsilon_\textrm{s}$ is 3.653 for BTBPD, which is in accord with the measured values (ranging from 2.0 to 5.0 [49, 50]) in most organic materials. As for PC $_{61}$ BM, the experimental $\varepsilon_\textrm{s}$ value of 3.9 [51] was used in this work. The total $\varepsilon_\textrm{s}$ of BTBPD-PC $_{61}$ BM system was taken as an average of their respective contributions. Our calculation showed that the $\Delta G_{\textrm{dis}}$ is about -0.316 eV, while the $\Delta G_{\textrm{rec}}$ is smaller (-0.640 eV). As seen, the $\Delta G$ $_{\textrm{dis}}$ and $\Delta G_{\textrm{rec}}$ are consistently calculated to be negative, denoting that the charge-dissociation and charge-recombination are always favorable in thermodynamics. In addition, the smaller $\Delta G_{\textrm{rec}}$ indicated a larger driving force in charge-recombination process.

2. Reorganization energies in charge dissociation and recombination

Generally, in organic solids the total reorganization energy ( $λ$ ) of electron transfer can be divided into two parts, namely the internal reorganization energy ( $λ_{\textrm{int}}$ ) and the external one ( $λ_{\textrm{ext}}$ ). The $λ_{\textrm{int}}$ term can be calculated with the adiabatic potential energy surface (PES) method [52, 53]. In the case of charge dissociation, the $λ_{\textrm{int}}$ is actually taken as an average of the following $λ_1$ and $λ_2$ [54],

 $\begin{eqnarray}λ _1 \hspace{-0.1cm}& =&\hspace{-0.1cm} (E_\textrm{D}^* (Q^ + ) + E_\textrm{A}^0 (Q^ {-} )) - (E_\textrm{D}^* (Q^* ) + E_\textrm{A}^0 (Q^0 ))\\ \end{eqnarray}$ (12)
 $\begin{eqnarray}λ _2 \hspace{-0.1cm}& =&\hspace{-0.1cm} (E_\textrm{D}^ + (Q^* ) + E_\textrm{A}^ - (Q^0 )) - (E_\textrm{D}^ + (Q^ + ) + E_\textrm{A}^ - (Q^\textrm{-})) \end{eqnarray}$ (13)

where $E_\textrm{D}^* (Q^+)$ and $E_\textrm{D}^* (Q^*)$ are the energies of donors in the lowest excited-state with the equilibrium geometries of cationic and excited state respectively, $E_\textrm{D}^ + (Q^*)$ and $E_\textrm{D}^ + (Q^ + )$ are the energies of donors in the cationic states with the equilibrium geometries of excited and cationic states respectively, $E_\textrm{A}^0 (Q^ {-} )$ and $E_\textrm{A}^0 (Q^0 )$ are the energies of acceptors in the neutral states with the equilibrium geometries of anion and neutral states, respectively, and $E_\textrm{A}^ {-} (Q^0)$ / $E_\textrm{A}^ 0 (Q^{-} )$ are the energies of acceptors in the anionic states with the equilibrium geometries of neutral and anionic states. As seen in Table Ⅲ, our calculation showed that the $λ_{\textrm{int}}$ ( $λ_{\textrm{dis}}$ ) is $\mbox{0.191 eV}$ in charge-dissociation process for BTBPD-PC $_{61}$ BM, which remarkably increases to 0.348 eV in the case of charge recombination. Compared with the $λ_{\textrm{int}}$ , the $λ_{\textrm{ext}}$ was difficult to be accurately calculated. Here, we used the classical dielectric continuum model initially developed by Marcus for the electron-transfer reaction between spherical ions in solution to estimate the $λ_{\textrm{ext}}$ . According to this model, the $λ_{\textrm{ext}}$ term is given by [55],

 $\begin{eqnarray}\hspace{-0.2cm}λ _{\textrm{ext}} \hspace{-0.1cm} = \hspace{-0.1cm}\frac{{(\Delta e)^2 }}{{8\pi \varepsilon _0 }}\left( {\frac{1}{{\varepsilon _{\textrm{op}} }} \hspace{-0.1cm}-\hspace{-0.1cm} \frac{1}{{\varepsilon _\textrm{s} }}} \right)\hspace{-0.2cm}\left( {\frac{1}{{R_\textrm{D} }}\hspace{-0.1cm} +\hspace{-0.1cm} \frac{1}{{R_\textrm{A} }} \hspace{-0.1cm}-\hspace{-0.1cm} 2\sum\limits_\textrm{D} {\sum\limits_\textrm{A} {\frac{{q_\textrm{D} q_\textrm{A} }}{{r_{\textrm{DA}} }}} } } \right)\end{eqnarray}$ (14)

where $\varepsilon_{\textrm{op}}$ is the optical dielectric constant of medium, $R_\textrm{D}$ (=6.41 Å for BTBPD) and $R_\textrm{A}$ (=6.50 Å for PC $_{61}$ BM) are the effective radii of donor and acceptor estimated as the radius of sphere having the same surface as the surface accessible area of molecule. The $q_\textrm{D}$ and $q_\textrm{A}$ denote the atomic charges on the ions. The $\varepsilon_{\textrm{op}}$ which can be estimated with the Lorentz-Lorenz equation [56, 57],

 $\begin{eqnarray}\varepsilon _{\textrm{op}} = n^2 = \frac{{V_\textrm{m} + 2R}}{{V_\textrm{m} - R}}\end{eqnarray}$ (15)

where $n$ is the refractive index, $V_\textrm{m}$ is the molar volume ( $V_\textrm{m}=M$ / $\rho$ , $M$ is the molar mass, and $\rho$ is the density of material), $R$ is the molar refraction. Here, the $\rho$ was estimated with the molecular dynamics method, and the simulated detail was shown in the supplementary materials. Our results showed the $\varepsilon_{\textrm{op}}$ and $\rho$ for BTBPD were equal to 2.960 and 1.312 g/cm $^3$ respectively. As for PC $_{61}$ BM, the experimental refractive index ( $n$ =1.866) is used to estimate to the $\varepsilon_{\textrm{op}}$ , which is equal to 3.482 according to our estimation. With the above-mentioned parameters, the $λ_{\textrm{ext}}$ can be conveniently obtained. At the case of BTBPD-PC $_{61}$ BM, the $λ_{\textrm{ext}}$ is equal to $\mbox{0.103 eV}$ . Summary, the $λ$ is $\mbox{0.294 eV}$ in charge-dissociation process for BTBPD-PC $_{61}$ BM system, which remarkably increases to 0.451 eV for the charge-recombination process.

Table Ⅲ Calculated $λ_{\textrm{dis}}$ , $λ_{\textrm{rec}}$ , and $λ_{\textrm{ext}}$ values in the gas and solid state for BTBPD with two different DFT methods.
3. Charge transfer integral in charge dissociation and recombination

As seen in Eq.(8), the $V_{\textrm{DA}}$ is an important parameter that determines the $k_{\textrm{DA}}$ , in this work, the direct-coupling (DC) method coupled with the PW91PW91/6-31G(d) scheme was used to estimate $V_{\textrm{DA}}$ [58, 59], which have been illustrated to present the most accurate $V_{\textrm{DA}}$ value at the DFT level [60, 61]. In terms of the DC scheme, the $V_{\textrm{DA}}$ value of charge transfer can be calculated by the following expression [62],

 $\begin{eqnarray}V_{\textrm{D}(i)\textrm{A}(j)} = \frac{{T_{\textrm{D}(i)\textrm{A}(j)} - 0.5(e_{\textrm{D}(i)} + e_{\textrm{A}(j)} )S_{\textrm{D}(i)\textrm{A}(j)} }}{{1 - S_{\textrm{D}(i)\textrm{A}(j)}^2 }}\end{eqnarray}$ (16)

where $T_{\textrm{D}(i)\textrm{A}(j)}$ is the charge transfer integral of the $i$ th molecular orbital of donor and the $j$ th molecular orbital of acceptor, $S_{\textrm{D}(i)\textrm{A}(j)}$ is the spatial overlap integral of the above two molecular orbitals, and $e_{\textrm{D}(i)}$ / $e_{\textrm{A}(j)}$ is the site energy. $T_{\textrm{D}(i)\textrm{A}(j)}$ , $S_{\textrm{D}(i)\textrm{A}(j)}$ , and $e_{\textrm{D}(i)}$ / $e_{\textrm{A}(j)}$ can be obtained from the $T_{\textrm{D}(i)\textrm{A}(j)}=\langle \psi_{\textrm{D}(i)}|F^{\textrm{KS}}|\psi _{A(j)}\rangle$ , $S_{\textrm{D}(i)\textrm{A}(j)}=\langle \psi_{\textrm{D}(i)}|\psi_{\textrm{A}(j)}\rangle$ , and $e_{\textrm{D}(i)}$ / $e_{\textrm{A}(j)}=\langle \psi_{\textrm{D}(i)}|\psi_{\textrm{A}(j)}|F^{\textrm{KS}}|\psi_{\textrm{D}(i)}|\psi_{\textrm{A}(j)}\rangle$ . Among them, $\psi_{\textrm{D}(i)}$ is the HOMO (for charge-recombination) or LUMO (for charge-dissociation) of donor, $\psi_{\textrm{A}(j)}$ is the LUMO of acceptor, and $F^{\textrm{KS}}$ is the Kohn-Sham matrix of donor-acceptor system. The $F^{\textrm{KS}}$ can be obtained from

 $\begin{eqnarray} F^{\textrm{KS}} = SC\varepsilon C^{ - 1}\end{eqnarray}$ (17)

where $S$ is the intermolecular overlap matrix, $C$ is the molecular orbital coefficient matrix from the isolated monomer, and $\varepsilon$ is the orbital energy from one-step diagonalization without iteration. Generally, the $V_{\textrm{DA}}$ in the charge-dissociation process is taken as the coupling between the LUMO of donor and acceptor. However, since the LUMO+1 and LUMO+2 in PC $_{61}$ BM are degenerated energetically with its LUMO [63], the $V_{\textrm{DA}}$ between the LUMO of BTBPD and the LUMO+1/LUMO+2 of PC $_{61}$ BM was also computed. Finally, the average $V_{\textrm{DA}}$ value (( $V_1V_2V_3$ ) $^{1/3}$ ) was viewed as the total $V_{\textrm{DA}}$ and then applied to estimate the charge-dissociation rate. As for the charge-recombination, the same treatment was done. Calculations show that the $V_{\textrm{DA}}$ in the charge dissociation for the BTBPD-PC $_{61}$ BM is -31.82 meV, which is equal to 20.37 meV for the charge-recombination process. Based on the calculated $λ$ and $V_{\textrm{DA}}$ values, the charge-dissociation ( $k_{\textrm{dis}}$ ) and charge-recombination ( $k_{\textrm{rec}}$ ) rate constants were estimated to be as high as 3.079 $×$ 10 $^{13}$ and 4.808 $×$ 10 $^{12}$ s $^{-1}$ respectively in BTBPD-PC $_{61}$ BM complex. Recent studies illustrated that the decay rate constant ( $k_\textrm{d}$ ) of excited organic molecules typically ranges 1.0 $×$ 10 $^8$ s $^{-1}$ to 1.0 $×$ 10 $^{10}$ s $^{-1}$ [64]. Our results showed that the $k_{\textrm{dis}}$ is larger than the $k_\textrm{d}$ 3-5 orders of magnitude, which indicates high charge-dissociation efficiency ( $\sim$ 100%) in the BTBPD-PC $_{61}$ BM system. In addition, although the $k_{\textrm{rec}}$ is relatively large, the charge-recombination efficiency is still very low. According to previous studies, the electron transferred onto PC $_{61}$ BM can be rapidly converted to the triplet state from the singlet state [63, 64], which remarkably hinders from the recombination of free carriers.

C. Hole transfer rate and hole mobility in BTBPD thin-film

As is known to all, the charge transport ability of donor remarkably affects the solar cell's performance. Thus, it is essential to discuss charge transport properties of BTBPD thin-film. Generally, the charge transport ability of organic materials can be chartered with its carrier mobility, $\mu$ , which can be calculated by means of the Einstein-Smoluchowski equation [67, 68],

 $\begin{eqnarray}\mu = \frac{{eD}}{{k_{\rm{B}} T}}\end{eqnarray}$ (18)

where $D$ is diffusion coefficient, $e$ is elementary charge, $k_\textrm{B}$ is Boltzmann constant, and $T$ is absolute temperature, respectively. The $D$ can be estimated by means of the following appropriate relation [69, 70]:

 $\begin{eqnarray}D \approx \frac{1}{{2n}}\sum\limits_i {d_i^2 k_i P_i }\end{eqnarray}$ (19)

where $n$ is the spatial dimensionality, which is 3 in organic solids, $d_i$ is the centroids distance of the $i$ th hopping dimer, $k_i$ is the charge transfer rate constant, and $P_i$ is the hopping probability, $P_i=k_i$ / $\sum_ik_i$ . In this work, the charge mobility of BTBPD thin-film was evaluated by means of an amorphous cell with 100 BTBPD molecules built by means of the molecular dynamics simulation. Table Ⅳ lists the calculated the $λ_{\textrm{int}}$ term with two different DFT methods. As seen, the CAM-B3LYP/6-31G(d, p) scheme presents quite large $λ_{\textrm{int}}$ values due to the long-range correlation effect. In addition, it can be noticed that the $λ_{\textrm{int}}$ in solid state is obviously smaller than that in gas phase, denoting that the solid stack to some extent, limits the structural relaxation of BTBPD molecule in charge transfer process. Since the donor materials in OSC devices usually keep in solid state under operating conditions, the $λ_{\textrm{int}}$ estimated in the solid state is more reasonable. To explore possible charge transfer dimers, 21 molecular pairs with the relatively large $V_{\textrm{DA}}$ values were abstracted from the optimized amorphous cell, and their geometries, centroids distances, as well as estimated $V_{\textrm{DA}}$ values were shown in Table S1 (supplementary materials). Based on the $λ_{\textrm{int}}$ the solid state and $V_{\textrm{DA}}$ values, the hole carrier mobility, $\mu_\textrm{h}$ , was estimated to be as high as 3.970 $×$ 10 $^{-3}$ cm $^2$ V $^{-1}$ s $^{-1}$ in the solid BTBPD, which is in excellent agreement with its measured value (3.0 $×$ 10 $^{-3}$ -8.4 $×$ 10 $^{-3}$ cm $^2$ V $^{-1}$ s $^{-1}$ [10]. According to the previous investigation, for high-performance OSC devices, the $\mu_\textrm{h}$ should be not less than 10 $^{-3}$ cm $^2$ V $^{-1}$ s $^{-1}$ [26]. Our estimation showed as a potential donor material of OSC devices, the BTBPD can rapidly transports holes.

Table Ⅳ Calculated $λ_{\textrm{int}}$ for BTBPD in solid and gas states with two different DFT methods.
Ⅳ. CONCLUSION

In summary, BTBPD-PC $_{61}$ BM as a promising OSC system was theoretically studied by means of quantum-chemical and molecular dynamics calculations. Results showed that BTBPD-PC $_{61}$ BM system possesses middle-sized open-circuit voltage (0.7 V), large short-circuit current density (16.874 mA/cm $^2$ ), high fill factor (0.846), and high PCE ( $>$ 10%). In addition, BTBPD was also revealed to possess the strong optical response, and suitable charge-binding energy (0.457 eV). Using the Marcus model, the $k_{\textrm{dis}}$ was estimated to be as large as 3.079 $×$ 10 $^{13}$ s $^{-1}$ in the BTBPD-PC $_{61}$ BM blend, which indicated very high charge-dissociation efficiency. Moreover, by means of an amorphous cell model with 100 BTBPD molecules, the hole carrier mobility of BTBPD thin film was predicted to be as high as 3.970 $×$ 10 $^{-3}$ cm $^2$ V $^{-1}$ s $^{-1}$ . In brief, our calculation shows that BTBPD is a very potential donor material, and the BTBPD-PC $_{61}$ BM system is a high-performance OSC candidate. However, these results need to be verified by experiments.

Supplementary material: Detailed potential-surface scan, optimized BTBPD geometry, calculated the lowest-excited energy for BTBPD, HOMO and LUMO of BTBPD, and detailed description for molecular dynamics simulation are shown.

Ⅴ. ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No.21373132, No.21502109, No.21603133), the Education Department of Shaanxi Provincial Government Research Projects (No.16JK1142, No.16JK1134), and the Scientific Research Foundation of Shaanxi University of Technology for Recruited Talents (No.SLGKYQD2-13, No.SLGKYQD2-10, No.SLGQD14-10).

 FIG. S1 Energy evolution as a function of the centroids distance between PC61BM and BTBPD
 FIG. S2 Energy evolution as a function of the orientation angle (θ) between PC61BM and BTBPD
 FIG. S3 Optimized BTBPD with the B3LYP/6-31G(d, p) method
 FIG. S4 HOMO (up) and LUMO (down) of BTBPD calculated with the B3LYP/6-31G(d, p) method
Table S1 Possible transfer dimers and estimated VDE values for BTBPD cell The detailed description for molecular dynamics simulation
The detailed description for molecular dynamics simulation

As seen in Eq(14), to estimate the external reorganization energy λext, the optical dielectric constant of medium, εop, need to be firstly calculated by means of Eq(15), that is to say, we firstly need to estimate the material's density ρ. Here, in order to simulate the ρ, Materials Studio 5.5 software package is applied. The whole simulation process can be described as the following steps.

Firstly, the object molecule (BTBPD) was optimized at the B3LYP/6-31G(d) level with the aid of the Gaussian 09 software (Fig.S4), and the stable molecular structure was introduced the Materials studio 5.5 software package. Then, keep the molecular rigid, and apply the Dreiding force field and the Amorphous Cell module to establish a cubic lattice amorphous cell containing 100 molecules (The starting cubic cell parameter was selected to make the density close to 1.0g· cm-3). In the practical calculation, the summation method of atomic charges was used in evaluating the electrostatic and van der Waals interactions, and the quality of energy calculation was selected as fine. And then, the Discover module and smart minimizer method were used to roughly optimize the amorphous cell more than 3000 iterations. In the smart process, all options were selected as the default values. Thus, a relatively reasonable amorphous cell (Ⅰ) was obtained (Fig.S5).

 FIG. S4 Optimized BTBPD molecule at the B3LYP/6-31G(d) level
 FIG. S5 Optimized amorphous cell including 100 BTBPD molecules

Secondly, the amorphous cell (Ⅰ) was introduced into the Forcite Module, and an anneal dynamics calculation was carried out. Here, the cell temperature was firstly increased from 300 K to 600 K, and then was decreased from 600K to 300 K. After the two-step annealing processing including 8 cycles, the unreasonable structures can be efficiently eliminated. Then, the stable cell (Ⅱ) was obtained, which provides a relatively balance geometry for the subsequent molecular dynamics simulation. Fig.S6 showed the energy variation of amorphous cell with the simulated time. Fig.S7 exhibited the temperature change with the simulated time.

 FIG. S6 Correlation of energies and times in annealing simulation
 FIG. S7 Correlation of temperature and times in annealing simulation

Thirdly, in order to achieve balance as soon as possible, a molecular dynamics simulation with 50ps under the condition of constant temperature and volume (NTV) was firstly carried out. In the practical simulation, the Dreiding force field, the current charge, and the Andersen method of temperature control were applied. In addition, the total simulated time was selected to 50ps. Finally, the amorphous cell (Ⅲ) was obtained. Fig.S8 showed the energy variation of amorphous cell with the simulated time. Fig.S9 exhibited the temperature change with the simulated time. As seen, when the simulated time was 5ps, the cell was close to the balance.

 FIG. S8 Energies of amorphous cell vs the simulated times
 FIG. S9 Temperature of amorphous cell vs the simulated times

Fourthly, based on the cell structure obtained by means of the step 3, then a molecular dynamics simulation for 100ps under the condition of constant temperature and pressure (NTP), after 50ps system was balanced, finally, the density of material was obtained (Fig.S11). In the practical simulation, the Dreiding force field, the current charge, and the Berendsen method of temperature-bath coupling were used. In addition, the total simulated time was selected to 100ps. Fig.S10 showed the temperature variation of amorphous cell with the simulated times. Fig.S11 exhibited the material's density change with the simulated times. As seen, when the simulated time was 10ps, the cell trended toward to be a balance structure.

 FIG. S10 Temperature of amorphous cell vs the simulated times
 FIG. S11 Density of amorphous cell vs the simulated times
Reference