Theoretical Study of the Charge Transfer Exciton Binding Energy in Semiconductor Materials for Polymer:Fullerene-Based Bulk Heterojunction Solar Cells.

Recent efforts and progress in polymer solar cell research have boosted the photovoltaic efficiency of the technology. This efficiency depends not only on the device architecture but also on the material properties. Thus, insight into the design of novel semiconductor materials is vital for the advancement of the field. This paper looks from a theoretical viewpoint into two of the factors for the design of semiconductor materials with applications to bulk heterojunction solar cells: the charge transfer exciton binding energy and the nanoscale arrangement of donor and acceptor molecules in blend systems. Being aware that the exciton dissociation of local excitons in charge transfer states initiates the charge generation process, the excited state properties of four oligomers (one donor-type: PEO-PPV; and three donor-acceptor-types: PTFB, PTB7, and PTB7-Th) and two fullerene derivatives ([60]PCBM and [70]PCBM), previously reported in the literature as having high electrical conductance, are studied. With such a study, the donor molecules, either of donor-type or donor-acceptor type, are screened as candidates for [60]PCBM- and/or [70]PCBM-based bulk heterojunctions. The charge transfer energy and charge transfer exciton binding energy of suitable donor:acceptor bulk heterojunctions, some of them not yet fabricated, are studied. Further, the charge transfer exciton binding energies of [60]PCBM- and [70]PCBM-based blends are compared. A combination of molecular dynamics simulations with calculations based on Kohn-Sham density functional theory (KS-DFT) and its time-dependent extension (KS-TDDFT) is used. An important feature of this work is that it incorporates the effect of the environment of the quantum chemical system in KS-DFT or KS-TDDFT calculations through a polarizable discrete reaction field (DRF). Our predictions in terms of the influence of the nanoscale arrangement of donor and acceptor molecules on the performance of organic solar cells indicate that bulk heterojunction morphologies for donor-acceptor-type oligomers lead to their lowest excited states having charge transfer character. Further, we find that in terms of favorable charge transfer exciton binding energy, the PTB7-Th:[70]PCBM blends outperform the other blends.

Introduction

Organic photovoltaics (OPVs) as inexpensive, flexible, and lightweight solar cells have become a promising energy source. There are, however, many issues related to their low efficiencies that have to be addressed before mass production/commercialization; therefore, OPVs constitute an active area of research. In OPVs, the charge generation process involves the formation of excitons created by sunlight absorption. A current is generated if the exciton can be split into a free electron and a free hole. However, exciton dissociation is not easy to achieve, and in many cases, losses occur.[1] The efforts to understand and control the operation of OPVs have led to many device architectures, ranging from a single conductive layer to donor:acceptor (D:A) bulk heterojunctions (BHJs) through D/A bi or multilayer systems. Single layers are the simplest but also the least efficient in separating the exciton.[2] Multilayer junctions, be it stacked D/A films or BHJs, instead combine molecules with different potentials, D (or hole transport) and A (or electron transport) molecules, to overcome the exciton binding energy (Eb). It is believed that in such device architectures, the charge generation occurs through charge transfer (CT) processes from D to A molecules that lead to charge-separated (CS) states. BHJs, as interpenetrating networks of D and A materials dispersed in the bulk, however, have more D/A interfaces, and consequently, have more sites for the CT exciton dissociation, making them more efficient devices.[3]

Conjugated materials with a small band gap, large induced dipole moments, and polarizable fragments are potential candidates for BHJ solar cells.[4] The combination of a semiconducting polymer with a fullerene derivative as an organic blend has up to now been the norm for BHJs. One of the most common BHJs is based on poly(3-hexylthiophene) (P3HT)[5] and the [6,6]-phenyl-C61-butyric acid methyl ester (PCBM)[6] as the D and A molecules, respectively. Nevertheless, morphology disorders mainly associated with P3HT have led to P3HT/PCBM blends yielding low efficiencies.[7] In view of this, many other D and A molecules have emerged. For instance, Torabi et al.[8] recently functionalized conventionally known photovoltaic materials to enhance their dielectric constants (which in principle would reduce both the Eb and losses due to recombination).[9] For that, Torabi et al.[8] attached triethylene glycol (TEG) side-chains to conventional polymers such as diketopyrrolopyrrole (DPP) to phenylenevinylene (PV)-based ones and to fullero-pyrrolidine derivatives. They found that TEG-functionalized polymers and fulleropyrrolidines (PTEG-1 and PTEG-2[10]) have considerably higher dielectric constants than their respective reference polymers and PCBM. PTEG-1 with its high dielectric constant has not yet outperformed PCBM, but the experimental conditions for PTEG-1-based blends have not yet been fully optimized.

Another polymer successfully used in polymer:fullerene solar cells is the poly [[4,8-bis[(2-ethylhexyl)oxy] benzo[1,2-b:4,5-b′] dithiophene-2,6-diyl] [3-fluoro-2-[(2-ethylhexyl)carbonyl] thieno[3,4-b] thiophenediyl]], more commonly known as PTB7. PTB7, in contrast to other donor polymers, broadly absorbs in the near-infrared and has a low optical band gap, which leads to high-performance PTB7-based OPVs.[11,12] The absorption spectrum of PTB7 is further red-shifted when it is thiophene-functionalized as PTB7–Th. The functionalization certainly leads to more efficient devices, but as recently reported by Doumon et al.,[13] it also leads to devices that are less photostable.

The design of novel semiconductor materials for BHJs may be a long and costly process that involves many experiments including synthesis and characterization of the materials with techniques such as scanning force microscopy (to investigate surface structures linked to their electrical properties), current–voltage measurements, and so on. Therefore, it is convenient to run computational simulations to explore the applicability of new materials before synthesizing them. Computational methodologies, such as Kohn–Sham density functional theory (KS-DFT)[14,15] and its time-dependent extension (KS-TDDFT),[16] have proven to be good candidates when studying the electronic structure of photovoltaic materials.[17,18] For instance, the computational work by Few et al.[19] in the modeling of CT state properties at the D/A interface of several thiophene-based polymer:PCBM blends, revealed the influence of the chemical structure on the excitation energies. Calculated spectra of excited states, using KS-TDDFT, showed that hole delocalization in high electronically excited CT states can result in a decreased charge transfer exciton binding energy, ECT-b. Further, they demonstrated that functionalized polymers have a large impact on the degree of CT. Moreover, the TDDFT work by Yi et al.[20] on the electronic couplings and rates of exciton dissociation and charge recombination of pentacene:fullerene heterojunctions (HJs) stressed the role of the intermolecular configurations in such competitive processes. There, the superior performance of bilayer HJs over BHJs was already anticipated. Of course, the reliability of KS-DFT or KS-TDDFT predictions depend on their approximate functionals. For instance, it is well-known that although the generalized gradient approximation (GGA)[21] and hybrid functionals[22] yield good energies and good densities, they have poorly behaved potentials; thus, they underestimate nonlocal contributions.[23] As a consequence, the long-range electron–hole interaction between D and A fragments is underestimated.[24,25] There is, however, a class of corrected functionals, for which the local character of conventional functionals is overcome, namely, the long-range corrected (LC) functionals.[26] LC functionals split the exchange interaction into a long-range part, usually treated with Hartree–Fock (HF), and a short-range part treated by an exchange-correlation functional, usually a GGA functional. In general, when LC functionals instead of conventional functionals are applied to KS-TDDFT, the excited state properties are improved. Therefore, for studying CT states within the framework of the KS-TDDFT, we employ LC functionals.[25]

In the simulations on model systems, the size of the system is another point to take into account. While it is indisputable that large systems imply expensive calculations, sometimes unachievable, it is also true that gas-phase calculations can be misleading. Especially in BHJs, where the domain sizes of the D and the A play a crucial role, ground and excited state properties may be very sensitive to the environment, and partial or total neglect of the environment may lead to different conclusions. Alternatively, multiscale methods, which combine quantum mechanics (QM) and classical mechanics at different levels, may be used.[27]

McMahon et al.[28] studied, through molecular dynamics (MD) simulations and QM calculations, the morphology and electronic structure of P3HT/PCBM blends. They computed the density of states for P3HT chains at different distances from the P3HT/PCBM interface. The results indicated that the quasi-free charge-separated species at the interface are a result of the changes in the electronic structure of P3HT at the P3HT/PCBM interface compared to the electronic structure in the P3HT bulk.

D’Avino et al.[29] studied the exciton dissociation in P3HT/PCBM heterojunctions by combining atomistic MD simulations with QM and classical microelectrostatic calculations, the latter describing the embedding molecules as permanent charges and induced dipoles. They evaluated the energy landscape explored by mobile charges in the vicinity of donor–acceptor interfaces with realistic morphologies. These studies revealed that the exciton binding energy may be overcome by a favorable electrostatic energy landscape of the P3HT/PCBM interface, electronic polarization due to the environment, and interface-induced torsional disorder in P3HT chains.

de Gier et al.[30] demonstrated, through MD simulations and TDDFT calculations within the framework of the discrete reaction field (DRF) method,[31] that the inclusion of side-chains with dipole moments in photovoltaic materials lowers the Eb. Electronic state diagrams, including local excitations as well as CT and CS states for oligothiophene:PCBM BHJs, suggested that the inclusion of polarizable chains is a promising route to improve the efficiency of OPVs. This was further supported by experimental and theoretical work on the influence of permanent dipoles in fullerene derivatives.[32] There, a PCBM analogue with a side-chain containing a permanent dipole, namely, PCBDN, was synthesized and characterized. Complementary TDDFT/DRF calculations predicted the embedding effects on the CT and charge separation processes in close agreement to experiments.

In the present work, KS-DFT and KS-TDDFT are used to study ground and excited state properties, respectively, in single and embedded D/A pairs selected from large D:A BHJs. The DRF method[31] is used to mimic the embedding D:A molecules for a given D/A pair in a given bulk. An important advantage of using DRF is that the properties obtained from a polarizable medium are close to those obtained with full KS-DFT or KS-TDDFT, while the computing time is not substantially increased compared to a vacuum calculation.[31]

Here, a combination of quantum mechanics and polarizable force fields is used to study the electronic structure of a few semiconducting materials, with potential applications in organic photovoltaics. The LC CAM-B3LYP functional is used to study the absorption properties of four oligomers, either D or D–A-type conjugated (poly(ethylene oxide)–polyphenylenevinylene, PEO–PPV (D-type, previously reported as having a high dielectric constant[8]), polythiophenefluorobenzotriazole PTFB (D–A-type, previously reported as a good candidate for nonfullerene-based BHJ solar cells[33]), and PTB7 and PTB7–Th (both D–A-type, previously reported as good candidates for fullerene-based bulk heterojunction solar cells[13])) and two fullerene derivatives with similar absorption properties ([60]PCBM and [70]PCBM) (see Figure ). Next, the excited state properties of all possible donor:acceptor combinations, based on the absorption properties of single films, are theoretically studied. In particular, the CT energy (ECT) and ECT-b of [60]PCBM- and of [70]PCBM-based BHJs are compared. Further, the importance of including the surroundings in the estimation of excited state properties, for which the DRF method has been successfully used, is highlighted. Similarly to D’Avino et al.,[29] the description of the embedding subsystem is given by atomic charges and atomic polarizabilities; although in our model, charges and polarizabilities are placed at all atoms, not only in heavy atoms. Furthermore, atomic polarizabilities are distributed over all the atoms of the molecular mechanics (MM) region, in contrast to a layer-like model where the inner layers are described by the atomic polarizabilities, and the outer layers are described by a single anisotropic polarizability at the center of each MM molecule.[34]

Figure 1

Simplified molecular structures of donor (D) and acceptor (A) materials under study. The C3H7 side-chain at the triazole group on PTFB as well as the C8H17 side-chain at the alkoxycarbonyl and alkoxy (thiophene) groups of PTB7 (PTB7–Th) have been reduced to methyl-chains.

Overall, we investigate the influence of the embedding on the CT and CS states and on the ECT-b, whether simulations predict how morphology might limit the CT, and how all of these factors may guide us in the design of more efficient polymer:fullerene materials for BHJ solar cells. If our model succeeds in the prediction of the charge transfer exciton binding energy and the nanoscale arrangement of donor and acceptor molecules in blend systems, then this scheme may be applied for other OPV materials.

The organization of this paper is as follows. Section summarizes the computational information. Here, the QM and QM/DRF calculations (KS-DFT and KS-TDDFT), including the molecular dynamics simulations, are described. Section discusses the theoretical/computational results. The excited state properties of both single oligomers and tetramer:fullerene derivative BHJ blends together with the effect of their morphology are explained. Extended data are provided in the Supporting Information (SI). Finally, Section presents the conclusions.

Methods

Eight BHJs built from four tetramers (PEO–PPV, PTFB, PTB7, and PTB7–Th) and two fullerene derivatives ([60]PCBM and [70]PCBM), were studied theoretically (see Figure ).

To obtain representative structures to be used in the QM/MM calculations, atomistic MD simulations for neutral tetramer:fullerene BHJs were carried out with the GROMACS package.[35] For both, tetramers and fullerene derivatives, all-atom GROMOS 53A6 topologies[36] were used to simulate the ground state of the D:A BHJ. The topologies of tetramers were generated using an automated topology builder (ATB)[37] (selected geometrical features of the force-field-optimized structures were compared to DFT-optimized structures; see Table S1). Since the GROMOS force fields haven been mainly developed and refined for biomolecules, which generally do not exhibit an extensive cross-conjugation, the obtained oligomer morphologies showed structural variations mainly due to intermonomer torsions. The topologies of the fullerene derivatives were built from an optimized fullerene topology[38] in combination with the topologies of the side-chains generated by ATB. Two different D:A ratios were simulated, 1:1 (20:20/30:30 molecules) and 1:1.5 (10:15/20:30 molecules). The 1:1.5 D:A ratio was only used to mimic PTB7:PCBM and PTB7–Th:PCBM blends, as commonly used in experiment.[13]

D:A BHJs were simulated as follows. To a 30 × 30 × 30 nm oligomer-only box, acceptor molecules were added. Then, the D:A molecules in the box were progressively compressed through a series of 10 short MD simulations in a NPT ensemble. Each MD simulation ran during 100 ps with a 0.001 ps time step, a temperature of 298 K, and pressure of 500 bar (the temperature and pressure were controlled via the Berendsen thermostat and the Berendsen barostat,[39] with relaxation times of 0.1 and 0.5 ps, respectively). The resulting compressed D:A blend was progressively energy equilibrated through a series of eight MD simulations in a NPT ensemble. The pressure in the series ranged from 500 bar, passing by 400, 300, 200, 100, 50, 10, 5, and to finally 1 bar. At this point, each MD simulation ran during 250 ps with a 0.002 ps time step and a temperature of 298 K (as before, the temperature and pressure were controlled via the Berendsen thermostat and the Berendsen barostat, respectively). Box sizes vary depending on the number of D:A molecules and their ratio. In general, equilibrated [70]PCBM-based blends lead to larger box sizes than [60]PCBM sizes. The box sizes range from 3.8 × 3.8 × 3.8 nm, for PEO–PPV:[60]PCBM, to 4.4 × 4.4 × 4.4 nm for PTB7–Th:[70]PCBM (for the energy equilibration plots, see Figures S2 and S3; for a validation of the time scale of the MD simulations, see Table S2 and Figures S4 to S7). It is worth mentioning that these MD simulations are intended to model the spin coating process, for which the time scale should be appropriate rather than optimize thermally equilibrated blends. D/A configurations for QM calculations were selected from the energy equilibrated blends.

For ground state properties such as optimal geometry, ionization potentials (IP, computed as the energy difference between the total energy of the positively charged system and the neutral system), and electron affinities (EA, computed as the energy difference between the total energy of the neutral system and the negatively charged system), KS-DFT was used. For excited states such as local excited states and CT states, KS-TDDFT was used. Both ground state and excited state properties were computed using the long-range corrected CAM-B3LYP functional (with 65% of HF exchange at long-range)[23] with the DZP basis set as implemented in the Amsterdam density functional (ADF) modeling suite.[40,41] Only singlet excited state energies were determined for local excitons (LE) and CT states. The CS energy (ECS) was determined as the difference between the IP and the EA of the D/A pair. Here, a periodic boundary conditions (PBC)-like scheme is used. In this simplified scheme, that mimics the initial stage of charge separation at the interface, it is assumed that a CS state evolves from a given CT state in such a way that electron and hole move away from the active D/A pair to distant D and A molecules with the same conformation in the heterogeneous blend as they have in the CT state. Then, ECT-b is estimated as the difference between the corresponding ECS and ECT.

Embedded clusters were constructed from D/A isolated pairs, taken from MD, with surrounding molecules in a sphere with a radius of 3 nm (see Figure ; for a validation of the MM embedding radius, see Table S3). Embedded calculations were performed by combining either KS-DFT or KS-TDDFT for the active D/A pair and the DRF method, also implemented in the ADF modeling suite. In TDDFT/DRF, linear response theory is used to obtain the first-order change in the density to a time-dependent perturbation. The effective potential is given by the self-consistent field (SCF) potential (formed by the Coulomb, exchange correlation, and DRF potentials) and the external potential. The DRF potential accounts for the QM/MM interactions; it describes the MM region through atomic charges and dipole polarizabilities. The DRF contribution arises from the induced dipoles in the MM part because of the first-order change in the QM charge distribution. Thus, the charges and induced dipoles are obtained self-consistently by solving the DRF linear equations at each SCF iteration.[42] Here, DRF parameters, like atomic charges and atomic polarizabilities for all MM atoms, were obtained from multipole derived charges (MDC-Q)[43] and Thole’s model,[44] respectively.

Figure 2

Illustrative representation of a BHJ blend, from which D/A pairs were selected. (a) Full blend: in blue, the D/A pair; in gray, the embedding D, A molecules. (b) Example of a selected D/A pair.

Results and Discussion

Absorption Properties of Photovoltaic Materials

The mechanism by which excitons dissociate is still unclear. However, it is clear that the nature of D and A materials plays a role in the charge generation process. The optical properties of the single D and A materials shown in Figure were explored. First, the HOMO and LUMO energy levels of each single molecule were calculated. Second, for each D and A molecule, the absorption spectrum was computed. In both cases, the structural dynamics of tetramers of PEO–PPV, PTFB, PTB7, and PTB7–Th and of [60]PCBM and [70]PCBM were simulated by classical trajectories, from which QM geometries were selected. For the orbital energy calculations, reported in Table , BLYP, B3LYP, and CAM-B3LYP functionals with the DZP basis set were used. The reason why three functionals rather than only CAM-B3LYP, used to compute the ECT, were used lies in the orbital energies of virtual orbitals. It is expected that B3LYP and CAM-B3LYP, as a result of the HF exchange contribution, lead to virtual orbitals shifted to higher energies.[45] Nevertheless, HOMO–LUMO trends, as shown in Table , remain valid for the selection of D and A in BHJs.

Table 1

HOMO (H) and LUMO (L) Energy in eV of Isolated Tetramersa and of Isolated Fullerene Derivativesb Calculated with Different Functionals and the DZP Basis Set

	BLYP			B3LYP			CAM-B3LYP
molecule	HOMO (H)	LUMO (L)	L – H	HOMO (H)	LUMO (L)	L – H	HOMO (H)	LUMO (L)	L – H
PEO–PVV	–4.70	–3.09	1.61	–5.60	–2.62	2.98	–7.00	–1.53	5.47
PTFB	–5.01	–3.68	1.33	–5.72	–3.37	2.35	–6.98	–2.46	4.52
PTB7	–4.90	–3.81	1.09	–5.57	–3.54	2.03	–6.75	–2.70	4.05
PTB7–Th	–4.93	–3.86	1.07	–5.59	–3.60	1.99	–6.91	–2.96	3.95
[60]PCBM	–5.97	–4.76	1.21	–6.69	–4.50	2.19	–7.60	–3.39	4.21
[70]PCBM	–6.02	–4.66	1.36	–6.70	–4.43	2.27	–7.48	–3.36	4.12

PEO–PPV, PTFB, PTB7, and PTB7–Th.

[60]PCBM and [70]PCBM.

PEO–PPV has the largest HOMO–LUMO gap, while the lowest is for PTB7–Th, which is very close to that of PTB7. In terms of energy of HOMO and LUMO on D and A, respectively, it can be seen that all D:A combinations seem to fit the requirements for energy level differences in OPVs, i.e., HOMO and LUMO levels on D must be at higher energies than HOMO and LUMO levels on A, respectively.

The KS-TDDFT absorption spectra of D tetramers (Figure ) show that for all the tetramers, the main absorption peaks lie in the visible region, between 2.4 and 2.7 eV (454–519 nm). However, there are clear differences between PEO–PPV, which is D-type, and the other tetramers, which are D–A-type. D–A-type tetramers absorb at lower energies than the PEO–PPV tetramer, which is consistent with their smaller HOMO–LUMO gap. In addition, the backbone of the tetramer determines the absorption more than the side-chains, as suggested by the electronic structure calculations at CAM-B3LYP/DZP level. When comparing PTB7 and PTB7–Th for instance, no significant differences, at least in terms of absorption energies, are found. Important differences might come from their morphology or photostability. However, simulations to investigate this, especially for stability, are outside the scope of this research. When enlarging PTB7 and PTB7–Th tetramers to hexamers or octamers, energy trends get closer to experiments as shown in Table , with measured absorption maximum peaks at 1.85 and 1.77 eV for PTB7 and PTB7–Th, respectively.[13] Computations on infinite chains would lead to improved agreement with experiments but would also require other DFT implementations like periodic DFT or the density-functional-based tight-binding (DFTB[46]) method.

Figure 3

KS-TDDFT (CAM-B3LYP/DZP) absorption spectra of isolated PEO–PPV, PTFB, PTB7, and PTB7–Th tetramers (line broadening as interpolation of excited states via Gaussian broadening; peak width = 0.086 eV).

Table 2

CAM-B3LYP/DZP Local Excitations (LE) (in eV) and Oscillator Strengths (f)a

molecule	LE	f
PEO–PPV	2.73	0.34
PTFB	2.39	1.08
PTB7	2.44t	1.24
	2.15h	0.65
	2.02°	0.61
PTB7–Th	2.45t	1.12
	2.08h	0.62
	2.03°	1.37

For isolated PEO–PPV, PTFB, PTB7, and PTB7–Th oligomers (superscript indices t, h, and °, on PTB7 and PTB7–Th refer to tetramer, hexamer, and octamer, respectively).

With regards to the absorption spectra of [60]PCBM and [70]PCBM, from Figure , it can be seen that absorptions with significant oscillator strengths appear from 3.40 eV onward. [60]PCBM has a peak around 3.77 eV (329 nm) and a broad absorption band between 3.87 and 4.02 eV (320–308 nm) with a maximum at 3.95 eV (314 nm). [70]PCBM has an increased optical absorption in the visible region compared to [60]PCBM. It has two absorption peaks, centered at 3.35 eV (370 nm) and at 3.53 eV (351 nm), in agreement with the experimental trends reported earlier in the literature: the UV–vis spectra of [60]PCBM and [70]PCBM in toluene present main peaks at ∼340 and ∼380 nm, respectively.[47] As a reference, the main absorption peak of [60]PCBM computed using KS-TDDFT (BHandH/DZP) is ∼315 nm.[32]

Figure 4

KS-TDDFT (CAM-B3LYP/DZP) absorption spectra of isolated [60]PCBM and [70]PCBM (line broadening as interpolation of excited states via Gaussian broadening, peak width 0.086 eV).

Charge Transfer Energy and Exciton Binding Energy in BHJs

A study of isolated D and A molecules in terms of excitation energies may guide us in the preselection of photovoltaic molecules; nevertheless, it does not guarantee good performance of BHJ solar cells. There are several parameters that determine the efficiency of BHJs, among those here, the ECT and the ECT-b are considered. In BHJs, a CT state can be the result of a local absorption on the D molecule (tetramer/polymer) followed by an electron transfer from the absorber molecule to a neighboring acceptor molecule at a D/A interface. The energy ECT-b needed to break the attraction between the so-formed electron–hole pair is indicative of efficiency.[9]

Conventionally, OPVs include hole and electron transport layers to drive the generated charges in the active layer toward their respective electrodes. Ideally, molecular dynamics and quantum mechanics simulations should include such transport layers, but in practice, that is computationally unfeasible. Given these difficulties, the molecular dynamics and quantum mechanics simulations were limited to only active layers consisting of D and A molecules. For each blend, from an equilibrated ensemble, several isolated and embedded D/A pairs were selected, for which calculations to determine their ECT, ECS, and ECT-b were performed. First, the influence of the environment on the properties of the tetramer:[60]PCBM blends was evaluated, and then, the performance of [60]PCBM was compared to that of [70]PCBM in their corresponding blends.

Isolated and embedded average properties from a set of 10 D/A pairs for each tetramer:[60]PCBM blend are given in Table . As illustrative examples, four D/A pairs of PTB7/[60]PCBM are shown in Figure . For simplicity, only D/A pairs in vacuum are shown; for ECT, ECS, and ECT-b of all embedded D/A samples, see Tables S4 to S7.

Table 3

CAM-B3LYP/DZP Lowest ECT, ECS, and ECT-b in eV of Different Isolated and Embeddeda Tetramer/PCBM Pairsb

	ECT		ECS		ECT-b
blend	x̅	SD	x̅	SD	x̅	SD
PEO–PPV:[60]PCBM	2.22	0.23	3.77	0.19	1.55	0.04
	2.29	0.37	2.97	0.40	0.68	0.19
PTFB:[60]PCBM	2.26	0.27	3.75	0.17	1.49	0.27
	2.03	0.34	2.87	0.47	0.84	0.22
PTB7:[60]PCBM	1.99	0.16	3.56	0.24	1.57	0.21
	1.88	0.12	2.80	0.24	0.92	0.21
PTB7–Th:[60]PCBM	1.88	0.07	3.39	0.10	1.51	0.06
	1.94	0.15	2.84	0.20	0.90	0.22

In bold.

x̅ stands for average values, and SD stands for standard deviation.

Figure 5

PTB7/[60]PCBM configurations as illustrative examples of D/A conformations, for which ECT, ECS, and Eb were computed.

From the vacuum calculations, it can be seen that PEO–PPV:[60]PCBM and PTFB:[60]PCBM blends have rather similar excited state properties. Likewise, the ECT, ECS, and therefore ECT-b of PTB7:[60]PCBM and PTB7–Th:[60]PCBM blends are rather close to each other. PEO–PPV:[60]PCBM and PTFB:[60]PCBM blends have larger ECT and ECS than PTB7:[60]PCBM and PTB7–Th:[60]PCBM blends; however, both have comparable ECT-b. When comparing results of vacuum and embedded calculations, it follows that in all the cases, the environment stabilizes the CS states much more than the CT states. This is as expected. The ΔECS for PEO–PPV:[60]PCBM is ∼0.8 eV, and for PTFB:[60]PCBM, it is ∼0.9 eV. For PTB7:[60]PCBM, it is ∼0.7 eV, and for PTB7–Th:[60]PCBM, it is ∼0.6 eV. The ΔECT for PEO–PPV:[60]PCBM is ∼0.07 eV, and for PTFB:[60]PCBM, it is ∼0.23 eV. For PTB7:[60]PCBM, it is ∼0.11 eV, and for PTB7–Th:[60]PCBM, it is ∼0.06 eV. Therefore, lower ECT-bs are obtained. In general, PTB7:[60]PCBM and PTB7–Th:[60]PCBM blends exhibit a lower average ECT than those of PEO–PPV:[60]PCBM and PTFB:[60]PCBM blends, while embedded PTB7–Th:[60]PCBM has a slightly larger ECT than PTB7:[60]PCBM.

The DRF energy stabilization to the ECS is further analyzed by a decomposition in contributions from the permanent charge distribution and induced atomic dipoles, in order to reveal the mechanism with which the environment influences the CS states (for selected PTB7/[60]PCBM pairs, see Table S11). In the selected PTB7/[60]PCBM pairs, the QM/MM interaction energy is largely dominated by the polarization energy and the charge-induced dipole interaction term, rather than the electrostatic energy. That is, the change in the charge distribution of the environment as a result of the interaction with the D/A pair (QM system) and other D and A embedding molecules contributes more to the ECS than the Coulombic interaction between the D/A pair and the permanent charge distribution of the environment. In absence of static charges in the DRF region, only QM charge-induced dipole interactions occur. In such cases, the polarization term is comparable to the QM/MM interaction energy because of the induced dipole interactions with the whole system, i.e., because of both QM and MM charges. DRF energies due to only the polarization contribution of selected PTB7/[60]PCBM pairs are listed in Table S12. These results demonstrate that accounting for electrostatic interactions alone omits almost half of the effects of the surroundings. Thus, a rigorous description of embedded excited state properties requires a polarizable force field such as the DRF model.

The CT energies depend also on the following factors. First, the ECT is very dependent on the relative position of the D molecule with respect to the A molecule (few instances of D/A pair configurations are depicted in Figure ; for single ECT, ECS, and ECT-b from different D/A pair configurations, see Tables S4 to S7). Second, the CT state energies depend on the interaction between D and A molecules in the BHJ blend. Even more crucial, the CT state is very sensitive to the proximity between the hole on the D molecule and the conjugated system on the [60]PCBM molecule. Configurations where the local exciton on the D molecule is next to the buckyball, as those where the D molecule wraps the A molecule, lead to lower ECT (for contour plots of the molecular orbitals involved in the lowest CT state of the D/A pair configurations shown in Figure , see Figure S8). Next, for most of the PEO–PPV:[60]PCBM blends, the lowest CT states are higher in energy than the lowest excited states on PEO–PPV (for the excited states manifold, see Tables S4 to S7). This means that excitons on PEO–PPV may decay to other low-lying excited states, such as local states, rather than being transferred to [60]PCBM. Under such conditions, losses due to recombination of electrons and holes are quite likely, a pattern that the polarizable side-chain cannot break. Indeed, it is to be expected that in PEO–PPV-based blends, the predicted CT state is difficult to access. Last, and in relation to the previous point, D–A-type conjugated tetramers combined with [60]PCBM lead to the lowest excited states with a strong CT character. Some blends in the vacuum have high CT states; however, when surrounding molecules are included, the CT states become the lowest states. This suggests that in blends, excitons profit from the environment to quickly reach the D/A interface. It implies, also, that the arrangements of D and A molecules in the blend play a crucial role in the CT process. This shows again that predictions based on vacuum calculations do not sufficiently reflect the physics behind the process of charge generation, and therefore, hereafter, only embedded systems are discussed.

For the D–A tetramers:[70]PCBM-based blends, the same procedure described for [60]PCBM blends was followed. PEO–PPV, despite having the lowest ECT-b across the series, was excluded due to its large HOMO–LUMO gap and high-lying excited states in both pristine states and in the blend, which, as shown in Table S1, indicate that the CT migration is energetically hardly feasible. In contrast to [60]PCBM blends, [70]PCBM blends evolved in MD simulations toward very heterogeneous D and A domains (see Figure ). The simulation revealed that the accessibility of the A molecules is limited by the side-chains of PTB7–Th. Thus, for PTB7–Th-based blends, because of steric effects induced by the thiophene side-chains, the A molecules were surrounded by fewer D molecules than those in [60]PCBM-based blends (as will be discussed later and shown in Figure ). Consequently, the [70]PCBM-based blends lead to more D/A interfaces. From the so-formed D/A interfaces, several embedded D/A pairs were selected, for which ECT and ECS were determined.

Figure 6

MD simulated nanoscale arrangements of PTB7 and PTB7–Th-based blends in a 1:1.5 D:A ratio.

Table lists the average values of ECT, ECS, and ECT-b from a set of 10 embedded D/A pairs (for single ECT, ECS, and ECT-b values from different D/A pair configurations, see Tables S8 to S10, and for contour plots of the molecular orbitals involved in the lowest CT state, see Figure S9).

Table 4

CAM-B3LYP/DZP Lowest ECT, ECS, and ECT-b in eV of Different Embedded Tetramer/[70]PCBM Pairsa

	ECT		ECS		ECT-b
blend	x̅	SD	x̅	SD	x̅	SD
PTFB:[70]PCBM	2.32	0.15	3.21	0.16	0.88	0.16
PTB7:[70]PCBM	2.03	0.15	2.82	0.40	0.79	0.37
PTB7–Th:[70]PCBM	2.02	0.18	2.65	0.33	0.63	0.31

x̅ stands for average values, and SD stands for standard deviation.

The comparison of CT and CS states of [70]PCBM blends (Table ) to those of [60]PCBM blends (Table ) shows that the [60]PCBM energies are slightly lower, especially for PTFB-based blends. The ECSs in particular are close to each other (except for PTFB-based blends). The ECS depends on the IP of the D molecule and the EA of the A molecule. In this case, the EAs of [60]PCBM and [70]PCBM (obtained at the CAM-B3LYP/DZP level) are 3.05 and 3.13 eV, respectively (compared with the gas-phase EA of [60]PCBM of 2.63 eV, measured by low-temperature photoelectron spectroscopy[48]). Furthermore, the ECT-b is expressed as the difference between the ECS and the ECT; thus, the actual difference between [60]PCBM and [70]PCBM blends lies in the CT states. The differences in ECT-b imply differences in the ease of exciton separation.[9] Ultimately, the performance of the blends, whether with [60]PCBM or [70]PCBM, will also depend on the charge diffusion barrier. The morphology of the active layer helps in the charge dissociation and charge transport.[9,49] However, charge diffusion is outside the scope of this work.

For comparison, Figure shows MD simulated arrangements of PTB7 or PTB7–Th blended with fullerene derivatives, [60]PCBM and [70]PCBM, which gives insight into the nanomorphology of the blend layers. The radial distribution function (RDF) of the D molecules with respect to the A molecules at their center of mass (COM) is shown in Figure . This gives an indication of the correlation between D and A domains in the blend. All the RDFs show a broad band between 0.4 and 2.0 nm. At 0.4 nm (4 Å), the density of D and A molecules in the PTB7–Th-based blends is lower than that for PTB7-based blends. At 6 Å, the RDF approaches unity, and the trend remains; however, there is a clear difference between PTB7–Th:[60]PCBM and PTB7–Th:[70]PCBM blends, with the former having a higher D/A density. At larger distances, above 6 Å, we can expect that the CT process is unlikely. PTB7–Th-based blends, which have the lowest D/A density, have more D:A domains with more possible D/A interfaces. 3D pictures would show that the PTB7–Th:[70]PCBM blend has more coupled D:A domains than the other D:A blends, which suggests that in such blends, the CT processes is favored, as is also suggested by experimental work.[13]

Figure 7

Radial distribution functions of the center of mass of the D molecules with respect to the A molecules in PTB7 and PTB7–Th-based blends in a 1:1.5 D:A ratio.

The distribution of ECT across the set of PTFB:[70]PCBM is more homogeneous than that of the set of PTFB:[60]PCBM blends, consistent with the standard deviation that drops from 0.34 eV for [60]PCBM to 0.15 eV for [70]PCBM. This suggests that [70]PCBM is a more favorable acceptor for PTFB than [60]PCBM (see Tables and 4). It would be interesting to set the experimental conditions for its fabrication and see if in agreement with a low ECT, high performance is obtained. Statistically, the ECT-b values of PTB7–Th and PTB7-based blends are close, with PTB7–Th having a lower ECT-b than PTB7, suggesting that the former would slightly outperform the latter in terms of efficiency if the morphology also favors the charge diffusion.

In general, CT states for [60]PCBM-based blends tend to have a complete CT from D to A (∼90% HOMO → LUMO, see Tables S4 to S7), while [70]PCBM-based blends tend to have partial CT character including also partially local excitations on D and A. The lowest-lying excited states of [70]PCBM-based blends are mainly because of transitions between the HOMO on the D molecule and the LUMO on the A molecule; although for some D/A configurations, there are contributions from other transitions between lower-occupied orbitals and higher-unoccupied orbitals (see Tables S8 to S10).

From simulations, it is found that the ECT values of [60]PCBM blends are lower in energy than those for [70]PCBM blends;, however, in terms of ECT-b, [70]PCBM blends lead to weaker electron–hole pairs. As the ECT-b determines the ease of exciton dissociation, one could conclude that [70]PCBM blends, having lower ECT-b, would be more efficient. However, as mentioned above, the estimation of the ECT-b was done through approximations. The large error margin in computed values is a result of the fact that (1) the ECT-b depends on the ECT and the ECS. The ECS in turn depends on the IP of the D molecule and the EA of the A molecule; thus, the errors in IP and EA are propagated. (2) The embedding varies from one D/A pair to another D/A pair. Depending on the configuration of the D/A pair in the bulk, CT and CS states are more or less favored. (3) The computed values might even be closer to experimental values if more repeating units of the polymers were used in the simulations. However, as indicated earlier, calculations on such long chains are computationally very expensive and are not expected to change the observed trends.

To verify our suggestions that (1) PTB7–Th-based blends are more efficient than the remaining tetramer-based blends and (2) [70]PCBM-based blends outperform [60]PCBM blends, the predictions should be complemented by experimental evidence. That evidence would include the determination of the local Eb values for the polymers themselves and the embedded ECT-b values for the blends obtained from measurements on real BHJ solar cells. Work in this direction is currently in progress in collaboration with the Photophysics and Optoelectronics group at the University of Groningen.

Conclusions

By KS-TDDFT studies, particularly with the CAM-B3LYP LC functional, we predicted the ECT-b of tetramer:fullerene derivative BHJ blends. Through QM calculations, we found that the D:A cluster arrangements in the blend influence the exciton dissociation. We demonstrated that the inclusion of many D:A molecules as embedding is fundamental to mimic experimental active layers. Further, we demonstrated that the stabilization of CT and CS states on D/A pairs induced by the embedding can be effectively taken into account by combining KS-DFT/KS-TDDFT with the DRF method. We confirmed that moving from isolated to embedded systems, CS states are much more stabilized than CT states. We observed that the ECT strongly depends on the configuration of the D/A pair, which in turn depends on the interactions between D and A molecules in the BHJ. We showed that the accessibility of the A molecules is limited by the side-chains of the oligomer, thus, influencing the morphology. We infer that despite the predicted ECT-b for PEO–PPV:[60]PCBM blends being the lowest across the series, the CT states of the blends are energetically inaccessible. Our predicted values for ECT and ECT-b values for PTFB with [60]PCBM or [70]PCBM indicate that PTFB:[60]PCBM blends would work better. However, the experimental conditions have neither been set nor optimized. PTB7 and PTB7–Th are structurally quite similar, and from simulations, we only observed more heterogeneous D and A domains for PTB7–Th, mainly because of the thiophene side-chains. Even so, we found that the predicted ECT-b values are lower when these tetramers are combined with [70]PCBM. These results suggest that our modeling of the CT process in BHJ blends may be used to scan the absorption and electrical conductance properties of (novel) semiconductors, being then a guide for further simulations or experiments on the performance of polymer:fullerene-based BHJ solar cells.

As a closing remark, we believe that in the quest of designing novel materials for organic solar cells, polarizable materials as conjugated donor–acceptor copolymers are crucially important not only for fullerene-based cells but also for small molecule acceptor-based devices. As a consequence, for the prediction of the microscopic behavior of organic photovoltaic materials, the inclusion of a polarizable embedding in the quantum mechanical calculations is decisive. A major remaining challenge is understanding the role of the molecular orientation in the charge separation.