Generic Model for Lamellar Self-Assembly in Conjugated Polymers: Linking Mesoscopic Morphology and Charge Transport in P3HT.

We develop a generic coarse-grained model of soluble conjugated polymers, capable of describing their self-assembly into a lamellar mesophase. Polymer chains are described by a hindered-rotation model, where interaction centers represent entire repeat units, including side chains. We introduce soft anisotropic nonbonded interactions to mimic the potential of mean force between atomistic repeat units. The functional form of this potential reflects the symmetry of the molecular order in a lamellar mesophase. The model can generate both nematic and lamellar (sanidic smectic) molecular arrangements. We parametrize this model for a soluble conjugated polymer poly(3-hexylthiophene) (P3HT) and demonstrate that the simulated lamellar mesophase matches morphologies of low molecular weight P3HT, experimentally observed at elevated temperatures. A qualitative charge-transport model allows us to link local chain conformations and mesoscale order to charge transport. In particular, it shows how coarsening of lamellar domains and chain extension increase the charge carrier mobility. By modeling large systems and long chains, we can capture transport between lamellar layers, which is due to rare, but thermodynamically allowed, backbone bridges between neighboring layers.

Introduction

Soluble semiconducting polymers are promising materials for manufacturing flexible, lightweight electronic devices using scalable and low-cost technologies, such as printing.[1−4] Polymer solubility and processability are achieved by mitigating the attraction of conjugated backbones with flexible alkyl side chains. The underlying molecular architecture, together with processing conditions, define chain conformations and packing in a film and, therefore, its electronic properties. In particular, charge mobility is intimately related to film morphology. Understanding the link between morphology and mobility is therefore essential for improving electronic properties of polymeric films.

Because of slow polymer dynamics, conjugated polymers seldom reach global thermodynamic equilibrium even after annealing.[5] The resulting thin films are normally heterogeneous[6] with kinetically trapped regions of varying molecular order.[5−8] In amorphous regions, for example, side chains and backbones are completely disordered. In crystalline lamellae, formed by layers of cofacially stacked backbones alternating with layers of side chains, both side chains and backbones are crystalline.[9,10] Finally, in partially ordered domains, backbones are stacked cofacially but the lamellae do not form three-dimensional crystals.[8]

Heterogeneity of a polymeric film complicates the analysis of its morphology, since a single spectroscopic technique normally targets only a limited spatial resolution. Optical spectroscopies, for example, use wavelengths that are significantly larger than the molecular scale.[11] X-ray scattering[12,13] resolves the molecular-scale structure but provides only area-averaged information.[14] Scanning-probe techniques provide real-space imaging of morphologies on a molecular scale[14,15] but have limited resolution and imaging depth.[7] In most cases morphological analysis is complemented by a phenomenological model.[16]

In principle, models of morphology can be devised using computer simulations. The complications here are the time and length scales involved: morphological features as large as hundreds of nanometers[12] prohibit the use of all-atom simulations. Coarse-grained (CG) models reduce the computational cost, but their development and parametrization are far from trivial since they should quantitatively capture various local enthalpic and entropic contributions from side chains as well as “π”-stacking of backbones.

In coarse-grained models of semiconducting polymers, where CG particles correspond to small groups of atoms, molecular features that drive structure formation can be resolved explicitly.[17−23] Bonded and nonbonded CG interactions are chosen to reproduce atomistic distribution functions and hand-picked thermodynamic properties, e.g., solvation free energies. In Martini-based CG models,[24] local chain planarity is additionally incorporated by resolving aromatic rings.[21] Models with explicit side chains are capable of driving the system into a lamellar arrangement[19] and are useful to study polymer aggregation in solutions[20] as well as effects of polymer architecture[19] and processing[21] on morphology of heterojunction blends.

Despite the reduction of the degrees of freedom, relaxation times in these CG models become prohibitively large once the polymers are long enough to be entangled.[25] Topological constraints appear because these potentials conserve the (hard) excluded volume, similarly to all-atom descriptions. Another fundamental issue is that any CG description averages over many underlying microscopic states. Accordingly, CG potentials represent an atomistic potential of mean force[26−28] (PMF). This PMF is a complex many-body function of translational and orientational degrees of freedom, particularly in structured phases, and is usually unknown.[28] Therefore, CG potentials are only approximations to the actual PMF: typical examples are isotropic and pairwise CG potentials. In fact, the choice of a suitable potential function is not trivial, especially if the CG potential cannot be directly related to the structure of the mesophase.

One way to increase efficiency even further is to use coarser models where a single CG particle represents an entire repeat unit, including all its cyclic moieties and side chains.[29] Because many atomistic states contribute to a single CG configuration, the long-range part of the underlying PMF is soft, i.e., comparable to the thermal energy.[26,27,30] Therefore, it becomes possible to relax local excluded volume constraints and focus on the long-range packing of the repeat units by introducing soft nonbonded potentials. Their softness facilitates efficient sampling and boosts computational efficiency.

In a lamellar phase the PMF between entire atomistic repeat units is expected to be anisotropic. Accordingly, nonbonded potentials in coarser models must be anisotropic as well. One can, for example, approximate these potentials by analytic functions, such as the Gay–Berne potential, and parametrize them using atomistic simulations. This approach, however, retains the hard excluded volume.[31,32] In soft models, a more general route is to construct the anisotropic potentials using the symmetry of the molecular order in different mesophases.[33] Using this approach, we devise a generic soft model that can simulate amorphous, nematic, and partially ordered lamellar mesophases, even though cyclic moieties and side chains are not explicitly resolved.

We then use this model to study self-assembly and charge transport in poly(3-hexylthiophene) (P3HT). This material has been a fruit-fly system of polymeric organic semiconductors for decades,[34−36] but many fundamental questions remain unanswered. Because of that, P3HT is a perfect case for new theoretical approaches. For example, the link between charge mobility and morphology of P3HT or, in general, semiflexible conjugated polymers is still under debate.[37−39] It is generally assumed that the value of the mobility depends on the amount of amorphous and semicrystalline material in a film. In P3HT, the mobility in fully amorphous regions, ∼10–5 cm2/(V s), is due to interchain hopping. In partially crystalline P3HT, the mobility is 2 orders of magnitude higher but never reaches the magnitude typical for a crystalline polymer, ∼10–2 cm2/(V s). This reduction is due to random orientations of crystallites and boundaries between them.[40] Chain rigidity plays a crucial role here: stiffer chains have less conformational defects and lead to a higher overall mobility.[38,39]

We first perform Monte Carlo simulations of lamellar mesophases of experimentally relevant[41−43] chain lengths. These mesophases are compared with experimentally reported structures of partially ordered lamellae (phase III),[42] which belong to the general class of sanidic liquid-crystalline mesophases.[44,45] They can be used to drive soluble polymeric semiconductors into solid-state morphologies with improved electronic properties.[5,45,46] Subsequently, we simulate large-scale morphologies with different degrees of lamellar order and heterogeneity and use a qualitative charge-transport model[37−39] to link charge mobility to mesoscopic molecular organization.

Coarse-Grained Model

Degrees of Freedom

We base our CG scheme on a model developed for uniaxial and biaxial nematic liquid-crystalline (LC) mesophases.[33] Each P3HT monomer, that is, a thiophene ring and attached alkyl side chain, is represented by a single interaction site, located at the intersection of two lines extending from the bonds that connect the thiophene rings (see Figure a). A P3HT chain contains N such CG monomers, and the entire system is composed of n polymer molecules. The position vector of a CG site is denoted by r(s), where i = 1, ..., n is the chain index and s = 1, ..., N is the monomer index.

Figure 1

(a) Coarse-graining of atomistic P3HT (side chains are omitted for clarity). Each thiophene monomer is mapped on a single CG interaction site. The site has boardlike symmetry and is assigned three unit vectors {n((s)} (k = 1, 2, 3, while i and s are chain and monomer indices). Angular and torsional degrees of freedom, θ and ϕ, of the coarse-grained model are also shown (cf. eq ). The length of the CG bond is fixed to b. (b) Sketch of partial lamellar order considered in this work. Polymer backbones form stacks alternating with “empty” layers. They represent molten (noncrystalline) layers of hexyl side chains that are implicitly described in our model. r(s,t) is the interparticle vector connecting two CG sites. To clarify better the meaning of occupied and “empty” alternating layers, an atomistic P3HT molecule underlying a CG chain is shown.

Anisotropic nonbonded interaction potentials between CG sites require orientational CG degrees of freedom.[47] We describe the orientation of thiophene rings by three orthonormal unit vectors {n((s)} (k = 1, 2, 3), as shown in Figure a. In fact, the orientation of these vectors is fully determined by the local chain conformation. Indeed, if u(s) = r(s+1) – r(s) is a vector along the bond connecting the sth and (s + 1)th CG sites, then n(1)(s)||u(s) + u(s – 1), n(2)(s)||u(s) – u(s – 1), and n(3)(s) = n(1)(s) × n(2)(s) (cf. Figure a). The unit vectors at the two end monomers of a chain are defined by adding to each of them a fictitious site (indexed by s = 0 or N + 1) which does not introduce nonbonded interactions. The fictitious site is linked to the respective chain end via a “ghost” bond. These “ghost” bonds are subjected to the same bonded interactions as the normal CG bonds (see section ). On the basis of the local molecular frame, we associate[48,49] with each CG monomer a symmetric tensor with biaxial symmetry:

Bonded Interactions

The bonded interactions consist of angular and torsional potentials:[33]where θ denotes a bond angle and ϕ is a dihedral angle, as illustrated in Figure a.

The length of the bonds linking the CG monomers into a chain is fixed to b = 4 Å. ϕ = 180° of the dihedral angle corresponds to the trans conformation. The constants θ0 = 147.46°, kθ = 462.653 kJ/mol, c0 = 2.75248, c1 = −1.37645, c2 = −5.29397, c3 = 3.19667, c4 = 3.12177, and c5 = −2.41059, all in kJ/mol, are obtained by Boltzmann inversion of angular and dihedral probability distributions extracted from atomistic simulations of a single P3HT 20-mer in conditions mimicking implicitly Θ-solvent. To mimic these conditions, nonbonded interactions in atomistic simulations are active only if the participating atom pairs belong to monomers in the 1–2, 1–3, or 1–4 position.[33]

Nonbonded Interactions

We construct the nonbonded potential using the symmetry of the lamellar mesophase, which is sketched in Figure b. In lamellae of P3HT, these stacks alternate with layers of side chains. In our model the stacks of the backbones must be separated by empty space with “virtual” side chains. In other words, our lamellar mesophases will be characterized by strong density modulation along the direction of lamellar stacking (cf. Figure b). Maximum density will be observed in regions of stacked backbones, while minimum (almost zero) density will be found in between. Because side chains are not explicitly resolved, our model is more suitable for simulations of lamellae with noncrystalline side chains.[42,46,50]

The phenomenological nonbonded potential, Vnb, promoting the lamellar order has several contributions:Vnb acts between all pairs of CG sites, unless they belong to the same chain and are less than four bonds apart. This exclusion is consistent with the bonded CG potential, which incorporates the effect of atomistic interactions for these intramolecular pairs. Vnb is a sum of three interaction terms; their strength is set by non-negative parameters κ, λ, and ζ. Each term has a specific function and is motivated by simple arguments, as follows.

The isotropic repulsive potentialprovides finite compressibility and hence prevents the collapse of the liquid. This interaction can be used to model an isotropic polymer melt and depends only on the interparticle distance r ≡ r(s,t) = |r(s,t)| via the soft core functionThe Heaviside function Θ(r) sets the interaction range to 2σ. The physical meaning of Viso has been discussed previously.[33] It is proportional to the overlap integral of two spherical “clouds” with radius σ and uniform density w(r) = 3/4πσ3, centered on each CG site.[27,30] In an isotropic melt, these clouds approximate the spatial distribution of the atomistic degrees of freedom that were coarse-grained out. The repulsion between CG sites should begin at length scales where the underlying side chains come into contact. Therefore, we set σ = 7.6 Å, which is the length of a hexyl chain in the all-trans conformation. ρ0 stands for[33] some characteristic, reference, density and is seen here as a prefactor rescaling the interaction strength parameters. We use ρ0 = 4 nm–3, which is about the number density of monomers in crystalline P3HT.[9,51] Drastic coarse-graining implies[26,27,30] substantial overlap between CG P3HT monomers.[33] For the repulsive interaction to be soft, κU(0) should be comparable to the thermal energy, kBT, at a representative temperature T. Throughout this study we set κ/kBT = 15, leading to κU(0)/kBT ≃ 2. The temperature is set to T = 500 K, which is close to the melting point of P3HT.

The anisotropic potential[33,48,49,52−54]is defined via the Frobenius product b(s):b(t) of the biaxial tensors of two interacting CG sites. It is designed to promote biaxial nematic order of polymer chains. This mesophase lacks density modulation but already reproduces one feature of lamellae–parallel arrangement of planes of chain backbones (see Figure b). That is, in the biaxial nematic the local axes of the CG particles, n((s) and n((t) (k = 1, 2, 3), tend to be mutually parallel or antiparallel (we consider only nonpolar biaxial phases). Although Vbiaxial is only a special case[48,54] of a general quadrupolar pair potential between objects with D2 symmetry, it is sufficient to obtain biaxial nematic mesophases in simulations of polymers.[33] For simplicity, the biaxial potential has the same distance dependence as the repulsive interactions, U(r(s,t)).

The potentials Viso and Vbiaxial have been already incorporated into the coarse-grained model of ref (33). These interactions alone are not sufficient for obtaining lamellar order. Here, to enable the formation of lamellae, we add to Vnb a new interaction—the “stacking ” potential:Here r̂(s,t) = r(s,t)/r(s,t) is the unit interparticle vector and P2 is the second-order Legendre polynomial. We introduce this interaction considering that a sanidic lamellar mesophase, apart from the isotropic repulsion and biaxial order, requires two additional features. First, as shown in Figure b, polymer backbones must stack on top of each other. In a stack, CG particles tend to arrange face-to-face, so that n(3)(s)·r̂(s,t) = ±1 and n(3)(t)·r̂(s,t) = ±1. Second, density modulation must accompany stacking: in a side-by-side arrangement, where n(3)(s)·r̂(s,t) = 0 and n(3)(t)·r̂(s,t) = 0, the CG particles must repel each other strongly. The product of the isotropic core U(r(s,t)) with an anisotropic term of P2 symmetry promotes these two features.

As an illustration, Figures a–c show Vnb between two particles having perfect biaxial alignment with respect to each other. Without loss of generality, the axes of the particles {n(1), n(2), n(3)} are chosen parallel to the {X, Y, Z} axes of the laboratory frame, as sketched in Figure d. Figures a and 2b present contour plots of Vnb as a function of the magnitude of the components r⊥, r of the interparticle vector, r; r⊥ is the projection of the vector on the XY-plane, and r is the component parallel to the Z-axis. In Figure a only repulsive and biaxial interactions are active, κ/kBT = 15, λ/kBT = 7, and ζ = 0, while in Figure b these interactions are augmented with the stacking potential, ζ/kBT = 3.5. For the given configuration of particles, Vnb has cylindrical symmetry around Z and is fully characterized by the two-dimensional contours of Figures a and 2b. For the chosen values of interaction parameters, Vnb ≥ 0 everywhere in both plots. However, the symmetry of Vnb is different in the two cases. With only repulsive and biaxial terms active, Vnb is insensitive to the orientation of the interparticle vector, r, and the isosurfaces are spheres. Activating Vstack breaks the spherical symmetry. For a given interparticle distance r, particles repel the least when they arrange face-to-face, r⊥ = 0, and the most when they pack side-by-side, r = 0, as manifested by the “peanut-like” shape of the potential isosurfaces. To illustrate better the strength of interactions for these two specific particle arrangements, Figure c presents[55] a one-dimensional plot of Vnb as a function of r.

Figure 2

(a, b) Contour plots of the nonbonded potential Vnb as a function of the components r and r⊥ of the interparticle vector r (cf. scheme (d)). The plots are obtained, assuming that the two particles are oriented in a perfectly biaxial configuration, with their axes n(1), n(2), n(3) parallel respectively to the X, Y, Z axes of the laboratory frame. The interaction parameters are κ/kBT = 15, λ/kBT = 7 and (a) ζ/kBT = 0, (b) ζ/kBT = 3.5. White lines highlight selected isosurfaces of Vnb. (c) One-dimensional profiles of Vnb as a function of the interparticle distance r, for two limiting particle configurations: side-by-side (r = |r⊥| and r = 0, dashed line) and face-to-face (r = |r| and r⊥ = 0, solid line).

We have constructed Vnb using symmetry arguments. To connect qualitatively to the underlying molecular picture, it is helpful to consider the combinationwhere ξ = ζ/κ. For ξ < 0.5 the right-hand side of eq is positive and ζVstack perturbs the isotropic potential κViso, generating anisotropic repulsion. The combined interaction in eq mimics the average effect of steric interactions between side chains in a lamellar mesophase; the conformations of the side chains are very anisotropic because they extend outside the volume occupied by the stacked backbones. An approach for qualitatively connecting ξ to shape anisotropy of P3HT monomers is discussed in the Supporting Information. In principle, the effect of side chains can be captured more accurately, introducing a soft repulsion without the cylindrical symmetry of Figure b. This asymmetry would take into account that side chains are attached only to positions 3 and 4 of a thiophene ring. Using the spherically symmetric Viso only is acceptable for approximating effective repulsive interactions in isotropic and nematic (uniaxial or biaxial) mesophases where the liquid is more or less homogeneous.

Potentials with anisotropic terms P2(r̂·û) (here û is a generic molecular axis), analogous to those in Vstack, have been employed in studies of low-molecular-weight liquid crystals.[56−59] These studies, however, considered only particles with uniaxial symmetry and focused on objects with prolate shape, which typically form nematic or smectic A mesophases. Hence, the P2(r̂·û) terms were constructed to promote side-by-side configurations of particles.

Monte Carlo Sampling

We sample the configuration space with Monte Carlo (MC) simulations. Depending on the objectives, they are performed either in a canonical or an isostress ensemble. In both ensembles, the temperature T and the density of the system ρ = nN/V (V is the volume of the simulation cell) are fixed. We use orthorhombic simulation cells with edge lengths Lα (where α = X, Y, Z) and periodic boundary conditions (PBC) in all directions.

In the canonical ensemble, Lα are fixed. The configuration space of the system (translational and internal degrees of freedom of the chains) is sampled using the standard[60,61] “slithering snake”, reptation, MC move. The move has been adjusted[33] to the current CG model to account for “ghost” bonds.

The simulations in the isostress ensemble optimize Lα, making them commensurate with the natural geometry of the lamella morphology at the prescribed density ρ. “Natural geometry” refers to lengths characterizing the periodicity of a morphology in the bulk, free of any bias from PBC. In a lamella with natural geometry the stress acting on the simulation cell is isotropic;[62] i.e., the diagonal elements of the stress tensor should be equal. We implement a variable-shape-constant-volume (VSCV) MC algorithm. The method has been discussed extensively in the literature[63−70] so we provide only a summary, clarifying aspects specific to polymers. We apply the technique to systems where the directions of lamellar and “π”-stacking are parallel to the Y- and Z-axes, respectively. Chain backbones are oriented, on the average, along the X-axis. In a VSCV move, new edge lengths are proposed according to Lnew = Lold + ΔL and Lnew = Lold + ΔL. The random increments ΔL and ΔL are distributed uniformly in the interval [−ΔLmax, ΔLmax]. The new L follows from the constraint of constant volume: Lnew = V/LnewLnew. Subsequently, the new coordinates of the CG monomers, rnew(s), are obtained by applying an affine transformation only to the first monomer of every chain: rnew(1) = T·rold(1), i = 1, ..., n. Here T is a diagonal matrix with elements Tαα = Lαnew/Lαold. The coordinates of the remaining s = 2, ..., N monomers of each chain follow from rnew(s) = rnew(1) + ∑uold(t). In this scheme, uold(t) are the vectors of the N – 1 bonds of the ith chain in the old configuration. The move is accepted with probability: pacc = min (1, exp[−βΔ]) where βΔH = βΔUnb + ln(LnewLnew/LoldLold). Because the internal degrees of freedom of the chains are not affected by the move, only the difference of nonbonded energies, ΔUnb, in the new and the old configuration enters ΔH. The logarithmic term stems from a Jacobian associated with the constraint of constant volume. A typical mix of MC moves contains only 0.1% VSCV moves; the rest is given to reptation moves. This small fraction is motivated by the computational cost of VSCV moves—whenever a change of Lα is attempted, the nonbonded energy of the system must be calculated from scratch. Setting ΔLmax = b leads to an acceptance rate of VSCV moves of ∼2%.

We consider three cases of monodisperse systems, composed of molecules with N = 16, 24, and 32 monomers. These degrees of polymerization are comparable to low-molecular-weight polyalkylthiophenes in experiments.[42,43] For a first generic study, we consider only one temperature, i.e., T = 500 K. The number of molecules and volume of modeled systems are chosen such that ρ = ρ0 (for the definitions of ρ0 see section ).

Model Parametrization

We first investigate the phase behavior as a function of λ and ζ. In general, correlating quantitatively λ and ζ with molecular features is formidable because Vnb approximates a PMF which is not a conventional potential but a free energy. The ingredients of this free energy, such as the coefficients κ, λ, and ζ, are determined by subtle, unknown, entropic and enthalpic contributions from side chains and thiophene rings, underlying single interaction centers in our CG model. The difficulties in molecular-based interpretations of free-energy-like parameters have been illustrated, for example, by studies[71] estimating Maier–Saupe constants in polymer nematics. Although we provide in the Supporting Information some (very) qualitative arguments linking λ and ζ to molecular properties, the main body of our study considers them as phenomenological input parameters. In this section we heuristically identify a region of λ and ζ values where lamellae are formed. In the following, a representative combination of λ and ζ from this region will be used to study molecular organization and charge transport in the lamellar mesophase.

During the phenomenological development of Vnb, the biaxial and stacking terms were conceived to promote lamellar order synergistically. Therefore, we initially characterize the phase behavior for ζ = 0 and identify the range of λ for which the system exhibits biaxial order. Focusing on this range of λ, we activate the stacking potential to locate the region where lamellae are formed. We scan the (ζ, λ) plane, shown in Figure , only in the region where Vnb ≥ 0, for any choice of arguments. For cases where Vnb < 0 our soft potential can lead to an instability:[72,73] the system may collapse, gathering the molecules in a small region of space. The secure region Vnb ≥ 0 is defined through the condition λ + 2ζ – κ ≤ 0. The constraint follows from eq , requiring that Vnb ≥ 0 even when λVbiaxial + ζVstack is the most negative. This situation happens when two interacting particles are found in a perfectly biaxial, face-to-face registration, irrespective of their distance. In Figure , the boundary of the stability region is marked by the red dashed line.

Figure 3

Phase behavior as a function of the interaction parameters λ and ζ, for N = 16 and κ/kBT = 15. The boundary of the biaxial nematic–lamellar (NB → L) transition is shown with solid circles (the black solid line is a guide to the eye). The width of the error bars corresponds to the step Δζ/kBT = 0.25 used to scan the values of ζ. The red dashed line marks the boundary of the stability region of the nonbonded potential; the part of the (ζ, λ) plane above the stability line is not considered (see text for details). The region where lamellar structures are unambiguously formed is shaded with lines. The snapshots provide 2D XY-views of a biaxial morphology (λ/kBT = 7, ζ/kBT = 0) and a lamellar morphology (λ/kBT = 7, ζ/kBT = 3).

To probe phase behavior, we consider chains with N = 16 monomers. These molecules are sufficiently long to exhibit a biaxial nematic phase[33] and are, at the same time, short enough to allow fast exploration of the phase diagram. Simulations are performed in the NVT ensemble, using a cubic box with n = 512 chains and edge length L = 16.72σ. Two types of initial configurations are employed: (i) chains in all-trans conformation with perfect biaxial orientational but no positional order; (ii) chains with conformations drawn from the bonded-potential distribution and arranged randomly in the box, without any orientational or positional order.

We quantify orientational order in a standard way[74,75] by computing three molecular ordering tensors Q( (k = 1, 2, 3):where I is the unit matrix. Diagonalization of Q( provides a set of nine eigenvalues and eigenvectors. The largest of these eigenvalues is the major order parameter S; its eigenvector defines the principal phase director. The biaxiality of the phase is quantified via an additional order parameter B which is calculated from the remaining eigenvalues, eigenvectors, and respective ordering tensors. Details are available in several publications[33,74,75] and the Supporting Information.

Configurations with confirmed biaxial order are screened for lamellar order by observing distinct density modulations and quantifying long-range positional order. For this purpose, we use an intermolecular pair-correlation function[76,77] probing positional correlations along the n(2) particles axes, i.e., along the molecular axis pointing into the lamellar stacking direction (cf. Figure b). This function is defined asHere δ(r) is the Dirac function. The argument r∥, is the projection of the distance vector between a test particle and one of the surrounding particles on the vector n(2) of the test particle. As clarified in the inset of Figure a, surrounding particles contribute to g(r∥,) only when found in a cylindrical sampling region, with radius Rc = σ. Long-range order is revealed by oscillations in g(r∥,), extending over many molecular distances. Resolving positional correlations along a molecular axis, instead of a laboratory axis, enables us to detect positional order not only in a lamellar monodomain but also in multidomain morphologies. Introducing the cutoff radius, Rc, when considering particles in the direction perpendicular to n(2), reduces the effects of deformations, e.g., buckling, or fluctuations of lamellar layers, which can wash out density modulations. In addition to lamellar order, we probe positional order along the “π”-stacking direction using the correlation function g(r∥,). It is defined through an expression identical to eq , replacing n(2) by n(3). The radius of the sampling region is again Rc = σ (cf. inset of Figure b).

Figure 4

Intermolecular positional correlations resolved along the (a) n(2) and (b) n(3) particle axis. Different values of ζ are considered for fixed λ/kBT = 7. The plots are obtained from the simulations used to construct the phase diagram in Figure . The cartoons in the insets clarify the geometrical construction used to compute these functions (cf. eq ).

Figure reports the phase diagram. Without stacking interactions, ζ = 0, the system develops biaxial order at about λ/kBT = 5.25. Interestingly, the roughness of the orientational energy landscape of CG monomers produced by λVbiaxial for λ/kBT ≈ 5.25–7 is qualitatively comparable with estimates from ab initio calculations[78] for pairs of 3-methylthiophenes in a vacuum. More details are provided in the Supporting Information. We remind, however, that this comparison should not be taken too literally: λVbiaxial is a free energy that encapsulates entropic and enthalpic contributions when coarse-graining an interacting liquid. The (ζ, λ) plane above λ/kBT = 5.25 and below the stability line is scanned with resolution Δλ/kBT = 0.5 and Δζ/kBT = 0.25. For low values of ζ the system remains in a homogeneous biaxial state. By further increasing ζ, a biaxial–lamellar transition (NB → L) takes place. Solid circles mark the (ζ, λ) points at which a well-developed lamellar morphology is formed. Clearly, this boundary is only an approximation of the true thermodynamic phase boundary. The accurate determination of the latter within NVT simulations is not straightforward. Because we are working with a compressible model, phase coexistence might occur, complicating the analysis of the phase transition.[79] Moreover, uncertainties in the location of the biaxial–lamellar transition can result from a mismatch between the natural lamellar periodicity and the length of the edges of the simulation box. This effect is known, for example, for the order–disorder transition in block copolymer systems.[66] Our work focuses on lamellar morphologies, rather than on the phase transition itself. For our purposes the approximate phase boundary drawn in Figure is sufficient to identify a region where well-developed lamellar structures are unambiguously formed. This region is shaded in Figure with lines.

The snapshot in Figure (right) illustrates a lamellar morphology obtained for λ/kBT = 7 and ζ/kBT = 3. It consists of a single lamellar monodomain, where the chain long axes lie on average in the XY-plane and the backbone planes are on average perpendicular to the Z-axis. For the system size examined here, such monodomains are obtained even in simulations started from a disordered initial configuration. The tilt of the lamellar stacking direction with respect to the X- and Y-axes indicates a mismatch between the size of the simulation box and the equilibrium lamellar periodicity. For comparison, Figure (left) shows a homogeneous biaxial morphology obtained for the same value of λ but ζ = 0.

To illustrate how the “phase boundary” is determined in practice, Figure a reports the correlation function g(r∥,) obtained for λ/kBT = 7 and different values of ζ. For ζ = 0, g(r∥,) exhibits only a hump at r∥, ≈ 1.5σ. For ζ ≠ 0, oscillations start to appear; their amplitude and range grow with increasing ζ. The data suggest that long-range order sets in for ζ/kBT ≈ 2.5–3. Inspecting visually the morphologies obtained within this range of ζ, we conclude that well-developed lamellar order is clearly observed at ζ/kBT = 2.75. This value of ζ is chosen to mark the NB → L transition boundary for λ/kBT = 7. A similar procedure is applied to identify the other boundary points. The spatial modulation of g(r∥,) for these points is very close to that obtained at the boundary point λ/kBT = 7, ζ/kBT = 2.75. Figure b shows the correlation function g(r∥,) for the systems considered in Figure a. For ζ = 0, the shape of g(r∥,) is analogous to that of g(r∥,), in agreement with the spherical symmetry of the nonbonded potential (cf. Figure a). When the stacking potential is activated, the first-neighbor peak shifts to smaller distances, consistently with the decreased repulsion along the n(3) axis (cf. Figure b). Correlations become more pronounced, although their rapid decay indicates that positional order along the “π”-stacking direction remains rather short range. In the Supporting Information, we attempt to link qualitatively ζ, via the asymmetry parameter ξ = ζ/κ, to shape asymmetry of P3HT monomers, assuming that their conformations in the lamellar phase can be approximated by disc-like objects.

We observe in Figure that the NB → L transition boundary shifts to lower values of ζ as λ increases. However, the relationship is not linear: the boundary line becomes steeper as ζ decreases, suggesting that a minimum value of ζ is necessary to induce lamellar order. Another interesting result is that for our soft interactions no lamellar morphologies are formed in the presence of the stacking potential only, i.e., for λ = 0, even when ζ is just below the stability line. These observations indicate that the biaxial and stacking potentials work cooperatively to promote lamellar order.

In this study, we investigate the properties of lamellar mesophases described by our model on generic level. Therefore, we choose a representative set of values λ/kBT = 7 and ζ/kBT = 3.5, where the lamellar order is well developed. Indeed, Figure demonstrates that for N = 16 this pair of parameters is located well inside the lamellar phase. The lamellar order becomes even stronger for the longer N = 24 and 32 chains (see section ). The structural characterization in section will demonstrate that for this representative parametrization the lamellar and “π”-stacking distances are fairly close to experimental data. Taking into account the generic character of the study, we do not perform any parameter optimization to achieve more quantitative agreement.

Properties of Lamellae

Mesophase Structure and Chain Conformations

Before investigating structural and conformational properties, we optimize the geometry of lamellar monodomains with VSCV simulations. The procedure is described in the Supporting Information and is performed in such a way that the optimized monodomains have their lamellar stacking direction along the Y-axis, the “π”-stacking direction is along the Z-axis, and the chain long axes are along the X-axis of the laboratory frame. For all considered N, the monodomains contain n = 500 chains and their L is approximately twice as large as the end-to-end distance of a chain in the all-trans conformation. This condition avoids artifacts due to interactions between chains and their periodic images. To illustrate the monodomain orientation, Figure a presents an XY and an YZ view of an optimized configuration for a lamellar system formed by chains with N = 16 monomers.

Figure 5

(a) The snapshots illustrate an optimized lamellar morphology formed by chains with N = 16 monomers. Left and right panels correspond respectively to XY and YZ views of the same system. The plots in (b) and (c) present intermolecular positional correlations resolved along the n(2) and n(3) particle axis. In contrast to Figure , the plots are obtained from lamellar monodomains with optimized geometry equilibrated at λ/kBT = 7 and ζ/kBT = 3.5 for chains with N = 16 (black solid line) and N = 32 (red dashed line) monomers.

First, the strength of orientational order in the optimized monodomains is quantified via the major and the biaxial order parameters, S and B. We find high values of S and B, ranging from S ≈ 0.86 and B ≈ 0.82 for N = 16 to S ≈ 0.88 and B ≈ 0.85 for N = 32.

Figures b and 5c quantify positional order by presenting, respectively, the correlation functions g(r∥,) and g(r∥,), for two chain lengths: N = 16 (black solid lines) and 32 (red dashed lines). The long-range oscillatory behavior of g(r∥,) confirms the existence of long-range positional order along the lamellar stacking direction of our monodomains. From the distance of the peaks we estimate the optimum lamellar spacing: dlam ≈ 1.67σ = 12.7 Å. Similarly to the nonoptimized lamellae (cf. Figure b), the optimized monodomains have only short-range positional order along the “π”-stacking direction: the function g(r∥,) in Figure c exhibits peaks only at short distances. From the distance of the first two peaks we estimate the optimum π–π packing distance for our model: dπ ≈ 0.68σ = 5.2 Å.

Experimental studies of P3HT (for example, see refs (42 and 80)) typically report for lamellar spacing and “π”-stacking distance ≈16 and ≈3.8 Å, respectively. Considering the simple and phenomenological way our model is developed, dlam ≈ 12.7 Å and dπ ≈ 5.2 Å are reasonably close to these experimental observations. We expect that the geometrical characteristics of the lamellae can be brought closer to values reported in experiments by tuning the various parameters entering Vnb (including the characteristic length scale, σ).

Interestingly, the positional order increases slightly as the chains become longer (for fixed κ, λ, and ζ). In Figure b, the amplitudes of oscillations in g(r∥,) are a bit larger for N = 32 compared to the N = 16 monodomain. The trends in Figure c are consistent: the range and the amplitudes of oscillations in g(r∥,) are somewhat larger for N = 32 than for the N = 16 system. A plausible explanation for the increase of positional order with N is cooperativity effects, i.e., correlations in the order of monomers originating from chain connectivity and stiffness. These correlations play a role, for example, in uniaxial[71,81−83] and biaxial[33] nematic liquid crystalline polymers, where they cause a shift of the isotropic–nematic transition with chain length.

The snapshot in Figure a illustrates the organization of polymers in one representative layer of a lamellar monodomain composed of N = 16 chains. The chains are more or less planar. They are, on average, biaxially aligned, though some variations in orientation along the contour of the chains are clearly visible. There is no long-range order along the “π”-stacking direction, in agreement with the behavior of g(r∥,). Moreover, it is evident that there is no mutual registration of chains along their backbones, i.e., the X direction. Overall, the structure of a single lamella can be described as a quasi-2D biaxial nematic. The chains are extended, almost without hairpins (U-turns along the polymer contour). In addition to visual inspection, we compute cos(n(1)(s)·X), for each chain i and each site s on this chain. A hairpin is found when this quantity changes sign.[84] According to this analysis, only 0.1% of chains contain at least one backfold. The percentage of backfolding in monodomains with N = 24 and 32 chains is similar. Still, backfolding might become more significant for chains with higher molecular weights. We do not find any long-range correlations between the positions of chains belonging to two neighboring stacks (i.e., stacks that are beside each other along Y).

Figure 6

(a) XZ view of one representative layer in a lamellar monodomain. (b) Enlarged local XY view of three lamellar layers. The red dashed oval highlights a bridge connecting two neighboring lamellae. Both (a) and (b) are extracted from a snapshot of an optimized monodomain, composed of N = 16 chains (λ/kBT = 7, ζ/kBT = 3.5).

On the basis of all morphological features, we conclude that our monodomains belong to the class of sanidic, liquid-crystalline, smectic mesophases.[44,45] These mesophases, denoted as Σr, Σo, and Σd, have been initially described[44] for nonconjugated polymers (polyesters and polyamides). In all these mesophases, stacks of backbones alternate with layers of disordered side chains. Therefore, neighboring stacks of backbones are stochastically displaced with respect to each other and are uncorrelated. However, the mesophases differ in structure of individual stacks. In Σr, each stack is a two-dimensional crystal: there is long-range cofacial registration along the backbone axis and long-range order along the π-direction. In Σo, the long-range cofacial registration within each stack is lost, whereas the long-range order along the π-direction is maintained. In Σd there is neither long-range cofacial registration nor long-range order along the π-direction. Our monodomains correspond to the least ordered Σd mesophase. Experimentally, morphologies reproducing the phenomenology of the Σd mesophase have been observed[42] near 170 °C for P3HT molecules with lengths comparable to those in our simulations. In those experimental studies these mesophases were termed “phase III”.

Interlamellar Bridging

Figure b presents a close view of a monodomain with N = 16 chains and illustrates an interesting effect: chains, bridging neighboring stacks through the lamellar stacking direction. In terms of atomistically resolved P3HT lamellae, this effect corresponds to aromatic backbones bridging neighboring stacks by crossing the intermediate layer of hexyl side chains. Such bridges can provide conducting pathways through the insulating aliphatic region, so it is useful to quantify their number. For this purpose, for each monodomain configuration, we compute the density profile ρ(Y) along the lamellar stacking direction. A monomer is considered to be inside a lamella if located in a region where ρ(Y) is larger than the average density of the system and part of a bridge otherwise. For monodomains with N = 16 the average number of bridges per chain is found to be nbridge/n = 0.014, and their average length is Nbridge = 3 monomers. Similar values are obtained for N = 24 and N = 32. These results demonstrate that bridging segments, though present, are quite rare and relatively short. To the best of our knowledge, bridging through the side-chain layer, along the lamellar stacking direction, has not been explicitly addressed in experiments. Currently, it is difficult to say whether the bridging events are specific to our model or whether they indeed occur also in the actual material. It is worth mentioning that bridges are observed in other coarse-grained models,[85] though their formation has not been explicitly discussed.

A rough estimate suggests that bridges may be at least thermodynamically possible in the actual P3HT. We quantify the free-energy cost of a bridge, simply as Ubridge ≈ NbridgekBTχ. Here Nbridge is the number of thiophene monomers in the bridge, and χ is the Flory–Huggins interaction parameter describing the incompatibility between the backbone and the side chains. We evaluate χ through[86] Hildebrand solubility parameters, considering thiophene (T) and hexane (H) as the mixing species: kBTχ = (δT – δH)2/(ρTρH)1/2. The Hildebrand solubility parameters are given by δα, while ρα is the number density of the α compound (α = T, H). Setting Nbridge = 3, and using experimental data for the Hildebrand solubility parameters and densities (δT = 20.22 J1/2 cm–3/2, δH = 14.988 J1/2 cm–3/2, ρT = 7.58 nm–3, ρH = 4.595 nm–3),[87] the free-energy cost of a bridge is estimated to be Ubridge ≈ 8 kJ/mol. Even for rather low temperatures, e.g., T = 50 °C, this free-energy cost is moderate: Ubridge is only about 3 kBT.

Scattering Patterns

We calculate scattering patterns that can be qualitatively compared with GIWAXS experiments. Namely, we consider a scattering function with the general formTwo scattering geometries are addressed where (i) q∥ is oriented along the lamella stacking direction Y, while q⊥ lies in the XZ-plane and (ii) q∥ is oriented along the “π”-stacking direction Z, while q⊥ lies in the XY-plane. Accordingly, r(t) and r(t) are the projections of the position vector of the monomer on the (i) lamella stacking direction Y and XZ-plane and (ii) “π”-stacking direction Z and XY-plane. The Cartesian components of the vector q = q∥ + q⊥ comply with the PBC, i.e., qα = 2πm/Lα, where m is an integer. Angular brackets denote the canonical average.

The first scattering geometry is analogous to GIWAXS experiments on films with edge-on orientation and random in-plane distribution of lamellar domains. The snapshot on the left of Figure a illustrates this situation. The scattering function S(q⊥,q) is presented on the right of Figure a. A series of bright and sharp diffraction spots are visible along the q-axis, in analogy to observations made in experiments for morphologies with well-defined lamellar stacking.[8] The position of the spots corresponds to a periodicity dlam = 12.4 Å, which is consistent with the values extracted from the pair correlation function and the box size. Two additional scattering features are distinguished. The broad feature around q⊥ ≈ 1.7 Å–1 (corresponding to a distance of ≈3.7 Å) results from intramolecular scattering by monomers along the chain. The halo at q⊥ ≈ 1.0–1.3 Å–1 corresponds to a d-spacing of 4.8–6.3 Å, which is consistent with the “π”-stacking distance dπ = 5.2 Å estimated from the correlation functions.

Figure 7

2D scattering patterns computed for lamellar monodomains with optimized geometry (N = 32, λ/kBT = 7, ζ/kBT = 3.5), assuming (a) an edge-on and (b) a face-on orientation of the lamellae. The cartoons on the left of panels a and b clarify the setup. To mimic scattering from systems with domains oriented randomly in plane, the scattered intensities are azimuthally averaged (see rotation arrow in the cartoons).

The second scattering geometry mimics GIWAXS experiments on films with face-on orientation and random in-plane distribution of lamellar domains. An illustration is provided by the snapshot on the left of Figure b. The scattering function S(q⊥,q) is presented on the right of Figure b and allows us to resolve more clearly the scattering features from “π”-stacking. We observe a distinct spot at q ≈ 1.2 Å–1. However, the signal is diffuse, as expected for a system with very short-range positional order in the “π”-stacking direction.

Charge Transport

To simulate charge transport in partially ordered morphologies, we consider three types of events: (i) fast intrachain transitions which represent charge motion in conjugated segments, (ii) slower intrachain transport between conjugated segments, and (iii) even slower interchain hopping. In our transport model the rates of these three events are determined by the value of the respective values of the charge transfer integral J in the Marcus rate[88,89]which we fix to 1, 10–1, and 10–2 eV. For T = 300 K and the reorganization energy of Eλ = 0.1 eV,[35] we then get 2 × 1016, 2 × 1014, and 2 × 1012 s–1 for intraconjugated, interconjugated, and intermolecular rates. The smallest coupling is chosen to reproduce the “π”-stacking mobility of P3HT, which is ∼0.1 cm2/(V s).

Note that ΔE = F·r(s,t), where F is the external field, i.e., the energetic disorder and the dependence of electronic couplings on the local structure are not taken into account. Hence, we only capture a qualitative link between the morphology and the topology of the charge percolating network. Also note that the intraconjugated transport is not in a hopping regime, since the excess charge is delocalized in a conjugated segment. We model this delocalization by an effective (large) intraconjugated rate and use the (not applicable in this case) Marcus rate expression only to set physically meaningful prefactors (units) in the rate expression.

Each CG bead represents a site of a charge percolating network. Two beads in a chain belong to a conjugated segment if the dihedral angle between them deviates from the planar cis or trans conformations by less than ±45°. As we will see, conjugation breaking occurs only during the equilibration, even for the longest chains studied here, N = 32. P3HT chains are mostly extended, and the number of slow intrachain transitions due to conjugated breaks is small. The same conclusion is reached when we compare to the chain length used in experiments, where much longer chains are required to observe the charge localization effects.[90−92]

The interchain hops occur between two beads belonging to different chains separated by less than σ, which corresponds to the first minimum of the distribution function g(r∥,) in lamellar morphologies (see Figure c). An example of charge transport network is visualized in Figure for one lamella.

Figure 8

Charge transport network for a lamella in a monodomain morphology with N = 16. Bonds are color-coded according to the rate constants used: fast rates within a conjugated fragment (black), medium rates between conjugated fragments (red), and slow intermolecular hops (cyan). The selected layer is identical to the one in Figure a.

Charge dynamics on this network is modeled by solving the corresponding master equation for state occupation probabilities pusing the kinetic Monte Carlo (KMC) algorithm.[93,94] KMC trajectories are used to calculate the components of the mobility tensor[95−97]where F is the external field and v is the charge velocity. Note that in the absence of energetic disorder mobility does not depend on the applied field. The averaging ⟨...⟩ is performed over 128 trajectories with different initial positions of a charge. Charge transport simulations are performed using the VOTCA package.[98,99]

Chain-Length Dependence

First, we consider three monodomain lamellae with chains of N = 16, 24, and 32 repeat units. The corresponding eigenvalues of the mobility tensors are shown in Figure . As expected, they differ by about 1 order of magnitude, since the transport along the chains, perpendicular to the chains, and between the lamellae is governed by inter- and intramolecular electronic couplings. The largest value, μ1 ≈ 1 cm2/(V s), corresponds to the intrachain transport with the fastest rates. The transport in the “π”-stacking direction is via the slowest, intermolecular rates and has mobility of μ2 ≈ 0.1 cm2/(V s). The lowest mobility is observed for the interlamellar transport, since here we have only a few chains bridging the lamellae, as shown in Figure . This is a remarkable result: in spite of the fact that the transport between the neighboring lamellae is mitigated by only a few bridges, the transport in the direction perpendicular to the lamellae is only an order of magnitude smaller than the transport along the stacks. This effect would not be possible to observe in atomistic simulations, where a relatively small number of lamellar stacks is preassembled. The stacks are well separated by insulating alkyl chains already in the starting configuration. This arrangement does not change on the time scales of molecular dynamics simulations.

Figure 9

Components of the mobility tensor, calculated in monodomains, as a function of chain length.

The mobility along the direction of chains increases with the chain length, which is due to a larger number of intramolecular charge transfers. Along the “π”-stacking direction and perpendicular to the lamellae it does not change, at least for the chain lengths studied here.

Coarsening of Domains

We now look how mobility behaves during the coarsening of lamellar domains. Eigenvalues of the mobility tensor for selected snapshots along the MC trajectory are shown in Figure a, together with the squared end-to-end distance, Re2. At t < 105 MC steps, chain extension leads to an increase of Re2. The mobility tensor, however, stays isotropic, indicating that the lamellar domains are randomly oriented and are smaller than the simulation box (see Figure b).

Figure 10

Properties of an evolving polydomain system with N = 16, obtained from simulations of n = 4096 chains in a cubic box of edge length L = 25.42 nm (≈33.45σ): (a) The square of the end-to-end distance R̃e2 normalized by its value in the all-trans conformation, order parameter of the end-to-end vector Se, and eigenvalues of the mobility tensor μ. (b) Morphologies obtained at various stages of evolution. Color coding reflects the orientation of the chain end-to-end vector.

After 105 MC steps, the end-to-end distance begins to converge, and the preferential (collective) chain alignment starts to dominate, as can be seen from the rapid increase of the orientational order parameter. This is reflected in the increase of the largest eigenvalue of the mobility tensor. In this regime, the well-formed lamellar domains coarsen, their average size approaches the box size, and a monodomain structure is formed in the simulation box at ∼106 MC steps. In this state, the end-to-end distance and mobilities saturate, matching the respective values of the monodomain as shown in Figure . Interestingly, the interlamellar component is still higher than that of a monodomain, implying that the equilibration of bridging chains requires even longer simulation times.

Summary and Outlook

We have developed a generic coarse-grained model that can be used to simulate amorphous, nematic, and partially ordered lamellar mesophases of polymeric semiconductors. The polymer architecture is described using a hindered-rotation chain model, where a single interaction site represents an entire atomistic repeat unit. The bonded potentials are defined such that coarse-grained and atomistic angular and dihedral distribution functions match in Θ-solvent conditions.[33] Soft anisotropic nonbonded interactions are introduced based on symmetries of molecular order in a lamellar mesophase.

We parametrize this model for the conjugated polymer P3HT and use Monte Carlo simulations to equilibrate monodomains of partially ordered lamellae. Structural analysis shows that these lamellae form smectic sanidic mesophases[44] in which disordered side chains alternate with stacks of backbones. In each stack, there is no long-range order along the “π”-stacking direction and no regular cofacial registration between backbones along their long axes. Therefore, these lamellae correspond to the phase III, reported experimentally in low-molecular-weight P3HT.[42]

We observe that the neighboring layers of stacked P3HT backbones can be bridged by chains. An estimate based on solubility parameters of thiophenes and hexane suggests that these sparse bridges are thermodynamically allowed. They connect adjacent P3HT layers of stacked backbones through the insulating hexyl side chains.

Using a simple charge-transport model, where three distinct charge transfer rates represent charge delocalization in conjugated segments, intrachain charge transfer between conjugated segments, and interchain charge hopping, we simulated charge carrier mobilities in lamellar monodomains. In agreement with the rates introduced in the model, systems with longer chains have higher mobilities along backbones. Mobilities in the orthogonal directions turn out to be independent of chain length. In real polymer systems, mobility along the chains saturates with the increase of the chain length. In our systems we do not reach this regime, even for chains of 32 repeat units, since the majority of chains are completely conjugated.

Subsequently, we modeled systems with evolving molecular order, from amorphous to lamellar, and observed two distinct regimes. First, chains extend and form small lamellar domains. The mobility in this regime stays isotropic and does not increase as the system order increases. The value of the mobility is comparable to the “π”-stack mobility in a lamellar arrangement; i.e., it is defined by the intermolecular charge transfer rates. The domains slowly coarsen and eventually reach the boundaries of the simulation box. The mobility tensor becomes anisotropic once there is only one lamellar domain in the box. In this lamellar domain the mobility perpendicular to the lamellar layers is only 2 orders of magnitude smaller than along the chains. This is remarkable since this type of transport is mitigated only by few chains bridging neighboring lamellar layers. Interestingly, the average mobility, which qualitatively describes an average over domain orientations, increases upon domain coarsening. This is due to the much faster transport along the chains, which dominates the average in the lamellar mesophase. This rationalizes the need of stiff chains: apart from better bridging of crystalline regions, stiffness also ensures faster average mobility in larger crystalline regions, since the mobility along the chains dominates the average. Note, however, that our polydomains are liquid-crystalline, and their structure is not exactly equivalent to semicrystalline morphologies.[6,45]

As an outlook, we would like to comment on the dynamical behavior of the system since experimentally the morphologies are kinetically trapped or are far from equilibrium. In our Monte Carlo simulations we use the reptation as well as a variable-shape constant-volume move. The reptation has the advantage of efficient sampling and is straightforward to implement. Replacing reptation by other ergodic moves will not change the equilibrium properties of lamellae. However, the molecular organization of nonequilibrium polydomain morphologies may change. It would be interesting to see how other types of MC moves, e.g., the crankshaft move,[61] which leads to Rouse-like pseudodynamics, will affect the evolution of the morphology.

Our charge transport model could be further improved by reintroducing the atomistic details. In an actual material, charge transport is not only determined by the mesoscopic arrangement of chains, but also by the local atomistic structure. This structure influences both the density of states and electronic couplings. Our current charge-transport model neglects these effects, in line with the simplified microstructure: in soft models, hard excluded volume constraints are relaxed and coarse-grained particles can overlap. By reinserting atomistic details into a morphology generated by a soft model,[100] one could achieve more rigorous description of charge-transport.[35,101]