Modelling the morphology and thermomechanical behaviour of low-bandgap conjugated polymers and bulk heterojunction films

1.2 Melt Phase System Initialization. Isolated chains in the fully extended conformation were first subjected to elevated temperature (800 K) simulations using Langevin dynamics (damping parameter = 3800 fs, time step = 2 fs). Random conformations were outputted on nanosecond intervals. This protocol resulted in 60 uncorrelated chain conformation for each polymer simulated. These chains were next assigned a randomly generated centre of mass and orientation with respect to an orthorhombic, periodic simulation box at a mass density of 0.01 g cm. Placement of chains was constrained such that no two atoms on different chains were within 20 Å of another. This was especially important for simulations containing PC71BM, as to avoid spearing the fullerene cage, which would result in unstable dynamics. The randomly generated configurations were next subjected to a conjugategradient energy minimization to further ensure that there were no unphysical interactions to make the simulation unstable in the first couple of time steps. Next, the temperature was ramped from 0 to 800 K over the course of 1 ns using Langevin dynamics (damping parameter = 3800 fs, dielectric constant = 9.8). Next, a Nosé-Hoover style barostat (time constant = 1000 fs) was used to relax the simulation box at 1 atmosphere of pressure until the density converged to the equilibrium melt-phase value at 800 K and 1 atm over the course of 5 ns. For all bulk heterojunction simulations, a mass fraction of 1:1.5 polymer: PCBM was used. This composition is consistent with experimentally optimized devices.


Introduction
Mechanical deformability underpins many of the advantages of organic semiconductors. 1 For example, use of low-cost, ultrathin substrates, fabrication of devices by roll-to-roll printing, stability of devices in portable applications, and biological integration require the active materials to accommodate levels of mechanical strain that are sometimes moderate (e.g., roll-to-roll printing) but often extreme (e.g., biological integration).The mechanical properties of an organic semiconductor depend not only on its molecular structure, but also on how the molecules pack in the solid state. 2 This morphology, however, is difficult to predict by simple inspection of the molecular structure. 3This paper describes the use of molecular dynamics (MD) to understand the nanoscale structural features and mechanical properties of low-bandgap polymers that exhibit the donor-acceptor (DA) motif. 4,5t has recently been demonstrated experimentally that many of the highest-performing organic semiconductors lack the mechanical stability required for flexible, stretchable, and mechanically robust applications. 3For example, all three of the high-performance polymers shown in Fig. 1a-c have been measured to crack at r3% strain using film-on-elastomer techniques. 3,6Furthermore, the addition of electron-accepting fullerene derivatives, such as PC 71 BM (Fig. 1d) to form an organic bulk heterojunction (BHJ) solar cell, has been shown experimentally 7 and computationally 8 to result in a significant increase in stiffness and brittleness of the composite film.
The science of organic materials has benefitted from the ability to modify molecular structure and measure the effects of these changes on function. 9This approach has helped to understand charge-transport in organic semiconductors, 10 and recently to understand their mechanical behaviour. 7,11Experimental correlations between molecular structure, mechanical, and optoelectronic properties of p-conjugated materials-even large libraries containing over 50 compounds 7 -have yielded insightful structure-property relationships.However, an understanding of how molecular structure of p-conjugated materials influences formation of films and dissipation of mechanical energy is far from complete.Computational resources and techniques have progressed to a level where it is now feasible to apply them to developing detailed structure-property relationships of polymers whose repeat units contain as many as 100 atoms, e.g., the DA polymers.
MD simulations are a useful tool for understanding both the conformational behaviour and the thermomechanical properties of organic molecular materials.The level of detail accessible by simulation is difficult or impossible to attain experimentally.Modelling of conjugated polymers entails a unique challenge for standard classical force fields.This difficulty arises in part because of the strong geometrical constraints imposed by p-conjugation along the backbone.Specifically, generic force fields are not directly transferrable to these heterocyclic systems, and thus fail to accurately predict the interactions that arise from torsional rotations between rings, as well as atomic partial charges.This shortcoming requires each polymer under study to be parameterized independently from electronic structure calculations. 12Density functional and wave-function based methods for building atomistic models for conjugated polymers have been established and applied to a variety of materials, such as the poly(3-alkylthiophene)s (P3ATs), 13,14 and recently to a library of low-bandgap polymers. 15While these atomistic models 14 -and detailed coarse-grained models derived from them 16 -have been used to investigate the mechanical properties of P3ATs, 8,17 they have not been applied to the more structurally complex low-bandgap polymers.Although P3HT is typically used as a model system, many results obtained from such case studies cannot be directly generalized to the structurally complex low-bandgap polymers. 18n this work, we studied the morphologies and mechanical properties of three low-bandgap polymers and their blends with PC 71 BM.The goal was to identify the structural determinants and molecular mechanisms for the accommodation of mechanical energy.Our choice of polymers (Fig. 1a-c) was influenced by two requirements: (1) accurate atomistic models exist for these materials 15 and (2) the mechanical properties have been measured experimentally. 3Moreover, each of these three polymers occupy specific points on the range of mesoscale ordering available to this class of materials, from the structurally disordered TQ1 polymer 19,20 to the locally ordered PTB7 21,22 and the semi crystalline PDTSTPD. 23,24Our efforts were guided by three questions: (1) Effect of molecular structure.How do molecular attributes such as the configuration of the side chain (i.e.branched vs. linear) and the connectivity of the conjugated backbone (i.e.fused vs. isolated rings vs. off-axis conjugated units) dictate the solidstate packing and mechanical behaviour of these materials?
(2) Effect of morphology.How do processing conditions affect the mechanical properties of these materials-specifically the ways in which solution-casting generates non-equilibrium kinetically trapped morphologies-and how can these conditions be used to increase the robustness of devices?
(3) Agreement with experiment.More generally, do the results of MD simulations match those of experiments?Can MD simulations be used to complement experimental research aimed at understanding and optimizing the morphology and thermomechanical properties of organic semiconductors?
To begin deconvoluting the effects of molecular structure on solid-state packing and bulk mechanical properties, we simulated the melt-phase equilibrium structure of each polymer, as well as their composites with PC 71 BM.We employed a coarse-grained

View Article Online
analysis to compare the equilibrium statistics of key intramolecular degrees of freedom, single-chain conformations, and multi-chain packing preferences.These melt-phase features provided an approximate representation of the thermodynamic forces driving the solidification of films cast from solution.
Although it was impractical to simulate the process of solvent evaporation in atomistic detail, we included effects of this process using two disparate simulation protocols which generated two different types of glassy structures.These two protocols approximately simulated the effects of casting films from solvents of different quality and evaporation rates.We simulated the mechanical behaviour of the resulting glassy structures by imposing a uniaxial tensile deformation and computing a stressstrain curve.Our results indicated that differences in processing, composition, and solution-phase behaviour play a significant role in determining the mechanical properties of these materials.These properties may in fact be influenced more strongly by differences between the ways in which the films are formed than by differences in molecular structure.That is, the extent to which the molecular structure dictates the mechanical properties of these materials is limited by the extent to which it influences its structure in solution, and thus in the solid state.

Simulation design
Choice of materials.The three polymers shown in Fig. 1a-c were chosen because they are among the highest performing materials in the OPV literature and have been extensively characterized.It was thus possible to compare predictions of tensile modulus, glass transition temperature, and p-stacking distance with experimental results.Additionally, these materials are structurally diverse among the backbones and side chains.This diversity allowed us to draw qualitative comparisons of properties resulting from chemical motifs including fused rings, off-axis conjugated units, and side-chains that are linear, branched, and orthogonally extended from the backbone.
We simulated a set of systems containing 60 monodisperse oligomers with 12 repeat units each.Experimental studies on the mechanical properties of these specific polymers used molecular weights on this order of magnitude (B10 kDa) and thus direct comparisons could be made. 3Moreover, this chain length marks the onset of self-folding and the formation of photophysical aggregates observable in solution-phase UV-Vis absorption and fluorescence spectra. 25Furthermore, this degree of polymerization has also proven sufficient for predicting the condensed phase behaviour and mechanical properties of P3HT (while allowing for computationally accessible simulations) and has served as the basis for detailed coarse-grained models capable of describing experimentally relevant chain dimensions and polydispersity. 16,17Finally, DA polymers can be difficult to synthesize at a high molecular weight and thus the properties of oligomers are also relevant. 26hoice of simulation protocol.The structures and properties of conjugated polymers thin films are known to depend critically on processing. 27Guo et al. have investigated the effect of solvents and solvent additives on the morphology of PTB7-F40 homopolymer films using X-ray scattering.These researchers have determined that the crystallinity of PTB7 films is not affected by the solvent used unless a high boiling point solvent additive, di-iodooctane, is added. 28This additive has been determined to act as a plasticizer and allows the surface of the polymer film to dynamically rearrange instead of being completely kinetically trapped. 29Instead of including these additives in our simulations explicitly, we have included the effect of film nanostructure on the mechanical properties by applying two different methods for generating glassy solid morphologies.
These morphologies represent opposing ends of a continuum of possible structures available to these chains in the condensed phase.In one method, we used a conventional MD approach for simulating polymer glasses: quenching from the equilibrated melt phase.This approach allowed for the characterization of equilibrium structure in the melt and estimation of the glass transition temperature.We termed the structures prepared in this way ''melt-quenched'' morphologies.In the other method, we performed simulated annealing on implicitly solvated individual chains to form self-aggregated structures and then packed them together at room temperature to mimic solution casting rapidly from a poor, low boiling point solvent (see Fig. S1, in the ESI, † for a visual comparison of the two protocols).The use of this procedure was influenced by experimental spectroscopic evidence suggesting these materials tend to form aggregated structures in solution. 25We expected that these pre-aggregated chains would be kinetically trapped in the solid state.Pre-aggregation would lead to a lower density of entanglements (and this might be the cause of experimentally observed brittleness in these materials).We termed the structures prepared in this way the ''self-aggregated'' morphologies.For PTB7 in particular, the structural effects of solution processing have been characterized rigorously by Al-Hussein et al.Their X-ray experiments have determined that spray-coated PTB7 films are kinetically-trapped in a morphology that is significantly more disordered than spin-coated films. 30he morphology of low-bandgap polymers is also known to depend strongly upon interactions with the surface to which they are deposited. 31For PTB7 in particular, the structure has been quantified rigorously by Hammond et al., who employed an array of sophisticated experimental techniques and models to determine the nanoscale morphology of homopolymer and bulk heterojunction films.These researchers determined that PTB7 has a moderate preference for p-stacking in the vertical direction (''face-on''), however, the coherence length for these orientational correlations was determined to be quite small (ca.1.7 nm). 32This experimental result (and computational constraints) lend credence to our decision to not include a substrate in our simulations and instead use a domain that is periodic in all three dimensions to mimic the mid-section of the thin film.
Choice of models.We chose classical atomistic MD models to simulate the structural properties of these materials because they can describe large systems containing greater than 100 000 atoms.We used custom force fields based on the all-atomistic optimized potential for liquid simulation 33 (OPLS-AA).Such models give excellent agreement with experimental values for density, heat capacities, and compressibility of organic liquids.Moreover, these force fields have previously been used to simulate the mechanical behaviour of P3HT and P3HT:PC 61 BM blends. 17,34A detailed description of the process used to build these models using quantum mechanical methods is given elsewhere. 15All simulations and visualizations were performed with LAMMPS 35 and OVITO 36 respectively.

Results and discussion
Intramolecular melt-phase structure Dihedral distributions.We began by untangling the contributions of various molecular features on the melt-phase structure of these polymers at 600 K and 1 atm.Under these conditions, the polymers are sufficiently mobile to allow for efficient computational sampling of equilibrium configurations.An important factor dictating the dynamics of conjugated polymers are the strong dihedral forces imposed by p-conjugation along the backbone of the polymer.These dihedral preferences play an important role in determining bulk photophysical and mechanical properties.Extended planar configurations facilitate intrachain charge transport, and the dihedral geometries play an important role in determining the overall conformational behavior of the polymer chains.
To analyse the equilibrium dihedral preferences, we defined an orthonormal basis set for each conjugated DA unit such as the one shown in Fig. 2a for PTB7.These geometrical definitions enabled the calculation of probability distributions for dihedral angles between adjacent conjugated rings shown in Fig. 2b-d.Due to the asymmetry of the acceptor unit of PTB7, we split the dihedral angle into two distinct distributions (the one labeled DA is shown in Fig. 2a).Comparing across the three materials we observed several notable features.We found that TQ1 had the highest fraction of syn conformations, while PDTSTPD and PTB7 had approximately equal preferences for syn and anti conformations.The preference of TQ1 for syn conformations is likely due to the small size and rotational flexibility of the donor unit, which contains no side chain.
By defining planar conformations as the regions within the green domains on Fig. 2b and c (AE401 from true planarity for optimal delocalization of p electrons, see the next paragraph) we found that PDTSTPD had the highest fraction of planar conformations (and thus the greatest energy barrier for rotation) followed by PTB7 and TQ1 (see Table S1 (ESI †) for a quantitative analysis of the distributions).We ascribe the planarity of PDTSTPD to the arrangement of the side chains on the donor unit (in addition to the intrinsic dihedral potential).Since the side chains are attached orthogonally to the conjugated unit, steric forces inhibit rotations from planar conformations.
Conjugation length.The broad absorption spectra in conjugated polymers are in part due to the distribution of effective conjugation lengths induced by the torsional disorder along the backbone.The semiconducting nature of these polymers thus also depends on the distribution of conjugation lengths (L c ).In order for p-electrons to delocalize, the p orbitals on adjacent atoms must be approximately coplanar: rotations from planarity greater than 401 have a detrimental effect on delocalization. 37sing this value as a cut-off, we grouped conjugated units into distinct chromophores and plotted the probability distribution P(L c ) in Fig. 3a.Taking the number-average (first moment) from these distributions, we found that TQ1 had the shortest conjugation length (2.7 units), compared to that of PTB7 and PDTSTPD (3.1 and 3.2 units), this result follows directly from previous results in which we found that PDTSPTD had the highest fraction of planar units.For PTB7 in particular, the fluorine atom plays a particularly important role.Guo and co-workers have determined experimentally that increasing the fraction of fluorinated acceptor units along the conjugated backbone results in a decreased conjugation length and subsequent blue-shift. 38Dubay et al. have computed the conjugation length distribution for P3HT using a similar approach.They found an average conjugation length of 1.06 monomers for a P3HT chain in vacuum.
Tangent correlations.Equilibrium chain structure can be quantified statistically using a tangent correlation function, ht ˆiÁt ˆi+j i.Here, t ˆi represents the unit tangent vector of the ith  View Article Online conjugated unit along the backbone, and the brackets denote an average over all separations j, polymer chains, and time.The computed tangent correlations, shown in Fig. 3b, exhibited an oscillatory decay, implying the existence of a helical structure.Such helical structure has previously been observing in quantum mechanical simulations of TQ1 oligomers, however, tangent correlation functions were not computed. 39his behaviour likely results from the ribbon-like shape of the backbone, which contains planar units of finite width.
The tangent correlations of ribbon-like objects have been described previously using the following functional form: Lp cos lj l . 40 Here, l is the length of the monomer unit taken as the average value of the donor and acceptor units.L p is the persistence length, and l represents the wavelength of helical structure.We found that exact fits to the simulation data required a series expansion of the analytical expression, implying multiple modes of periodicity and tangent decay in these complex polymer chains.A detailed description of the fitting procedure, parameters obtained, and their significance is given in Section S2.2 of the ESI.† Numerical values obtained from the first-order fit are shown in Table 1.The dimensionless quantity, L p l , describes the helical nature of the chain; in the limiting case of L p { l, the worm-like chain model is valid.In all cases, we found that L p l $ 1, suggesting that a worm-like chain is not the best model for interpreting experimental data obtained from neutron or light scattering experiments aimed at probing the conformational structure of low band-gap polymers. 41For both PTB7 and TQ1 we found that L p l o 1, while for PDTSTPD 1.This analysis suggests the PDTSTPD is the most ribbon-like of the three polymers, which agrees with the experimentally observed order in PDTSTPD when compared to PTB7 and TQ1.An important question that arises is whether the helical structures have a handedness.To determine this, we generated Ramachandran plots 42 which show the correlations between successive dihedral angles (see Fig. S5, ESI †).We found them to be symmetric in all cases, implying that there is no overall handedness to the helical structures.We expect the existence of a helical structure to have important mechanistic implications for charge-transport 43 and the entropic elasticity 44 in this class of materials.
Intermolecular melt-phase structure Donor-acceptor stacking.Intermolecular packing of the conjugated DA units dictates both the interchain charge transport and the cohesive energy of these materials.Stacking preferences can be quantified using a pair distribution function.The results from our simulations of the pure polymers are shown in Fig. 4a-d.Overall, we found that the p-stacking preferences were predominately controlled by the topology of the conjugated unit and the chemical nature of the side chains.The general shapes of the pair distribution functions compare favourably to those obtained from analogous simulations of P3HT (both atomistic and coarse-grained), in that the primary p-stacking peak does not exceed a value of 1. 16,17,45 Focusing on the results for PTB7, we found that although the benzodithiophene (donor) unit has two branched ethylhexyloxy side chains, it has a stronger p-stacking preference than the fluorinated thienothiophene (acceptor) unit that contains a single branched side-chain.This packing feature could be the result of the enhanced flexibility of the side-chains that arises from the higher flexibility of the ether linkage on the benzodithiophene as compared to the ester linkage on the fluorinated thienothiophene.Additionally, the symmetric shape of the benzodithiophene unit allows its side chains to align in the same plane, increasing their overall attraction due to favourable vdW interactions. 46Another distinctive feature in the plots is a broad peak in the A-A interaction at E14 Å which seems to indicate a tendency for interdigitation of side chains as evidenced geometrically in the Fig. S6 (ESI †).
The importance of molecular topology on these interactions is further exemplified by the pair distribution functions obtained for TQ1 in Fig. 4c.We observed that the bulky off-axis units on the benzopyrazine (acceptor) unit serve to inhibit efficient p-stacking in comparison to the thiophene (donor) unit, which has no side chain.The direction that the side chain emanates from the conjugated ring also plays an important role in how the polymers pack.This effect can be seen in the plots for PDTSTPD in Fig. 4d where the D-D peak is one third as prominent as the A-A peak due to the two branched side chains extending orthogonally outward from the dithienosilole (donor) unit.We also note that similar to PTB7, a broad peak at E12 Å can be seen for the A-A interactions, which is not surprising given the geometric similarity of the acceptor units of PTB7 and PDTSTPD.
For all three polymers, we plotted a vertical green line showing the experimentally determined p-stacking distances.The p-stacking View Article Online distance of the three polymers have been reported as 4.2 Å, 24 3.8-4.0Å, 32,47 and 4.3 Å 48 for PDTSTPD, PTB7 and TQ1 respectively.Comparing between simulation and experiment, it can be seen that the experimental stacking distance is slightly smaller than the simulated peak.We attributed this difference to two possible effects.The first is that these curves are obtained from melt-phase simulations at 600 K whereas the experimental results are obtained at room temperature.The second is that the inadequacy of treating the conjugated carbon atoms with a standard Lennard-Jones potential due to extended p-conjugation.Computational studies suggest that a buffered 14-7 potential might provide a more accurate description of excluded volume interactions of these atoms at small interatomic separations. 12olymer-fullerene stacking.The performance of OPV material systems depends strongly upon the intermolecular arrangements in the mixed regions between polymer and fullerene where charge separation takes place.Graham and co-workers have observed a general trend that higher-performing DA polymers generally have acceptor units that are more sterically accessible to the fullerene molecules than the donor units. 49That is, the acceptor units bear linear side chains as opposed to branched ones.The authors provided experimental evidence to show that this trend was due to enhanced coupling of fullerene molecules with the acceptor moiety.Recent computational studies by Wang and co-workers combined classical atomistic MD with density functional theory to provide theoretical evidence in support of this hypothesis for the PBDTTPD family of polymers. 50The authors systematically investigated the impact of side chain structure on packing and conformational arrangement in the mixed region.To further test this trend, we computed pair distribution functions between PC 71 BM (shortened to PCBM in the text and figures) and the DA units.The results of these calculations are shown in Fig. 4f-h.
For TQ1, shown in Fig. 4g, it can be seen that there is a clear preference (sharper peak) for PCBM-A as opposed to PCBM-D stacking.This preference could be the result of the additional van der Waals forces arising from the off-axis conjugated units forming p-stacks with the fullerene cage.For PDTSTPD, shown in Fig. 4h, we observed that the primary PCBM-A peak was sharper and higher amplitude than the corresponding PCBM-D peak.In addition to this peak, there was a clear secondary peak indicating a preference for bilayer stacking of this moiety with the fullerene cage.Finally, for PTB7, shown in Fig. 4f, we observed that the primary peaks were similar in amplitude, distance, and width.The only observable difference was the secondary peak in the PCBM-A curve indicating the presence of bilayer stacks.Overall, the results of our simulations supported the claim that branched side chains inhibit efficient packing with the fullerene unit.

Thermal properties
The kinetics of molecular organization during the annealing and solidification of polymer films is important for determining both the solid-state structure and thermal stability of solid films comprising DA polymers.The glass transition temperature (T g ), which marks the onset of thermally activated relaxation processes of the main chain, is thus a critical processing parameter for printing and coating techniques. 51To estimate T g , we performed thermal quenching simulations while monitoring the density to obtain the thermal curves shown in Fig. 5a-c.The glass-transition region was determined by the narrow range in temperature in which the thermal expansion coefficient underwent a change, and

View Article
T g was thus estimated as the intersection of linear fits to the glassy and melt regions.Glassy configurations obtained from these quenching simulations were subsequently used in uniaxial tensile loading simulations where they are denoted as the melt-quenched morphology.
Pure polymers.Beginning our discussion with the pure polymer systems, we found that all three materials underwent a glass transition well above room temperature, which has been observed experimentally in low-bandgap polymers and is the result of the strong inter-and intramolecular forces arising from extended p-conjugation.This is in contrast to P3HT, which has a glass transition that has been measured experimentally 51 and predicted by coarse-grained MD 8 to be just below room temperature.We found that PTB7 had the highest T g of E400 K followed by PDTSTPD at E380 K and TQ1 at E375 K (the error of these values is approximately AE10 K).These results are in excellent agreement with experimentally determined values of E382 K for PDTSTPD 24 and E373 K for TQ1. 19We note that while there is a disparity in molecular weight between the experimental and simulated systems, it has been observed experimentally that the glass transition was insensitive to changes in molecular weight above 10 kDa for other DA polymers. 52Interestingly, we found a direct correlation between T g and density in these materials; this finding is consistent with the free-volume model for glass transitions. 53ulk heterojunctions.T g is an important parameter for the BHJ because it represents the upper bound for the operating temperature at which a device can function without undergoing detrimental morphological rearrangements. 51Solar modules must be capable of withstanding significant temperature fluctuations (À40 to 85 1C) without coarsening of the morphology or degrading. 54For all three materials (Fig. 5a-c), we observed a shift in the glass transition temperatures and the densities with the addition of PC 71 BM; this finding was consistent with previous experimental 51 and computational studies. 8The shift can ultimately be attributed to the high dispersion forces arising from the fullerene cage and the spherical shape, which enables efficient packing and causes it to act as an anti-plasticizing agent.With glass transitions well above 100 1C we conclude that BHJs based on these DA polymers should be suitably stable to thermal cycling as part of an integrated solar module.

Solution-phase conformations
Strong inter-and intramolecular forces that arise from extensive p-conjugation make low bandgap polymers especially difficult to solubilize.These materials exhibit a wide range of solutionphase conformational behaviour that have a profound effect on the morphology of films cast from solution.To take the effects of solvation into account without explicitly simulating the evaporation process, we performed simulated annealing on implicitly solvated individual chains using Langevin dynamics, following an approach originally described by Jackson et al. 15 This protocol was used to generate minimum energy self-folded chains that have been demonstrated to be present in DA polymer solutions on the basis of spectroscopic evidence. 25he annealed structures can be classified according to the four general polymer conformations: globular, toroidal, folded, and extended (Fig. 6).For each polymer, we performed 60 independent annealing simulations and then classified the resulting structures according to a Q-tensor analysis, based on the Landau-de Gennes order parameter. 55The Q-tensor is a generalized order parameter whose maximum eigenvalue describes the amount of order present in the polymer backbone.A more detailed description of this analysis is given in Section S2.4 of the ESI.† A useful visualization of order in these systems can be obtained by computing an orientational contact

Energy & Environmental Science
Paper map, M ij = t ˆiÁt ˆj, shown in Fig. 6.The results of our classification scheme for each polymer are shown in Fig. 7a.
We found that the conformational preferences had strong agreement with the experimentally determined ordering in these materials.TQ1, which is known to be one of the most structurally disordered low-bandgap polymers, had the strongest tendency to form globular, disordered conformations.PDTSTPD, which exhibits a distinct lamellar stacking peak in X-ray diffraction experiments, had the strongest tendency to form ordered, folded conformations.PTB7, which has significant short-range order but limited long-range order, had the broadest distribution of conformational classes.Interestingly, we found that the annealing simulations of PDTSTPD resulted in no conformations in an extended state.This behaviour exemplifies the strong thermodynamic driving force for folding of this polymer in solution.We performed the same analysis on the meltquenched morphology obtained for each of the three polymers and found all the chains to be in an extended conformation according to this classification scheme.This finding was consistent with the absence of an enthalpic driving force for chain folding in the melt.Similar simulations have been performed for P3HT, using an implicitly solvated coarse-grained model. 56t was found that P3HT has a strong preference for adopting toroidal and folded conformations, however, the structures were not characterized using the same order parameter analysis, so a direct comparison cannot be made.
The increased propensity for chain folding in solution also results in a change in the tangent correlation functions when compared to the melt phase.As shown in Fig. 7b, we found that the tangent correlation functions of the isolated polymer chains (Fig. 3b) all exhibited a more pronounced oscillatory shape than chains in the melt phase, but could still be fit to an analytical expression, as described in Section S2.2 of the ESI.† The tangent correlation function is intrinsically related to the conformational preferences of the chain.The folded and toroidal structure result in more oscillatory correlations while the extended and globular conformations result in more of a worm-like structure.The toroidal structure is the only one that is truly a helix; its pitch is approximately equal to the p-stacking distance of the conjugated units.

Solid morphologies
Melt-quenched morphology.The glassy configurations (T = 300 K, p = 1 atm) obtained from the thermal quenching simulations of pure polymers represent one possible solid-state amorphous morphology accessible to the low-bandgap polymers.This structure is best compared to the morphological arrangement that might be obtained experimentally by spin-coating from a theta solvent, or by subjecting the thin film to thermal annealing followed by rapid quenching.The slowest quenching rates possible with MD are still orders of magnitude faster than the experimental timescale (see Section S1.8 of the ESI †).Analysis of the simulation trajectories reveals that the chains exist in an extended, amorphous state conducive to the formation of entanglements.
Self-aggregated morphology.Films of conjugated polymer that are cast from solution are known to form kinetically trapped morphologies.The quality of the solvent determines the free energy of polymer interactions in solution and thus dictates their morphological arrangement in the solid film cast from solution.To simulate a glassy morphology that might result from spin coating from a poor, low-boiling point solvent, we randomly packed the folded chain conformations described previously into a simulation box and allowed the density to converge at constant temperature and pressure (T = 300 K, p = 1 atm).See Fig. S1 in the ESI, † for a visual representation of this simulation process.This protocol led to a solid structure containing voids and self-folded chains (Fig. 8).This morphology is farther away from equilibrium than the melt-quenched morphology.We expected that the presence of voids would lead to a lower elastic modulus, and the self-folded chains would reduce the number of entanglements and make the selfaggregated morphology softer and more brittle than the meltquenched morphology.
Bulk heterojunction morphology.The OPV bulk heterojunction generally consists of a ternary system containing phase segregated regions of pure polymer, pure fullerene, and a mixed region containing both materials. 57][60] This morphology will be referred to as the bulk heterojunction morphology for comparative purposes (see Fig. 8).2][63] Due to computational constraints, phenomena that occur on such length and time scales can't be captured adequately using atomistic models.For all bulk heterojunction simulations, a mass fraction of 1 : 1.5 polymer to PCBM was used, this is a typical composition used in devices.
Entanglements.The topological constraint that chains can slide past but not pass through one another gives rise to a number of rheological and mechanical properties unique to polymers, including ultimate strain at fracture, tensile strength, and toughness. 646][67][68] We have applied this algorithm to our simulations in order to calculate the average number of interior kinks per chain, which is proportional to the entanglement density.The results of our analysis are shown in Fig. 9.As expected, we found that the melt-quenched morphology contained more entanglements than either the self-aggregated or bulk heterojunction morphologies.Across the three different polymers we found that PTB7 had the most entanglements in all cases, followed by PDTSTPD and finally TQ1.We attribute this finding to the fact that TQ1 had the shortest contour length and the longest persistence length of the three polymers.It is important to note that on average there is only B1 interior kink per chain, therefore these polymers do not form an entangled network at this molecular weight.Accurate coarse-grained models will be necessary for simulating high molecular weight entanglement properties. 8,69This is the subject of on-going work.
Uniaxial tensile response.The capacity for semiconducting materials to accommodate tensile strain is required for applications in stretchable and wearable devices.To simulate the tensile behaviour of the glassy configurations, we imposed a uniaxial deformation and computed the stress response shown in Fig. 10a-c.The tensile modulus (E) was estimated by fitting a straight line to the stress at less than 2% strain.We found that differences in composition and morphology played a more important role in determining the tensile response than differences in molecular structure.As expected, the melt-quenched morphology had a higher tensile modulus and overall toughness than the self-aggregated morphology.Additionally, the composites with PC 71 BM had increased tensile moduli.Increased modulus of the bulk heterojunction morphology is consistent with the increased T g and density of this structure compared to the other two.For all three polymers, we found that the computed tensile modulus obtained from the melt-quenched morphology was comparable to that obtained from similar atomistic simulations of P3HT. 17

View Article
A comparison of the results from the pure polymer simulations to the experiment reveals an interesting finding.As shown in Fig. 11, we observed that the experimental value for the tensile modulus agrees well with the melt-quenched morphology for PTB7, the self-aggregated morphology for PDTSTPD, and is intermediate between the two morphologies for TQ1.Given that PDTSTPD had the strongest tendency to form folded aggregates, it makes sense that the self-aggregated morphology is a more accurate representation of PDTSTPD films cast from solution than is the melt-quenched morphology.The way in which this line of reasoning can be applied to TQ1 and PTB7 is less clear, most likely due to explicit solvation effects that were not taken into account.Nonetheless, the fact that the experimental values for tensile moduli of the three materials span the range bound by the calculated self-aggregated and melt-quenched morphologies suggest a strong role for the different solubility properties of the molecules, and that explicit atomistic simulations of solution-phase structure will be necessary in future work.
In general, these results suggest that differences in the morphology-as opposed to molecular structure-produced the wide range of tensile moduli measured in a library of low-bandgap polymers. 3This conclusion indicates that solution-processing represents a viable route for optimizing the mechanical robustness of these materials.One specific testable prediction is that films cast at high temperature may have improved mechanical stability compared to those cast at room temperature due to an increased extension of polymer chains and corresponding density of entanglements.
It is important to note that we did not observe cracking in any of the tensile loading simulations.This finding, coupled with the fact that all three of these materials exhibited brittle fracture at r3% strain experimentally 3 led us to conclusion that the dominant fracture mechanism in these materials is likely through chain scission and the presence of macroscopic defects that localize stress, neither of which are present in these simulations.Although the experimental crack-on-set strain was small for all three of these materials, there still were minor differences.Of the three materials TQ1 was able to withstand the most strain (2.7%), followed by PTB7 (2.3%) and finally PDTSTPD (1.2%). 3These values correlate well with the glass transition temperature predicted and measured for these materials giving further evidence that the brittle fracture observed experimentally is the result of bond scission.Since the models employed in these simulations rely on a harmonic bond approximation, they are unsuitable for predicting such fracture behaviour.In future mechanical deformation simulations, bond potentials that allow for cleavage will be parameterized for this class of materials.It would also be instructive to perform similar simulations on materials that are known experimentally to exhibit ductile modes of fracture. 3Such an approach would allow for computational identification of the molecular-scale attributes that enable stress relaxation mechanisms other than bond scission.

Conclusions
In this paper we have demonstrated the utility of MD simulations to elucidate the structural and mechanical behaviour of low-bandgap polymers and composites.A summary of important material parameter predictions is given in Table 1.Through comparison to experimental results for p-stacking distances,  View Article Online glass transition temperatures, and tensile moduli, we have identified the strengths and limitations of this approach and provided suggestions for further refinement of custom models parameterized from electronic structure calculations (e.g., the use of a buffered 14-7 Lennard-Jones potential, and the allowance for bond scission).Through the computation of various distribution functions, we have provided insights that could be of use to synthetic chemists designing new materials.Furthermore, our simulations have revealed several new findings previously unconsidered in the literature that warrant further investigation.
We found that the rigid, planar, fused rings along the backbone give rise to the formation of helical structures in the melt-phase that are characteristic of ribbon-like chains.This finding suggests that the classic worm-like chain model might not be appropriate for the interpretation of scattering experiments investigating the conformational behaviour of these polymers.Consideration of DA polymers as ribbon-like chains should also have important implications for the statistical mechanical basis of entropic elasticity 44 and charge transport properties 43 of these chains.
We found that the extent to which the molecular structure dictates the mechanical behaviour of these materials is limited by the extent to which it contributes to the solution-phase behavior and resulting thin-film morphology.These results have led us to the hypothesis that spin-coating from heated solutions of optimized solvent quality should result in more entangled, mechanically robust films.This hypothesis will be the subject of future experimental and computational investigations.
To simulate the solvent evaporation process of a printing or coating procedure accurately would require detailed coarsegrained models such as those that exist for P3HT 16,70 and other complex polymers. 71In principle, the statistics gathered from our melt-phase simulations contain most of the necessary information to build stylistically similar models for these polymers.Additional information, however, must be gathered from explicitly solvated atomistic simulations to compute concentration dependent free energies of interactions capable of accurately modelling the solidification process.This is the subject of ongoing work and will allow for the simulation of an experimentally relevant morphology of a solution-cast thin film containing realistic molecular weight distributions and interfacial effects.

Fig. 2
Fig. 2 Inter-ring dihedral distributions.(a) Diagram showing definitions of basis sets defining ring orientations for PTB7 monomer.Dihedral distributions, P(f), obtained from simulations at 600 K and 1 atm for (b) PDTSTPD (c) PTB7 (d) TQ1.The syn conformations shown in Fig. 1 represent the planar arrangement that results from f = 01.

Fig. 3
Fig. 3 Single-chain statistics.(a) Plot showing probability distribution of conjugation length.(b) Plot showing the decay of the tangent correlation function (dotted line) as well as the fits to the analytical model (solid line).

Table 1
Summary of material property predictions reported.Values displayed are all from the melt-quenched simulations of pure polymers.For each material property, values are colored according to relative magnitudes of the three polymers simulated Energy & Environmental Science Paper Open Access Article.Published on 19 December 2016.Downloaded on 13/01/2017 22:48:59.This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Fig. 4
Fig. 4 Morphological characterization of pure polymers and BHJs.(a) Snapshot of condensed phase morphology of pure PTB7.Molecules are coloured by indexed values.(b-d) Pair-distribution functions showing intermolecular donor-acceptor stacking preferences in the melt (green lines indicate experimental p-p stacking distances 24,47,48 ).(e) Snapshot of condensed phase morphology of PTB7 BHJ.Polymer molecules are coloured by indexed values.(f-h) Pair-distribution functions showing intermolecular coupling between PCBM and donor/acceptor units.Note: the pair distributions are computed using the centres of masses of the conjugated units, indicated by the shading in the figures (i.e. the centre of mass of the fullerene cage for PCBM).

Fig. 6
Fig. 6 Conformational classes.Representative chains and their associated orientational contact maps illustrating the conformational classes available to single chains and the distinctive patterns used for classifying the chain conformations.

Fig. 8
Fig. 8 Comparison of the three glassy, solid morphologies (T = 300 K and p = 1 atm) for PTB7.Snapshots showing 2 nm slices of the simulation trajectory.Polymers are coloured separately.Hydrogen atoms and side chains are omitted for clarity.

Fig. 9
Fig. 9 Entanglement analysis.(a) Plot showing the number of interior kinks per chain for each system simulated.Interior kinks define the primitive path of the polymer chain and are proportional to the number entanglements per chain.Mathematically they are defined as nodes that represent the limit to which the contour can be reduced to a minimum length without violating the topological constraint imposed by excluded volume interactions.(b) Visual representation of two chain contours (side chains and hydrogen atoms removed for clarity) as well as their respective primitive paths demonstrating the concept of an interior kink.

Fig. 10
Fig. 10 Mechanical response.Stress-strain curves showing the responses of (a) the self-aggregated morphology, (b) the melt-quenched morphology, and (c) the bulk heterojunction morphology of each polymer.The solid lines overlayed on the plots are the linear fits used to calculate the tensile moduli.It can be seen in the plots that different morphologies and compositions play a more important role in determining the stress response than differences in molecular structure.

Fig. 11
Fig. 11 Summary of tensile modulus predictions and comparison to experiment.