Phase diagrams of conformationally asymmetric pentablock copolymer melts: a theory and simulation study

Abstract

Conformational asymmetry, a difference in Kuhn length between two monomers, is known to shift the order–order transition (OOT) boundaries of block copolymer (BCP) melts and stabilize complex morphologies (e.g., Frank–Kasper phases) with enticing optical and transport properties. In this work, we investigate the influence of conformational asymmetry on the self-assembly of A_xB_yA_zB_yA_x pentablock copolymers (pentaBCPs) by adopting our previously developed workflow involving self-consistent field theory (SCFT) calculations and coarse-grained (CG) molecular dynamics (MD) simulations aided by the RAPSIDY protocol. For varying conformational asymmetry ratios, CAR ≡ β_A/β_B, where β_i is related to the Kuhn length of each block, SCFT shows shifts in the OOTs as CAR changes. For example, as CAR increases, the stability window widths of the hexagonally closed-packed cylinder and body-centered cubic spherical phases expand. To further understand the effect of CAR on chain conformations, we conduct MD simulations using two CG models – the ‘semiflexible’ chain model and the ‘unequal-bead-diameter’ (UBD) model. To vary CAR, in the ‘semiflexible’ chain model we vary the stiffness of A and B blocks, and in the UBD chain model we vary the A bead diameter with respect to the B bead diameter. The simulated phase behavior is qualitatively similar with both models and agrees to a large extent with the predicted phase diagram from SCFT. The normalized chain end-to-end distance and A–B interface width are consistent between both models and generally insensitive to CAR. However, the normalized lamellae periodicity expands with increasing conformational asymmetry (i.e., as CAR moves away from 1) for the semiflexible chain model, and the opposite is observed with the UBD chain model.

---

1. Introduction

Block copolymers (BCPs) are macromolecules composed of blocks of two or more distinct monomers. Due to the thermodynamic incompatibility between these monomers, BCPs can undergo microphase separation and form ordered nanostructures including the canonical lamellae, hexagonally packed cylinders, and double gyroid morphologies. The mechanical and transport properties exhibited by these nanostructures combined with the precise control over copolymer architecture (i.e., the size and connectivity of blocks) afforded by modern synthesis techniques have led to BCPs being used for lithography, biotech, and chemical separation applications^1–7 among others. Engineering fit-for-purpose BCP materials requires a thorough understanding of how their self-assembly into ordered morphologies and resulting properties depend on BCP design parameters. The three design parameters that describe the self-assembly behavior of the simplest case of BCPs, the A–B diblock copolymers (diBCPs), are χN, a measure of the A–B monomer segregation strength, where χ is the Flory–Huggins parameter and N is the degree of polymerization, the A monomer volume fraction, f_A, and the conformational asymmetry between A and B blocks, typically parameterized as the ratio of the A and B monomer Kuhn (statistical segment) lengths.

The segregation strength χN can be tuned by the chain length (N) and the choice of monomer chemistries which dictate the effective interactions (χ). In some cases, monomer chemistries also bring conformational asymmetry, which relates to the difference in conformational and volume-filling behavior of the monomers. Conformationally asymmetric diBCPs have been studied experimentally^8–17 and computationally.^16,18–25 Early theory-based work^10,21–24 predicted that conformational asymmetry would shift the balance of enthalpic and entropic forces that drive BCP self-assembly by lowering the entropic stretching penalty for the block with the larger Kuhn length. Since blocks within minority domains must stretch to fill the domain core, this reduced penalty drives the stiffer blocks to occupy the interior of curved domains. For instance, if A has a higher Kuhn length than B, this promotes curvature towards the A domains and stabilizes morphologies with minority A domains surrounded by a matrix domain of the B monomers. These theoretical predictions have been experimentally validated for diBCPs of varying chemistries including mixed polyolefins,^10,11,17 polyethers,^11,14 poly(isoprene)–(poly)styrene,^10,15 poly(dimethylsiloxane),^11,12 and blends thereof.

Besides theory, molecular simulations have also been employed to study the phase behavior of conformationally asymmetric diBCPs. Beránek and Posel²⁶ used dissipative particle dynamics (DPD) simulations of conformationally asymmetric diBCPs and found shifts in phase boundaries towards higher f_A as well as stabilization of perforated lamellae over double gyroid as they increased the A segment's Kuhn length while keeping the B segment's Kuhn length constant. Nikoubashman et al.¹⁸ also used DPD simulations to study the lamellar alignment of confined thin films of conformationally asymmetric diBCPs. They found that the lamellar alignment changed from being perpendicular to parallel with respect to the confining walls due to entropic effects influenced by conformational asymmetry.

In addition to shifting phase boundaries and/or changing the stability of morphologies as compared to conformationally symmetric diBCPs, conformational asymmetry has been shown to unlock novel morphologies. The Frank–Kasper σ phase was first observed in diBCP melts by Lee, Bluemle, and Bates²⁷ in 2010 and the role that conformational asymmetry played in its formation was shown in the analysis of Xie et al.²⁸ using self-consistent field theory (SCFT). It has since been repeatedly demonstrated that conformational asymmetry can stabilize the Frank–Kasper family of low-symmetry micellar morphologies in diBCP systems.^13,29,30 There is keen interest in the Frank–Kasper phases as they exhibit properties^31–35 suitable for optoelectronics. Similarly, theory^36–39 suggests that the sought after double diamond (DD) morphology^40–42 is not a metastable precursor to the double gyroid (DG) morphology for linear conformationally symmetric diBCPs; however, conformational asymmetry renders DD as a metastable alternative to DG in linear diBCPs.⁴³

Moving from diBCPs to multiblock copolymers (multiBCPs), the behavior of conformationally asymmetric BCPs with 3+ blocks has received less attention as compared to diBCPs. To address this knowledge gap, this study is focused on elucidating the effect of conformational asymmetry in linear pentablock copolymers (pentaBCPs). As in our previous investigation,⁴⁴ we study pentaBCPs with midpoint symmetry (where the first and last blocks and the second and fourth blocks are identical) to enable access to self-assembly behavior unavailable with simpler BCP designs, including different modes of chain conformations. The increased design space, however, still poses a challenge for an exhaustive experimental exploration, motivating our use of ‘in silico’ techniques. We use a combination of self-consistent field theory (SCFT) and molecular dynamics (MD) simulations in this work. We are motivated to use SCFT because of its established applicability to study BCP melt phase behavior and its computational efficiency (as compared to MD simulations). In contrast, particle-based MD simulations capture the trajectory (i.e., positions of particles as a function of time) of every component of the BCP chains, allowing for detailed analyses of blocks’ and chains’ conformations that are difficult to obtain with SCFT. Such MD simulations are, however, computationally intensive especially as the BCP design parameter space becomes large. One could reduce the computational cost of simulations to some extent by pre-screening the BCP design parameter space using faster theory calculations⁴⁴ or by biasing the melts directly into ordered morphologies of interest to identify optimal unit cell dimensions and test their (meta)stability using methods like RAPSIDY (R̲apid A̲nalysis of P̲olymer S̲tructure and I̲nverse D̲esign strategY̲).^45,46 One could also use intermediate methods like field-theoretic simulations^47–51 that incorporate density and concentration fluctuation effects, unlike SCFT, while remaining relatively inexpensive as compared to particle-based simulations.

In this paper, we combine SCFT calculations, coarse-grained (CG) MD simulations, and RAPSIDY^45,46 to assess the equilibrium morphology and chain conformations of pentaBCP melts as a function of f_A, χN, and conformational asymmetry. We first perform SCFT calculations to identify equilibrium morphologies across a wide design space followed by CG MD simulations for a select few designs. For these CG MD simulations, inspired by previous work,^{16,18,20,26,52} we first develop CG models of pentaBCPs with conformational asymmetry arising from differences in the stiffness of A and B segments or differences in the diameter of A and B CG beads. We use RAPSIDY to identify the correct simulation box dimensions for simulating the predicted ordered morphologies. We compare phase diagrams from SCFT calculations and MD simulations with both CG chain models and extract representative chain conformation distributions from each ordered morphology as a function of conformational asymmetry.

2. Approach

2.1 Conformationally asymmetric ABABA pentaBCPs

We study linear pentaBCPs with three A monomer blocks and two B monomer blocks arranged in the sequence A₁B₁A₂B₂A₃; the subscripts in this notation denote the order of A and B segments, rather than their degree of polymerization. These A₁B₁A₂B₂A₃ pentaBCPs have midpoint symmetry, where the first and last blocks (A₁ and A₃) and the second and fourth blocks (B₁ and B₂) have the same number of repeat units.

The design parameters that define an A₁B₁A₂B₂A₃ pentaBCP are as follows:

• A–B segregation strength (χN),

• the volume fraction of monomer A in the pentaBCP (f_A),

• the fraction of the middle A₂ block among all A blocks, , where N_i is the degree of polymerization of block i and N_Ai = {N_A1,N_A2,N_A3}; in this work, we only consider τ_A2 = 0.5, which corresponds to a design where the middle A₂ block contains half of the total A monomers in the pentaBCP, and

• the conformational asymmetry ratio CAR ≡ β_A/β_B where with Kuhn (statistical segment) length b and reference Kuhn segment volume v.^10,53

2.2 Self-consistent field theory (SCFT) for pentaBCP melts

We modify the canonical ensemble SCFT formalism presented in our previous paper on conformationally symmetric pentaBCPs⁴⁴ to incorporate the conformational asymmetry we focus on in this paper. We direct the reader to the SCFT formalism for the A₁B₁A₂B₂A₃ pentaBCP melt presented in our previous paper⁴⁴ or to the work of Arora et al.⁵⁴ for a more generalized SCFT formalism for multiBCP melts. Here, we briefly introduce the overall SCFT workflow and describe how we incorporate conformational asymmetry.

We consider a pentaBCP melt system with n identical A₁B₁A₂B₂A₃ pentaBCP chains. The pentaBCP chain in SCFT is modeled as a continuous Gaussian chain with a total degree of polymerization N, where the chain contour parameter s ∈ [0, N] is used to describe the forward propagator q(r, s) and the backward propagator q^†(r, s) of the polymer chain. The propagator functions are solved by the following modified diffusion equations:

with the initial conditions q(r, 0) = 1 for the free end of the A₁ block and q^†(r, N) = 1 for the free end of the A₃ block. Unlike the SCFT formalism in our previous paper⁴⁴ for the conformationally symmetric A₁B₁A₂B₂A₃ pentaBCP melts, in this work the statistical segment length b_γ(s) of a coarse-grained monomer of type γ(s) is a function of s which determines the chemical type of blocks γ = A or B for the A₁B₁A₂B₂A₃ pentaBCPs.

In our numerical computations for eqn (1) and (2), we set the statistical segment length for the B blocks as a length unit, b_B = 1, and vary the statistical segment length for the A blocks, b_A = 0.6–1.4, to incorporate the conformational asymmetry in the SCFT calculations. As noted before, the conformational asymmetry ratio between A and B is defined as CAR = β_A/β_B, where β_i = b_i/(6v_i)^1/2 – this β_i is a model-invariant parameter defined by the statistical segment length b_i and the reference monomer volume v_i.^10,53 Since the A and B monomer volumes, v_A and v_B, are equal in the SCFT formalism, the CAR parameter reduces to CAR = b_A/b_B. The rest of the formalism where the propagator functions are used to compute the A and B monomer density fields, ϕ_A(r) and ϕ_B(r), and the chemical potential fields, ω_A(r) and ω_B(r), in eqn (1) and (2), is the same as the one we used in our previous paper.⁴⁴

Once all of the self-consistent mean field solutions, ϕ_A(r), ϕ_B(r), ω_A(r), and ω_B(r), are solved, the Helmholtz free energy of the pentaBCP melt system is calculated by the following expression:

where χ is the Flory–Huggins interaction parameter, k_B is the Boltzmann constant, and T and V are the temperature and volume of the pentaBCP melt system. The total partition function Q for the entire polymer chain is obtained by volume-averaging the propagator function as .

The Helmholtz free energy is computed for each candidate phase and the homogeneous disordered phase (Dis) for identifying the equilibrium phase at each chain design; we define each chain design by its conformational asymmetry parameter, CAR = b_A/b_B, A block volume fraction, f_A, and segregation strength, χN. For the candidate phases, we include the canonical phases observed in BCP melts: body-centered cubic spheres (BCC), hexagonally packed cylinders (C₆), double gyroid (DG), Fddd orthorhombic network (O⁷⁰), and lamellae (L). Detailed information on the space group symmetry of each candidate phase structure and the spatial resolutions in the SCFT calculations are provided in Table S1 (SI). We used the C++-based open-source polymer self-consistent field (PSCF) software⁵⁴ to calculate the self-consistent mean field solutions and Helmholtz free energies and follow the same SCFT calculation details presented in our previous paper.⁴⁴

2.3 Coarse-grained (CG) molecular dynamics (MD) simulations

2.3.1 CG models of conformationally asymmetric pentaBCPs. To enable simulations of conformationally asymmetric pentaBCPs, we need to develop models that capture the molecular characteristics that give rise to the conformational asymmetry, namely differences in stiffness [e.g., iso- vs. syndiotactic poly(methylmethacrylate)^55,56] or monomer excluded volumes or “bulkiness” [e.g., BCPs of mixed polyolefins^8–11,17] between two polymer species. The stiffness and bulkiness of the monomers affect the orientational correlation along a polymer chain and the polymer Kuhn length. To capture these molecular characteristics, we use two CG models to represent conformationally asymmetric pentaBCPs. In the first CG model, which we call the ‘semiflexible chain model’, we include an angle potential to induce higher stiffness in one block as compared to the other (Fig. 1a). In the second CG model, which we call the ‘unequal-bead-diameter (UBD) chain model’, we choose different CG bead sizes for the A and B beads (Fig. 1b). For both models, the volume fraction of the A monomer , where V_bead,j is the volume of a type j bead and V_chain is the chain volume . The fraction of the block A₂ among all A blocks in the chain , and the conformational asymmetry ratio CAR = β_A/β_B, where with R_g being the radius of gyration.

Fig. 1 The two coarse-grained (CG) chain models we use to represent conformationally asymmetric A₁B₁A₂B₂A₃ pentaBCPs. N_i is the degree of polymerization for block i of the pentaBCP where i ∈ {A₁,B₁,A₂,B₂,A₃} and N_A1 = N_A3 and N_B1 = N_B2. (a) Semiflexible chain model depicting A blocks and B blocks with equal bead diameters (D_A = D_B) but different flexibility modeled using the unequal angle potential constants k_angle. (b) Unequal-bead-diameter (UBD) chain model depicting A and B blocks with different bead diameters (D_A ≠ D_B) and the same flexibility (k_angle = 0). The Kuhn lengths are defined with the reference volume (ν) equal to the B bead volume.

Next, we present the details (bonded and nonbonded potentials) for both models. After that, we describe how we calibrate our two CG models to compare the phase behavior of pentaBCPs at equivalent CAR values.

We model each pentaBCP as a chain of CG beads of the A and B monomers, hereafter referred to as A- and B-type beads, with diameters D_A and D_B, respectively. For the ‘semiflexible model’, D_A = D_B = 1d, while for the ‘UBD chain model’ D_A may be greater or less than 1d and D_B = 1d. In both models, the bonded interaction between two adjacent beads of a pentaBCP chain is captured using

U_bond(r) = k_bond(r − r₀)² (4)

where k_bond is the bond force constant and is equal to 50k_BT/d² for all bonds, r is the bead center-to-center distance, and r₀ is the bond equilibrium length, which we define as for bonded beads of types i and j.

In the case of ‘UBD chain model’ pentaBCPs, we keep the entire chain flexible and do not use an angle potential. In the ‘semiflexible chain model’ we define an angle potential between triplets of bonded beads as

U_angle,ijk = k_angle(1 + cos(θ_ijk)) (5)

where k_angle,ijk is the angle potential coefficient for beads i, j, and k, and θ_ijk is the angle formed by the vectors between beads i and j and between j and k. We find k_angle,ijk as the average of k_angle assigned to beads i, j, and k, and we outline how we choose k_angle for A and B beads for specific CAR in the next section.

In both models, the interaction between nonbonded beads is modeled using the cut-and-shifted Lennard-Jones (LJ) potential

where ε_ij and σ_ij are the LJ well-depth and particle diameter parameters between beads of types i and j and r_cut is the potential cutoff distance. We also define σ_ij and r_cut,ij using the bead diameters: and r_cut,ij = σ_ij + 1d.

We relate χ to the LJ well-depth parameters of interactions between beads of the same (ε_AA, ε_BB) and different (ε_AB) types with the following linear expression:

where z is the ensemble-average coordination number, which we take to be 6 as noted in the literature^57–59 and as confirmed in our simulations for this work. We set ε_AA = ε_BB ≥ ε_AB and adjust ε_AB to model microphase-separation induced by thermodynamically incompatible A and B monomers. We note that this mapping of χ to ε is convenient and a variety of other rigorous methods exist for relating χ to model parameters.^60,61 However, these methods were developed for diBCP systems and may not be accurate for our pentaBCP melts. Furthermore, if we were to do the rigorous mapping we would need map χ for each value of f_A, CAR, and the LJ ε parameters making it an extensive task.

2.3.2 Calibrating the conformational asymmetry ratio using stiffness or bead diameter. Inspired by previous simulation studies of the influence of semiflexibility^{20,52,62–65} and bead size^62,63,65 on polymer conformations, we adjust the angle potential strength (k_angle) and bead diameters (D), respectively, to tune the A and B Kuhn lengths and reach a target CAR. To calibrate how these parameters relate to CAR, we performed homopolymer melt simulations in which we varied either k_angle or D and sampled the resulting conformational asymmetry parameter . We present the resulting data graphically in Fig. 2 and in a tabular form along with the homopolymer chain designs in Table S2 (SI).

Fig. 2 Conformational asymmetry ratios β/β_ref for homopolymer melts for (a) varying angle potential coefficient k_angle in the semiflexible chain model (β/β_kangle=0kBT) or (b) varying bead diameter D in the UBD chain model (β/β_D=1d). (a) All semiflexible homopolymer chains have a degree of polymerization N = 50 with beads with D = 1d. (b) For homopolymer chains with unequal bead diameters, we specify N and D such that V_chain is equal to that of a chain with N = 50 and D_bead = 1d. Homopolymer chain designs and β/β_ref are tabulated in Table S2 (SI). We simulate each homopolymer melt for 2 × 10⁶ timesteps in the NPT ensemble (T* = 1, P* = 1) and sample each quantity every 5 × 10⁴ timesteps over the final 5 × 10⁵ timesteps of this simulation. All error bars are ±1 standard deviation of the sampled quantity.

To model semiflexible pentaBCPs, we keep D_A = D_B = 1d and apply an angle potential to only one bead type while leaving the other flexible. We use our β/β_kangle=0kBT vs. k_angle correlation in Fig. 2a to find k_angle for the stiffened bead type to achieve a desired CAR. Specifically, for CAR > 1, A beads are stiffened (β_A > β_kangle=0kBT) and B beads are left flexible (β_B = β_kangle=0kBT) and therefore, CAR ≡ β_A/β_B = β/β_kangle=0kBT. Similarly, for CAR < 1, CAR ≡ β_A/β_B = (β/β_kangle=0kBT)⁻¹.

To model UBD pentaBCPs, we use the β/β_D=1d vs. D relation shown in Fig. 2b to find the appropriate A bead diameter D_A for the desired CAR. As we maintain D_B = 1d, CAR ≡ β_A/β_B = β/β_D=1d. We then set r_0,A and σ_A equal to D_A.

Using our calibrated ‘semiflexible chain model’ and ‘UBD chain model’, we perform simulations of conformationally asymmetric pentaBCPs to sample equilibrium morphologies and chain conformations and compare the phase behavior at equivalent CAR. We list the pentaBCP designs for the semiflexible chain model in Table S3 (SI) and for the UBD chain model in Table S4 (SI). Briefly, we choose pentaBCP designs with a total number of beads N ≈ 50 and τ_A2 = 0.5 for each combination of the parameters varied on the phase diagrams: f_A = 0.32, 0.38, and 0.46, and CAR = 0.6, 0.8, 1.0, 1.2, and 1.4.

Our CAR vs. D correlation (Fig. 2b) indicated that a UBD pentaBCP with CAR = 0.6 would have large A beads (D_A ≅ 2d) and would likely need to have end (A₁ and A₃) and middle (A₂) blocks of only 1 and 2 beads, respectively, to ensure that N ≅ 50 and τ_A2 ≅ 0.5. For simplicity, we chose UBD pentaBCP designs with D_A = 0.6, 0.8, 1, 1.2, and 1.4d as these nearly cover the target CAR range (≅0.8–1.4) and avoid any pentaBCP designs with single-bead blocks. We chose UBD pentaBCP designs that matched our target N, f_A, and τ_A2, and we added extra designs with D_A = 0.6 and 0.8d to increase the resolution of certain phase boundaries in the UBD phase diagram.

2.3.3 Traditional MD simulation protocol. We simulate melts of pentaBCPs with the LAMMPS software package⁶⁶ with a simulated annealing protocol as described in our previous work.⁴⁴ To initialize well-mixed melts, we relax the chains with the harmonic potential for bonded interactions and the cut-and-shifted LJ (with ε_AA = ε_AB = ε_BB = 0.5k_BT) or soft push-off potentials for nonbonded interactions to increase chain mobility. We shrink the simulation box to reach a cubic box with a side length of 30d and then equilibrate in the NVT ensemble (T* = 1) using the Nosé–Hoover thermostat^67,68 to finish initialization. The number of chains in each melt is chosen so that the bead volume fraction ∼0.45, consistent with previous work,^57,59,69 in the cubic box with side lengths = 30d. We use a simulation timestep of Δt = 0.01τ, where τ is the LJ unit of time (not to be confused with our τ_A2 pentaBCP design parameter) for all pentaBCP designs except those with the smallest A beads (i.e., D_A = 0.6d, CAR ≅ 1.4) for which we use a smaller time step of Δt = 0.005τ for stability.

To begin the simulated annealing protocol, we switch to the NPT ensemble by adding the Nosé–Hoover barostat to maintain T* = 1 and P* = 1 and increase the LJ interactions to full strength. We keep ε_AA = ε_BB = 1k_BT and change ε_AB according to eqn (7) to model a specific χN. We simulate each melt sequentially at χN = 90, 150, and 200 to mimic a cooling (annealing) protocol that would allow morphologies to develop over time in real-world melts. We simulate at each χN for a total of 8 × 10⁵τ and define the last quarter of this period as the production stage, from which we sample morphologies and chain conformations. The time required for equilibration, analyzed using the profiles of box side length and averaged chain end-to-end distance, R_ee, is sensitive to pentaBCP semiflexibility and bead diameter disparity. We find that our simulation time is sufficient to equilibrate all systems except the stiffest semiflexible pentaBCPs with CAR = 0.6. For this case of semiflexible pentaBCPs with CAR = 0.6, we simulate each melt for 8 × 10⁵τ at χN from 40 to 90 in increments of Δ(χN) = 10 before continuing to χN = 150 and 200.

We sample melt configurations every 5 × 10⁴τ, which we identified as sufficient time to sample decorrelated configurations for all simulated pentaBCPs; the decorrelation time was calculated using the autocorrelation of the chain end-to-end distance for flexible chains or the block end-to-end distance for the flexible blocks in the semiflexible chains. We determine the morphology in each sampled frame using an isosurface drawn at the A–B interface and from the characteristics of the R_ee distribution, as outlined in Section 2.3.5.

We run three replicate trials for each simulated pentaBCP design to ensure our results are consistent and not kinetically trapped.

2.3.4 RAPSIDY-informed simulation box size. A known challenge in simulating BCP self-assembly is the potential mismatch between an ordered morphology's unit cell and the simulation box dimensions.^48,70–72 In short, the simulation box is required to accommodate an integer multiple of an ordered morphology's repeating unit for that morphology to be stable in the simulation. Some phases with one or two periodic dimensions (e.g., lamellae, cylinders) can overcome this mismatch by “tilting” such that their periodic axes are not parallel to the box axes. Triply periodic morphologies such as the double gyroid (DG) or Fddd (O⁷⁰) networks, however, cannot adapt in the same fashion and tend to form frustrated morphologies; therefore, there is a need to efficiently determine suitable box dimensions to alleviate frustration of morphology self-assembly. The R̲apid A̲nalysis of P̲olymer S̲tructure and I̲nverse D̲esign Strategy̲ (RAPSIDY)^45,46 workflow is used to identify the suitable box dimensions for simulating BCPs in any ordered morphology. While we choose to use RAPSIDY, we note that there are other approaches to address incommensurability, including independently varying all box dimensions^73,74 and calculating the correct unit cell dimensions for “tilted” morphologies.⁷⁵ We utilized RAPSIDY whenever we found “flawed” or “frustrated” versions of conventional ordered morphologies; the details of the RAPSIDY workflow are presented in the “RAPSIDY workflow” section of the SI. In short, we use RAPSIDY to first bias pentaBCP melts into a target morphology, then allow the melt to relax without biasing, and finally, evaluate the biased structure's stability in the target morphology after relaxation. In this work, we use this procedure to identify the unit cell sizes for the pentaBCP designs for target ordered phases, when the standard simulation results in “flawed” ordered structures. In the next section, we describe how we identify morphologies as conventional or frustrated/flawed.

2.3.5 Identifying morphologies with A–B interface and chain conformations. We identify ordered morphologies using A–B interface visualization and quantification of chain conformations. The ordered morphologies we observe from CG MD simulations include hexagonally packed cylinders (C₆), network (N), double gyroid (DG), and lamellae (L). We identify the disordered microphase (DM) as one where there are distinct microphase-separated domains that appear disordered as they fluctuate in shape and location.

We visualize the A–B interface using the isosurface function implemented in the OVITO software⁷⁶ to help us identify ordered morphologies. We use the Gaussian density method⁷⁷ to calculate a voxel grid of bead density, where the particle-to-mesh transfer is performed using a Gaussian function with a standard deviation equal to the bead diameter. We specify a density threshold such that the isosurface is drawn where the A bead densities are equal to ∼45–55% of the bulk density.

While ordered phases such as L or C₆ can be easily identified with the isosurface visualization, morphologies with multiple interpenetrating networks like N and DG require more detailed analysis. We label a melt's morphology as DG only if the ordered networks we observe have the true DG structure, including two gyroid structures with three-fold nodes; otherwise we classify the morphology N. Previously, we attempted to distinguish N and DG using their structure factors,⁴⁴ but we could not differentiate between these similar structures in this study due to fluctuations in the configurations and the liquid state of the melt.

We quantify the distribution of chain end-to-end distances R_ee and find that the distributions of ordered phases including L, C₆, DG, and N have two or more distinct peaks (modes), which correspond to the categorically different conformations adopted by their constituent chains;⁴⁴ in contrast the DM phase has a single, rounded peak. The R_ee distributions we sample for the semiflexible and UBD chain models are shown with the rest of our supplementary simulation results in Fig. S10 and S11 (SI).

We identify the morphologies for each simulation run for a given pentaBCP design based on the isosurface visualization and the distribution of chain end-to-end distances. In cases where we identify degenerate morphologies across separate simulation trials of one pentaBCP design, we report as equilibrium morphology the one we observe in the majority of trials.

2.3.6 A–B interface width and domain sizes. To analyze the effect of conformational asymmetry on the A–B interface width and domain sizes of CG pentaBCPs, we examine the A and B composition profiles of several melts of similar or identical pentaBCP designs and varying CAR but exhibiting the L morphology. We define the A–B interface width and A and B domains using the volume fraction of A beads Φ_A and B beads Φ_B along the normal axis of the lamellae, as shown in Fig. 3a. We show representative Φ_A and Φ_B profiles in Fig. 3b and highlight the intervals corresponding to the A–B interface and the lamellae unit cell, which includes a pair of A and B domains. We define the interface width L_int as the distance along the normal axis between the points where Φ_A = 0.2 and Φ_A = 0.8. We use the points where Φ_A = Φ_B = 0.5 to define the boundaries of A and B domains, and we report the lamellae periodicity L_LAM as the combined width of pairs of adjacent A and B domains. For comparison, we sample L_LAM,q* = 2π/q*, where q* is the microphase peak's q-value in the structure factor S(q) for each melt. We sample each quantity every 5 × 10⁴τ in the final 2 × 10⁵τ of simulation in one trial of each selected pentaBCP design. We exclude lamellae in the first and last 10% of the normal axis from our analysis to avoid edge effects.

Fig. 3 (a) Snapshot of semiflexible chain model pentaBCPs with f_A = 0.46 and CAR = 1 at χN = 150 in the lamellae ‘L’ phase, with the normal axis labeled with a black arrow. (b) Volume fractions of A (red line) and B beads (blue line) along the lamellae normal axis, with an example interfacial region and lamellae unit cell highlighted.

3. Results

3.1 SCFT phase diagrams

In Fig. 4, we present the equilibrium phase diagrams obtained from SCFT calculations for the A₁B₁A₂B₂A₃ pentaBCP melts with τ_A2 = 0.5 at χN = 40 (Fig. 4a) and χN = 60 (Fig. 4b). First, we see that for both χN, the OOT boundary shifts to higher f_A values as CAR increases from 0.6 to 1.4. One can explain this shift in OOT by extending the intuitive explanation provided by Matsen using strong stretching theory.⁷⁸ The formation of ordered structures in BCP melts is explained by the spontaneous curvature adopted by BCPs as a result of a balance of A and B block stretching energies. For compositionally asymmetric BCPs (e.g., f_A < 0.5), the minority block tends to be easily stretched by curving the A–B interface toward the minority domains so that the majority blocks can be relaxed and the total stretching energy is minimized. However, the selection of an ordered morphology is also governed by the incompressibility constraint of melts that requires space to be filled at uniform density. Therefore, the position of OOT boundaries is determined according to which space-filling structure (e.g., L, DG, C₆, and BCC) provides the average interfacial curvature that most closely matches the preferred curvature of BCPs at given f_A and CAR. The degree of spontaneous curvature increases as the A to B composition becomes more asymmetric and favors more curved structures. This explains the conventionally observed phase sequence L → DG → C₆ → BCC as f_A deviates from f_A = 0.5. In Fig. 4, we observe the conventional sequence of ordered phases (BCC → C₆ → networks (DG or O⁷⁰) → L) as f_A increases for all CAR values except at CAR = 0.6, where an inverted DG is formed at f_A = 0.5 and χN = 40 (Fig. 4a). The shift of OOT boundaries as CAR varies indicates that as the Kuhn length of the A block increases (i.e., as CAR increases), curved structures are more thermodynamically favored in balancing the stretching energy between A and B blocks. This is because A blocks with longer Kuhn lengths suffer a lower stretching energy penalty when stretched to form curved domains, which is explained by the polymer theory with Gaussian chain models for which the stretching energy is inversely proportional to b².⁷⁸

Fig. 4 Phase diagrams of A₁B₁A₂B₂A₃ pentaBCP melts of τ_A2 = 0.5 from the SCFT calculations at (a) χN = 40 and (b) χN = 60. Each symbol represents the equilibrium phase identified by SCFT, including homogeneous disorder (Dis), body-centered cubic spheres (BCC), hexagonally packed cylinders (C₆), double gyroid (DG), Fddd (O⁷⁰), and lamellae (L) phases.

The next observation in Fig. 4a is that as CAR is increased, the stability window widths of the C₆ and BCC phases expand significantly while the ODT boundary positions change insignificantly with f_A. In addition, increasing CAR significantly alters the position of the stability window for the network phases, including DG and O⁷⁰, but the phase window width remains relatively small compared to that for the C₆ phase; this can be attributed to a higher degree of packing frustration for BCPs when forming complex network phases such as DG and O⁷⁰ than the classical phases such as BCC, C₆, and L.⁷⁸ Interestingly, as CAR decreases to 0.6, the DG window shrinks while the O⁷⁰ window does not change, suggesting that the O⁷⁰ morphology provides less packing frustration than the DG morphology when the structure forming blocks (A blocks) have a shorter Kuhn length than the matrix forming blocks (B blocks).

At χN = 60 (Fig. 4b), we see the same trend of shifting OOT boundaries with CAR as we see for χN = 40 (Fig. 4a), but with a larger DG-forming stability window and the L stability window extended to low f_A ≈ 0.32 at CAR = 0.6. The disordered phases (Dis) at χN = 40 and f_A < 0.34 in Fig. 4a change to the C₆ and DG phases at χN = 60 in Fig. 4b. Specifically, at the same f_A values, as χN increases from 40 to 60, the Dis phases at CAR > 1 change to C₆ phases while the Dis phases at CAR < 1 change to the C₆ or DG phases. This again reflects the tendency to form morphologies with high interfacial curvature when A blocks have longer Kuhn lengths than B blocks (CAR > 1).

A SCFT phase diagram of diBCPs plotted with respect to CAR and f_A, and at χN = 40 was reported by Bates et al.³⁰ Compared to Fig. 4a, the ODT boundaries are positioned at much lower values of f_A ≈ 0.1 in contrast to f_A ≈ 0.34 in Fig. 4a. This is because of the dispersed sequence of A and B blocks in pentaBCPs, so a higher χN is required for microphase separation of pentaBCPs as compared to the diBCPs. The dramatic shift of ODT boundaries from χN = 40 to χN = 60 in Fig. 4 indicates that the phase diagrams are in a weak segregation regime, while diBCPs above χN = 40 correspond to a relatively strong segregation regime.

3.2 CG MD simulation phase diagrams

To confirm if some or all of the SCFT predictions with varying CAR are observed in the presence of excluded volume and concentration fluctuations, we conduct CG MD simulations for select pentaBCP designs, first using the ‘semiflexible chain model’. In Fig. 5, we plot the phase diagram constructed by CG MD simulations using the semiflexible chain model; the pentaBCP designs studied are listed in Table S3 (SI). We note that the χN values we study in our MD simulations are higher than in our SCFT calculations because of the increase in (χN)_ODT from compositional fluctuations, which are present in simulation but not represented in the mean-field technique we use. Because of this intrinsically different behavior, we focus only on the qualitative correspondence between SCFT and MD, specifically on how OOT boundaries vary with CAR and the resulting phase sequences.

Fig. 5 Phase diagrams of A₁B₁A₂B₂A₃ pentaBCP melts of τ_A2 = 0.5 from the MD simulations using the ‘semiflexible chain model’ for pentaBCPs at (a) χN = 90 and (b) χN = 150. Each symbol represents the self-assembled morphologies, including hexagonally packed cylinders (C₆), double gyroid (DG), lamellae (L), network (N), and disordered microphase (DM). The asterisk on the DG symbol on the chain design with f_A = 0.38 and CAR = 1.4 indicates that the equilibrium DG morphologies were formed in RAPSIDY-informed simulation boxes.

Generally, in CG MD simulations using the ‘semiflexible chain model’ for A₁B₁A₂B₂A₃ pentaBCP melts of τ_A2 = 0.5, we see that as CAR increases, the OOT boundaries shift to higher values of f_A. In agreement with the SCFT phase diagrams, as CAR increases, the morphology with higher interfacial curvature is preferred. For example, the ordered phase sequence L → N → C₆ as CAR is increased at fixed f_A = 0.32 (Fig. 5b) is consistent with the ordered phase transitions in the SCFT phase diagram (see Fig. S6 (SI) for the constant f_A lines overlaid on Fig. 4).

In the CG MD simulations at χN = 90 (Fig. 5a), we observe the DG, N, and L ordered morphologies. Interestingly, at CAR = 0.6 for two chain designs (f_A = 0.38 and f_A = 0.46) we observe the inverted N morphology, where the B blocks form the network structure despite being the majority blocks, and the minority A blocks form the surrounding matrix. To confirm that the inverted N morphologies are not kinetically trapped morphologies which would form L at equilibrium, we performed additional MD simulations with the initial chain configurations biased into the L phase with the two most stable lamellae periodicities (12d and 14d) as predicted by RAPSIDY in Fig. S1a (SI). Except for the simulations with f_A = 0.38 and lamellae periodicity 12d, all of these additional MD simulations transitioned from the L phase to the inverted N morphology during the simulated annealing procedure, confirming that the inverted N phases at CAR = 0.6 in Fig. 5a are not kinetically trapped phases.

At low f_A = 0.32, we observe the L phase for CAR = 0.6 and the disordered microphase (DM) morphology at other CAR values (Fig. 5a). At f_A = 0.38, however, we only observe DM phases at CAR = 1.0 and 1.2. We attribute the decrease in (χN)_ODT at high and low CAR to the higher degree of compositional fluctuations for the more conformationally symmetric chain designs (i.e., CAR ≈ 1.0) than for the conformationally asymmetric chain designs (i.e., CAR = 1.4 and 0.6). As either A or B type blocks become stiffer than the other type (i.e., as one approaches CAR = 1.4 or 0.6), the number of neighboring chains interacting with a single polymer chain increases. In fluctuation theory, (χN)_ODT is inversely proportional to the total number of interacting neighboring polymers in a polymer's pervaded volume, which is associated with the invariant degree of polymerization and degree of compositional fluctuations.^79,80 In Fig. S2 (SI), we quantify the number of interacting neighbor polymers, interpreted as the chain coordination number, and confirm that the chain designs with CAR = 1.0 and 1.2 have lower coordination numbers than the chain designs with other CAR values, supporting the observation of higher (χN)_ODT. On the other hand, for SCFT phase diagrams in Fig. 4a, we observe relatively consistent ODT boundary positions as CAR is varied. Considering that the mean field approximation in SCFT ignores the effect of compositional fluctuations, the difference in (χN)_ODT across CAR in SCFT calculations is expected to be smaller than in MD simulations where the relative strength of fluctuations varies significantly with CAR.

At χN = 150 (Fig. 5b), all the chain designs form ordered phases, and for most of the CAR values, we observe the phase sequence C₆ → N or DG → L as f_A is increased, consistent with the SCFT phase diagrams (Fig. 4). Exceptionally, for CAR = 0.6, we observe a transition from a flat curvature structure (L) to a moderate curvature structure (inverted N) as f_A increases. This L → inverted N transition for a semiflexible chain model pentaBCP is analogous to the L → inverted DG transition observed as f_A increases in SCFT at the same CAR value.

At χN = 150 (Fig. 5b), f_A = 0.38, and CAR = 1.4, the original three simulation trials gave three different ordered morphologies (perforated lamellae, N, and C₆), and an additional three simulation trials also did not show any consistent morphology (Fig. S3, SI). We assume that a simulation box size mismatch with the unit cell likely led to these frustrated or inconsistent morphologies. This motivated us to use the RAPSIDY workflow in Section 2.3.4 to identify the correct simulation box. Using the simulation box size (26d) informed by the RAPSIDY stability profile (Fig. S1b, SI) we observe well-ordered DG morphologies from two out of three trials. The DG morphologies from SCFT calculations and from MD simulations are compared in Fig. S4 (SI). The characteristic features of the DG structure from SCFT calculations (double networks in Ia3̄d space group symmetry and the three-fold connectors) are also observed in the DG morphology from MD simulations.

At χN = 200 (Fig. S5c, SI), the phase diagram remains quite similar to χN = 150, suggesting that in this high segregation regime the effect of increasing χN on self-assembled morphology is negligible.

Next, in Fig. 6, we present the phase diagrams obtained from CG MD simulations using the ‘UBD chain model’ for the pentaBCP designs listed in Table S4 (SI). In the UBD model, we change CAR by changing the A bead diameter D_A, so the same pentaBCP designs that we used for semiflexible pentaBCPs could not be used for UBD pentaBCPs for all CAR. Additionally, CAR = 0.8 is the lowest CAR value we tested with UBD pentaBCPs due to the challenges of having large A beads for CAR = 0.6, as noted in Section 2.3.2.

Fig. 6 Phase diagrams of A₁B₁A₂B₂A₃ pentaBCP melts of τ_A2 = 0.5 from the MD simulations of unequal-bead-diameter (UBD) chain model pentaBCPs at (a) χN = 90 and (b) χN = 150. The symbols represent the same morphologies as in Fig. 4 and 5. Note that the range in CAR plotted along the vertical axis is different from the previous phase diagrams. Equilibrium morphologies which were formed in RAPSIDY-informed simulation boxes are marked with an asterisk.

In the χN = 90 phase diagram (Fig. 6a) we observe similar trends in morphology across the phase diagram (varying CAR or f_A) to those from SCFT calculations (Fig. 4a) and MD simulation of CG semiflexible pentaBCPs (Fig. 5a). In Fig. 6a, we observe the ordered morphologies C₆, N, and L, as well as the DM morphology. Increasing CAR shifts OOT boundaries to stabilize greater curvature towards A domains. For instance, the equilibrium morphology transitions from L → C₆ as CAR = 1.2 → 1.4 at f_A ≅ 0.38 and from L → N as CAR = 1.0 → 1.4 at f_A ≅ 0.46. We note that the equilibrium morphology (C₆) of the pentaBCP design marked with an asterisk (f_A = 0.39, CAR = 1.4) is determined after additional trials with RAPSIDY-informed simulation box dimensions (see details below and in Fig. S8 (SI)). As previously explained, we attribute the variation in (χN)_ODT to a reduction in compositional fluctuations resulting from conformational asymmetry. Interestingly, while increasing the conformational asymmetry of semiflexible pentaBCPs decreases their (χN)_ODT, for UBD pentaBCPs, increasing CAR > 1 (decreasing D_A) lowers (χN)_ODT but decreasing CAR < 1 does not reduce (χN)_ODT to the same degree.

At χN = 90, the semiflexible pentaBCP with f_A = 0.38 and CAR = 0.8 forms L, while UBD pentaBCPs with CAR = 0.8 do not order at any f_A we tested. The lack of ordering at low CAR is partially explained by the low average chain coordination number (Fig. S7d, SI) for the UBD pentaBCPs with CAR = 0.8. We observe that the chain coordination number changes non-monotonically with CAR; for instance, the UBD pentaBCP with CAR = 1.0 and f_A = 0.46 has a higher coordination compared to those with CAR = 0.8 or 1.4 and comparable f_A values of 0.44 and 0.46, respectively.

At χN = 150 (Fig. 6b), the fully ordered UBD pentaBCP phase diagram bears greater resemblance to the corresponding phase diagrams from SCFT (χN = 60, Fig. 4b) and MD simulation of semiflexible pentaBCPs (χN = 150, Fig. 5b). All UBD pentaBCPs with CAR ≤ 1.0 formed L at χN = 150, while most of these pentaBCPs formed DM at χN = 90. Conversely, curvature towards A domains becomes more favorable as CAR increases: N is first observed at f_A = 0.32 with CAR = 1.0, at f_A = 0.37 with CAR = 1.2, and at f_A = 0.42 with CAR = 1.4. Likewise, the highest f_A at which the C₆ phase is observed is f_A = 0.34 with CAR = 1.2 and f_A = 0.39 with CAR = 1.4.

At both χN = 90 and 150, for the UBD pentaBCP with CAR = 1.4 and f_A = 0.39 (denoted with an asterisk in Fig. 6) the true equilibrium morphology is identified using the RAPSIDY workflow as outlined in Section 2.3.4 to screen for the optimal unit cell dimensions of the C₆ and DG phases. After screening unit cell dimensions using RAPSIDY, we perform three trials using the box dimensions optimized for either morphology following our simulated annealing protocol (Section 2.3.4). We report the stability profiles for the screened C₆ and DG lattice constants in Fig. S8 (SI).

At χN = 200, as shown in Fig. S7c (SI), the diagram is mostly unchanged from the one at χN = 150 (Fig. S7b, SI), similar to what we saw with the semiflexible pentaBCPs. The UBD pentaBCP with CAR = 1.2 and f_A = 0.34 is the only design that changes phase, going from C₆ → N as χN increases from 150 to 200.

In summary, we can conclude the following from the SCFT and CG MD simulations’ phase diagrams: the influence of conformational asymmetry on pentaBCP equilibrium morphology seen in the SCFT phase diagrams is qualitatively consistent with that predicted in past literature for the phase behavior of diBCPs;^23,78 as CAR increases, the OOT boundaries shift to higher f_A values as a result of the increased spontaneous curvature toward the A domains. A similar phase behavior is observed in the phase diagrams of the MD simulations for both CG models, particularly at high χN = 150, where the morphology transitions and overall shifts of OOT boundaries are consistent with SCFT predictions. At lower segregation strength (χN = 90), however, the ODT boundary positions are appreciably different between both CG models and vary significantly as CAR changes; in contrast, the ODT as predicted by SCFT is relatively insensitive to CAR, which we attribute to the lack of compositional fluctuations in the mean-field method.

3.3 CG MD simulations – chain conformations, A–B interface width, and domain sizes

Next, we seek to understand how CAR affects chain conformations, A–B interface width, and microphase-separated domain sizes of systems with the same morphology and similar chain designs.

From the CG MD simulations of both the semiflexible chain model and UBD chain model we identify several systems that exhibit the same morphology for similar pentaBCP designs. In Fig. 7, we present the A–B interface width L_int, lamellae periodicity L_LAM, and chain end-to-end distance R_ee for four semiflexible (abbreviated in Fig. 7 as ‘SF’) and three UBD pentaBCP designs in the lamellae morphology, with the design parameters listed in the caption. To account for the effect of chain length on conformational properties, we normalize each quantity by contour length L_c, which is the R_ee of a fully extended chain. We provide the original L_int, L_LAM, and R_ee values as well as the A and B domain widths for the semiflexible pentaBCPs in Fig. S12 (SI) and for UBD pentaBCPs in Fig. S13 (SI).

Fig. 7 (a) A–B interface widths L_int, (b) lamellae periodicities L_LAM, and (c) chain end-to-end distance R_ee normalized by contour length L_c, sampled at χN = 150 for semiflexible (‘SF’) pentaBCPs (black circles with solid lines) and UBD pentaBCPs (green squares with dashed lines). The semiflexible pentaBCPs have CAR = 0.8, 1, 1.2, and 1.4, and all have the same chain design with f_A = 0.45 and L_c = 52d. The UBD pentaBCPs have different chain designs with f_A ≅ 0.46 and L_c = 58.8, 53.2, and 44.8d for chains with CAR = 0.8, 0.9, and 1.2, respectively. Each quantity is sampled as described in Section 2.3.6. Error bars are ±1 standard deviation.

In Fig. 7a we observe that the normalized A–B interface widths L_int/L_c of all seven pentaBCPs lie close together, nearly all within each other's uncertainty intervals, and do not significantly vary with CAR. The discrepancy in L_int/L_c between semiflexible and UBD pentaBCPs is greatest at CAR = 0.8 but decreases as CAR increases, such that both L_int/L_c values at CAR = 1 are nearly identical. Similarly, R_ee/L_c (Fig. 7c) changes little with CAR, and the average R_ee/L_c values of all semiflexible and UBD pentaBCPs that we examine here are within each other's confidence intervals. In contrast, in Fig. 7b we observe opposing trends between L_LAM/L_c and CAR between both chain models. For the semiflexible pentaBCPs, L_LAM/L_c is minimized for conformationally symmetric chains and increases as either bead type is stiffened. For the UBD pentaBCPs, L_LAM/L_c decreases as the A bead diameter D_A is either increased or decreased from 1d to adjust CAR. We do not believe that this disagreement between the two models is due to the use of Φ_A and Φ_B profiles to capture these domain characteristics as we observe good agreement between L_LAM,volfrac and L_LAM,q* values of each chain model in Fig. S12 and S13 (SI). We are unable to explain this disagreement between the two models regarding the influence of CAR on the ratio of domain size to contour length.

4. Conclusions

In this work, we studied the self-assembly behavior of melts of A₁B₁A₂B₂A₃ pentaBCPs with conformational asymmetry between A and B monomers using SCFT calculations and CG MD simulations. We conducted CG MD simulations using one of two models – semiflexible vs. unequal bead diameter – to capture conformational asymmetry in monomers due to differences in either stiffness or bulkiness in monomers, respectively. In the semiflexible chain model, we modeled a flexible A block (β_A = β_kangle=0kBT) and a semiflexible B block (β_B > β_kangle=0kBT) for CAR ≡ β_A/β_B < 1, and vice versa for CAR ≡ β_A/β_B > 1. In the UBD chain model, we always keep β_B fixed and only increase or decrease β_A through the bead diameter of A. We compared the phase diagrams of self-assembled morphologies for A₁B₁A₂B₂A₃ pentaBCPs with the middle A₂ block representing half of all A blocks in the pentaBCP (τ_A2 = 0.5) from the SCFT calculations and the MD simulations using the two models. We also compared the chain conformations, A–B interface widths, and domain sizes for the MD simulations using the two models.

Both CG models’ phase diagrams exhibit stabilization of curved phases at high CAR and the shift of OOT boundary locations. The agreement between the phase diagrams of both models and the SCFT phase diagram (Fig. 4) is best at higher χN. To facilitate this direct comparison for the reader, in Fig. S14 (SI) we present these three phase diagrams side by side. To calculate f_A using the actual bead volume in simulation, which may vary with chain composition and CAR, we replot the phase diagrams of semiflexible and UBD pentaBCPs in Fig. S15 and S16 (SI), respectively, using bead diameters sampled using radial distribution functions (RDFs). We find that sampling f_A with RDFs does not change the trends between CAR and morphology we observed in the original phase diagrams (Fig. 5 and 6).

Beyond phase behavior, we analyzed the effect of CAR on A–B interface width, domain size, and chain conformations for the two models as chain design and morphology are held constant. For both the semiflexible model and the UBD model, the value of CAR has little effect on the A–B interface width or chain extension (R_ee). However, the lamellae periodicity normalized by contour length, L_LAM/L_c, shows opposing trends for the two models with varying CAR; we are unable to explain this disagreement between the two models.

We find that the effects of CAR on pentaBCP morphology and chain conformations are consistent with previous work on conformationally asymmetric diBCPs and multiBCPs. In our SCFT calculations and MD simulations with both CG models, we see an increase in the stability of high-curvature morphologies (e.g., C₆) as CAR increases, thereby increasing the minority block (A block) Kuhn length. As a result, phase boundaries shift to higher f_A as CAR increases, as is seen in many theory-based and experimental studies of conformationally asymmetric BCP behavior.

For computational researchers who may want to choose between the two coarse-grained models, besides the trends in phase behavior and chain conformations, the computational cost is a key consideration. The computational cost of simulating semiflexible pentaBCPs is correlated with their conformational asymmetry; e.g., stiffening either bead type results in slower equilibration and the need for longer simulations. Increasing bead size and thereby lowering CAR for the UBD model results in fewer beads being needed to fill a given volume, which can increase simulation efficiency. Conversely, this increases the computational cost of simulating UBD pentaBCPs with high CAR as we achieve this by decreasing D_A. This could instead be addressed by increasing D_B to increase CAR (changing either β_A or β_B, as in the semiflexible chain model), although one should anticipate less ideal (i.e., less universal) self-assembly behavior for chains composed of fewer CG beads.

This study also provides some practical information for choosing conformationally asymmetric pentaBCPs in experiments. For example, the morphology transitions observed as CAR is changed suggest that if one has a flexible block and another semiflexible block, pentaBCPs with a semiflexible minority block and a flexible majority block (CAR > 1) will have a higher tendency to form spontaneous curvature towards the minority domains. Our results will guide those looking to engineer nanostructures with high- and low-curvature domain shapes of varying A–B interface width and domain sizes with a single conformationally asymmetric BCP chemistry by selecting which component forms the minority block.

Author contributions

Tristan Myers: conceptualization, data curation, formal analysis, investigation, methodology, software, visualization, writing – original draft, writing – review & editing; So Jung Park: conceptualization, data curation, formal analysis, investigation, methodology, software, visualization, writing – original draft, writing – review & editing; Arthi Jayaraman: conceptualization, funding acquisition, methodology, project administration, resources, supervision, writing – review & editing.

Conflicts of interest

The authors declare that there are no known conflicts of interest.

Data availability

The supplementary information (SI) includes: the list of SCFT candidate morphologies; the designs of all semiflexible and UBD pentablock chains; a detailed outline and results of our RAPSIDY workflow; all diagrams of equilibrium phase and chain coordination number for both semiflexible and UBD pentaBCPs; a comparison of the double gyroid (DG) morphology in SCFT and in MD simulation; chain conformation distributions for semiflexible and UBD pentaBCP chains; data on the sensitivity of chain end-to-end distance and lamellae domain size and interface width to conformational asymmetry; and MD phase diagrams adjusted for bead partial molar volume. See DOI: https://doi.org/10.1039/d5sm01223f.

The input files of the PSCF program for the SCFT calculations and the final configurations of one trial of CG MD simulation for each pentablock copolymer design are available in a Zenodo repository at https://doi.org/10.5281/zenodo.17137334.

Acknowledgements

The authors acknowledge financial support from Multi-University Research Initiative (MURI) grant from the Army Research Office, Award Number W911NF2310260. The SCFT calculations and CG MD simulations were run on the Caviness supercomputing cluster at the University of Delaware (UD).

References

1. F. S. Bates and G. H. Fredrickson, Phys. Today, 1999, 52, 32–38.
2. Y. Mai and A. Eisenberg, Chem. Soc. Rev., 2012, 41, 5969–5985.
3. C. M. Bates, M. J. Maher, D. W. Janes, C. J. Ellison and C. G. Willson, Macromolecules, 2014, 47, 2–12.
4. M. Karayianni and S. Pispas, J. Polym. Sci., 2021, 59, 1874–1898.
5. T. Lohmuller, D. Aydin, M. Schwieder, C. Morhard, I. Louban, C. Pacholski and J. P. Spatz, Biointerphases, 2011, 6, MR1-12.
6. C. M. Bates and F. S. Bates, Macromolecules, 2017, 50, 3–22.
7. D. Woo, H. Yoon and J. K. Kim, J. Polym. Sci., 2024, 62, 639–661.
8. M. Sikka, N. Singh, F. S. Bates, A. Karim, S. Satija and C. F. Majkrzak, J. Phys. II, 1994, 4, 2231–2248.
9. M. Sikka, N. Singh, A. Karim, F. S. Bates, S. K. Satija and C. F. Majkrzak, Phys. Rev. Lett., 1993, 70, 307–310.
10. F. S. Bates, M. F. Schulz, A. K. Khandpur, S. Forster, J. H. Rosedale, K. Almdal and K. Mortensen, Faraday Discuss., 1994, 98, 7–18.
11. K. Almdal, M. A. Hillmyer and F. S. Bates, Macromolecules, 2002, 35, 7685–7691.
12. Z. M. Gdowski, S. M. Swartzendruber, M. K. Mahanthappa and F. S. Bates, Macromolecules, 2024, 57, 7926–7940.
13. M. W. Schulze, R. M. Lewis, J. H. Lettow, R. J. Hickey, T. M. Gillard, M. A. Hillmyer and F. S. Bates, Phys. Rev. Lett., 2017, 118, 207801.
14. N. Zhou, T. P. Lodge and F. S. Bates, J. Phys. Chem. B, 2006, 110, 3979–3989.
15. C. Lai, W. B. Russel, R. A. Register, G. R. Marchand and D. H. Adamson, Macromolecules, 2000, 33, 3461–3466.
16. R. J. Sanchez-Leija, J. A. Mysona, J. J. de Pablo and P. F. Nealey, Macromolecules, 2024, 57, 2019–2029.
17. E. W. Cochran and F. S. Bates, Macromolecules, 2002, 35, 7368–7374.
18. A. Nikoubashman, R. A. Register and A. Z. Panagiotopoulos, Macromolecules, 2014, 47, 1193–1198.
19. M. A. Carignano and I. Szleifer, Europhys. Lett., 1995, 30, 525–530.
20. K. Hagita and T. Murashima, Comput. Mater. Sci., 2024, 243, 113105.
21. J. D. Vavasour and M. D. Whitmore, Macromolecules, 1993, 26, 7070–7075.
22. M. W. Matsen and M. Schick, Macromolecules, 1994, 27, 4014–4015.
23. M. W. Matsen and F. S. Bates, J. Polym. Sci., Part B:Polym. Phys., 1997, 35, 945–952.
24. E. Helfand and Z. R. Wasserman, Developments in Block Copolymers, Elsevier, New York, 1982.
25. E. Reister, M. Müller and S. K. Kumar, Macromolecules, 2005, 38, 5158–5169.
26. P. Beránek and Z. Posel, J. Nanosci. Nanotechnol., 2016, 16, 7832–7835.
27. S. Lee, M. J. Bluemle and F. S. Bates, Science, 2010, 330, 349–353.
28. N. Xie, W. Li, F. Qiu and A.-C. Shi, ACS Macro Lett., 2014, 3, 906–910.
29. S. Lee, C. Leighton and F. S. Bates, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 17723–17731.
30. M. W. Bates, J. Lequieu, S. M. Barbon, R. M. Lewis, 3rd, K. T. Delaney, A. Anastasaki, C. J. Hawker, G. H. Fredrickson and C. M. Bates, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 13194–13199.
31. W. Man, M. Megens, P. J. Steinhardt and P. M. Chaikin, Nature, 2005, 436, 993–996.
32. M. Florescu, S. Torquato and P. J. Steinhardt, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 80, 155112.
33. M. C. Rechtsman, H. C. Jeong, P. M. Chaikin, S. Torquato and P. J. Steinhardt, Phys. Rev. Lett., 2008, 101, 073902.
34. A. P. Hynninen, J. H. Thijssen, E. C. Vermolen, M. Dijkstra and A. van Blaaderen, Nat. Mater., 2007, 6, 202–205.
35. M. A. Klatt, P. J. Steinhardt and S. Torquato, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 23480–23486.
36. M. W. Matsen and M. Schick, Phys. Rev. Lett., 1994, 72, 2660–2663.
37. M. W. Matsen and F. S. Bates, Macromolecules, 1996, 29, 7641–7644.
38. A. Reddy, M. S. Dimitriyev and G. M. Grason, Nat. Commun., 2022, 13, 2629.
39. M. S. Dimitriyev, X. Feng, E. L. Thomas and G. M. Grason, Phys. Rev. Lett., 2024, 132, 218201.
40. C. Y. Chang, G. M. Manesi, C. Y. Yang, Y. C. Hung, K. C. Yang, P. T. Chiu, A. Avgeropoulos and R. M. Ho, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2022275118.
41. X. Feng, M. S. Dimitriyev and E. L. Thomas, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2213441120.
42. W. Shan and E. L. Thomas, ACS Nano, 2024, 18, 9443–9450.
43. M. S. Dimitriyev, B. R. Greenvall, R. Mathew and G. M. Grason, ACS Macro Lett., 2025, 14, 1291–1298.
44. S. J. Park, T. Myers, V. Liao and A. Jayaraman, Mol. Syst. Des. Eng., 2024, 9, 1235–1253.
45. V. Liao and A. Jayaraman, JACS Au, 2025, 5, 2810–2824.
46. V. Liao, T. Myers and A. Jayaraman, Soft Matter, 2024, 20, 8246–8259.
47. J. Lequieu, J. Chem. Phys., 2023, 158, 244902.
48. J. Huh, C. Park and Y. K. Kwon, J. Chem. Phys., 2010, 133, 114903.
49. G. H. Fredrickson and K. T. Delaney, Field-Theoretic Simulations in Soft Matter and Quantum Fluids, Oxford University Press, Oxford, 2023.
50. T. M. Beardsley and M. W. Matsen, J. Chem. Phys., 2021, 154, 124902.
51. M. W. Matsen, T. M. Beardsley and J. D. Willis, Phys. Rev. Lett., 2023, 130, 248101.
52. T. Y. Pan, B. B. Patel, D. J. Walsh, S. Dutta, D. Guironnet, Y. Diao and C. E. Sing, Macromolecules, 2021, 54, 3620–3633.
53. F. S. Bates and G. H. Fredrickson, Macromolecules, 1994, 27, 1065–1067.
54. A. Arora, J. Qin, D. C. Morse, K. T. Delaney, G. H. Fredrickson, F. S. Bates and K. D. Dorfman, Macromolecules, 2016, 49, 4675–4690.
55. O. N. Tretinnikov and K. Ohta, Langmuir, 1998, 14, 915–920.
56. O. N. Tretinnikov, Macromol. Symp., 2000, 149, 269–275.
57. G. S. Grest, M. D. Lacasse, K. Kremer and A. M. Gupta, J. Chem. Phys., 1996, 105, 10583–10594.
58. M. Murat, G. S. Grest and K. Kremer, Macromolecules, 1999, 32, 595–609.
59. H. Guo and K. Kremer, J. Chem. Phys., 2003, 119, 9308–9320.
60. A. Petrov, H. Huang and A. Alexander-Katz, Macromolecules, 2024, 57, 8212–8222.
61. J. D. Willis, T. M. Beardsley and M. W. Matsen, Macromolecules, 2020, 53, 9973–9982.
62. C. Svaneborg and R. Everaers, Macromolecules, 2020, 53, 1917–1941.
63. R. Everaers, H. A. Karimi-Varzaneh, F. Fleck, N. Hojdis and C. Svaneborg, Macromolecules, 2020, 53, 1901–1916.
64. R. Faller, A. Kolb and F. Müller-Plathe, Phys. Chem. Chem. Phys., 1999, 1, 2071–2076.
65. R. K. W. Spencer and M. W. Matsen, Macromolecules, 2022, 55, 1120–1126.
66. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
67. S. Nose, Mol. Phys., 1984, 52, 255–268.
68. W. G. Hoover, Phys. Rev. A:At., Mol., Opt. Phys., 1985, 31, 1695–1697.
69. G. S. Grest and K. Kremer, Phys. Rev. A:At., Mol., Opt. Phys., 1986, 33, 3628–3631.
70. F. J. Martinez-Veracoechea and F. A. Escobedo, J. Chem. Phys., 2006, 125, 104907.
71. J. Škvor and Z. Posel, Macromol. Theory Simul., 2015, 24, 141–151.
72. A. Arora, D. C. Morse, F. S. Bates and K. D. Dorfman, Soft Matter, 2015, 11, 4862–4867.
73. A. J. Schultz, C. K. Hall and J. Genzer, J. Chem. Phys., 2004, 120, 2049–2055.
74. Y. Seo, J. R. Brown and L. M. Hall, Macromolecules, 2015, 48, 4974–4982.
75. Y. Feng, B. Li and Q. Wang, Soft Matter, 2022, 18, 4923–4929.
76. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
77. M. Krone, J. E. Stone, T. Ertl and K. Schulten, Eurographics Conference on Visualization, 2012.
78. M. W. Matsen, J. Phys.:Condens. Matter, 2002, 14, R21–R47.
79. J. Glaser, P. Medapuram, T. M. Beardsley, M. W. Matsen and D. C. Morse, Phys. Rev. Lett., 2014, 113, 068302.
80. P. Medapuram, J. Glaser and D. C. Morse, Macromolecules, 2015, 48, 819–839.

---

Footnote

† T. M. and S. J. P. contributed equally to the article.

---