Densely Packed Semiflexible Macromolecules in a Rigid Spherical

Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

Densely Packed Semiflexible Macromolecules in a Rigid Spherical Capsule Andrey Milchev,† Sergei A. Egorov,§ Daniel A. Vega,∥ Kurt Binder,‡ and Arash Nikoubashman*,‡ †

Institute for Physical Chemistry, Bulgarian Academia of Sciences, 1113 Sofia, Bulgaria Institute of Physics, Johannes Gutenberg University Mainz, Staudingerweg 7, 55128 Mainz, Germany § Department of Chemistry, University of Virginia, Charlottesville, Virginia 22901, United States ∥ Department of Physics, Universidad Nacional del Sur-IFISUR-CONICET, 8000 Bahía Blanca, Argentina ‡

S Supporting Information *

ABSTRACT: The ordering of semiflexible polymers with persistence length lp and contour length L confined in a sphere of radius R is studied by molecular dynamics simulations of a coarse-grained model. Monomer densities are chosen where the corresponding bulk lyotropic solution or melt is a well-ordered nematic, and purely repulsive walls of the rigid confining sphere are considered. It is found that polymers close to the walls are bent according to the curvature of the confining spheres with all their monomers in a few layers parallel to the sphere surface, whereas the remaining macromolecules closer to the sphere center have one chain end and their center of mass far from the surface. The latter chains are responsible for the average nematic order of the system for sufficiently stiff chains. If L exceeds R, a single transition from isotropic to nematic states is found when lp increases. However, if R is close to an even multiple of L, a second transition occurs where chain ends are enriched near the equatorial plane (and two further planes parallel to it, if R ≈ 2L); i.e., the system develops a distorted smectic structure. The order in the surface shell then is characterized by topological defects.

1. INTRODUCTION Semiflexible macromolecules find various applications as versatile materials, in particular due to their possible liquid crystalline order,1,2 and are also important constituents of living matter.3,4 On a coarse-grained level, the most important characteristic of a semiflexible polymer is its persistence length lp,5−8 characterizing its stiffness (although we note that the precise definition and characterization of lp is problematic9,10). In particular for biopolymers, lp is typically much larger than the size of monomeric units, ranging, for example, from about 50 nm for double-stranded (ds) DNA to about 10 μm for filamentous (F) actin; these linear dimensions are comparable to cell sizes or sizes of nanofluidic/microfluidic devices.11 Hence, confinement of semiflexible macromolecules in various geometries has found widespread attention, particularly for the case of cylindrical11,12 or slit pores.13−24 Also, the confinement of single semiflexible polymers inside spheres has been extensively studied,25−33 e.g., with the motivation of understanding the packaging of very long dsDNA in bacteriophage capsids. However, much less attention has been paid to the confinement of many (shorter) semiflexible macromolecules inside of spherical containers.34,35 These systems may be of interest for problems such as the self-organization of actin filaments in cell-size confinement,36 drug delivery from spherical vesicles,37,38 etc. Indeed, experiments on the organization of DNA strands or DNA−actin filament mixtures inside of vesicles have been made.39,40 In this work, we address the challenging theoretical problem to understand the conformation of the semiflexible macro© XXXX American Chemical Society

molecules in spherical confinement, when the persistence length lp, the contour length L, and the sphere radius R are all of the same order. We recall that linear dimensions of bacteriophage capsids are in the range from 20 to 120 nm,28,30 while vesicles and cells are in the micrometer size range. To capture the general underlying physics of this problem, we employ a generic model instead of describing a specific experiment. Here, we model the confinement through a completely rigid sphere, neglecting the possibility that the rather stiff polymers may deform the confining membranes.41 We further ignore the possibility of attractive interactions between the confining surface and the monomeric units, although for slit pores it is known that in this case interesting ordering phenomena occur at the confining walls.23,24 We assume also that the confined polymers are dissolved in a good (implicitly modeled) solvent, as in our previous work concerned with the properties of such solutions in the bulk;20,42,43 in this description, the effective interactions between monomeric units are purely repulsive excluded volume type forces, and hence conformations of the chains and their possible orientational ordering are entirely controlled by entropy. Previously, we have studied this problem for scaled number densities of monomeric units ρ ≤ 0.2, where for the chosen parameters the bulk solution still is in the isotropic phase.34 It Received: December 13, 2017 Revised: February 15, 2018

A

DOI: 10.1021/acs.macromol.7b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

2. MODEL Following previous work,20−24,34,35,42,43 the macromolecules are described by a coarse-grained model, where each chain is a sequence of N spherical beads connected by springs. All pairs of beads (i, j) at distance rij interact through the purely repulsive Weeks−Chandler−Andersen (WCA) potential54

was shown that depending on chain stiffness, bond orientations perpendicular to the surface are significantly suppressed, and rather nonuniform distributions of the center of mass (CM) of the chains result. In the present paper, we provide a complementary study for a densely packed system, i.e., ρ = 0.7, for which the bulk solution is deep in the nematic phase.42,43 However, the constraining hard spherical surface creates a conflict with global nematic order which requires uniform orientation of the director field.44 For short rigid molecules (e.g., modeled by hard rods or spherocylinders45) at average densities close to the onset of nematic order in the bulk, tactoids were observed, i.e., spindle-shaped nematic droplets46 surrounded by the disordered phase; at high densities, a strongly frustrated nematic order with topological defects at the poles (bipolar orientational order) was found.45 These findings, as well as earlier studies47 of a lattice model for liquid crystals inside a spherical cavity, can be discussed in the framework of the Frank−Oseen continuum theory for the elasticity of liquid crystals,44,48,49 since the radius of the confining sphere is much larger than the length of the rods. For semiflexible chains, however, a description solely in terms of the nonuniform director field n(r) cannot capture the entire physical picture, as the contour length L and persistence length lp of the chains are significantly larger than any linear dimension of small molecules. Rather, the connectivity of the macromolecules and the arrangement of their chain ends introduce new facets into the problem. For instance, now also the case that R and lp are of the same order is physically meaningful. In this work, we shall present results of molecular dynamics (MD) simulations,50−52 combined with density functional theory (DFT).21,34,43 While in principle MD yields numerically exact results (apart from statistical errors) for the chosen model, in practice one must face the problem that for certain conditions the system may get trapped in long-lived metastable states, which are (almost) indistinguishable from the equilibrium states, as indeed noted in our preliminary communication on this subject.35 In general, it is possible to distinguish metastable states from stable ones by identifying the configuration which minimizes the free energy, but estimation of the free energy within a MD framework requires cumbersome thermodynamic integration procedures.50,53 In contrast, DFT is based on formulating the free energy as a functional of appropriate inputs and therefore would be well suited for this problem. However, minimizing the free energy functional using the inhomogeneous distribution of the molecules ρp(r,ω), where ω describes the local orientation of the bond vectors at r, is not practically feasible here. Thus, we have circumvented this problem by using the necessary information on ρp(r,ω) from the MD simulations as an input for computing the free energy using DFT. The technical details on this DFT-MD combination will be discussed in the Appendix, while we briefly summarize the model under study in section 2. Section 3 then describes results for the case when the radius of the confining sphere, R, is comparable to the chain contour length L. In section 4 we then discuss how the chain arrangement in the densely packed spheres changes when R significantly differs from L. Finally, we present our conclusions and outlook in section 5. Additional details on the equilibration procedure and data analysis are provided in the Supporting Information.

⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ 1/6 ⎪ ⎪ 4ε⎢⎜⎝ ⎟⎠ − ⎜⎝ ⎟⎠ ⎥ + ε , r ≤ 2 σ r r ⎦ UWCA(r ) = ⎨ ⎣ ⎪ ⎪ 0, r > 21/6σ ⎩

(1)

where the strength ε = 1 sets the energy scale and the bead diameter σ = 1 sets the length scale. The spring potential then is taken as the finitely extensible nonlinear elastic (FENE) potential55 ⎧ 2 ⎤ ⎡ ⎪− 1 kr0 2 ln⎢1 − r ⎥ , r ≤ r0 UFENE(r ) = ⎨ 2 r0 2 ⎦ ⎣ ⎪ r > r0 ⎩∞ ,

(2)

with k = 30 ε/σ and r0 = 1.5σ. In our model, chain stiffness is controlled by the bending potential 2

Ubend(Θijk ) = κ(1 − cos Θijk )

(3)

where Θijk is the bond angle formed between two subsequent bond vectors ai = rj − ri and aj = rk − rj with j = i + 1 and k = j + 1. The bending stiffness is controlled through the parameter κ (fully flexible chain for κ = 0). At a temperature T = 1, this model yields an effective bond length lb ≈ 0.97σ. The contour length of the chain is simply given by L = (N − 1)lb, while its persistence length lp is defined in terms of ⟨cos Θijk⟩ as10 lp lb

=

−1 ln⟨cos(Θijk )⟩

(4)

which simplifies to lp/lb ≈ κ/kBT for κ ≳ 2kBT. While eqs 1−4 are convenient in the framework of MD simulations, for the formulation of DFT it is more convenient to use a tangent hard sphere chain model.56 Then the contour length is slightly larger, L = (N − 1)σ, since lb = σ = 1, but eqs 3 and 4 still apply. For the MD simulations, we used the HOOMD-blue software package57,58 on graphics processing units (GPUs), allowing the study of systems containing on the order of 1 million beads. Here, the standard velocity-Verlet algorithm is employed to integrate the equations of motion.50−52 We employed a time step of Δt = 0.01τMD, with intrinsic MD time unit τMD = mσ 2/ε = 1, m = 1 being the mass of the monomers. We carried out runs up to an observation time of 107τMD and verified that the systems reached equilibrium before we took our measurements. The temperature was kept constant using a Langevin thermostat with a friction coefficient of ξ = 0.25. The wall−particle interaction was modeled in analogy to eq 1, but the interparticle distance r is replaced with r′ = R − r in spherical coordinates, where the origin of the system coincides with the center of the sphere. The scaled average density of the beads in the system is ρ = 35σ 3/(4πR3)

(5)

where 5 is the total number of beads inside the sphere. B


Article

Macromolecules In previous work,42,43 we have studied the phase diagram of this model in the bulk for various choices of N and lp. Varying the density ρ, the isotropic−nematic phase transition has been located. For the choices of N and lp considered here, the chosen density ρ = 0.7 would always correspond to a well-ordered nematic phase in the bulk (see Figure S1 in the Supporting Information). In contrast, in spherical confinement no welldeveloped phase transition and no simple order can any longer be expected.

3. CONFINEMENT OF CHAINS WITH N = 32 IN SPHERES OF RADIUS R = 35 When the spherical cavity is filled with fully flexible chains (κ ≤ 4) with N = 32 monomers, then the sphere interior is found to be homogeneously occupied with polymers. Further, the conformation of the polymers differs only marginally from the bulk situation, since their radius of gyration, Rg ≈ 2.9, is much smaller than the sphere radius, R = 35, and thus the majority of the enclosed chains do not feel the confinement. However, the situation changes drastically as the chain stiffness κ is increased: in the central region of the sphere, the ordering of the chains transforms from isotropic (κ ≤ 4) to biaxial nematic (8 ≤ κ ≤ 20) to uniaxial nematic (κ ≥ 24). At the same time, the surface region of the sphere goes from a disordered phase to a bipolar structure and then to a quadruploar “tennisball” like texture. In the ordered phases, chain ends get strongly enriched in the equatorial plane of the spherical container. Figure 1 shows exemplary simulation snapshots taken for a system with κ = 96, highlighting the apparent separation between the polymers in the central bulk-like region of the sphere and the chains close to the confining surface. Further, we can clearly identify the equatorial plane occupied by the chain ends in the sphere center. In the following, we will provide a detailed discussion of these effects and relate the macroscopic ordering to the microscopic conformations of the individual polymers. In order to quantify the arrangement of the (semiflexible) macromolecules, we computed the bond orientational order in the system. Denoting the unit vector from monomer i to monomer i + 1 in the nth chain as uni, the local bond orientational order Qαβ ni (r) is 1 Q niαβ(r) = ⟨3uniαuniβ − δαβ⟩ (6) 2 where r = (ri + ri+1)/2 is the center-of-mass position of the considered bond, and α, β = x, y, z. When we average Qαβ ni over all the bonds in the system, we obtain some average characteristic Q̅ αβ of the nematic order in the system via its eigenvalues λ3 > λ2 > λ1. The magnitude of the largest eigenvalue (λ3) characterizes the alignment of the bonds, and the corresponding eigenvector gives the director of the nematic phase. For standard nematic order, λ1 = λ2 = −λ3/2, and therefore the single order parameter S ≡ λ3 suffices to completely characterize nematic order in the system. This behavior is indeed the case for semiflexible chains in the bulk, as shown in Figure 2a. For the chosen conditions (ρ = 0.7, N = 32), one finds λ3 ≥ 0.9 for κ ≥ 32 in the bulk system, and even for moderate stiffnesses κ = 8 the system is already deep in the nematic phase. (The Supporting Information contains an additional phase diagram for various chain lengths N, where the monomer density ρ was varied at fixed stiffness κ = 32.) The situation becomes more complex in the spherical capsule because global nematic ordering of the chains is suppressed

Figure 1. (a) Cross-sectional view on the equatorial plane in a system with N = 32, κ = 96, ρ = 0.7, and R = 35. The inset on the right-hand side shows a zoomed-in representation of the “crack” formed by the chain ends. The inset on the left shows a zoomed-in representation of the separation between the “bulk” (blue) and “surface-attached” (red) chains. (b) Snapshot picture of all the chain ends in the system, for the same case as (a). Three different colors are used, green for chains whose center of mass is in the surface shell so that both of their chain ends are near the surface, blue for all chain ends near the equatorial plane region, and red for the remaining ends of those chains which have their other end near the equatorial plane. Snapshots rendered using Visual Molecular Dynamics 1.9.2.60

through the confinement. Figure 2b shows the eigenvalues λ1, λ2, and λ3 as a function of chain stiffness κ, and it is clear that for κ < 8 the system is still in the isotropic phase, i.e., λ1 = λ2 = λ3 = 0. For moderate stiffness, 8 ≤ κ ≤ 20, we have some order, since λ3 is strictly nonzero. However, this order is quite distinct from standard nematic order in the bulk: first of all, the largest eigenvalue stays rather small, i.e., λ3 ≈ 0.1, rather than approaching a value of order unity with increasing κ in this regime; second, the symmetry relation valid in the uniaxial nematic phase λ1 = λ2 = −λ3/2 is dramatically violated. Rather, the eigenvalue spectrum of Q̅ αβ is characteristic of a biaxial nematic,61−65 where the order parameter B ≡ (λ2 − λ1)/2 is a measure of the biaxial character. Note that the tensor in eq 6 is averaged over all bonds in the system, irrespective of their position r in the system and position i along the chains. In the general case, a full characterization of order in bulk biaxial nematics requires a more detailed representation through the matrix QAB ab = ⟨(3laAlbB − δabδAB)/2⟩, where the subscripts a, b denote the molecular axes, A, B indicate the laboratory axes, and laA is the direction cosine between the molecular axis a and the laboratory axis A, etc.65 This formulation can then be reduced to four independent components of QAB ab (or four order parameters, respectively). However, since for our model in the bulk no biaxial order is possible, the above description in terms of the two order parameters S ≡ λ3 and B ≡ (λ2 − λ1)/2 suffices here. C


Article

Macromolecules

“biaxial phases”,18 this behavior clearly is a finite size effect, unrelated to the intrinsic biaxiality resulting from special molecules whose structure singles out two axes, e.g., V-shaped or banana-shaped molecules.63−65 In a bulk system of our semiflexible linear molecules, no other direction is singled out apart from the director orientation, and hence no biaxiality occurs. Although spherical confinement does not single out any direction, as stressed above, it causes a more or less uniform curvature of surface-attached stiff chains, reminiscent of the curvature of banana-shaped molecules. Of course, in our case this surface-induced curvature disappears as R → ∞, and thus the biaxiality observed here basically also results as a surface effect, though having a rather different character from the biaxiality found for planar walls. Figure 2c now shows the order parameter B plotted vs κ, with its value extracted from the tensor obtained from an average over all bonds, and also for the case where only bonds with r > 0.85R and r < 0.85R (i.e., in the surface and bulk region, respectively) are included. Interestingly, biaxiality in the surface shell is largest for 8 ≤ κ ≤ 20 and decreases when biaxility in the bulk region reaches its maximum. At the same time, λ3 in the surface region increases monotonically until κ = 24 and is also consistently larger than λ3 in the bulk region of the sphere (see Figure 2b). (Note that the tensor defined in eq 6 is additive with respect to all the bonds, but there is no additivity with respect to its eigenvalues: thus, the eigenvalues of the average tensor of all the bonds is not the sum of the eigenvalues of the surface region and of the interior.) This behavior is likely due to surface-induced ordering of the chains at high densities (see, e.g., the snapshot shown in Figure 3a). As κ is increased, however, λ3 in the surface region decreases monotonically and is well below the value in the central region. Visual inspection suggests that this transition is related to a change in character of the topological defects on the sphere surface from bipolar to quadrupolar (“tennis ball”), as shown in Figure 3a−c. To better understand this surprising behavior, we introduce an order parameter suitable for the distinction of topological defects in the surface layer.66 We define a tensor

Figure 2. Eigenvalues λ3 > λ2 > λ1 of the average nematic order tensor, eq 6, of the bonds plotted vs κ at fixed monomer density ρ = 0.7 and N = 32 in (a) the bulk and (b) in a sphere with radius R = 35. Solid lines were computed for all bonds in the system, while dashed and dotted curves in (b) show λ3 of bonds with r > 0.85R (i.e., from the surface region) and r < 0.85R, respectively. (c) Biaxiality parameter B = (λ2 − λ1)/2 plotted vs chain stiffness κ for the same case as in (b). The blue curve shows the data averaged over all bonds in the system, while the red and green curves show the bonds in the surface and bulk region, respectively.

αβ

Ω

1 = M−1

M−1

∑ i=1

⎛ u i × u i + 1 ⎞α ⎛ u i × u i + 1 ⎞ β ⎜ ⎟ ⎜ ⎟ ⎝ |u i × u i + 1| ⎠ ⎝ |u i × u i + 1| ⎠

(7)

for all unit vectors ui(r) lying in the shell r > 0.85R. Here, M denotes the total number of bonds in this shell, and ui × ui+1 characterizes the orientation of the plane formed by two successive bonds. The eigenvalues Ω1, Ω2, and Ω3 and their asymmetry αΩ are plotted in Figure 4 and give clear hints that two transitions occur where the orientational ordering changes: in the disordered phase (κ ≤ 4), the orientation of these planes is random, and hence the eigenvalues of the tensor Ωαβ are roughly equal to each other, namely Ω3 ≈ Ω2 ≈ Ω1 ≈ 1/3. Near κ = 8, the system transforms to the bipolar order in the surface layer (see the snapshot in Figure 3a), and we find Ω3 > Ω2 > Ω1, so that the asymmetry αΩ (defined in analogy to αΛ, eq 8) clearly is nonzero and positive. Further, the eigenvector npol belonging to the smallest eigenvalue (Ω1) provides a good description for the location of the poles. [Note that when we talk about “bipolar” structures, we have simplified matters somewhat, since actually in many cases there is a pair of topological defects close to each other in the vicinity of each pole (see Figure 3b).] Finally, we can identify a further transition to the quadrupolar “tennis-ball” configuration near κ

In confined systems, biaxial order in standard nematics arises due to the orientational ordering field caused by planar walls (see, e.g., ref 18 and references therein); however, while planar walls also single out a particular direction in space (normal to it), the walls of a perfect sphere do not, and in this respect the present situation is novel. Note that for all uniaxial molecules that exhibit standard nematic order (S > 0) in the bulk, a local biaxiality B(z) exists across a nematic−isotropic interface, with B(z) → 0 at large normal distances z from this interface. Similarly, a nonzero B(z) may occur also close to a wall15 confining a lyotropic solution, already under conditions where the bulk is still in the isotropic phase. As a consequence, in a slit pore confined by two planar walls a distance Lz apart, a small L average biaxiality B̅ = Lz−1∫ 0 zdz B(z) of order Lz−1 results. Although such thin films therefore are sometimes denoted as D


Article

Macromolecules

Figure 4 is similar to the biaxiality in the surface region, Bshell, plotted in Figure 2c. We emphasize these special orderings with topological defects, since clearly related structures occur for other parameter choices as well (e.g., the case κ = 144, N = 48, R = 35, and ρ = 0.7 shown in Figure 3d−f). We believe that these phenomena are general when conflicts appear between the long-range order that would arise in the bulk and geometries incompatible with such orderings enforced by confinement. To better understand the origin of the observed structures and their transitions, it is instructive to study the distribution and conformation of the confined polymers. The quantities that are most straightforward to study in our geometry are the radial density profiles of the chains and their constituent monomers. Figure 5 presents the radial density profile of all monomers, ρ(r), as well as of chain ends, ρE(r), and of mid-monomers, ρM(r), and the radial density profile of the chain’s center of mass, ρCM(r). As expected for a dense system with a hard repulsive surface, there is a pronounced layering of the total density near the hard confining sphere surface, irrespective of chain stiffness. However, the systematic decrease of total density near the sphere center for large κ is unexpected, and even more surprising is the behavior of end-monomers and mid-monomers: for stiff chains, mid-monomers are almost completely expelled from the center of the sphere, while endmonomers are enhanced in the center, and both types of monomers are enhanced near the sphere surface. Both of them exhibit a pronounced layering effect, and four layers can be distinctly recognized. The nontrivial distribution of the midmonomers in fact resembles the distribution of the chain CM (Figure 5d) since stiff chains behave as somewhat flexible rods, and then the chain CM must be close to its mid-monomer position. However, one important distinction is the fact that due to chains attached to the surface of the spherical wall, midmonomers can occur at a distance of ≈1.2σ from the surface, while the CM coordinate must be somewhat further inside (cf. Figure 5c,d) due to geometric considerations. The fact that the densities of both end-monomers, midmonomers, and the CM are enhanced in four layers adjacent to the confining sphere indicates that there a shell of stiff bent chains occurs. Assuming that wall-attached rather stiff chains have their monomers all at a distance ≈1.2σ from the wall (cf. Figure 5a) and that such chains take essentially the conformation of a bent rod (see inset of Figure 5d), one can predict the distance of the polymer’s CM from the sphere surface. Within this idealized picture, we find a value of R − r = 2.37σ for the CM position of a chain with N = 32 (indicated by the arrow in Figure 5d), which is in excellent agreement with the simulation data. Further, we note that the contour length of a chain in the last layer adjacent to the sphere surface (R′ ≈ 34) occupies only a fraction of the circumference, L/Lcircle = L/ (2πR′) ≈ 0.14. Thus, the angle between the last and the first bond, i.e., a1 and aN−1, along the chain is about 50.7°, but the angle Θijk between subsequent bond vectors only is about 1.63°. Thus, the bending energy per bond κΘijk2/2 due to this systematic chain bending of surface-attached chains still is rather small in comparison with the thermal energy; the chains hence significantly differ from rigid rods, even though the polymers are rather stiff. In this context, it is of interest to clarify to what extent the packing of the monomeric units is similar to the behavior of a crystalline solid phase, that is expected to form in the considered model system in the bulk for large enough densities,

Figure 3. (a) Snapshot of a system with N = 32, κ = 16, and ρ = 0.7 showing bipolar nematic order in a spherical cavity with radius R = 35. (b) Same system as in panel (a), but from the top perspective, highlighting the two defects which lie rather close to each other. (c) “Tennis-ball” order at N = 32, κ = 32, and ρ = 0.7. (d) A snapshot of a phase of densely packed semiflexible chains of length N = 48, for the case κ = 144, ρ = 0.7, and R = 35, showing only a view of chain end positions. (e) Side view of the chains for the case of panel (d) from the same perspective. Along the lines where chain ends accumulate clear discontinuities of orientational order can be recognized. (f) Shown is the system with R = 35, N = 48, κ = 144, and ρ = 0.7, viewed from outside showing also the nematic director. Because of the high density of the beads, only chains in the vicinity of the sphere surface are visible in all panels except (d).

Figure 4. Eigenvalues Ωi of the tensor Ω (left axis) and their asymmetry αΩ (right axis) plotted vs κ at N = 32 and fixed monomer density ρ = 0.7.

= 24 (Figure 3c), where Ω3 ≈ Ω2 > Ω1 and 0 < αΩ < 0.1. In this regime, the bonds in the bulk region (r < 0.85R) still exhibit nematic order (cf. Figure 2b) and the ordering of the surface region with its “tennis-ball” texture has decoupled somewhat. Note that the evolution of αΩ with growing stiffness κ shown in E


Article

Macromolecules

Figure 5. Radial distribution of (a) the total monomer density ρ(r), (b) the end-monomer density ρE(r), (c) the mid-monomer density ρM(r), and (d) the center-of-mass (CM) density ρCM(r) plotted vs the distance from the walls, R − r, for N = 32, R = 35, ρ = 0.7, and several choices of κ, as indicated in (a). A logarithmic scale for the abscissa is chosen in order to display the oscillatory structure near the confining sphere surface more clearly. The arrow in (d) indicates the estimated CM distance when the polymer is fully aligned to the curved sphere surface, as schematically shown in the same panel (mid- and end-monomers are colored in green and orange, respectively; CM position indicated through cross).

= 0) and stiff polymers (κ = 32) at a high density ρ = 0.7, and stiff polymers (κ = 32) at a low density (ρ = 0.1). Note that in the region 0 < r ≤ R/2 the monomer density is essentially constant and hence g(rij) is independent of the choice of origin for rij = rj − ri. We can identify a very sharp peak at rij = lb for all systems, which results from the bonded monomers along a chain. The nearest-neighbor distance along a chain can fluctuate only very little, irrespective of stiffness. For rather stiff polymers, locally the chains behave like rigid rods, and therefore g(rij) exhibits rather well-defined peaks also at twice and three times the nearest-neighbor distance. These intrachain peaks are the only structure visible at low densities but are present at high densities as well. For κ = 32 and ρ = 0.7, we also recognize a relatively pronounced peak at rij ≈ 1.3lb, which is related to the closest monomer of neighboring chains. However, more distant peaks due to monomers of neighboring chains are rather broad and smeared out. If the monomer positions of neighboring chains were strongly correlated with the monomer positions of the considered chains, a much more pronounced structure in g(rij) would occur. Thus, we can conclude that in spite of the strong orientational order that is present in the system with ρ = 0.7 and κ = 32, the system behaves like a disordered fluid with respect to positional order. The fact that most of the stiff chains have their midmonomers and CM in the range 13 ≤ r ≤ 26 with a clear peak near r ≈ 17 (see Figure 5c,d), while the end-monomers have a peak near the sphere origin, indicates that they form a “chain population” with very different behavior. This peculiar chain distribution reflects the formation of a bipolar structure in the sphere center, with an equatorial plane perpendicular to the axis connecting the poles, z. From the snapshot pictures shown in Figure 1, we can conclude that chains in the central bulk-like

amounting to a regular periodic packing of the chains (that then are rod-like stretched out). Such structures most easily are identified by studying the nearest-neighbor pair correlation function of monomeric units (not distinguishing whether the monomers belong to the same chain or different chains) and the Steinhardt bond orientational order parameters (based on the same definition what the nearest neighbors are).59 In Figure 6, we compare the radial pair distribution function g(rij) between monomers i and j located in the central region of the sphere, i.e., |ri | ≤ R/2 and |rj | ≤ R/2, for three cases: flexible (κ

Figure 6. Radial pair distribution function g(rij) between monomers i and j in the central region of the confining sphere (R = 35) with their positions |ri | ≤ R/2 and |rj| ≤ R/2. The schematic illustrates the structure when chains are completely stretched out and monomers of neighboring chains are perfectly correlated (which is not the case here, see text). F


Article

Macromolecules region of the sphere have their CM position roughly at a distance Δz = L/2 above or below the equatorial plane. Chains which have their ends close to the sphere center can orient straight along the z-axis toward either the “north” or “south pole” of the sphere, respectively; however, chains with chain ends near the boundary of the equatorial plane cannot point along the z-axis but must bend gradually (in general having a nonplanar conformation) so that their other chain end can be close to the sphere surface. The snapshot picture shown in Figure 1b clearly indicates that the resulting density distribution of chain ends near the sphere surface is strongly nonuniform: chain ends of both the chains with their CM in the sphere interior and of chains with their CM near the surface are enriched in a broad region near the poles and also near the equator. The intermediate regime of the near-surface region of the sphere is taken predominantly by the inner monomers of the surface-attached chains. As a result, we can conclude that the chains in the wall-attached region arrange by no means independently of the chains filling the bulk region of the sphere. Rather, the geometric arrangement of the chains both in the bulk and the surface region are strongly coupled. To identify more quantitatively this spontaneous breaking of the spherical symmetry, we computed the gyration tensor of all chains ends with r < 0.5 R and its eigenvalues (Λ3 > Λ2 > Λ1) and corresponding eigenvectors. The choice r < 0.5R is made to safely exclude the surface shell, but since ρE(r) is rather small and only weakly dependent on r near r = 0.5R, the precise value of this cutoff does not matter. Alternatively, the equatorial plane can also be determined as the plane which has the minimum distance with respect to all chain ends within r < 0.5R. If there were no surplus of end-monomers in the central region, then this procedure will still yield a plane normal, resulting from statistical fluctuations in the end-monomer distributions. To clarify whether the obtained plane is real (as shown in Figure 1) or an artifact of the algorithm, we estimate the asymmetry of the distribution of end-monomers in the entire system, r < R, defined as αΛ = Λ3 − (Λ 2 + Λ1)/2

(8)

noting that αΛ = 0 for a perfectly isotropic system. In Figure S2, we have plotted this asymmetry parameter as a function of κ for ρ = 0.7 as well as a function of ρ for fixed κ = 32. There, one can see that strong asymmetry (αΛ > 50) occurs for κ ≳ 20. We have also checked for κ = 32 how sensitive this parameter is to the specific choice of density. Indeed, for ρ ≳ 0.4 the asymmetry is large and depends on ρ only weakly.35 Once the equatorial plane has been identified, one can also study some of the quantities presented in Figure 5 as a function of the distance z from the plane (see Figure 7). As suspected from the snapshot pictures, chain ends are enriched in the region near the equator (z ≤ 2) and in a broad regime near the poles (25 ≤ z ≤ 32), while the CM for stiff chains is enriched near z ≈ L/2 ≈ 15. All these phenomena happen only for rather stiff chains, κ ≥ 32, which behave as somewhat flexible rods already. Interestingly, confined chains with κ = 8 and κ = 16 still have almost flat distributions of ρE(z) and ρCM(z), although strong nematic order was found in the bulk for the same parameters. This discrepancy is likely due to confinementinduced bending of the surface-attached chains, which overrules the emerging nematic ordering in the center of the sphere (see also Figure 2). Figure 7c shows the parallel and the perpendicular components of the squared end-to-end vector, Ree,∥2 and

Figure 7. Density distribution of (a) chain ends, ρE, and (b) of the center of mass (CM), ρCM, plotted vs their distance z from the equatorial plane for 8 choices of κ at ρ = 0.7, as indicated. Panel c shows the variation of the squared end-to-end length Ree2(z), Ree,∥2(z), and Ree,⊥2(z) vs the CM position z of the chains for the particular case ρ = 0.7 and κ = 96. Note that irrespective of their CM position, chains are almost fully stretched out (which would mean Ree2/L2 = 1).

Ree,⊥2, respectively. Near the sphere center, both components are close to the theoretical value expected for randomly oriented chains, i.e., Ree,∥2 ≈ L2/3 and Ree,⊥2 ≈ 2L2/3. For z = L/2 ≈ 15 (where chain ends are near the equatorial plane z = 0), the parallel component of Ree2 reaches a maximum close to the limit of rigid rods, i.e., Ree,∥2 ≈ 0.9L2. At distances z ≳16, the end-to-end vector of the polymers becomes more and more aligned perpendicular to the CM vector, so that Ree,∥2 → 0 and Ree,⊥2 maximized. At the end of this section, we briefly comment on the way the systems have been prepared, since equilibration of such dense and strongly confined systems as studied here is notoriously difficult. For the data shown in Figure 2, we started with κ = 96 and reduced κ gradually. However, using the inverse process starting at κ = 0 (i.e., fully flexible chains), some G


Article

Macromolecules hysteresis near κ = 24 was found, which coincides with the location of the isotropic−bipolar transition in this system. Such hysteresis associated with phase transitions is a rather common experience from simulations50,53 and also from experiments. In order to judge which state is the stable phase and which is only metastable, a knowledge of the free energy is required. By the combination of MD and DFT described in the Appendix, we have attempted to compute the free energies associated with the two paths. Figure 8 shows the total free energy per volume,

orientational parts and on the MD procedure for reaching equilibrium can be found in the Supporting Information.

4. SPHERICAL CONFINEMENT OF SEMIFLEXIBLE POLYMERS WHEN THEIR CONTOUR LENGTH SIGNIFICANTLY DIFFERS FROM THE SPHERE RADIUS 4.1. Chains That Are Too Long to Allow Formation of an Equatorial Plane. In Figure 3d−f we have already drawn attention to the special ordering of the chain ends at the surface as well as the surface-attached chains as a whole when N = 48, and thus L distinctly exceeds R. It is clear that then the structure with an equatorial plane formed by end-monomers, depicted in Figure 1, cannot fit into the sphere, and so also the chain arrangement must be strongly affected. Thus, we study again the radial variation of end-monomers, mid-monomers, and mean-squared end-to-end vector components parallel and perpendicular to the radius vector for this case (see Figure 9). It is remarkable that for flexible chains the end-monomers are distributed almost homogeneously inside the sphere (note that this homogeneous density would be 2ρ/N ≈ 0.0292) apart from two rather weak layering peaks adjacent to the sphere surface. When κ increases, the end-monomer density near the surface strongly increases, and up to six layering peaks can be recognized. Near the center of the sphere, the end-monomer density strongly decreases, completely different from the behavior for N = 32. Mid-monomers, on the other hand, also show pronounced layering at the surface, with a range that increases with the stiffness κ, but at the same time their density is decreased in the region where the layering occurs. Midmonomers for stiff chains have always a density maximum near the sphere center. Only for the stiffest chains (κ = 64, 96, and 144) a second maximum appears, corresponding to a distance of R − r ≈ L/2 ≈ 24 from the sphere surface. A rather nontrivial behavior is also seen for the length of the end-to-end vector, Ree2, shown in Figure 9c. Near the sphere center, both the parallel and the perpendicular components are close to the theoretical value for isotropic chain orientation, i.e., Ree,∥2 ≈ L2/3 and Ree,⊥2 ≈ 2L2/3, respectively, for all investigated κ. However, chains with their CM beyond ≈R/2 are almost fully aligned perpendicular with respect to the radial vector due to confinement, so that Ree,⊥2 ≈ L2. Note that the degree of alignment in this region also increases with increasing chain stiffness. The dependence of Ree2 on the chain stiffness can be understood when the nematic order in the system is considered (see Figure 10a). It is seen that for κ = 12 the system is still essentially disordered, unlike a corresponding bulk system.42,43 The nematic order develops gradually with increasing κ until around κ ≈ 48, where lp ≈ L, and then stays more or less independent of κ. This variation with κ is very different from the case of N = 32 (with the same sphere radius R = 35), where there was evidence for two successive transitions (at κ ≈ 8 and κ ≈ 24), and the value of the nematic order parameter S ≡ λ3 in the region of large κ was distinctly smaller (about 0.4, while here it is close to 0.6). This behavior is somewhat counterintuitive: one would think that the longer the chains become, the larger the frustration of the nematic order is due to the confinement, and hence the smaller S ≡ λ3 should be. In contrast, λ2 (which is exactly equal to λ1 in the bulk nematic phases) is here around λ2 ≈ −0.2, as previously (cf. Figures 10a

Figure 8. Free energy per volume F/V plotted vs stiffness κ at fixed N = 32, R = 35, and ρ = 0.7. Inset: eigenvalues λ3 > λ2 > λ1 of the average nematic order parameter of the bonds plotted vs κ (same conditions as in main panel). In both panels, open symbols connected by dashed lines show the data initialized with stiff chains (κ = 96) and gradually decreasing κ, whereas filled symbols connected by solid lines correspond to simulations initialized with fully flexible chains (κ = 0) and gradually increasing κ.

F/V, as a function of chain stiffness, κ, for fixed values of the chain length, monomer density, and sphere radius. For both paths, F/V decreases nearly monotonically with increasing κ. In order to understand this behavior, one needs to consider the decomposition of F into its orientational and translational contributions (see Figure S4). The translational term shows only weak dependence on the chain stiffness. In contrast, the orientational term exhibits a much stronger variation, with its ideal part increasing and the excess part decreasing with κ. The former effect is due to the increasing difficulty of packing stiffer chains in a sphere, while the latter one is due to the decrease of the chain’s effective thickness with increasing κ.56 Overall, the excess term dominates over the ideal one, which results in the decrease of the orientational (and, hence, the total) free energy with increasing chain stiffness. For κ much smaller than the value of the expected isotropic− bipolar transition both paths yield the same free energy. In contrast, for κ = 24 the path coming from the ordered side seems to overestimate the free energy, while for κ ≥ 28 the free energy along the path from the ordered side seems to be distinctly lower. These findings are qualitatively as expected and put the transition in the vicinity of κ ≈ 26, i.e., inside the region where the hysteresis of the eigenvalues occurs (see inset of Figure 8). We conclude from these findings that the results which we report here are not strongly affected by such nonequilibrium effects, apart from the regions where the character of the ordering changes. More details about the decomposition of the free energy into translational and H


Article

Macromolecules

Figure 9. (a) End-monomer density and (b) mid-monomer density plotted vs the radial distance from the sphere surface for the case N = 48, ρ = 0.7, and R = 35. Eight choices of κ from κ = 1 to κ = 144 are included, as indicated. (c) Ree,∥2(r) (solid lines) and Ree,⊥2(r) (dashed lines) vs the radial distance from the origin for the same choice of N, ρ, and R, and the three values of κ = 12, 48, and 96, as indicated.

Figure 10. (a) Eigenvalues λ3 > λ2 > λ1 plotted vs κ for the case N = 48, ρ = 0.7, and R = 35. (b) Density ρE(z) of end-monomers plotted against the distance from the sphere center (z) along the director for the case R = 35, N = 48, and ρ = 0.7 for eight choices of κ, as indicated. (c) Same as (b), but for the mid-monomers. Note that the high densities of both end-monomers and mid-monomers for 32 ≤ z ≤ 34 show that some chains are attached to the surface of the sphere orthogonal to the director.

and 2b). Thus, a slight biaxiality B develops for κ ≈ 12 and changes relatively little when κ is increased further. Having found the eigenvalues of the nematic order, we also know the orientation of the director and can ask how are chain ends, mid-monomers, etc., distributed along this direction (see Figure 10b,c). While for the case where L and R are comparable (Figure 7) chain ends in the ordered phase are accumulated near the plane z = 0; now a depression of chain ends near z = 0 occurs, and rather a peak at large z develops (first near z = 25), which gradually moves toward z = 32 for κ = 144. The midmonomers, however, first are enhanced near z = 0 when nematic ordering has started, but then this peak near z = 0 diminishes in height, and rather another peak near z = 10 develops for the stiffest chains. At larger values of z, a minimum of ρM(z) develops, which for intermediate κ is rather shallow. For very stiff chains, κ ≥ 96, we observe an almost zero mid-

monomer density in 25 ≤ z ≤ 29, while then another maximum is reached in the layers adjacent to the sphere surface. Thus, the optimization of the chain conformations of densely packed rather stiff polymers in spherical confinement clearly is a delicate problem, and a comparison of the cases N = 32 and N = 48 in spheres with R = 35 reveals that there are subtle commensurability issues. 4.2. Short Chains: Formation of Several Planes Orthogonal to the Director Where Chain Ends Get Enriched. In this section, we first focus on the case N = 16, R = 35, and ρ = 0.7, since due to the rather short chain length and not too large total number of monomers (5 mon = 125 696) equilibration of this system is fairly easy. Figure 11 presents typical snapshot pictures, which clearly reveal that now three parallel planes have formed where chain ends get enriched. The I


Article

Macromolecules

poles of the sphere (we consider here the situation where the director is oriented along the axis through the poles). Figure 12a then shows the eigenvalues λ3 > λ2 > λ1 vs κ, indicating two transitions at κ ≈ 8 and κ ≈ 40. It is also remarkable that λ1 ≈ λ2 ≈ −λ3/2, indicating uniaxial ordering. Thus, the biaxiality parameter B is rather small for these conditions, in contrast to our previous findings for the systems with N = 32 and R = 35 (cf. Figure 2c). The corresponding distributions of end-monomers and midmonomers are plotted in Figure 12b,c. The first two maxima of ρM(z) fall just halfway in between the corresponding maxima of ρE(z), as would be true for a perfect smectic arrangement. However, this “rule” breaks down for the outmost maxima, as expected, since the sphere surface massively disturbs the quasismectic arrangement. It is also seen that this quasi-smectic distribution already develops gradually (but weakly) in the regime 16 ≤ κ ≤ 32, where λ3 still increases and the nematic order improves (see Figure 12a); however, for κ ≥ 48 the modulation is very strong and no longer looks like a superimposed sine function, and this strong modulation goes along with a reduction of the average order parameter (cf. Figure 12a). Figure 12d shows the second Legendre polynomial, P2(z) = (3 cos2 Θ − 1)/2 of the angle Θ between a polymer bond and the director as a function of the distance z from the sphere center. Here, a value of P2 = 1 indicates parallel alignment, whereas a value of P2 = −1/2 corresponds to a perpendicular orientation. For P2 = 0, the bonds are oriented isotropically with respect to the director. It is clear from the data plotted in Figure 12d that the reduction of the average order parameter happens predominantly in the outer region of the sphere, where increasing chain stiffness enforces that more beads follow the

Figure 11. Simulation snapshots for a system with N = 16, R = 35, ρ = 0.7, and κ = 96. (a) General side view and (b) distribution of chain ends. (c) General side view from the same viewing position as in (a), but including only the monomers inside of r ≤ 30. (d) View through the sphere center along the director. All snapshots are drawn at the same scale.

snapshot in Figure 11c looks like a section cut out from a bulk smectic C phase. Indeed, for high enough density in the bulk, semiflexible polymers can experience a nematic−smectic transition, as will be discussed in a future publication. However, the snapshot in Figure 11d reveals a whirl in the chain orientation, projected on the plane orthogonal to the director. This behavior demonstrates that spherical confinement also involves topological defects in the order associated with the

Figure 12. (a) Eigenvalues λ3 > λ2 > λ1 plotted vs κ for the case N = 16, R = 35, and ρ = 0.7. (b) Density of end-monomers, ρE(z), plotted against the distance z from the sphere center along the director for the same case as in (a) for eight choices of κ, as indicated. (c) Same as (b), but for midmonomers. Note that the maxima occur for z ≈ 6.5, 22, and 32, while the end-monomer distribution maxima occur for z ≈ 0, 13, and 30. (d) Local bond vector orientation P2(z) vs z, for eight choices of κ, as indicated. J


Article

Macromolecules

chain linear dimensions (not shown here for the sake of brevity). Finally, we turn to the case N = 32, R = 70, and ρ = 0.7, where hence essentially the mesoscopic lengths (L, R) have been scaled up by a factor of 2 in comparison with the previous case, while otherwise the situation should resemble it. From a simulation perspective, this is the most challenging case, involving in total 1 005 216 monomeric units, and hence we can make only much less definitive statements about the quality of equilibration. Figure 14 shows a few snapshot pictures. Naively,

orientation perpendicular to the director and hence parallel to the sphere surface. This conclusion is substantiated by examining radial distribution functions of end-monomers, mid-monomers, and bond orientations (see Figure 13). The

Figure 14. (a, b) Snapshots of chains with N = 32 and κ = 48 in the sphere with radius R = 70 and density ρ = 0.7. Panel (a) shows two topological defects on the sphere surface, while panel (b) shows a view on the inner part of the sphere for r ≤ 60, from the same perspective. (c, d) Snapshots of chains with N = 32 and κ = 128 in the sphere with radius R = 70 and density ρ = 0.7. Panel (c) shows the distribution of chain ends only, while panel (d) shows a view on the inner part of the sphere for r ≤ 60, from the same perspective. All snapshots are drawn at the same scale.

one would expect that the situation is qualitatively the same as for the case N = 16 and R = 35, but this is not always the case: comparison of Figures 11 and 14c,d shows that indeed structures with planes where chain ends are enriched do occur again, but in the present case these structures do not fill the whole sphere but are rather restricted to a single hemisphere only. Further, the axis normal to these planes in each hemisphere is rotated with respect to each other. Similarly (not shown here) the two hemispheres exhibit also different orientation of the director. It is not obvious whether the structures in Figure 14 represent the true equilibrium behavior or should be considered as a frozen-in domain pattern. Finally, Figure 15 shows one more time radial profiles of end-monomers and midmonomers. One recognizes again that for large κ the density of end-monomers is enhanced and the density of midmonomers is depressed near the sphere center, over a range up to about r = 10. Both types of monomers are enhanced near the sphere surface and exhibit strong layering there. Surprisingly, the range over which these distributions are enhanced or depressed near the confining center has about the same size as for N = 16 in the sphere with R = 35 (cf. Figure 13). Thus, there clearly is not a simple scaling property which would suggest that systems with the same ratio R/L behave similarly for stiff enough chains.

Figure 13. (a) Radial distributions of end-monomer density ρE(r), (b) mid-monomer density ρM(r), and (c) local bond orientation P2(r) plotted vs the distance from the walls, R − r, for the case N = 16, R = 35, ρ = 0.7, and eight values of κ, as indicated. A logarithmic scale for the abscissa is chosen in order to display the oscillatory structure near the confining sphere surface more clearly.

enhancement of end-monomer and mid-monomer densities near the surface of the sphere increases with increasing stiffness. Also, the layering of P2(r) (defined in an analogous way as P2(z), but with respect to the radial distance r) becomes more pronounced. At the same time, P2(r) becomes more and more negative, over a wider region of r adjacent to the sphere surface, indicating that for large enough κ the need of the bonds to orient parallel to the sphere surface (perfect parallel orientation would mean P2(r) = −0.5) wins in comparison with the tendency of the system to develop uniform long-range order. These conclusions are compatible with the analysis of local K


Article

Macromolecules

be more or less linearly oriented perpendicular to the equatorial plane; the larger the distance of the chain end from the sphere center, the stronger the chain must bend. In order to avoid crowding of chain ends near the poles, the nematic order develops special topological defects (see the whirl-type structure in Figure 3). The ordering of the surface-attached chains then is strongly coupled with the ordering in the sphere interior; while for intermediate stiffness it is also bipolar, for very large stiffness the quadrupolar “tennis-ball” texture appears. This bipolar−“tennis ball” transition goes along with a strong reduction in the value of the average nematic order parameter. For such systems, where nematic order is strongly frustrated by the spherical confinement, we observe a very nontrivial behavior of all the eigenvalues λ3, λ2, and λ1 of the tensor describing orientational order. Unlike the case where L exceeds R, mid-monomers are repelled from the region near the sphere center, while the density of end-monomers is enhanced, simply because of the many ends in the equatorial plane. Finally, when R is much larger than L, such as the case for N = 16 (L ≈ 15), R = 35, three parallel planes where chain ends are enriched form, so the structure in the sphere interior looks like a smectic liquid crystalline order. However, the more one proceeds toward the surface of the confining sphere, the more the ordering gets distorted, and topological defects characterize the chain arrangement in the shell containing the wall-attached polymers as well. On the methodological front, we have demonstrated in the present work the usefulness of combining MD and density functional theory (DFT) to compute free energy differences between different states. In DFT, the free energy is expressed as a function of the local density which depends on bond orientation as well in this case. Because of the nontrivial inhomogeneity of the system in these situations (and the spontaneous symmetry breaking of orientational order with confinement-induced topological defects), calculation of this orientation-dependent local density by minimizing the free energy functional would be a premature task. However, since this local density is computed by MD more readily, one can use these MD results as input into the DFT free energy expression. In this way, we succeeded in dealing with the hysteresis observed in MD for the isotropic−bipolar transition.

Figure 15. (a) End-monomer density ρE(r) and (b) mid-monomer density ρM(r) plotted vs the distance from the walls, R − r, for the case N = 32, R = 70, ρ = 0.7, and eight values of κ, as indicated. A logarithmic scale for the abscissa is chosen in order to display the oscillatory structure near the confining sphere surface more clearly.

5. DISCUSSION In this paper, we studied spherical containers rather densely filled with semiflexible polymers using molecular dynamics (MD) simulations of a coarse-grained bead−spring model. We have restricted attention to the case where both monomer− monomer and monomer−wall interactions are purely repulsive and focus on situations where contour length L, persistence length lp, and sphere radius R are all comparable. We find that close to the confining walls the chains are bent due to their residual flexibility and that their monomers form a layering structure of a few shells of beads. Further inside the spherical container, the chains have (at least) one chain end far away from the surface of the confining sphere, and their precise arrangement and hence also the degree of nematic order inside the sphere depend very much on commensurability issues between L and R. When L exceeds R, the mid-monomers are enhanced near the sphere center for stiff enough chains, and almost ten shells in their layered arrangement near the sphere surface can be recognized. End-monomers, however, are effectively attracted to the sphere surface, their density near the center is suppressed, and a rather large value of the nematic order parameter (S almost 0.6) is found for sufficiently stiff chains (lp ≳ L). The surface region of the system exhibits an interesting arrangement of topological defects in the orientational order, the so-called “tennis-ball” texture. When L is comparable to R, a bipolar structure can form in the sphere interior, since now two stretched-out chains can fit along the axis connecting the poles. As a consequence, an equatorial plane can form, where chain ends get enriched. Of course, only chains with one end close to the sphere center can

■

APPENDIX. FREE ENERGY CALCULATION: HYBRID MD-DFT APPROACH In the anisotropic regime, the nonuniform molecular density of semiflexible polymers under spherical confinement ρp(r,ω) depends on both the location r inside the sphere and on the two polar angles (θ,ϕ) ≡ ω specifying the angular orientation of the chain. Hence, the minimization of the Helmholtz free energy F with respect to ρp(r,ω), which must be performed in order to obtain the equilibrium distribution, presents a formidable numerical challenge. Instead, we adopt a hybrid MD-DFT method for calculating free energy, which was previously used to study the interaction between two spherical polymer brushes.67 Namely, we take the required nonuniform density distributions ρp(r,ω) directly from MD simulations and substitute them into the DFT expression for the free energy, without performing any further optimization of the density field. It should be noted that this approach leads to an approximation for F: while the MD distributions are exact for the given model (within inevitable statistical noise), the DFT free energy functional is necessarily approximate, and its minimization with respect to ρp(r,ω) (if it could be carried L


Article

Macromolecules

sampling of the local density and of the orientational distributions. We found that Δa = 3 σ gives a good compromise between these factors for the sphere with R = 35. For the free energy calculation, the distributions in each cell were averaged over all individual snapshots. Error bars were estimated by subdividing the simulation data into three subsets and computing the standard deviation between these smaller sample sizes. At the end of this section, we would like to rationalize our choice for the excess free energy terms, for which exact expressions are not available, and approximations are thus unavoidable. The excess orientational term is the key contribution to the free energy in the present context because it largely governs the isotropic−nematic transition. To calculate this term, we have used a semiphenomenological approach due to Fynewever and Yethiraj,56 whereby the excluded volume for a pair of semiflexible chains is obtained as a function of their mutual orientation from Monte Carlo simulations and fit to a simple three-parameter analytical function. This approach is not devoid of some limitations. In particular, it would predict an isotropic−nematic transition even for fully flexible chains with κ = 0. Furthermore, for semiflexible chains with κ = 32 (which are an important part of the present study) the approach predicts that the isotropic−nematic transition density initially decreases with the contour chain length (as expected and seen in MD), but for N > 64 it starts increasing with N; the latter (spurious) prediction of the theory has not been confirmed by MD.43 Finally, the Fynewever−Yethiraj expression for the excess orientational term involves a phenomenological prefactor which approximately accounts for hard-sphere-like excluded interactions at a given packing fraction of monomers. This prefactor is based on the standard CS expression for the free energy of (unconnected) hard spheres.69 Accordingly, in the translational part of our excess free energy expression, we also use the simple CS form. The primary rational for this choice has to do with avoiding the double-counting of excluded volume effects between rotational and translational terms. More specifically, we modify the rotational part by subtracting an offset value which is obtained by computing the Fynewever−Yethiraj expression assuming uniform orientational distribution functions for both molecules. With the CS prefactor, this offset exactly cancels the excess translational term in the isotropic phase, so that the double-counting is avoided. The disadvantage of this approximation is that it completely neglects the connectivity of the chains. However, one can point out that several (approximate) excess translational free energy functionals for polymeric molecules available in the literature67 are all limited to f ully f lexible chains and therefore do not take into account the effect of chain stiffness on the excess translational free energy. Given the importance of the stiffness parameter for the present study and the fact that the Fynewever−Yethiraj expression for the excess translational term effectively accounts for both stiffness and connectivity (via two-chain simulations), we feel that neglecting connectivity in the translational part is a reasonable approximation.

out) would have yielded distributions somewhat different from the ones obtained in the MD simulation. However, our main interest here is not in the absolute values of the free energy but rather in the trends they display, e.g., as a function of the chain stiffness κ (at a given monomer density ρ) or as a function of ρ (at a given κ). The DFT-based calculation of F(r) proceeds as follows. First, following the standard practice, F(r) is split into ideal and excess terms, with the former known exactly:67 Fid(r) = kBT

∫ dωρp(r, ω)(ln[4πρp(r, ω)] − 1)

(9)

The excess term is decomposed into translational and tr or orientational components Fexc and Fexc , respectively. The 34 translational contribution is given by tr Fexc (r) = ρ(r)fexc (η ̅ (r)) kBT

(10)

where η(r) = ρ(r)π/6 is the monomer packing fraction, η(r) ̅ is its corresponding weighted packing fraction,68 and fexc(η(r)) is ̅ the excess translational free energy density from the Carnahan− Starling (CS) equation of state:69 fexc (η ̅ (r)) =

η ̅ (r)(4 − 3η ̅ (r)) (1 − η ̅ (r))2

(11)

Finally, the excess orientational term has the following form:34 or Fexc (r) = kBT

∫ dω ∫ dω′ 8(14 −− 3ηη(̅ (rr))) 2

̅ iso ′ × (Vexc(ω , ω ) − ⟨Vexc⟩)ρp (r, ω)ρp (r, ω′)

(12)

where Vexc(ω,ω′) is the excluded volume for two semiflexible polymers with angular orientations ω and ω′, respectively, for which we take a phenomenological form obtained via fitting the simulation data.56 (Note that we have subtracted its spherical average ⟨Viso exc⟩ in order to avoid the double-counting of the isotropic contribution to the excess free energy, which is already taken into account via eq 10.) In order to perform the free energy calculation of the MD simulation data, we first split the entire simulation volume into small cubic integration cells, each with an edge length of Δa and centered around location r inside the sphere. Combining the contributions given from eqs 9, 10, and 12, one obtains the total free energy for a given cell, and summing over all cells yields the total free energy in the cavity. In order to generate the distribution ρp(r,ω) for each cell, we first read in the position data and then rotated all snapshots so that the nematic director in each frame coincides with the z-axis. After this transformation, we computed the local polymer (ρp) and monomer density (ρ). At this stage, we also determined the orientational distribution in each cell by recording the polar (θ) and azimuthal (ϕ) angles of the corresponding bond vectors. Here, we used the approximation f(r,ω) ≈ f(r,θ) f(r,ϕ) to improve the sampling of our data. The range of the angles was chosen 0 ≤ ϕ < 2π and 0 ≤ θ < π, and we made the distributions f(r,θ) and f(r,ϕ) symmetric around their midpoints because the bonds have neither “heads” nor “tails”. For the numerical integration over the spherical volume, the size of the cubic cells has to be chosen small enough to resolve fluctuations and so that the edges of the sphere are sufficiently smooth. However, a too small value of Δa will lead to poor

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.7b02643. M


Article

Macromolecules

■

semiflexible chains at nematic ordering transitions in thin films: a monte carlo simulation. Macromolecules 2014, 47, 1206−1220. (18) Ye, S.; Zhang, P.; Chen, J. Z. Y. Surface-induced phase transitions of wormlike chains in slit confinement. Soft Matter 2016, 12, 2948−2959. (19) Luzhbin, D. A.; Chen, Y.-L. Shifting the isotropic-nematic transition in very strongly confined semiflexible polymer solutions. Macromolecules 2016, 49, 6139−6147. (20) Egorov, S. A.; Milchev, A.; Binder, K. Semiflexible polymers in the bulk and confined by planar walls. Polymers 2016, 8, 296. (21) Egorov, S. A.; Milchev, A.; Virnau, P.; Binder, K. Semiflexible polymers under good solvent conditions interacting with repulsive walls. J. Chem. Phys. 2016, 144, 174902. (22) Egorov, S. A.; Milchev, A.; Binder, K. Capillary nematization of semiflexible polymers. Macromol. Theory Simul. 2017, 26, 1600036. (23) Milchev, A.; Egorov, S. A.; Binder, K. Semiflexible polymers confined in a slit pore with attractive walls: two-dimensional liquid crystalline order versus capillary nematization. Soft Matter 2017, 13, 1888−1903. (24) Milchev, A.; Binder, K. Smectic C and nematic phases in strongly adsorbed layers of semiflexible polymers. Nano Lett. 2017, 17, 4924. (25) Sakaue, T. Semiflexible polymer confined in closed spaces. Macromolecules 2007, 40, 5206−5211. (26) Morrison, G.; Thirumalai, D. Semiflexible chains in confined spaces. Phys. Rev. E 2009, 79, 011924. (27) Ostermeir, K.; Alim, K.; Frey, E. Buckling of stiff polymer rings in weak spherical confinement. Phys. Rev. E 2010, 81, 061802. (28) Micheletti, C.; Marenduzzo, D.; Orlandini, E. Polymers with spatial or topological constraints: Theoretical and computational results. Phys. Rep. 2011, 504, 1−73. (29) Cifra, P.; Bleha, T. Free energy of polymers confined in open and closed cavities. Macromol. Theory Simul. 2012, 21, 15−23. (30) Reith, D.; Cifra, P.; Stasiak, A.; Virnau, P. Effective stiffening of DNA due to nematic ordering causes DNA molecules packed in phage capsids to preferentially form torus knots. Nucleic Acids Res. 2012, 40, 5129−5137. (31) Marenduzzo, D.; Micheletti, C.; Orlandini, E.; Sumners, D. Topological friction strongly affects viral DNA ejection. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 20081. (32) Fošnarič, M.; Iglič, A.; Kroll, D. M.; May, S. Monte Carlo simulations of a polymer confined within a fluid vesicle. Soft Matter 2013, 9, 3976−3984. (33) Gao, J.; Tang, P.; Yang, Y.; Chen, J. Z. Y. Free energy of a long semiflexible polymer confined in a spherical cavity. Soft Matter 2014, 10, 4674. (34) Milchev, A.; Egorov, S. A.; Nikoubashman, A.; Binder, K. Conformations and orientational ordering of semiflexible polymers in spherical confinement. J. Chem. Phys. 2017, 146, 194907. (35) Nikoubashman, A.; Vega, D. A.; Binder, K.; Milchev, A. Semiflexible polymers in spherical confinement: Bipolar orientational order versus tennis ball states. Phys. Rev. Lett. 2017, 118, 217803. (36) Soares e Silva, M.; Alvarado, J.; Nguyen, J.; Georgoulia, N.; Mulder, B.; Koenderink, G. Self-organized patterns of actin filaments in cell-sized confinement. Soft Matter 2011, 7, 10631. (37) Qi, W.; Yan, X.; Fei, J.; Wang, A.; Cui, Y.; Li, J. Triggered release of insulin from glucose-sensitive enzyme multilayer shells. Biomaterials 2009, 30, 2799. (38) Wang, Y.; Wang, L.; Li, B.; Cheng, Y.; Zhou, D.; Chen, X.; Jing, X.; Huang, Y. Compact vesicles self-assembled from binary graft copolymers with high hydrophilic fraction for potential drug/protein delivery. ACS Macro Lett. 2017, 6, 1186. (39) Kato, A.; Shindo, E.; Sakaue, T.; Tsuji, A.; Yoshikawa, K. Conformational transition of giant dna in a confined space surrounded by a phospholipid membrane. Biophys. J. 2009, 97, 1678. (40) Negishi, M.; Sakaue, T.; Takiguchi, K.; Yoshikawa, K. Cooperation between giant DNA molecules and actin filaments in a microsphere. Phys. Rev. E 2010, 81, 051921.

Additional information on the simulation procedure and data analysis (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (A.N.). ORCID

Andrey Milchev: 0000-0003-1991-0648 Arash Nikoubashman: 0000-0003-0563-825X Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful to the German Research Foundation (DFG) for partial support under grants NI 1487/2-1 (A.N.) and BI 314/24-1 (A.M.). Further, D.A.V. acknowledges support from The National Scientific and Technical Research Council (CONICET) and the Universidad Nacional del Sur (UNS). Computing time was granted on the supercomputer Mogon at Johannes Gutenberg University Mainz.

■

REFERENCES

(1) Ciferri, A., Ed.; Liquid Crystallinity in Polymers: Principles and Fundamental Properties; VCH Publishers: New York, 1983. (2) Donald, A. M.; Windle, A. H.; Hanna, S. Liquid Crystalline Polymers; Cambridge University Press: Cambridge, 2006. (3) Mofrad, M. R. K., Kamm, R. D., Eds.; Cytoskeletal Mechanics; Cambridge University Press: Cambridge, 2006. (4) Alberts, B.; Johnson, A.; Lewis, J.; Raff, M.; Roberts, K.; Walter, P. Molecular Biology of the Cell; Garland Science: New York, 2007. (5) Kratky, O.; Porod, G. Röntgenuntersuchung gelöster Fadenmoleküle. Recl. Trav. Chim. Pays-Bas 1949, 68, 1106−1122. (6) Saitô, N.; Takahashi, K.; Yunoki, Y. The statistical mechanical theory of stiff chains. J. Phys. Soc. Jpn. 1967, 22, 219−226. (7) Grosberg, A. Y.; Khokhlov, A. R. Statistical Physics of Macromolecules; American Institute of Physics: Woodbury, NY, 1994. (8) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: Oxford, 2003. (9) Brelsford, G.; Krigbaum, W. Liquid Crystallinity in Polymers: Principles and Fundamental Properties; VCH Publishers: New York, 1983; Chapter Experimental evaluation of the persistence length for mesogenic polymers, pp 61−94. (10) Hsu, H.-P.; Paul, W.; Binder, K. Standard definitions of persistence length do not describe the local “intrinsic” stiffness of real polymer chains. Macromolecules 2010, 43, 3094−3102. (11) Reisner, W.; Pedersen, J. N.; Austin, R. H. DNA confinement in nanochannels: physics and biological applications. Rep. Prog. Phys. 2012, 75, 106601. (12) Chen, J. Z. Y. Theory of wormlike polymer chains in confinement. Prog. Polym. Sci. 2016, 54-55, 3−46. (13) Balducci, A.; Mao, P.; Han, J.; Doyle, P. S. Double-stranded DNA diffusion in slitlike nanochannels. Macromolecules 2006, 39, 6273−6281. (14) Chen, J. Z. Y.; Sullivan, D. E. Free energy of a wormlike polymer chain confined in a slit: Crossover between two scaling regimes. Macromolecules 2006, 39, 7769−7773. (15) Ivanov, V. A.; Rodionova, A. S.; Martemyanova, J. A.; Stukan, M. R.; Müller, M.; Paul, W.; Binder, K. Orientational ordering transitions of semiflexible polymers in thin films: A Monte Carlo simulation. Phys. Rev. E 2011, 84, 041810. (16) Ivanov, V. A.; Rodionova, A. S.; Martemyanova, J. A.; Stukan, M. R.; Müller, M.; Paul, W.; Binder, K. Wall-induced orientational order in athermal semidilute solutions of semiflexible polymers: Monte Carlo simulations of a lattice model. J. Chem. Phys. 2013, 138, 234903. (17) Ivanov, V. A.; Rodionova, A. S.; Martemyanova, J. A.; Stukan, M. R.; Müller, M.; Paul, W.; Binder, K. Conformational properties of N


Article

Macromolecules (41) Nakaya, K.; Imai, M.; Komura, S.; Kawakatsu, T.; Urakami, N. Polymer-confinement-induced nematic transition of microemulsion droplets. Europhys. Lett. 2005, 71, 494. (42) Egorov, S. A.; Milchev, A.; Binder, K. Anomalous fluctuations of nematic order in solutions of semiflexible polymers. Phys. Rev. Lett. 2016, 116, 187801. (43) Egorov, S. A.; Milchev, A.; Virnau, P.; Binder, K. A new insight into the isotropic - Nematic phase transition in lyotropic solutions of semiflexible polymers: Density-functional theory tested by molecular dynamics. Soft Matter 2016, 12, 4944. (44) de Gennes, P. G.; Prost, J. The Physics of Liquid Crystals, 2nd ed.; Clarendon: Oxford, 1992. (45) Trukhina, Y.; Schilling, T. Computer simulation study of a liquid crystal confined to a spherical cavity. Phys. Rev. E 2008, 77, 011701. (46) Prinsen, P.; van der Schoot, P. Shape and director-field transformation of tactoids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 021701. (47) Berggren, E.; Zannoni, C.; Chiccoli, C.; Pasini, P.; Semeria, F. Computer simulations of nematic droplets with bipolar boundary conditions. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1994, 50, 2929. (48) Oseen, C. W. The theory of liquid crystals. Trans. Faraday Soc. 1933, 29, 883−899. (49) Frank, F. C. I. Liquid crystals. On the theory of liquid crystals. Discuss. Faraday Soc. 1958, 25, 19−28. (50) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, CA, 2002. (51) Rapaport, D. C. The Art of Molecular Dynamics Simulation, 2nd ed.; Cambridge University Press: Cambridge, 2004. (52) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 2nd ed.; Oxford University Press: Oxford, 2017. (53) Landau, D. P.; Binder, K. A Guide to Monte Carlo Simulation in Statistical Physics, 4th ed.; Cambridge University Press: 2015. (54) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237−5247. (55) Bishop, M.; Kalos, M. H.; Frisch, H. L. Molecular dynamics of polymeric systems. J. Chem. Phys. 1979, 70, 1299−1304. (56) Fynewever, H.; Yethiraj, A. Phase behavior of semiflexible tangent hard sphere chains. J. Chem. Phys. 1998, 108, 1636. (57) Anderson, J. A.; Lorenz, C. D.; Travesset, A. General purpose molecular dynamics simulations fully implemented on graphics processing units. J. Comput. Phys. 2008, 227, 5342−5359. (58) Glaser, J.; Nguyen, T. D.; Anderson, J. A.; Lui, P.; Spiga, F.; Millan, J. A.; Morse, D. C.; Glotzer, S. C. Strong scaling of generalpurpose molecular dynamics simulations on GPUs. Comput. Phys. Commun. 2015, 192, 97−107. (59) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bondorientational order in liquids and glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 28, 784−805. (60) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (61) Freiser, M. J. Ordered states of a nematic liquid. Phys. Rev. Lett. 1970, 24, 1041−1043. (62) Galerne, Y.; Marcerou, J. P. Temperature behavior of the orderparameter invariants in the uniaxial and biaxial nematic phases of a lyotropic liquid crystal. Phys. Rev. Lett. 1983, 51, 2109−2111. (63) Madsen, L. A.; Dingemans, T. J.; Samulski, E. T. Thermotropic biaxial nematic liquid crystals. Phys. Rev. Lett. 2004, 92, 145505. (64) Görtz, V.; Goodby, J. W. Enantioselective segregation in achiral nematic liquid crystals. Chem. Commun. 2005, 3262−3264. (65) Bates, M. E.; Luckhurst, G. R. Biaxial nematic phases and Vshaped molecules: A Monte Carlo simulation study. Phys. Rev. E 2005, 72, 051702. (66) Zhang, W.-Y.; Chen, J. Tennis-ball state of a self-avoiding wormlike polymer on a spherical surface. Europhys. Lett. 2011, 94, 43001.

(67) Lo Verso, F.; Yelash, L.; Egorov, S. A.; Binder, K. Interactions between polymer brush-coated spherical nanoparticles: The good solvent case. J. Chem. Phys. 2011, 135, 214902. (68) Xu, H.; Honglai, L.; Ying, H. Dynamic density functional theory based on equation of state. Chem. Eng. Sci. 2007, 62, 3494. (69) Hansen, J.-P.; McDonald, I. R. In Theory of Simple Liquids, 3rd ed.; Academic Press: New York, 2006.

O


Densely Packed Semiflexible Macromolecules in a Rigid Spherical

Recommend Documents