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

Abstract

This paper describes the use of molecular dynamics (MD) to predict the nanoscale morphology and thermomechanical behavior of three low-bandgap semiconducting polymers and their blends with PC₇₁BM. While the three polymers modeled in this study—PTB7, PDTSTPD, and TQ1—all exhibit the donor–acceptor motif characteristic of high-performance donor materials in organic solar cells, they exemplify different morphologies in the solid state. Predictions from the atomistic simulations presented here include the average conjugation length of the polymers, the structural arrangement of conjugated donor and acceptor units in neat and bulk heterojunction (BHJ) films, as well as the glass transition temperature and tensile modulus of neat and BHJ polymer films. Calculated tangent correlation functions exhibit oscillatory decay. This finding suggests that DA polymers are more appropriately modeled as ribbon-like chains as opposed to worm-like chains. To account for the range of morphologies accessible by processing manipulations, both a melt-quenched and a self-aggregated morphology are prepared. Owing to the greater free volume of the self-aggregated morphology, these solid structures are found to be softer and weaker than the melt-quenched morphologies. The experimental modulus measured previously for PDTSTPD is similar to the predicted self-aggregated morphology, while the experimental modulus of PTB7 is similar to the predicted melt-quenched modulus. Our comparisons with experiment suggest that solution-processing plays a critical role in optimizing the mechanical properties of conjugated polymeric materials. Overall, the results of this study suggest the promise of MD simulations in determining the ways in which molecular structure influences the morphology and mechanical properties of bulk heterojunction films for solar cells and other organic electronic devices.

---

Broader context

Semiconducting polymers can form the active layers of organic solar cells that have achieved power conversion efficiencies of over ten percent. The highest efficiency materials generally contain electron-donating units alternating with electron accepting units along the molecular backbone (the “DA polymers”). When these polymers (and composites when mixed with fullerenes) are processed from solution, the morphology that forms upon evaporation of the solvent can be difficult to predict. The morphology, in turn, influences not only the efficiency with which charges are generated and separated, but also the lifetime and stability against fluctuations in temperature and mechanical insults (i.e., in portable applications and in the outdoor environment). This paper describes a computational method to predict the morphology and mechanical properties of DA polymers. It is based on fully atomistic molecular dynamics and can predict the glass transition temperature (i.e., the transition from brittle, glassy behavior to deformable, rubber-like behavior), the tensile modulus, and the statistical distribution of molecular conformation in the solid state. The study highlights the important role of the morphology—as opposed to the molecular structure—in determining the thermal and mechanical properties of DA polymers and their blends with fullerenes in solar cells.

Introduction

Mechanical deformability underpins many of the advantages of organic semiconductors.¹ 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.² This morphology, however, is difficult to predict by simple inspection of the molecular structure.³ This 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,5

It 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.³ For example, all three of the high-performance polymers shown in Fig. 1a–c have been measured to crack at ≤3% strain using film-on-elastomer techniques.^3,6 Furthermore, the addition of electron-accepting fullerene derivatives, such as PC₇₁BM (Fig. 1d) to form an organic bulk heterojunction (BHJ) solar cell, has been shown experimentally⁷ and computationally⁸ to result in a significant increase in stiffness and brittleness of the composite film.

Fig. 1 Structures of the polymers and molecules simulated in this work, along with their common names. All structures are shown in a syn conformation. (a) PDTSTPD, poly[(5,6-dihydro-5-octyl-4,6-dioxo-4H-thieno[3,4-c]pyrrole-1,3-diyl)[4,4-bis(2-ethylhexyl)-4H-silolo[3,2-b:4,5-b′]dithiophene-2,6-diyl]]. (b) PTB7, poly({4,8-bis[(2-ethylhexyl)oxy]benzo[1,2-b:4,5-b′]dithiophene-2,6-diyl}{3-fluoro-2-[(2-ethylhexyl)carbonyl]thieno[3,4-b]thiophenediyl}). (c) TQ1, poly[[2,3-bis(3-octyloxyphenyl)-5,8-quinoxalinediyl]-2,5-thiophenediyl]. (d) PC₇₁BM, [6,6]-phenyl C₇₁ butyric acid methyl ester.

The science of organic materials has benefitted from the ability to modify molecular structure and measure the effects of these changes on function.⁹ This approach has helped to understand charge-transport in organic semiconductors,¹⁰ and recently to understand their mechanical behaviour.^7,11 Experimental correlations between molecular structure, mechanical, and optoelectronic properties of π-conjugated materials—even large libraries containing over 50 compounds⁷—have yielded insightful structure–property relationships. However, an understanding of how molecular structure of π-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 π-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.¹² Density 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.¹⁵ While these atomistic models¹⁴—and detailed coarse-grained models derived from them¹⁶—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.¹⁸

In this work, we studied the morphologies and mechanical properties of three low-bandgap polymers and their blends with PC₇₁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¹⁵ and (2) the mechanical properties have been measured experimentally.³ Moreover, 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,24 Our 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 solid-state 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₇₁BM. We employed a coarse-grained 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 stress–strain 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 π-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 (∼10 kDa) and thus direct comparisons could be made.³ Moreover, 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.²⁵ Furthermore, 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,17 Finally, DA polymers can be difficult to synthesize at a high molecular weight and thus the properties of oligomers are also relevant.²⁶

Choice of simulation protocol. The structures and properties of conjugated polymers thin films are known to depend critically on processing.²⁷ Guo 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.²⁸ This 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.²⁹ Instead 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.²⁵ We 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.³⁰

The morphology of low-bandgap polymers is also known to depend strongly upon interactions with the surface to which they are deposited.³¹ For 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 π-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).³² This 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³³ (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₆₁BM blends.^17,34 A detailed description of the process used to build these models using quantum mechanical methods is given elsewhere.¹⁵ All simulations and visualizations were performed with LAMMPS³⁵ and OVITO³⁶ 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 π-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.

Fig. 2 Inter-ring dihedral distributions. (a) Diagram showing definitions of basis sets defining ring orientations for PTB7 monomer. Dihedral distributions, P(ϕ), 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 ϕ = 0°.

By defining planar conformations as the regions within the green domains on Fig. 2b and c (±40° from true planarity for optimal delocalization of π 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 π-electrons to delocalize, the p orbitals on adjacent atoms must be approximately coplanar: rotations from planarity greater than 40° have a detrimental effect on delocalization.³⁷ Using 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.³⁸ Dubay 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.

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).

Tangent correlations. Equilibrium chain structure can be quantified statistically using a tangent correlation function, 〈t̂_i·t̂_i+j〉. Here, t̂_i represents the unit tangent vector of the ith 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.³⁹ This 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: .⁴⁰ 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 λ 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, , describes the helical nature of the chain; in the limiting case of L_p ≪ λ, the worm-like chain model is valid. In all cases, we found that , 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.⁴¹ For both PTB7 and TQ1 we found that , while for PDTSTPD . 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.

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

---

An important question that arises is whether the helical structures have a handedness. To determine this, we generated Ramachandran plots⁴² 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⁴³ and the entropic elasticity⁴⁴ 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 π-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 π-stacking peak does not exceed a value of 1.^16,17,45

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 π–π 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).

Focusing on the results for PTB7, we found that although the benzodithiophene (donor) unit has two branched ethylhexyloxy side chains, it has a stronger π-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.⁴⁶ Another distinctive feature in the plots is a broad peak in the A–A interaction at ≈14 Å 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 π-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 ≈12 Å 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 π-stacking distances. The π-stacking distance of the three polymers have been reported as 4.2 Å,²⁴ 3.8–4.0 Å,^32,47 and 4.3 Å⁴⁸ 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 π-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.¹²

Polymer–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.⁴⁹ That 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.⁵⁰ The 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₇₁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 π-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.⁵¹ To 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 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.

Fig. 5 Thermal properties, density vs. temperature curves showing glass-transitions of polymer: PC₇₁BM bulk heterojunctions (top curves) and pure polymers (bottom curves) for (a) PDTSTPD, (b) PTB7, and (c) TQ1.

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 π-conjugation. This is in contrast to P3HT, which has a glass transition that has been measured experimentally⁵¹ and predicted by coarse-grained MD⁸ to be just below room temperature. We found that PTB7 had the highest T_g of ≈400 K followed by PDTSTPD at ≈380 K and TQ1 at ≈375 K (the error of these values is approximately ±10 K). These results are in excellent agreement with experimentally determined values of ≈382 K for PDTSTPD²⁴ and ≈373 K for TQ1.¹⁹ We 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.⁵² Interestingly, 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.⁵³

Bulk 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.⁵¹ Solar modules must be capable of withstanding significant temperature fluctuations (−40 to 85 °C) without coarsening of the morphology or degrading.⁵⁴ For all three materials (Fig. 5a–c), we observed a shift in the glass transition temperatures and the densities with the addition of PC₇₁BM; this finding was consistent with previous experimental⁵¹ and computational studies.⁸ The 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 °C 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 π-conjugation make low bandgap polymers especially difficult to solubilize. These materials exhibit a wide range of solution-phase 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.¹⁵ 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.²⁵

The 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.⁵⁵ The 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 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.

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. 7 Characterization of conformational preferences. (a) Chart showing fraction of each conformation present in the self-aggregated morphology for each type of polymer. (b) Tangent correlation function showing a large tendency for folded chains.

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 melt-quenched 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.⁵⁶ It 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 π-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 self-aggregated morphology softer and more brittle than the melt-quenched morphology.

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.

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.⁵⁷ The glassy configurations obtained from the thermal quenching simulations of the composite system (analogous to the melt-quenched morphology of the pure polymers) provide a representation of the mixed region, which is of utmost importance for charge separation in OPVs.^49,58–60 This morphology will be referred to as the bulk heterojunction morphology for comparative purposes (see Fig. 8). In reality, it is well known that fullerenes have a thermodynamic preference to phase segregate in the bulk heterojunction.^61–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.⁶⁴ Entanglements can be quantified in MD simulations using a primitive path analysis such as the Z1 algorithm of Kröger and coworkers.^65–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 ∼1 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,69 This is the subject of on-going work.

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.

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₇₁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.¹⁷

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.

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.

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

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.³ This 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 ≤3% strain experimentally³ 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%).³ These 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.³ Such 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 π-stacking distances, 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⁴⁴ and charge transport properties⁴³ 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 coarse-grained models such as those that exist for P3HT^16,70 and other complex polymers.⁷¹ In 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.

Acknowledgements

This work was supported by the Air Force Office of Scientific Research (AFOSR) Grant Number FA9550-16-1-0220. Additional support was provided by the Hellman Fellowship awarded to D. J. L. and the ARCS Fellowship awarded to S. E. R. Computational resources to support this work were provided by the Extreme Science and Engineering Discovery Environment (XSEDE) Program through the National Science Foundation grant number ACI-1053575.⁷² Additionally, N. E. J. and S. E. R. would like to thank Professor Mark Ratner for facilitating this collaboration.

References

1. A. Facchetti, Chem. Mater., 2011, 23, 733–758.
2. A. D. Printz and D. J. Lipomi, Appl. Phys. Rev., 2016, 3, 021302.
3. B. Roth, S. Savagatrup, N. De Los Santos, O. Hagemann, J. E. Carle, M. Helgesen, F. Livi, E. Bundgaard, R. R. Søndergaard, F. C. Krebs and D. J. Lipomi, Chem. Mater., 2016, 7, 2363–2373.
4. M. C. Scharber, D. Muehlbacher, M. Koppe, P. Denk, C. Waldauf, A. J. Heeger and C. J. Brabec, Adv. Mater., 2006, 18, 789–794.
5. K. Müllen and W. Pisula, J. Am. Chem. Soc., 2015, 137, 9503–9505.
6. S. Savagatrup, A. D. Printz, T. F. O’Connor, A. V. Zaretski, D. Rodriquez, E. J. Sawyer, K. M. Rajan, R. I. Acosta, S. E. Root and D. J. Lipomi, Energy Environ. Sci., 2015, 8, 55–80.
7. S. Savagatrup, A. S. Makaram, D. J. Burke and D. J. Lipomi, Adv. Funct. Mater., 2014, 24, 1169–1181.
8. S. E. Root, S. Savagatrup, C. J. Pais, G. Arya and D. J. Lipomi, Macromolecules, 2016, 49, 2886–2894.
9. G. M. Whitesides, Isr. J. Chem., 2016, 56, 66–82.
10. J. L. Brédas, D. Beljonne, V. Coropceanu and J. Cornil, Chem. Rev., 2004, 104, 4971–5003.
11. A. D. Printz, S. Savagatrup, D. J. Burke, T. N. Purdy and D. J. Lipomi, RSC Adv., 2014, 4, 13635–13643.
12. K. H. Dubay, M. L. Hall, T. F. Hughes, C. Wu, D. R. Reichman and R. a. Friesner, J. Chem. Theory Comput., 2012, 8, 4556–4569.
13. V. Marcon and G. Raos, J. Am. Chem. Soc., 2006, 128, 1408–1409.
14. V. Marcon and G. Raos, J. Phys. Chem. B, 2004, 108, 18053–18064.
15. N. E. Jackson, K. L. Kohlstedt, B. M. Savoie, M. Olvera de la Cruz, G. C. Schatz, L. X. Chen and M. a. Ratner, J. Am. Chem. Soc., 2015, 137, 6254–6262.
16. D. M. Huang, R. Faller, K. Do and A. J. Moul, J. Chem. Theory Comput., 2010, 6, 526–537.
17. N. R. Tummala, C. Risko, C. Bruner, R. H. Dauskardt and J.-L. Brédas, J. Polym. Sci., Part B: Polym. Phys., 2015, 1–9.
18. A. T. Kleinschmidt, S. E. Root and D. J. Lipomi, J. Mater. Chem. A, 2016 10.1039/C6TA08317J.
19. R. Kroon, R. Gehlhaar, T. T. Steckler, P. Henriksson, C. Müller, J. Bergqvist, A. Hadipour, P. Heremans and M. R. Andersson, Sol. Energy Mater. Sol. Cells, 2012, 105, 280–286.
20. E. Wang, L. Hou, Z. Wang, S. Hellström, F. Zhang, O. Inganäs and M. R. Andersson, Adv. Mater., 2010, 22, 5240–5244.
21. Y. Liang, Z. Xu, J. Xia, S. T. Tsai, Y. Wu, G. Li, C. Ray and L. Yu, Adv. Mater., 2010, 22, 135–138.
22. L. Lu and L. Yu, Adv. Mater., 2014, 26, 4413–4430.
23. T.-Y. Chu, J. Lu, S. Beaupre, Y. Zhang, J.-R. Pouliot, S. Wakim, J. Zhou, M. Leclerc, Z. Li, J. Ding and Y. Tao, J. Am. Chem. Soc., 2011, 133, 4250–4253.
24. T. Y. Chu, J. Lu, S. Beaupré, Y. Zhang, J. R. Pouliot, J. Zhou, A. Najari, M. Leclerc and Y. Tao, Adv. Funct. Mater., 2012, 22, 2345–2351.
25. T. J. Fauvell, T. Zheng, N. E. Jackson, M. A. Ratner, L. Yu and L. X. Chen, Chem. Mater., 2016, 8, 2814–2822.
26. F. Brouwer, J. Alma, H. Valkenier, T. P. Voortman, J. Hillebrand, R. C. Chiechi and J. C. Hummelen, J. Mater. Chem., 2011, 21, 1582.
27. Y. Diao, L. Shaw, Z. Bao and S. C. B. Mannsfeld, Energy Environ. Sci., 2014, 7, 2145–2159.
28. S. Guo, E. M. Herzig, A. Naumann, G. Tainter, J. Perlich and P. Müller-Buschbaum, J. Phys. Chem. B, 2014, 118, 344–350.
29. L. Ye, Y. Jing, X. Guo, H. Sun, S. Zhang, M. Zhang, L. Huo and J. Hou, J. Phys. Chem. C, 2013, 117, 14920–14928.
30. M. Al-Hussein, E. M. Herzig, M. Schindler, F. Löhrer, C. M. Palumbiny, W. Wang, S. V. Roth and P. Müller-Buschbaum, Polym. Eng. Sci., 2016, 889–894.
31. A. Salleo, R. J. Kline, D. M. DeLongchamp and M. L. Chabinyc, Adv. Mater., 2010, 22, 3812–3838.
32. M. R. Hammond, R. J. Kline, A. a. Herzing, L. J. Richter, D. S. Germack, H. Ro, C. L. Soles, D. a. Fischer, T. Xu, L. Yu, M. F. Toney and D. M. Delongchamp, ACS Nano, 2011, 8248–8257.
33. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
34. N. R. Tummala, C. Bruner, C. Risko, J.-L. Brédas and R. H. Dauskardt, ACS Appl. Mater. Interfaces, 2015, 18, 9957–9964.
35. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
36. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 1–7.
37. G. B. Street, J. L. Bredas, B. Themans and J. M. Andre, J. Chem. Phys., 1985, 83, 1323–1329.
38. S. Guo, J. Ning, V. Körstgens, Y. Yao, E. M. Herzig, S. V. Roth and P. Müller-Buschbaum, Adv. Energy Mater., 2015, 5, 1–11.
39. S. Hedstrom and P. Persson, J. Phys. Chem. C, 2012, 116, 26700–26706.
40. L. Giomi and L. Mahadevan, Phys. Rev. Lett., 2010, 104, 2–5.
41. B. Kuei and E. D. Gomez, Soft Matter, 2016 10.1039/c6sm00979d.
42. V. Ramachandran, G. N. Ramakrishnan and C. Sasisekharan, J. Mol. Biol., 1963, 7, 95–99.
43. P. Carbone and A. Troisi, J. Phys. Chem. Lett., 2014, 5, 2637–2641.
44. J. F. Marko and E. D. Siggia, Macromolecules, 1995, 28, 8759–8770.
45. K. Do, D. M. Huang, R. Faller and A. J. Moulé, Phys. Chem. Chem. Phys., 2010, 12, 14735–14739.
46. H. Yao, L. Ye, H. Zhang, S. Li, S. Zhang and J. Hou, Chem. Rev., 2016, 116, 7397–7457.
47. H. J. Son, W. Wang, T. Xu, Y. Liang, Y. Wu, G. Li and L. Yu, J. Am. Chem. Soc., 2011, 133, 1885–1894.
48. T. Yamamoto, B. L. Lee, H. Kokubo, H. Kishida, K. Hirota, T. Wakabayashi and H. Okamoto, Macromol. Rapid Commun., 2003, 24, 440–443.
49. K. R. Graham, C. Cabanetos, J. P. Jahnke, M. N. Idso, A. El Labban, G. O. N. Ndjawa, T. Heumueller, K. Vandewal, A. Salleo, B. F. Chmelka, A. Amassian, P. M. Beaujuge and M. D. Mcgehee, J. Am. Chem. Soc., 2014, 136, 9608–9618.
50. T. Wang, M. K. Ravva and J.-L. Brédas, Adv. Funct. Mater., 2016 DOI:10.1002/adfm.201601134.
51. C. Müller, Chem. Mater., 2015, 27, 2740–2754.
52. C. Müller, E. Wang, L. M. Andersson, K. Tvingstedt, Y. Zhou, M. R. Andersson and O. S. Inganä, Adv. Funct. Mater., 2010, 20, 2124–2131.
53. D. Turnbull and M. H. Cohen, J. Chem. Phys., 1961, 34, 120.
54. ASTM Stand. E 1171S, 2001, i, 1–5.
55. J. P. Pierre-Gilles de Gennes, The Physics of Liquid Crystals, Clarendon, Oxford, 1974.
56. K. N. Schwarz, T. W. Kee and D. M. Huang, Nanoscale, 2013, 5, 2017–2027.
57. N. D. Treat and M. L. Chabinyc, Annu. Rev. Phys. Chem., 2014, 65, 59–81.
58. J. A. Bartelt, Z. M. Beiley, E. T. Hoke, W. R. Mateker, J. D. Douglas, B. A. Collins, J. R. Tumbleston, K. R. Graham, A. Amassian, H. Ade, J. M. J. Fréchet, M. F. Toney and M. D. Mcgehee, Adv. Energy Mater., 2013, 3, 364–374.
59. K. Vandewal, S. Albrecht, E. T. Hoke, K. R. Graham, J. Widmer, J. D. Douglas, M. Schubert, W. R. Mateker, J. T. Bloking, G. F. Burkhard, A. Sellinger, J. M. J. Fréchet, A. Amassian, M. K. Riede, M. D. McGehee, D. Neher and A. Salleo, Nat. Mater., 2014, 13, 63–68.
60. A. J. Moulé and K. Meerholz, Adv. Mater., 2008, 20, 240–245.
61. D. T. Duong, B. Walker, J. Lin, C. Kim, J. Love, B. Purushothaman, J. E. Anthony and T. Q. Nguyen, J. Polym. Sci., Part B: Polym. Phys., 2012, 50, 1405–1413.
62. O. Oklobia, T. S. Shafai, B. Friedel, T. J. K. Brenner, C. R. McNeill, U. Steiner, N. C. Greenham, G. Li, G. Lu, S. Li, H. Tang, X. Yang, G. Lu, L. Li, J. Li, Y. Wang, X. Yang, S. J. Lou, J. M. Szarko, T. Xu, L. Yu, T. J. Marks, L. X. Chen, W.-H. Baek, T.-S. Yoon, H. H. Lee and Y.-S. Kim, J. Mater. Chem., 2010, 11, 1736–1745.
63. X. Fan, G. J. Fang, P. L. Qin, F. Cheng and X. Z. Zhao, Appl. Phys. A: Mater. Sci. Process., 2011, 105, 1003–1009.
64. R. Everaers, Science, 2004, 303, 823–826.
65. S. Shanbhag and M. Kröger, Macromolecules, 2007, 40, 2897–2903.
66. R. S. Hoy, K. Foteinopoulou and M. Kröger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 14–16.
67. N. C. Karayiannis and M. Kröger, Int. J. Mol. Sci., 2009, 10, 5054–5089.
68. M. Kröger, Comput. Phys. Commun., 2005, 168, 209–232.
69. M. L. Jones, D. M. Huang, B. Chakrabarti and C. Groves, J. Phys. Chem. C, 2016, 120, 4240–4250.
70. P. C. Tapping, S. N. Clafton, K. N. Schwarz, T. W. Kee and D. M. Huang, J. Phys. Chem. C, 2015, 119, 7047–7059.
71. V. Agrawal, G. Arya and J. Oswald, Macromolecules, 2014, 47, 3378–3389.
72. T. John, T. Cockerill, I. Foster and K. Gaither, Comput. Sci. Eng., 2014, 16, 62–74.

---

Footnote

† Electronic supplementary information (ESI) available: (1) Simulation details. (2) In-depth analysis of results. (3) LAMMPS data files³⁵ containing atomistic model parameters. See DOI: 10.1039/c6ee03456j

---