Computational Methodologies for Developing Structure–Morphology

Sep 8, 2016 - We outline a step-by-step protocol that incorporates a number of ... the performance properties of bulk-heterojunction organic solar cel...
0 downloads 0 Views 822KB Size
Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Computational Methodologies for Developing Structure-MorphologyPerformance Relationships in Organic Solar Cells: A Protocol Review Khanh Do, Mahesh Kumar Ravva, Tonghui Wang, and Jean-Luc Bredas Chem. Mater., Just Accepted Manuscript • DOI: 10.1021/acs.chemmater.6b03111 • Publication Date (Web): 08 Sep 2016 Downloaded from http://pubs.acs.org on September 8, 2016

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

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

Page 1 of 31

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

Chemistry of Materials

Computational Methodologies for Developing Structure-Morphology-Performance Relationships in Organic Solar Cells: A Protocol Review

Khanh Do,‡ Mahesh Kumar Ravva, Tonghui Wang, and Jean-Luc Brédas*

Laboratory for Computational and Theoretical Chemistry of Advanced Materials Physical Science and Engineering Division King Abdullah University of Science and Technology Thuwal 23955-6900, Kingdom of Saudi Arabia



On leave from: School of Chemistry and Biochemistry & Center of Organic Photonics and Electronics (COPE), Georgia Institute of Technology, Atlanta, Georgia 30332-0400 *

Corresponding author: [email protected]

ACS Paragon Plus Environment

Chemistry of Materials

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

Abstract We outline a step-by-step protocol that incorporates a number of theoretical and computational methodologies to evaluate the structural and electronic properties of π-conjugated semiconducting materials in the condensed phase. Our focus is on methodologies appropriate for the characterization, at the molecular level, of the morphology in blend systems consisting of an electron donor and electron acceptor, of importance for understanding the performance properties of bulk-heterojunction organic solar cells. The protocol is formulated as an introductory manual for investigators who aim to study the bulk-heterojunction morphology in molecular details, thereby facilitating the development of structure-morphology-property relationships when used in tandem with experimental results.

2 ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

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

Chemistry of Materials

I. Introduction Organic solar cells (OSCs) based on solution-processable π-conjugated semiconductors have garnered tremendous interest as a renewable energy-harvesting technology.1 The devices that have achieved the highest power conversion efficiencies (PCEs) exploit a bulk-heterojunction (BHJ) architecture, where the active layer is formed from a blend solution containing an electron donor and an electron acceptor material.2,3 While the adoption of the BHJ architecture has led to much higher efficiencies compared to those of bilayer devices, the task of controlling the active layer morphology has also become significantly more challenging due to the large parameter space involved in the device fabrication process. Many efforts have focused on improving the efficiency of BHJ solar cells by optimizing4 a number of crucial parameters – most of which are aimed at controlling the morphology – including, for example, the functionality and molecular weight5 of the active materials, blend ratio, processing solvent and additives,6 deposition method, and thermal and solvent annealing methods. Morphology control techniques in addition to molecular design strategies have currently led to record efficiencies of 11.7% and 7.7% for polymer:fullerene7 and all-polymer8 solar cells, respectively. A key to further systematic improvements in performance is the ability to develop a molecular understanding of the major factors that affect the morphology and in turn of how the resulting morphology impacts the electronic properties and processes relevant to solar-cell operation. The structure-property relationships should be elucidated for a wide-ranging number of materials such that general understanding and specific guidelines can be deduced for the rational design of new active materials and their processing conditions. The fact that the same processing conditions can lead to stark differences in the morphology – and consequently performance – for independent systems employing only slightly different active materials, underlines the need to 3 ACS Paragon Plus Environment

Chemistry of Materials

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

understand the intricate molecular interactions taking place during the formation of the BHJ active layer. To this end, theoretical and computational methodologies can be very effective given that all atomic coordinates take known values in a calculation or simulation. This helps to minimize ambiguity when making structure-property correlations, especially when used in tandem with experimental morphology and performance characterizations. Moreover, understanding the precise molecular packing behavior in the pure electron donor and electron acceptor phases and their mixed phases is of the utmost importance given that the transfer, dissociation, separation, and transport of excitations and/or charge carriers are all predicated on the nature of these phases and their prevalence in the active layer. Here, we lay out a holistic computational protocol for studying the properties of organic π-conjugated materials at varying length scales ranging from single molecules and molecular complexes to molecular aggregates and nanoscale/mesoscale domains. While it is by no means a comprehensive overview of the advances made in the area of organic solar cells, we hope that this contribution, built essentially on our experience, will serve as a pedagogical manual for those interested in investigating the BHJ morphology of organic solar cells in molecular details, via computational means.

II. Computational Methodologies In the following Sections, we will cover a variety of computational methodologies (Figure 1) for studying the structural and electronic properties of π-conjugated materials in the condensed phase, in particular, as related to bulk heterojunctions. These methodologies include: (I) quantum mechanical (QM) calculations; (II) molecular dynamics (MD) and (III) coarse-grained (CG) 4 ACS Paragon Plus Environment

Page 4 of 31

Page 5 of 31

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

Chemistry of Materials

simulations; and (IV) kinetic Monte Carlo (KMC) simulations. Specifically, we will discuss how these methodologies can be integrated to build a multiscale picture of the morphology and performance properties of the active materials in BHJ solar cells. For this purpose, we have organized the methodologies into a step-by-step protocol. Figure 1 shows the ordering of the methodologies, including the intermediate steps to process the inputs and outputs for each of them. Our goal here is not to give detailed technical instructions on how to perform such calculations or simulations. Instead, we aim to give a “hitchhiker’s guide” to the computational techniques available for studying the materials and systems of interest to BHJ solar cells. In each step, if applicable, we describe our rationale, personal insights, and recommendations for carrying out the procedure. We will primarily present the methodologies in the order they are depicted in Figure 1, starting with QM calculations. For convenience, we will generally use the term “molecule” to refer to the individual unit of the active materials regardless of whether they are actually small molecules, oligomers, or polymers.

5 ACS Paragon Plus Environment

Chemistry of Materials

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

Figure 1. Flow chart depicting the computational methodologies discussed in this contribution. The numbering of the methodologies indicates the possible order of implementation when they are integrated to form a comprehensive, multiscale study of the BHJ morphology. The notations correspond to: QM, quantum mechanical calculations; MD and CG, molecular-dynamics and coarse-grained simulations; and KMC, kinetic Monte Carlo simulations.

Step 1. Quantum mechanical (QM) calculations The most basic properties to evaluate for the active materials are those pertaining to their electronic structures, which can be determined by performing QM calculations, the results of which include the wavefunctions of relevant molecular orbitals such as HOMO and LUMO 6 ACS Paragon Plus Environment

Page 6 of 31

Page 7 of 31

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

Chemistry of Materials

(highest occupied and lowest unoccupied molecular orbital, respectively), the ionization potential (IP) and electron affinity (EA) values, and excited-state energies. These properties are of interest, because the energy gap between the donor IP (often incorrectly referred to in the literature as the donor “HOMO energy”) and acceptor EA (“LUMO energy”) has been shown empirically to correlate with the device open-circuit voltage (VOC).9 This correlation has been made more precisely for the charge-transfer (CT) state energy,1 which corresponds to an excited state at the interface between the donor and acceptor components, with a hole on a donor molecule and an electron on an adjacent acceptor.10 The difference in donor IP and acceptor EA can give an estimation of the energy of the final charge-separated (CS) state, where the hole resides in a donor domain and the electron in an acceptor domain without being any longer Coulombically bound.10 Moreover, given that the solar cell operation depends largely on an energetic cascade within the various materials used in the device (the mixed donor-acceptor phases, the pure hole and electron transport phases, and the electrodes), it is desirable to gain a molecular understanding of the complete energetic landscape. When performing a QM calculation, the chosen methods generally fall into three classes, which we list in descending order of the accuracy that can typically be expected: (i) high-level (correlated) wavefunction-based methods, (ii) density functional theory (DFT), and (iii) semiempirical methods. The high-level methods include post-Hartree-Fock methodologies such as configuration interaction (CI) and Møller-Plesset perturbation theory (MPPT), which are expected to have high accuracy due to their exact treatment of exchange and extensive evaluation of electron correlation. In DFT, exchange and correlation effects are taken into account approximately as a functional of the electron density; modern DFT functionals can capture exchange-correlation effects at a level of accuracy comparable to that of the high-level 7 ACS Paragon Plus Environment

Chemistry of Materials

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

wavefunction methods, generally at a fraction of the computational cost. The semi-empirical methods are even more computationally inexpensive, because they largely simplify the evaluation of the electron-electron interactions while recovering some level of accuracy and including correlation effects via incorporation of empirical parameters. In terms of computational requirements, the classes of QM methods as listed above thus come in ascending order of speed. With optimal combination of computational cost and accuracy in mind, along with the prospect of studying the electronic-structure properties of π-conjugated systems that are a few hundred atoms large, DFT has very much become the QM method of choice. A word of caution is, however, in order at this stage. While the B3LYP functional has been commonly chosen, it has critical shortcomings that should prevent its use to describe π-conjugated systems especially in the context of organic solar cells. Indeed hybrid functionals such as B3LYP tend to overdelocalize the electron density and, as a result, have been shown to overestimate torsion potential barriers in extended π-conjugated systems by over-stabilizing planar conformations;11–14 we note that this delocalization error stems from the electron self-interaction error.13 This same issue should prevent the use of B3LYP and similar functionals when evaluating the charge-transfer states we discuss below. For these reasons, it is strongly recommended using a long-range corrected (LC) functional such as ωB97XD15 for studying systems with extended π-conjugation pathways. In using a LC functional, the delocalization error is minimized, for instance, via “IPtuning” of the range-separation parameter ω, such that the HOMO energy of the molecular system matches its IP, which is a property obeyed by the exact functional.16 Analogously, the ωvalue can also be optimized by “gap-tuning” where in addition the LUMO energy is parameterized to the EA.

8 ACS Paragon Plus Environment

Page 8 of 31

Page 9 of 31

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

Chemistry of Materials

When the total number of systems to be evaluated becomes excessive, semi-empirical methods can become an appropriate choice to make the task computationally feasible. In the context of πconjugated systems, the class of methods involving the intermediate neglect of differential overlap (INDO) Hamiltonian and its variants, ZINDO and ZINDO/S, has proven to be useful. In particular, these methods give reasonable values for the electronic couplings between adjacent molecules, in comparison with DFT-based methods,17–19 and have been successfully exploited to calculate the rates of exciton dissociation and charge recombination in prototypical donoracceptor systems.20,21 Regardless of which methodology one chooses, what is needed to perform a QM calculation is simply the coordinates of the molecular system, which can be generated manually with generic or specific values for the bond lengths, angles, and dihedrals or obtained from a crystal structure, when available. The equilibrium geometry of the system can be predicted by performing geometry optimizations. The geometry-optimized structure can then be used to determine the electronic-structure properties of the molecular system. In the next step, we will show how to use the results of a QM calculation for a single molecule in vacuo to parameterize the force field for performing MD simulations. Beyond a single molecule, the electronic properties of dimers or molecular complexes are often relevant (we note terminology-wise that in a dimer the two molecules are the same, i.e., donordonor or acceptor-acceptor, while in a complex, they are different, i.e., donor-acceptor). The molecular configurations chosen to study dimers or complexes can in fact come, in the context of OSC characterization, from the results of an atomistic MD simulation or a CG simulation followed by a reverse-mapping procedure (Figure 1), which we will discuss below. The configurations can also be manually assigned with particular coordinates such that systematic 9 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 10 of 31

effects of a particular structural feature, e.g. molecular backbone stacking distance, can be studied.22 In any case, an atomistic representation of the molecular packing configuration is necessary. With that in hand, it becomes possible, for instance, to evaluate the charge-transport properties by calculating the electronic couplings for holes and electrons (th and te, respectively) and the reorganization energies (λ); these parameters are inputs to the semiclassical Marcus equation for determining the rates of charge hopping between neighboring molecules (discussed in Step 7):  + Δ° 2  1 || = exp −  ℏ 4 B 

4 B 

1

Here, ∆G° is the change in free energy, kB is the Boltzmann constant, and T is the absolute temperature. While the electronic couplings between HOMOs [LUMOs] on adjacent identical molecules are relevant to hole [electron] transport via hopping,23 other electronic couplings are useful to evaluate as well. For instance, the electronic coupling between the electron-donor LUMO [HOMO] and the electron-acceptor LUMO is relevant to exciton dissociation [charge recombination].20,24 To determine electronic couplings, the accurate method of Valeev et al. is highly recommended,25 as it accounts for polarization effects and the possibility that the molecular orbitals are non-orthogonal. The reorganization energy λ is typically divided into two parts, the external (outer) reorganization energy (λext) due to polarization effects of the medium and the internal (inner) reorganization energy (λint) due to structural relaxation upon charge transfer between the molecules involved. For the internal reorganization energies, λh [λe] for holes [electrons] can be calculated by comparing the energies of the neutral and cationic [anionic] states at the optimized geometries of the neutral and cationic [anionic] states:26 10 ACS Paragon Plus Environment

Page 11 of 31

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

Chemistry of Materials

 =  −   +  −   = ! − !!  + ! −  

2 3

where  [! ] is the energy of the cation [anion] calculated at the optimized geometry of the neutral molecule; likewise,  [!! ] is the energy of the cation [anion] calculated at the

optimized cation [anion] geometry;  [! ] denotes the energy of the neutral molecule

calculated at the cation [anion] geometry, while  is the energy of the neutral molecule at the ground-state geometry.

We now turn our attention specifically to a donor-acceptor molecular complex. Here, it is especially important to investigate the characteristics of the charge-transfer (CT) states, as the formation of separated (free) holes and electrons in OSCs are understood to proceed through these intermediate states.27 While CT state energies can be accurately evaluated using a number of high-level computational methodologies, including correlated excited state (ADC(n) family), coupled-cluster (CC) or CI, these methodologies become quickly computationally prohibitive due to the large system sizes that need to be considered (it must be borne in mind that the most common electron acceptor, a fullerene derivative, contains over sixty carbon atoms). Various other techniques such as the Maximum Overlap Method (MOM),28 Restricted Open-Shell KohnSham (ROKS),29,30 and a modified linear response theory for time-dependent density functional theory31 (TD-DFT) can also give good estimations of CT state energies. Among these techniques, TD-DFT based on a long-range corrected functional (and not B3LYP or similar hybrid functionals, see our caveat above) can be recommended, as this method has been successfully used to study CT-state properties and comes standard in most quantum chemical software packages.13 The nature of the CT states can then be characterized using natural transition orbitals (NTOs) derived from the excited states.32 In a typical CT state, the hole-NTO 11 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 12 of 31

is essentially localized on the donor and the electron-NTO is essentially on the acceptor.33 The spatial overlap between the centroid of the hole-NTO and electron-NTO is often used as a metric to describe the Coulombic interaction strength between the hole and electron.34,35 We note that CT state energies are largely dependent on the molecular packing configuration, i.e., the intermolecular distances, and relative orientations of the donor and acceptor in the molecular complex.33,36 Another word of caution here: In the case of extended π-conjugated systems, we strongly advise against the simple consideration of HOMO and LUMO wavefunctions to provide (even a crude) idea of the nature of the lowest excitation, be it intramolecular or intermolecular; the consideration of NTOs is necessary if a reliable picture is to be gained.37 In many instances, the stabilization (binding) energy due to the packing (interaction) of two or several molecules is a useful parameter to determine. Higher stabilization energies for specific packing configurations indicate a higher prevalence for these to exist in the active layer and point to higher thermal stability for the phase. Generally, interaction energies are often evaluated as the difference between the energies of the isolated fragments and their molecular complex. For example, in the case of a donor-acceptor complex, the interaction energy can be calculated using the following equation: int = DA − D + A 

4

where ED and EA are the total energies of the isolated donor and acceptor, respectively, and EDA is the total energy of the complex. As the stabilization energy due to molecular packing in πconjugated systems is governed by weak van der Waals interactions, it is necessary to use a density functional that can treat dispersion interactions accurately when performing interaction energy calculations. For that purpose, it can be recommended once again to use the LC functional, ωB97XD,15 which includes empirical corrections to account for dispersion 12 ACS Paragon Plus Environment

Page 13 of 31

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

Chemistry of Materials

interactions. Also, it is important to eliminate basis set superposition errors (BSSEs), for instance, by applying the counterpoise correction method proposed by Boys and Bernardi.38 Without such correction, BSSEs lead to an artificial over-stabilization of the binding energy, which is a result of each molecule in the complex having access to the additional basis functions of its neighbor; thus, its effect varies as a function of intermolecular distance and basis set quality. Obtaining accurate binding energies therefore requires to eliminate the inconsistent BSSEs. In addition to interaction energy calculations, an energy decomposition analysis can be performed using symmetry-adapted perturbation theory (SAPT) in order to obtain a breakdown of the various types of interactions into physically meaningful components, that is, the electrostatic, induction, dispersion, and exchange-repulsion terms that contribute to the overall interaction energy.39 Several levels of SAPT methodologies are available, with SAPT0 being the simplest and least computationally expensive truncation of SAPT.40,41 The implementation of SAPT0 is included in the PSI4 package,42 where system sizes of over 200 carbon-like atoms can be reasonably investigated. We note that the interaction energies obtained from an energy decomposition analysis do not suffer from basis set superposition errors.

Step 2. Force field parameterization for MD simulations Here, we first show how to modify the parameters of a well-established force field to account for the specific molecular properties of the materials of interest. Since our aim is to study the molecular packing configurations of organic materials in the bulk phase, the OPLS-AA (all-atom optimized potentials for liquid simulations) force field developed by Jorgensen and co-workers 13 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 14 of 31

appears as a method of choice.43 This force field has shown good reproducibility for bulk-phase properties such as density, compressibility, and heat capacity for a variety of small organic molecules. However, the atomic partial charges from this force field are empirically derived and generalized, as are the bond lengths and angles. The description of dihedrals is from lower-level Hartree-Fock methods or DFT calculations at the B3LYP level. In light of these features, we recommend using the OPLS-AA force field as a starting point and update as many parameters as practically feasible using state-of-the-art DFT methods, in the way we detail below. We note that, in principle, any force field could work as a starting point so long as the resultant force field after parameterization can reproduce the relevant condensed-phase properties – e.g., the material density is a simple parameter to reproduce. In the procedure to update the OPLS-AA parameters for π-conjugated systems, we can take the values for bonds and angles from the geometry-optimized molecules as determined from longrange corrected DFT calculations. For simplicity, the stock values of the force constants can be retained or derived by performing constrained geometry optimizations with fixed bond lengths and angles at the DFT level; those energies can then be fit to a curve to obtain a potential. Likewise, deriving new dihedral parameters requires performing constrained geometry optimizations with fixed dihedral angles. It is highly recommended to obtain new parameters for the dihedrals corresponding to inter-ring or inter-monomer bonds, as they dictate the conformation and thus the degree of conjugation along the π-conjugated backbones. Again, we emphasize that such constrained geometry optimizations should be performed using a LC functional such as ωB97XD in order to generate accurate torsion barriers and local minima, which can be completely missed when using B3LYP.44 The dihedral parameters are ultimately obtained by performing an analogous sampling of the dihedral angles at the force-field level 14 ACS Paragon Plus Environment

Page 15 of 31

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

Chemistry of Materials

while tuning the parameters to reproduce the LC-DFT torsion profile. This is usually done by performing molecular mechanics (MM) calculations using the same constrained-optimized geometries at the various dihedral angles from the DFT calculations; however, this method is lacking in two ways. First, the dihedral parameters obtained in this way correspond to reproducing the DFT torsion at absolute zero only, whereas the interest is in studying the system at elevated temperatures (room and melt temperatures). Secondly, using only one configuration for each dihedral angle neglects other possible molecular configurations with the same angle. To address these two issues, it is recommended to derive the dihedral parameters by performing an MD simulation of a single molecule or polymer chain in vacuo over a range of temperatures, e.g., 300-600 K. The parameters should be tuned so that the potential of mean force ()* +

corresponding to the dihedral angle distribution ,-. - + matches the LC-DFT torsion profile /*01 + within the temperature range:

()* + ≈ −3 ln5,-. - +6 ≈ /*01 +

5

An example describing the implementation of this method is shown in Ref. 45.45 Finally, the partial atomic charges should be updated to reflect the unique electron density distribution of the molecule. Charges derived from the Mulliken population analysis have been shown to be sensitive to the basis set used. Therefore, we recommend deriving charges based on reproducing the molecular electrostatic potential, which is the basis of the ESP46,47 and RESP48,49 methods, among others. If the material of interest is a polymer, the charges on the terminal monomers will have to be modified slightly to ensure total charge neutrality for each chain. Similarly, the charges may need to be modified to enforce symmetry conditions by taking the average of two similar values.

15 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 16 of 31

With the updated force field in hand, we are now ready to study the local structure in the bulk phase by performing atomistic MD simulations. We describe how to do this in the next Step.

Step 3. Atomistic molecular dynamics (MD) simulations As we mentioned earlier, MD simulations are necessary for studying the molecular packing behavior in large systems that are beyond the feasibility of QM calculations. These systems can include molecules or polymer chains in vacuum,50 solution,50–52 crystal,53–55 bilayer,56 or melt of the pure19,55,57 or blend58,59 active materials. Atomistic MD simulations can typically handle tens of thousands of atoms. The total number of molecules to include in a simulation should reflect this capacity. For condensed-phase simulations, one should build a system that is sufficiently large to minimize the artificial effect of having periodic boundary conditions. For instance, an appropriate criterion for choosing the simulation box size would be to ensure that intermolecular structural correlations do not extend further than half the box length, i.e., choose a box length that is at least twice the correlation length. For bulk systems, the correlation length can be determined by measuring the radial distribution function, g(r), and identifying where g(r) approaches unity. For polymers in dilute solution, the correlation length is roughly equal to the radius of gyration. To begin a simulation, an initial configuration must be generated, which can be random if studying the amorphous or solution phase or can correspond to a super cell of the crystal structure if studying the crystalline phase. In either case, periodic boundary conditions should be enforced when considering the condensed phase. For simulations in vacuo, a simulation box much larger than the system must be used and periodic boundary conditions turned off. For an 16 ACS Paragon Plus Environment

Page 17 of 31

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

Chemistry of Materials

amorphous/solution system starting from a random initial configuration, it may be necessary to perform a short simulation with a “soft” nonbonded interaction potential to displace overlapping atoms that stem from placing them randomly inside the simulation box. A “soft” potential is one that reaches a small finite value as the distance between two atoms approaches zero. In building a crystalline system, there are usually no overlapping atoms due to the high order of the system, so this step is not required. Next, the system needs to adjust to the temperature at which it is simulated. This is accomplished by performing a constant volume and temperature (NVT ensemble) simulation using a small timestep (0.1 fs) for about 50,000 time steps (such a very small timestep allows the system to slowly enter the ensemble). Then, the procedure must be repeated using the timestep that will be used for subsequent simulations (1-2 fs). We are now ready to perform the constant pressure and temperature (NPT ensemble) simulation to sample the packing configurations, for which we recommend the Nosé-Hoover barostat60 and Nosé-Hoover thermostat, which can produce the correct fluctuations and trajectories consistent with the canonical ensemble.61 Note that the Parinello-Rahman barostat62 will have to be used for simulating a crystalline system that has a triclinic simulation box. How long the simulation should run depends on the type of system being studied. For crystalline systems at low or room temperatures, a few hundred picoseconds should suffice to completely sample that configuration space. For amorphous systems in the melt, on the other hand, longer simulations are usually required. The guiding figure we use to determine how long to run the simulation is the relaxation time 8 of the backbone of the molecule or polymer. The relaxation time is obtained by fitting the autocorrelation function of a unit vector u (t ) between the ends of a chain or molecular backbone to the equation:

17 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 18 of 31

, 59 ⋅ 906 ≈ exp−/8  where , = ≡ 3=  − 1/2 is the second-order Legendre polynomial. Higher temperatures lead to shorter relaxation times. Simulations should run for around 6-10 times 8 and statistics should

be sampled from the latter part of that period where the structural distributions no longer vary systematically with time. Some software packages recommended for performing the MD simulations discussed here are LAMMPS,63 GROMACS,64 and Materials Studio.65 The output of an MD simulation is a trajectory of all the atoms. With these atomic coordinates, the dynamic and structural properties of the system can be analyzed. A list of the most basic values and functions that can be computed from the trajectory (along with references showing examples for π-conjugated systems) includes: density45; mean square displacement (MSD)58; radial distribution function (RDF)58; orientational correlation function (OCF)66; structural distribution functions45,56,57; or mechanical properties.67,68 In the instances where the sizes of the systems are very large, such as those containing more than 50,000 atoms, it can be desirable to reduce the degrees of freedom in the system by moving to coarse-grained simulations. In the next Step, we show how to parameterize the CG potentials using the structural distributions of the atomistic MD simulations.

Step 4. Parameterization of the CG force field In coarse-grained simulations, we simplify the system by turning groups of atoms into “super atoms,” taking their centers-of-mass as the new coordinates. Then, we derive a set of potentials such that when we apply them to the super atoms, the resulting structural distributions match

18 ACS Paragon Plus Environment

Page 19 of 31

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

Chemistry of Materials

those from the atomistic MD simulations. The following instructions are based on the coarsegraining procedure detailed in the work of Huang et al.45 The CG potentials are obtained iteratively following Equation 6: ,? = /?@ = = /? = + A? 3 ln B G ,CDEF C =

6

Here, /? = is the potential (bond, angle, dihedral, nonbonded, etc.) of the current iteration

resulting in the structural distribution ,? =. We obtain an updated potential /?@ = by comparing ,? = to the target atomistic distribution ,CDEF C =. Convergence is reached when the

two distributions match. The parameter A? is used to control the rate of convergence, where a

value between 0 and 1 is reasonable. The initial function to start the iteration cycle is shown in Equation 7: / = = −3 lnI,CDEF C =J

7

The optimization process for the CG potentials should be carried out in the NVT ensemble where the volume of the system is fixed. This allows the potentials to converge relative to each other without the need to reproduce the volume. Once all potentials have converged, the nonbonded interactions between CG sites j and k are scaled linearly following Equation 8 so the volume (and hence density) of the system is reproduced when performing the simulation in the NPT ensemble: LMN O = PMN QR −

O

OSTU

V,

O ≤ OSTU

8

The coarse-graining procedure described above is based on reproducing the atomistic structural distributions, which can come from various sources such as a system containing a molecule or 19 ACS Paragon Plus Environment

Chemistry of Materials

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

polymer chain in vacuum,69,70 solution,71,72 or the melt.45,58 Similar coarse-graining procedures are reported, for example, in Refs. 69 and 70. In contrast, the methodology of Voth et al. derives the CG potentials by matching the force data from the atomistic simulations.73,74 When constructing the CG sites, it is imperative not to integrate out the relevant degrees of freedom that are involved in the structural features of interest. For example, the degree of interdigitation of acceptor molecules between the alkyl side chains of donor molecules in the BHJ morphology is a crucial parameter influencing solar cell performance. Therefore, sufficient representation of the donor alkyl chains must be retained in the selection of CG sites. Lastly, we would like to mention the effect of coarse-graining on the simulation timescale. CG simulations typically lead to significantly faster dynamics as compared to the corresponding atomistic MD simulations for the same system, which is a result of having drastically reduced degrees of freedom. One way to relate the CG and atomistic timescales is to measure the mean square displacement and obtain the diffusion constant over a temperature range.

Step 5. Coarse-grained (CG) simulations The execution of CG simulations is analogous to that of atomistic MD simulations described above – the only difference is now we are using CG particles whose motion is governed by CG potentials. The same systems that were used during the parameterization procedure can now be studied at much larger length and time scales (tens of nanometers and hundreds of nanoseconds, respectively), which allows for the observation of nanoscale behavior such as the phase separation of donor and acceptor domains,58,69,70,75 self-assembly of polymer nanostructures in

20 ACS Paragon Plus Environment

Page 20 of 31

Page 21 of 31

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

Chemistry of Materials

solution,71,72,76 and short-range polymer crystallization.77 Additionally, CG simulations can be used to predict the mechanical properties of organic semiconductors.78

Step 6. Perform reverse-mapping to obtain the atomistic system The results of CG simulations can yield realistic nanoscale morphologies for which the processes of charge-carriers dynamics including generation, transport, and recombination can be assessed. However, a reverse-mapping procedure is first required to relate the CG morphology back to its atomistic representation, where the resulting molecular configurations can be evaluated via QM calculations to obtain the necessary electronic properties that can be fed into charge-transport models (Step 7). While there is no rigorously unique way to fine grain the CG system, there may be constraints that need to be taken into consideration, such as energy minimization of the atomistic system via the Metropolis algorithm69 or maintaining regioregularity of the polymer through geometrical constraints.77 In any case, the resulting atomistic representations will serve as realistic systems, having been derived from long-range correlations in CG simulations, to relate morphology to performance properties relevant to photovoltaic operation.

Step 7. Perform kinetic Monte Carlo (KMC) simulations Kinetic Monte Carlo (KMC) simulations allow the investigation of the dynamics of excitons and charge carriers in the simulated BHJ morphology, which can consist of a Cartesian lattice or an atomistic-resolution morphology, such as one obtained from reverse-mapping the CG nanoscale morphology. In the latter case, the relevant processes in solar cell operation can be directly

21 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 22 of 31

linked back to the morphology, and in turn, to its contributing factors such as chemical functionalities, molecular weights, blending ratio of the active materials, and thermal annealing temperatures. In KMC simulations, the dynamics of the electronic species proceed through a queue of processes or events such as recombination and hopping between adjacent molecular sites, each of which has a waiting time 8 that is determined by the associated rate k and a random number X between 0 and 1:79 8=−

lnZ 

9

The behavior of the particles of interest then follows the event with the shortest waiting time 8. The KMC simulation integrates over a long chain of events for many particles similarly to MD simulations where the integration is over a large number of time steps. Average quantities of interest can then be calculated from the equilibrated or stabilized trajectory of particles. The processes that are most relevant to organic solar cell operation and that can be readily modelled with KMC simulations are exciton transport, exciton dissociation, charge recombination, charge transport, and charge collection, each with its own rate. We refer to the recent review by Groves for a detailed description of how to implement a KMC description of these processes.79 While previous KMC simulation studies have investigated geminate and bimolecular recombination and charge transport largely using lattice-based morphologies,79–83 it is ultimately necessary to consider molecular morphologies where the precise molecular packing configurations dictate the rates involved in the electronic events. To the best of our knowledge, the recent work of Jones et al.77 is a first demonstration relating the molecular morphology of a -conjugated polymer system [poly(3-hexylthiophene) (P3HT)] to the hole mobility, which was 22 ACS Paragon Plus Environment

Page 23 of 31

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

Chemistry of Materials

approximated through semi-empirical quantum-chemical evaluations of the charge-transfer rates between molecular sites following the semiclassical Marcus equation (Equation 1). We note that a similar procedure has also been performed to calculate the electron mobilities in amorphous and crystalline systems of the fullerene acceptor, [6,6]-phenyl-C61-butyric acid methyl ester (PCBM).84 However, in that case, the morphologies were derived from systems which were constructed to be ordered (crystal structure configuration) or disordered (random initial configuration) rather than coming from self-assembly (crystallization) due to long-range correlations as considered by Jones et al.77 More importantly, the work of Jones et al. was able to elucidate the molecular mechanism responsible for the observed higher hole mobility in P3HT as a function of increased annealing temperature and molecular weight. The next step forward would be to study the charge dissociation, recombination, and transport processes in donoracceptor systems that are derived from solution. In this way, the blend morphologies can be as realistic as possible.

III. Conclusion In this contribution, we have surveyed a number of theoretical and computational methodologies for evaluating the structural and electronic properties of π-conjugated semiconductors in the condensed phase. The methodologies include: density functional theory (DFT), atomistic and coarse-grained (CG) molecular dynamics (MD), and kinetic Monte Carlo (KMC), which range over a multitude of length and time scales. We described how these methodologies can be combined to form a step-by-step protocol for studying the molecular morphology of various thinfilm phases and its impact on relevant parameters. In particular, we focused on blend systems

23 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 24 of 31

consisting of an electron donor and an acceptor, which is important for understanding the performance properties of bulk-heterojunction (BHJ) organic solar cells (OSCs). We hope this overview can serve as an introductory manual for performing simulations to elucidate the molecular packing configurations among the component materials in BHJ OSCs, thereby facilitating the development of structure-morphology-property relationships when used in tandem with experimental results. We expect that in the foreseeable future, systematic relationships between the chemical functionality of π-conjugated materials and their photovoltaic performance properties can be developed directly from theoretical and computational methodologies.

Acknowledgements This work has been supported by King Abdullah University of Science and Technology (KAUST), the KAUST Competitive Research Grant Program, and the Office of Naval Research Global (Award N62909-15-1-2003). The authors acknowledge the IT Research Computing Team and Supercomputing Laboratory at KAUST for providing computational and storage resources and thank Dr. Sean M. Ryno for stimulating discussions.

24 ACS Paragon Plus Environment

Page 25 of 31

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

Chemistry of Materials

References (1) (2)

(3)

(4)

(5)

(6)

(7) (8) (9)

(10) (11)

(12)

(13)

(14)

(15)

Tress, W. Organic Solar Cells; Springer Series in Materials Science; Springer International Publishing: Cham, 2014; Vol. 208. Halls, J. J. M.; Walsh, C. A.; Greenham, N. C.; Marseglia, E. A.; Friend, R. H.; Moratti, S. C.; Holmes, A. B. Efficient Photodiodes from Interpenetrating Polymer Networks. Nature 1995, 376 (6540), 498–500. Yu, G.; Gao, J.; Hummelen, J. C.; Wudl, F.; Heeger, A. J. Polymer Photovoltaic Cells: Enhanced Efficiencies via a Network of Internal Donor-Acceptor Heterojunctions. Science 1995, 270 (5243), 1789–1791. Neil D. Treat; Paul Westacott; Natalie Stingelin. Organic Semiconductors: Manipulation and Control of the Microstructure of Active Layers. In The WSPC Reference on Organic Electronics: Organic Semiconductors; Materials and Energy; WORLD SCIENTIFIC, 2016; pp 159–193. Brinkmann, M.; Rannou, P. Molecular Weight Dependence of Chain Packing and Semicrystalline Structure in Oriented Films of Regioregular Poly(3-Hexylthiophene) Revealed by High-Resolution Transmission Electron Microscopy. Macromolecules 2009, 42 (4), 1125–1130. Moon, J. S.; Takacs, C. J.; Cho, S.; Coffin, R. C.; Kim, H.; Bazan, G. C.; Heeger, A. J. Effect of Processing Additive on the Nanomorphology of a Bulk Heterojunction Material. Nano Lett. 2010, 10 (10), 4005–4008. Zhao, J.; Li, Y.; Yang, G.; Jiang, K.; Lin, H.; Ade, H.; Ma, W.; Yan, H. Efficient Organic Solar Cells Processed from Hydrocarbon Solvents. Nat. Energy 2016, 1 (2), 15027. Hwang, Y.-J.; Courtright, B. A. E.; Ferreira, A. S.; Tolbert, S. H.; Jenekhe, S. A. 7.7% Efficient All-Polymer Solar Cells. Adv. Mater. 2015, 27 (31), 4578–4584. Scharber, M. C.; Mühlbacher, D.; Koppe, M.; Denk, P.; Waldauf, C.; Heeger, A. J.; Brabec, C. J. Design Rules for Donors in Bulk-Heterojunction Solar Cells—Towards 10 % Energy-Conversion Efficiency. Adv. Mater. 2006, 18 (6), 789–794. Brédas, J.-L.; Norton, J. E.; Cornil, J.; Coropceanu, V. Molecular Understanding of Organic Solar Cells: The Challenges. Acc. Chem. Res. 2009, 42 (11), 1691–1699. Karpfen, A.; Choi, C. H.; Kertesz, M. Single-Bond Torsional Potentials in Conjugated Systems:  A Comparison of Ab Initio and Density Functional Results. J. Phys. Chem. A 1997, 101 (40), 7426–7433. Sutton, C.; Körzdörfer, T.; Gray, M. T.; Brunsfeld, M.; Parrish, R. M.; Sherrill, C. D.; Sears, J. S.; Brédas, J.-L. Accurate Description of Torsion Potentials in Conjugated Polymers Using Density Functionals with Reduced Self-Interaction Error. J. Chem. Phys. 2014, 140 (5), 54310. Körzdörfer, T.; Brédas, J.-L. Organic Electronic Materials: Recent Advances in the DFT Description of the Ground and Excited States Using Tuned Range-Separated Hybrid Functionals. Acc. Chem. Res. 2014, 47 (11), 3284–3291. Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. Excitation Gaps of Finite-Sized Systems from Optimally Tuned Range-Separated Hybrid Functionals. J. Chem. Theory Comput. 2012, 8 (5), 1515–1531. Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom–atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10 (44), 6615–6620. 25 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 26 of 31

(16) Kohn, W. Nobel Lecture: Electronic Structure of Matter--Wave Functions and Density Functionals. Rev. Mod. Phys. 1999, 71 (5), 1253–1266. (17) Brédas, J.-L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Charge-Transfer and EnergyTransfer Processes in π-Conjugated Oligomers and Polymers:  A Molecular Picture. Chem. Rev. 2004, 104 (11), 4971–5004. (18) Lemaur, V.; da Silva Filho, D. A.; Coropceanu, V.; Lehmann, M.; Geerts, Y.; Piris, J.; Debije, M. G.; van de Craats, A. M.; Senthilkumar, K.; Siebbeles, L. D. A.; Warman, J. M.; Brédas, J.-L.; Cornil, J. Charge Transport Properties in Discotic Liquid Crystals:  A Quantum-Chemical Insight into Structure−Property Relationships. J. Am. Chem. Soc. 2004, 126 (10), 3271–3279. (19) Cheung, D. L.; Troisi, A. Theoretical Study of the Organic Photovoltaic Electron Acceptor PCBM: Morphology, Electronic Structure, and Charge Localization. J. Phys. Chem. C 2010, 114 (48), 20479–20488. (20) Yi, Y.; Coropceanu, V.; Brédas, J.-L. Exciton-Dissociation and Charge-Recombination Processes in Pentacene/C60 Solar Cells: Theoretical Insight into the Impact of Interface Geometry. J. Am. Chem. Soc. 2009, 131 (43), 15777–15783. (21) Yi, Y.; Coropceanu, V.; Brédas, J.-L. A Comparative Theoretical Study of ExcitonDissociation and Charge-Recombination Processes in Oligothiophene/Fullerene and Oligothiophene/Perylenediimide Complexes for Organic Solar Cells. J. Mater. Chem. 2011, 21 (5), 1479–1486. (22) Ryno, S. M.; Fu, Y.-T.; Risko, C.; Bredas, J.-L. Polarization Energies at Organic-Organic Interfaces: Impact on the Charge Separation Barrier at Donor-Acceptor Interfaces in Organic Solar Cells. ACS Appl. Mater. Interfaces 2016, 8 (24), 15524–15534. (23) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.; Silbey, R.; Brédas, J.-L. Charge Transport in Organic Semiconductors. Chem. Rev. 2007, 107 (4), 926–952. (24) Fu, Y.-T.; da Silva Filho, D. A.; Sini, G.; Asiri, A. M.; Aziz, S. G.; Risko, C.; Brédas, J.L. Structure and Disorder in Squaraine–C60 Organic Solar Cells: A Theoretical Description of Molecular Packing and Electronic Coupling at the Donor–Acceptor Interface. Adv. Funct. Mater. 2014, 24 (24), 3790–3798. (25) Valeev, E. F.; Coropceanu, V.; da Silva Filho, D. A.; Salman, S.; Brédas, J.-L. Effect of Electronic Polarization on Charge-Transport Parameters in Molecular Organic Semiconductors. J. Am. Chem. Soc. 2006, 128 (30), 9882–9886. (26) Sun, F.; Jin, R. DFT and TD-DFT Study on the Optical and Electronic Properties of Derivatives of 1,4-bis(2-Substituted-1,3,4-Oxadiazole)benzene. Arab. J. Chem. 2013. doi: 10.1016/j.arabjc.2013.11.037 (27) Burke, T. M.; Sweetnam, S.; Vandewal, K.; McGehee, M. D. Beyond Langevin Recombination: How Equilibrium Between Free Carriers and Charge Transfer States Determines the Open-Circuit Voltage of Organic Solar Cells. Adv. Energy Mater. 2015, 5 (11). doi: 10.1002/aenm.201500123 (28) Gilbert, A. T. B.; Besley, N. A.; Gill, P. M. W. Self-Consistent Field Calculations of Excited States Using the Maximum Overlap Method (MOM). J. Phys. Chem. A 2008, 112 (50), 13164–13171. (29) Filatov, M.; Shaik, S. A Spin-Restricted Ensemble-Referenced Kohn–Sham Method and Its Application to Diradicaloid Situations. Chem. Phys. Lett. 1999, 304 (5–6), 429–437.

26 ACS Paragon Plus Environment

Page 27 of 31

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

Chemistry of Materials

(30) Kowalczyk, T.; Tsuchimochi, T.; Chen, P.-T.; Top, L.; Voorhis, T. V. Excitation Energies and Stokes Shifts from a Restricted Open-Shell Kohn-Sham Approach. J. Chem. Phys. 2013, 138 (16), 164101. (31) Hu, C.; Sugino, O.; Miyamoto, Y. Modified Linear Response for Time-Dependent Density-Functional Theory: Application to Rydberg and Charge-Transfer Excitations. Phys. Rev. A 2006, 74 (3), 32508. (32) Martin, R. L. Natural Transition Orbitals. J. Chem. Phys. 2003, 118 (11), 4775–4777. (33) Zhang, C.-R.; Sears, J. S.; Yang, B.; Aziz, S. G.; Coropceanu, V.; Brédas, J.-L. Theoretical Study of the Local and Charge-Transfer Excitations in Model Complexes of Pentacene-C60 Using Tuned Range-Separated Hybrid Functionals. J. Chem. Theory Comput. 2014, 10 (6), 2379–2388. (34) Lu, T.; Chen, F. Multiwfn: A Multifunctional Wavefunction Analyzer. J. Comput. Chem. 2012, 33 (5), 580–592. (35) Guido, C. A.; Cortona, P.; Mennucci, B.; Adamo, C. On the Metric of Charge Transfer Molecular Excitations: A Simple Chemical Descriptor. J. Chem. Theory Comput. 2013, 9 (7), 3118–3126. (36) Yang, B.; Yi, Y.; Zhang, C.-R.; Aziz, S. G.; Coropceanu, V.; Brédas, J.-L. Impact of Electron Delocalization on the Nature of the Charge-Transfer States in Model Pentacene/C60 Interfaces: A Density Functional Theory Study. J. Phys. Chem. C 2014, 118 (48), 27648–27656. (37) Pandey, L.; Doiron, C.; Sears, J. S.; Brédas, J.-L. Lowest Excited States and Optical Absorption Spectra of Donor–acceptor Copolymers for Organic Photovoltaics: A New Picture Emerging from Tuned Long-Range Corrected Density Functionals. Phys. Chem. Chem. Phys. 2012, 14 (41), 14243–14248. (38) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19 (4), 553–566. (39) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van Der Waals Complexes. Chem. Rev. 1994, 94 (7), 1887–1930. (40) Hohenstein, E. G.; Sherrill, C. D. Wavefunction Methods for Noncovalent Interactions. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (2), 304–326. (41) Parker, T. M.; Burns, L. A.; Parrish, R. M.; Ryno, A. G.; Sherrill, C. D. Levels of Symmetry Adapted Perturbation Theory (SAPT). I. Efficiency and Performance for Interaction Energies. J. Chem. Phys. 2014, 140 (9), 94106. (42) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; Russ, N. J.; Leininger, M. L.; Janssen, C. L.; Seidl, E. T.; Allen, W. D.; Schaefer, H. F.; King, R. A.; Valeev, E. F.; Sherrill, C. D.; Crawford, T. D. Psi4: An Open-Source Ab Initio Electronic Structure Program. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (4), 556–565. (43) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225–11236. (44) Wolf, J.; Cruciani, F.; El Labban, A.; Beaujuge, P. M. Wide Band-Gap 3,4Difluorothiophene-Based Polymer with 7% Solar Cell Efficiency: An Alternative to P3HT. Chem. Mater. 2015, 27 (12), 4184–4187. 27 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 28 of 31

(45) Huang, D. M.; Faller, R.; Do, K.; Moulé, A. J. Coarse-Grained Computer Simulations of Polymer/Fullerene Bulk Heterojunctions for Organic Photovoltaic Applications. J. Chem. Theory Comput. 2010, 6 (2), 526–537. (46) Singh, U. C.; Kollman, P. A. An Approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5 (2), 129–145. (47) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11 (4), 431–439. (48) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97 (40), 10269–10280. (49) Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P. The R.E.D. Tools: Advances in RESP and ESP Charge Derivation and Force Field Library Building. Phys. Chem. Chem. Phys. 2010, 12 (28), 7821–7839. (50) Jackson, N. E.; Kohlstedt, K. L.; Savoie, B. M.; Olvera de la Cruz, M.; Schatz, G. C.; Chen, L. X.; Ratner, M. A. Conformational Order in Aggregates of Conjugated Polymers. J. Am. Chem. Soc. 2015, 137 (19), 6254–6262. (51) Caddeo, C.; Mattoni, A. Atomistic Investigation of the Solubility of 3-Alkylthiophene Polymers in Tetrahydrofuran Solvent. Macromolecules 2013, 46 (19), 8003–8008. (52) Newbloom, G. M.; Hoffmann, S. M.; West, A. F.; Gile, M. C.; Sista, P.; Cheung, H.-K. C.; Luscombe, C. K.; Pfaendtner, J.; Pozzo, L. D. Solvatochromism and Conformational Changes in Fully Dissolved Poly(3-Alkylthiophene)s. Langmuir 2015, 31 (1), 458–468. (53) Cho, E.; Risko, C.; Kim, D.; Gysel, R.; Cates Miller, N.; Breiby, D. W.; McGehee, M. D.; Toney, M. F.; Kline, R. J.; Bredas, J.-L. Three-Dimensional Packing Structure and Electronic Properties of Biaxially Oriented Poly(2,5-bis(3-Alkylthiophene-2yl)thieno[3,2-B]thiophene) Films. J. Am. Chem. Soc. 2012, 134 (14), 6177–6190. (54) Poelking, C.; Andrienko, D. Effect of Polymorphism, Regioregularity and Paracrystallinity on Charge Transport in Poly(3-Hexylthiophene) [P3HT] Nanofibers. Macromolecules 2013, 46 (22), 8941–8956. (55) Alexiadis, O.; Mavrantzas, V. G. All-Atom Molecular Dynamics Simulation of Temperature Effects on the Structural, Thermodynamic, and Packing Properties of the Pure Amorphous and Pure Crystalline Phases of Regioregular P3HT. Macromolecules 2013, 46 (6), 2450–2467. (56) McMahon, D. P.; Cheung, D. L.; Troisi, A. Why Holes and Electrons Separate So Well in Polymer/Fullerene Photovoltaic Cells. J. Phys. Chem. Lett. 2011, 2 (21), 2737–2741. (57) Do, K.; Huang, D. M.; Faller, R.; Moulé, A. J. A Comparative MD Study of the Local Structure of Polymer Semiconductors P3HT and PBTTT. Phys. Chem. Chem. Phys. 2010, 12 (44), 14735–14739. (58) Do, K.; Risko, C.; Anthony, J. E.; Amassian, A.; Brédas, J.-L. Dynamics, Miscibility, and Morphology in Polymer:Molecule Blends: The Impact of Chemical Functionality. Chem. Mater. 2015, 27 (22), 7643–7651. (59) Huang, D. M. Computational Study of P3HT/C60-Fullerene Miscibility. Aust. J. Chem. 2014, 67 (4), 585–591. (60) Hoover, W. G. Constant-Pressure Equations of Motion. Phys. Rev. A 1986, 34 (3), 2499– 2500.

28 ACS Paragon Plus Environment

Page 29 of 31

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

Chemistry of Materials

(61) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31 (3), 1695–1697. (62) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182–7190. (63) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117 (1), 1–19. (64) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91 (1–3), 43–56. (65) Accelrys Software, Inc. Materials Studio; San Diego, 2007. (66) Do, K.; Saleem, Q.; Ravva, M. K.; Cruciani, F.; Kan, Z.; Wolf, J.; Hansen, M. R.; Beaujuge, P. M.; Brédas, J.-L. Impact of Fluorine Substituents on π-Conjugated Polymer Main-Chain Conformations, Packing, and Electronic Couplings. Adv. Mater. 2016. doi: 10.1002/adma.201601282 (67) Tummala, N. R.; Bruner, C.; Risko, C.; Brédas, J.-L.; Dauskardt, R. H. Molecular-Scale Understanding of Cohesion and Fracture in P3HT:Fullerene Blends. ACS Appl. Mater. Interfaces 2015, 7 (18), 9957–9964. (68) Tummala, N. R.; Risko, C.; Bruner, C.; Dauskardt, R. H.; Brédas, J.-L. Entanglements in P3HT and Their Influence on Thin-Film Mechanical Properties: Insights from Molecular Dynamics Simulations. J. Polym. Sci. Part B Polym. Phys. 2015, 53 (13), 934–942. (69) Lee, C.-K.; Pao, C.-W.; Chu, C.-W. Multiscale Molecular Simulations of the Nanoscale Morphologies of P3HT:PCBM Blends for Bulk Heterojunction Organic Photovoltaic Cells. Energy Environ. Sci. 2011, 4 (10), 4124–4132. (70) Lee, C.-K.; Pao, C.-W. Solubility of [6,6]-Phenyl-C61-Butyric Acid Methyl Ester and Optimal Blending Ratio of Bulk Heterojunction Polymer Solar Cells. J. Phys. Chem. C 2012, 116 (23), 12455–12461. (71) Schwarz, K. N.; Kee, T. W.; Huang, D. M. Coarse-Grained Simulations of the SolutionPhase Self-Assembly of poly(3-Hexylthiophene) Nanostructures. Nanoscale 2013, 5 (5), 2017–2027. (72) Chiu, M.; Kee, T. W.; Huang, D. M. Coarse-Grained Simulations of the Effects of Chain Length, Solvent Quality, and Chemical Defects on the Solution-Phase Morphology of MEH-PPV Conjugated Polymers. Aust. J. Chem. 2012, 65 (5), 463–471. (73) Izvekov, S.; Voth, G. A. Multiscale Coarse Graining of Liquid-State Systems. J. Chem. Phys. 2005, 123 (13), 134105. (74) Izvekov, S.; Voth, G. A. A Multiscale Coarse-Graining Method for Biomolecular Systems. J. Phys. Chem. B 2005, 109 (7), 2469–2473. (75) Huang, D. M.; Moule, A. J.; Faller, R. Characterization of Polymer–fullerene Mixtures for Organic Photovoltaics by Systematically Coarse-Grained Molecular Simulations. Fluid Phase Equilib. 2011, 302 (1–2), 21–25. (76) Lee, C.-K.; Pao, C.-W. Nanomorphology Evolution of P3HT/PCBM Blends during Solution-Processing from Coarse-Grained Molecular Simulations. J. Phys. Chem. C 2014, 118 (21), 11224–11233. (77) Jones, M. L.; Huang, D. M.; Chakrabarti, B.; Groves, C. Relating Molecular Morphology to Charge Mobility in Semicrystalline Conjugated Polymers. J. Phys. Chem. C 2016, 120 (8), 4240–4250.

29 ACS Paragon Plus Environment

Chemistry of Materials

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

Page 30 of 31

(78) Root, S. E.; Savagatrup, S.; Pais, C. J.; Arya, G.; Lipomi, D. J. Predicting the Mechanical Properties of Organic Semiconductors Using Coarse-Grained Molecular Dynamics Simulations. Macromolecules 2016. (79) Groves, C. Developing Understanding of Organic Photovoltaic Devices: Kinetic Monte Carlo Models of Geminate and Non-Geminate Recombination, Charge Transport and Charge Extraction. Energy Environ. Sci. 2013, 6 (11), 3202–3217. (80) Lyons, B. P.; Clarke, N.; Groves, C. The Relative Importance of Domain Size, Domain Purity and Domain Interfaces to the Performance of Bulk-Heterojunction Organic Photovoltaics. Energy Environ. Sci. 2012, 5 (6), 7657–7663. (81) Groves, C. Suppression of Geminate Charge Recombination in Organic Photovoltaic Devices with a Cascaded Energy Heterojunction. Energy Environ. Sci. 2013, 6 (5), 1546– 1551. (82) Jones, M. L.; Chakrabarti, B.; Groves, C. Monte Carlo Simulation of Geminate Pair Recombination Dynamics in Organic Photovoltaic Devices: Multi-Exponential, FieldDependent Kinetics and Its Interpretation. J. Phys. Chem. C 2014, 118 (1), 85–91. (83) Jones, M. L.; Dyer, R.; Clarke, N.; Groves, C. Are Hot Charge Transfer States the Primary Cause of Efficient Free-Charge Generation in Polymer:fullerene Organic Photovoltaic Devices? A Kinetic Monte Carlo Study. Phys. Chem. Chem. Phys. 2014, 16 (38), 20310– 20320. (84) Tummala, N. R.; Mehraeen, S.; Fu, Y.-T.; Risko, C.; Brédas, J.-L. Materials-Scale Implications of Solvent and Temperature on [6,6]-Phenyl-C61-Butyric Acid Methyl Ester (PCBM): A Theoretical Perspective. Adv. Funct. Mater. 2013, 23 (46), 5800–5813.

30 ACS Paragon Plus Environment

Page 31 of 31

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

Chemistry of Materials

TOC Figure:

31 ACS Paragon Plus Environment