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 threedimensional 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][27][28]. 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 cm 2 /(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 cm 2 /(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 largescale 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
2.1. Degrees of Freedom. We base our CG scheme on a model developed for uniaxial and biaxial nematic liquidcrystalline (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 1a). 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 i (s), where i = 1, ..., n is the chain index and s = 1, ..., N is the monomer index.
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 i (k) (s)} (k = 1, 2, 3), as shown in (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 2). 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 ij (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.  Figure 1a. In fact, the orientation of these vectors is fully determined by the local chain conformation. Indeed, if u i (s) = r i (s+1) − r i (s) is a vector along the bond connecting the sth and (s + 1)th CG sites, then n i Figure 1a). 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 2.2). 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 1a.
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, c 0 = 2.75248, c 1 = −1.37645, c 2 = −5.29397, c 3 = 3.19667, c 4 = 3.12177, and c 5 = −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 2.3. Nonbonded Interactions. We construct the nonbonded potential using the symmetry of the lamellar mesophase, which is sketched in Figure 1b. 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 1b). 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, V nb , promoting the lamellar order has several contributions: V nb 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. V nb 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 potential ij iso (4) provides 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 ij (s,t) = |r ij (s,t)| via the soft core function The Heaviside function Θ(r) sets the interaction range to 2σ. The physical meaning of V iso 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 alltrans 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, k B T, at a representative temperature T. Throughout this study we set κ/k B T = 15, leading to κU(0)/ k B T ≃ 2. The temperature is set to T = 500 K, which is close to the melting point of P3HT. The anisotropic potential 33 is defined via the Frobenius product b i (s):b j (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 1b). That is, in the biaxial nematic the local axes of the CG particles, n i (k) (s) and n j (k) (t) (k = 1, 2, 3), tend to be mutually parallel or antiparallel (we consider only nonpolar biaxial phases). Although V biaxial is only a special case 48,54 of a general quadrupolar pair potential between objects with D 2h 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 ij (s,t)).
The potentials V iso and V biaxial 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 V nb a new interactionthe "stacking " potential: Ur s t P s s t P t s t n r n r ( ( , )) ( ( ) ( , )) ( ( ) ( , )) Here rî j (s,t) = r ij (s,t)/r ij (s,t) is the unit interparticle vector and P 2 is the second-order Legendre polynomial. We introduce this interaction considering that a sanidic lamellar mesophase, apart

Macromolecules
Article from the isotropic repulsion and biaxial order, requires two additional features. First, as shown in Figure 1b, polymer backbones must stack on top of each other. In a stack, CG particles tend to arrange face-to-face, so that n i (3) (s)·rî j (s,t) = ±1 and n j (3) (t)·rî j (s,t) = ±1. Second, density modulation must accompany stacking: in a side-by-side arrangement, where n i (3) (s)·r̂i j (s,t) = 0 and n j (3) (t)·r̂i j (s,t) = 0, the CG particles must repel each other strongly. The product of the isotropic core U(r ij (s,t)) with an anisotropic term of P 2 symmetry promotes these two features.
As However, the symmetry of V nb is different in the two cases. With only repulsive and biaxial terms active, V nb is insensitive to the orientation of the interparticle vector, r, and the isosurfaces are spheres. Activating V stack breaks the spherical symmetry. For a given interparticle distance r, particles repel the least when they arrange face-toface, r ⊥ = 0, and the most when they pack side-by-side, r Z = 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 2c presents 55 a one-dimensional plot of V nb as a function of r.
We have constructed V nb using symmetry arguments. To connect qualitatively to the underlying molecular picture, it is helpful to consider the combination where ξ = ζ/κ. For ξ < 0.5 the right-hand side of eq 8 is positive and ζV stack perturbs the isotropic potential κV iso , generating anisotropic repulsion. The combined interaction in eq 8 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 2b. 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 V iso 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 P 2 (r·u) (here ûis a generic molecular axis), analogous to those in V stack , 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 P 2 (r·u) terms were constructed to promote side-byside 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- Figure 2. (a, b) Contour plots of the nonbonded potential V nb as a function of the components r Z 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 κ/k B T = 15, λ/k B T = 7 and (a) ζ/k B T = 0, (b) ζ/k B T = 3.5. White lines highlight selected isosurfaces of V nb . (c) One-dimensional profiles of V nb as a function of the interparticle distance r, for two limiting particle configurations: sideby-side (r = |r ⊥ | and r Z = 0, dashed line) and face-to-face (r = |r Z | and r ⊥ = 0, solid line).

Macromolecules
Article 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 The new L X follows from the constraint of constant volume: . Because the internal degrees of freedom of the chains are not affected by the move, only the difference of nonbonded energies, ΔU nb , 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 ΔL max = 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 lowmolecular-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 2.3).

MODEL PARAMETRIZATION
We first investigate the phase behavior as a function of λ and ζ. In general, correlating quantitatively λ and ζ with molecular features is formidable because V nb 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-energylike 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 V nb , 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  3, only in the region where V nb ≥ 0, for any choice of arguments. For cases where V nb < 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 V nb ≥ 0 is defined through the condition λ + 2ζ − κ ≤ 0. The constraint follows from eq 3, requiring that V nb ≥ 0 even when λV biaxial + ζV stack 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  3, the boundary of the stability region is marked by the red dashed line.
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 bondedpotential 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) (k = 1, 2, 3): where I is the unit matrix. Diagonalization of Q (k) 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

Macromolecules
Article 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 1b). This function is defined as , c , , (2) Here δ(r) is the Dirac function. The argument r ∥,n (2) 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 4a, surrounding particles contribute to g(r ∥,n (2) ) only when found in a cylindrical sampling region, with radius R c = σ. Long-range order is revealed by oscillations in g(r ∥,n (2) ), 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, R c , 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 ∥,n (3) ). It is defined through an expression identical to eq 10, replacing n (2) by n (3) . The radius of the sampling region is again R c = σ (cf. inset of Figure 4b). Figure 3 reports the phase diagram. Without stacking interactions, ζ = 0, the system develops biaxial order at about λ/k B T = 5.25. Interestingly, the roughness of the orientational energy landscape of CG monomers produced by λV biaxial for λ/ k B T ≈ 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: λV biaxial is a free energy that encapsulates entropic and enthalpic contributions when coarse-graining an interacting liquid. The (ζ, λ) plane above λ/k B T = 5.25 and below the stability line is scanned with resolution Δλ/k B T = 0.5 and Δζ/k B T = 0.25. For low values of ζ the system remains in a homogeneous biaxial state. By further increasing ζ, a biaxial−lamellar transition (N B → L) takes place. Solid circles mark the (ζ, λ) points at which a welldeveloped 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 3 is sufficient to identify a region where well-developed lamellar structures are unambiguously formed. This region is shaded in Figure 3 with lines.
The snapshot in Figure 3 (right) illustrates a lamellar morphology obtained for λ/k B T = 7 and ζ/k B T = 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 3 (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 4a reports the correlation function g(r ∥,n (2) ) obtained for λ/k B T = 7 and different values of ζ. For ζ = 0, g(r ∥,n (2) ) exhibits only a hump at r ∥,n (2) ≈ 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 ζ/k B T ≈ 2.5−3. Inspecting visually the morphologies obtained within this range of ζ, we conclude that welldeveloped lamellar order is clearly observed at ζ/k B T = 2.75. This value of ζ is chosen to mark the N B → L transition boundary for λ/k B T = 7. A similar procedure is applied to identify the other boundary points. The spatial modulation of g(r ∥,n (2) ) for these points is very close to that obtained at the boundary point λ/k B T = 7, ζ/k B T = 2.75. Figure 4b shows the correlation function g(r ∥,n (3) ) for the systems considered in Figure 4a. For ζ = 0, the shape of g(r ∥,n (3) ) is analogous to that of g(r ∥,n (2) ), in agreement with the spherical symmetry of the nonbonded potential (cf. Figure 2a). 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 2b). 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 3 that the N B → 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 λ/k B T = 7 and ζ/k B T = 3.5, where the lamellar order is well developed. Indeed, Figure 3 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 5.1). The structural characterization in section 5.1 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.

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 X 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 5a presents an XY and an YZ view of an optimized configuration for a lamellar system formed by chains with N = 16 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 Figures 5b and 5c quantify positional order by presenting, respectively, the correlation functions g(r ∥,n (2) ) and g(r ∥,n (3) ), for two chain lengths: N = 16 (black solid lines) and 32 (red dashed lines). The long-range oscillatory behavior of g(r ∥,n (2) ) 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: d lam ≈ 1.67σ = 12.7 Å. Similarly to the nonoptimized lamellae (cf. Figure 4b), the optimized monodomains have only short-range positional order along the "π"-stacking direction: the function g(r ∥,n (3) ) in Figure 5c exhibits peaks only at short distances. From the distance of the first two peaks we estimate the optimum π−π packing distance for our model: 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, d lam ≈ 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 V nb (including the characteristic length scale, σ).
Interestingly, the positional order increases slightly as the chains become longer (for fixed κ, λ, and ζ). In Figure 5b, the amplitudes of oscillations in g(r ∥,n (2) ) are a bit larger for N = 32 compared to the N = 16 monodomain. The trends in Figure 5c are consistent: the range and the amplitudes of oscillations in g(r ∥,n (3) ) 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 6a 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 ∥,n (3) ). 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 i (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 longrange correlations between the positions of chains belonging to two neighboring stacks (i.e., stacks that are beside each other along Y).
On the basis of all morphological features, we conclude that our monodomains belong to the class of sanidic, liquidcrystalline, 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 twodimensional 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". 5.2. Interlamellar Bridging. Figure 6b 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 n bridge /n = 0.014, and their average length is N bridge = 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 U bridge ≈ N bridge k B Tχ.
Here N bridge 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: k B Tχ = (δ 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 N bridge = 3, and using experimental data for the Hildebrand solubility parameters and densities (δ T = 20.22 J 1/2 cm −3/2 , δ H = 14.988 J 1/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 U bridge ≈ 8 kJ/ mol. Even for rather low temperatures, e.g., T = 50°C, this free-energy cost is moderate: U bridge is only about 3 k B T.
5.3. Scattering Patterns. We calculate scattering patterns that can be qualitatively compared with GIWAXS experiments. Namely, we consider a scattering function with the general form Two scattering geometries are addressed where (i) q ∥ is oriented along the lamella stacking direction Y, while q ⊥ lies in

Macromolecules
Article the XZ-plane and (ii) q ∥ is oriented along the "π"-stacking direction Z, while q ⊥ lies in the XY-plane. Accordingly, r j∥ (t) and r j⊥ (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 inplane distribution of lamellar domains. The snapshot on the left of Figure 7a illustrates this situation. The scattering function S(q ⊥ ,q Y ) is presented on the right of Figure 7a. A series of bright and sharp diffraction spots are visible along the q Y -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 d lam = 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.
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 7b. The scattering function S(q ⊥ ,q Z ) is presented on the right of Figure 7b and allows us to resolve more clearly the scattering features from "π"-stacking. We observe a distinct spot at q Z ≈ 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 × 10 16 , 2 × 10 14 , and 2 × 10 12 s −1 for intraconjugated, interconjugated, and intermolecular rates. The smallest coupling is chosen to reproduce the "π"-stacking mobility of P3HT, which is ∼0.1 cm 2 /(V s).
Note that ΔE = F·r ij (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 ∥,n (3) ) in lamellar morphologies (see Figure 5c). An example of charge transport network is visualized in Figure 8 for one lamella.
Charge dynamics on this network is modeled by solving the corresponding master equation for state occupation probabilities p i p k pk i j j ji i ij (13) using the kinetic Monte Carlo (KMC) algorithm. 93

Macromolecules
Article 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 6.1. 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 9. 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 cm 2 /(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 cm 2 /(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 6. 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.
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.
6.2. 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 10a, together with the squared end-to-end distance, R e 2 . At t < 10 5 MC steps, chain extension leads to an increase of R e 2 . The mobility tensor, however, stays isotropic, indicating that the lamellar domains are randomly oriented and are smaller than the simulation box (see Figure  10b).
After 10 5 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 ∼10 6 MC steps. In this state, the end-to-end distance and mobilities saturate, matching the respective values of the monodomain as shown in Figure 9. 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

Macromolecules
Article 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.

Macromolecules
Article 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 ■ ASSOCIATED CONTENT

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.8b01863.
Qualitative connection of parameters λ and ζ to molecular features; procedure used to optimize the geometry of lamellar monodomains with VSCV simulations; procedure for determining orientational order parameters (PDF)