Self-consistent field theory and coarse-grained molecular dynamics simulations of pentablock copolymer melt phase behavior

Abstract

Block copolymer (BCP) self-assembly leads to nanostructured materials with diverse ordered morphologies, some of which are attractive for transport applications. Multiblock AB copolymers are of interest as they offer a larger design parameter space than diblock copolymers allowing researchers to tailor their self-assembly to achieve target morphologies. In this study, we investigate the phase behavior of symmetric A_xB_yA_zB_yA_x and B_xA_yB_zA_yB_x pentablock copolymers (pentaBCPs) where A and B monomers have the same statistical segment length. We use a combination of self-consistent field theory (SCFT) calculations and molecular dynamics (MD) simulations to link the polymer design parameters, namely the fraction of middle block volume to the volume of all blocks of same type, τ, overall volume fraction of A block, f_A, and segregation strength, χN, to the equilibrium morphologies and the distributions of chain conformations in these morphologies. In the phase diagrams calculated using SCFT, we observe broader double gyroid windows and the existence of lamellar morphologies even at small values f_A in contrast to what has been seen for diblock copolymers. We also see a reentrant phase sequence of double gyroid → cylinder → lamellae → cylinder → double gyroid with increasing τ at fixed f_A. The chain conformations adopted in these morphologies are sampled in coarse-grained MD simulations and quantified with distributions of the chain end-to-end distance and fractions of chains whose middle (A or B) and end (A or B) blocks remain within domains of same chemistry (A or B). These analyses show that the pentaBCP chains adopt “looping”, “bridging”, and “hybrid” (both looping and bridging) conformations, with a majority of the chains adopting the hybrid conformation. The spatial distributions for each of the blocks in the pentaBCPs show that blocks of the same type in a chain locally segregate within the same domains, with shorter blocks segregating towards the domain boundaries and longer blocks filling the domain interior. This combined SCFT-MD approach enables us to rapidly screen the extensive pentaBCP design space to identify design rules for transport-favorable morphologies as well as verify the chain conformations and spatial arrangements associated with the theory predicted reentrant phase behavior.

---

Design, System, Application

Block copolymers (BCPs) self-assemble into a variety of nanostructures, such as lamellae, hexagonal-packed cylinders, and double gyroid, which in turn enable engineering of materials with desired transport and mechanical properties. The morphology formed by a given BCP is highly dependent on its design in terms of the monomer chemistry, number of blocks (diblock to multiblock), block lengths, and block sequence. As compared to diblock copolymers, multiblock copolymers have been studied to a smaller extent due to their larger design parameter space. However, it is noteworthy that a handful of computational studies of multiBCPs have uncovered novel nanostructures and phase behavior not seen in the well-studied diblock copolymers. In this study, we focus on linking the design of pentablock copolymers (pentaBCPs) to their morphology in the melt state using self-consistent field theory (SCFT) and molecular dynamics (MD) simulations. Our goal is to identify design rules for experimentalists looking to broaden phase windows of transport-friendly morphologies, such as double gyroid. The combination of theory and simulation allows for faster screening of design parameter space using SCFT as compared to MD simulations and quantification of chain conformations using MD simulations especially when getting distribution of chain conformations from SCFT is non-trivial.

I. Introduction

Block copolymers (BCPs) self-assemble into ordered periodic nanostructures through microphase separation arising from the chemical incompatibility and connectivity between the blocks.^1,2 The rich phase diagram makes BCPs ideal to use as materials for applications in next-generation memory devices, membranes, sensors, photonic crystals, and solar cells.^3–8 For the simplest case of linear AB diblock copolymer (in short, diBCP), the equilibrium morphology depends on three parameters – χN, f_A, and ε, where χ is Flory–Huggins parameter quantifying the segregation strength between A and B monomers, N is degree of polymerization, f_A is A block composition, and ε is the conformational asymmetry between two blocks.^1,2,9 According to mean field theory, the order-to-disorder transition (ODT) for a linear AB diBCP occurs at a critical value of (χN)_ODT = 10.5 with the volume fractions of A and B component equal (f_A = 0.5).¹⁰ Above the ODT, the ordered (canonical) morphologies including lamellae (L), hexagonal-packed cylinders (C₆), double gyroid (DG), and body-centered cubic (BCC) spheres are accessed in order of increasing asymmetry in copolymer composition. The chains self-assemble into these morphologies to minimize free energy by minimizing entropy loss associated with chain stretching and maximizing energetically favorable A–A and B–B contacts (enthalpic gain). The curvature of the A domain–B domain interface is determined by balance between stretching energies of two blocks, under the constraints associated with filling space at uniform density and minimizing packing frustration.^11,12

While the simplest case of linear diBCP has been studied extensively with experiments, theory, and simulations, there is a need to establish the above concepts – ODT, assembled equilibrium morphologies, order-to-order transition (OOT), adopted chain conformations within these morphologies – for other copolymer designs. With significant advances in polymer synthesis in the past couple of decades, researchers have been able to go beyond the simple linear diBCPs to complex multi-block and tapered sequences^13,14 as well as non-linear polymer architectures (e.g., star, cyclic, coil-brush, bottlebrush).^15–19 As these new polymers are created, it has also motivated researchers to understand how changing sequence and architecture impacts the assembled morphologies at equilibrium or during processing. Depending on the number of blocks (e.g., triblock, tetrablock, pentablock) and the sequence of A and B monomers (e.g., BAB vs. ABA for triblock), the phase diagrams for linear multiblock copolymers (in short multiBCPs) exhibit subtle or major differences from the linear AB diBCP phase diagram. For example, for symmetric ABA triBCPs, the overall topology of phase diagram remains similar to AB diBCP because an ABA triBCP if snipped in the middle forms two diBCPs.^20,21 There are, however, some clear differences between the ABA triBCP and AB diBCP: (i) at the same composition, f_A, the ODT is shifted to lower temperatures for ABA triBCPs as compared to AB diBCPs with a higher critical value of (χN)_ODT ≈ 18 compared to that of diBCPs, (χN)_ODT = 10.5,²² and (ii) while the morphologies accessed by ABA triBCPs are same as those of AB diBCPs, unlike the AB diBCPs, the ABA triBCPs adopt either looping (both end A blocks in the same A domains) or bridging (each A block in different A domains) conformations. Two different shifts of ODT to lower temperatures were found via predictions from self-consistent field theory (SCFT)²³ and experiments where a relatively short A end block, when synthetically added to the B block of the AB diBCP, causes the ODT to be lower than the parent AB diBCP.²⁴ For asymmetric (A)_N1(B)_NB(A)_N2 copolymers, an interesting phase behavior was predicted based on the SCFT calculations by Matsen;²³ he found that with varying parameter τ = N₁/(N₁ + N₂), where N₁ and N₂ is degree of polymerization of A₁ block and A₂ block, the OOT phase boundary shifts nonmonotonically. This nonmonotonic shift of phase boundaries results in reentrant morphological transitions (e.g., lamella → double gyroid → cylinder → double gyroid → lamella) with decreasing τ at fixed f_A. This reentrant morphological change (cylinder → lamella → cylinder) was also experimentally observed with asymmetric polystyrene-block-polybutadiene-block-polystyrene triBCP.²⁵

In the case of ABAB tetraBCPs, specifically styrene-b-isoprene-b-styrene-b-isoprene tetraBCPs, experiments showed that the ODT temperature is considerably lower than that of comparable AB diBCPs.²⁶ The phase diagrams of ABAB tetraBCPs constructed by SCFT theory^27,28 and molecular dynamics (MD) simulation²⁹ include canonical morphologies (spheres, cylinders, and lamellae) as observed in simpler AB diBCPs. However, compared to ABA triBCPs, existence of an additional tail B block in ABAB tetraBCP leads to the emergence of morphologies not seen with ABA triBCPs. For example, a hybrid lamella-sphere structure phase, where lamellar and spherical domains of the same component are separated by domains of the other component,²⁷ and Frank–Kasper (FK) phases, a group of complex spherical packing phases,²⁸ are predicted by the SCFT. The OOT phase behavior of ABAB tetraBCPs is also significantly different from that of ABA triBCPs. The reentrant phase transition (lamellae → double gyroid → hexagonal-packed cylinders → lamellae) was predicted by SCFT calculations of ABAB tetraBCPs with increasing ratio of the end A block length to total length of both A blocks, at a symmetric A–B volume fraction (f_A = f_B = 0.5).²⁷ The same reentrant morphological transition was experimentally observed with polystyrene-b-polyisoprene-b-polystyrene-b-polyisoprene (SISI) tetraBCPs.³⁰ It is notable that in these tetraBCPs, the cylinder and double gyroid phases are stable even at symmetric composition (A volume fraction = B volume fraction) of the copolymer unlike the simpler AB linear BCPs which only form lamellar phases at symmetric A and B volume fraction. In another experiment of SISI tetraBCPs,³¹ an unusual transition from cylindrical to lamellar morphology was observed upon increasing temperature; this is opposite to the OOT from lamella to cylinder morphology observed in AB diBCPs or ABA triBCPs. All of the above work demonstrates that presence of an additional B block in ABAB tetraBCP as compared to the ABA triBCP produces different phase behavior.

Compared to the ABA triBCPs and the ABAB tetraBCPs, experimental and computational studies of A and B containing pentaBCPs are sparse. For example, Zhao et al. calculated the SCFT phase diagrams of symmetric BABAB pentaBCPs, where the lengths of the two end B blocks are the same and the lengths of two A blocks are the same.³² The phase diagram with respect to χN and f_A exhibits (χN)_ODT > 30, at symmetric composition, f_A = 0.5, in contrast to diBCP which has an (χN)_ODT = 10.5, and symmetric triBCP, (χN)_ODT ≈ 18. For these BABAB pentaBCPs, the same order of canonical morphologies (L → DG → C₆ → BCC) was observed with decreasing f_A as in the AB diBCPs and ABA triBCPs. Interestingly, for these symmetric BABAB pentaBCPs, the phase diagram with respect to f_A and τ (where τ is the volume fraction of two end B blocks among all B blocks) exhibits, with decreasing τ, a reentrant phase behavior of body-centered cubic (BCC) and face-centered cubic (FCC) spherical phases (BCC → FCC → BCC → FCC → BCC).³² The reentrant phase behavior was supported by the SCFT density distributions of blocks, but the explicit chain conformations responsible for the reentrant phase behavior have not been verified yet.

In another recent SCFT study, it was shown that by optimizing design parameters in pentaBCP, i.e. individual length ratios of constituent blocks, one can achieve nontraditional phases that are difficult to access through simpler multiBCPs.^33,34 SCFT calculations from Xie et al. proposed a strategy to stabilize single gyroid in BABAB pentaBCPs where the gyroid network structure is formed by the two A blocks.³³ The constituent block lengths in pentaBCP was precisely tuned in the way that the stability of single gyroid is induced by the synergetic effect of lowered packing frustration and stretching energy of the bridging B blocks. The packing frustration in B domains of single gyroid is relieved through the local segregation of the two tail B blocks of asymmetric length, filling near and far space from the A–B interface. Further, the double gyroid which competes thermodynamically with single gyroid is destabilized by the stretched bridging B block connecting the two networks. Moreover, the equilibrium phase diagram of the BABAB pentaBCP with respect to the two design parameters – degree of asymmetry in length of two tail B blocks and length ratio of bridging B block to total B blocks – shows a stable phase of a square array of cylinders as compared to the canonical hexagonal-packed cylinders.^33,34

Morphologies where the domains have a lower coordination number (e.g., square vs. hexagonal-packed cylinders or single gyroid vs. double gyroid) have attracted attention due to potential applications in photonic crystals⁵ and next-generation integrated circuits,³ and have been attempted to be accessed using more complex branched BCPs^16,35,36 or blending different types of linear BCPs.^3,37 The prior pentaBCP studies suggest that linear AB-type multiBCP designs have the potential of offering noncanonical structures with fewer molecular design parameters to consider, i.e. block sequence and relative length ratios between blocks, as compared to the branched BCP or the multicomponent BCP systems. Moreover, with experimentally tractable design parameters, systematic analysis of pentaBCP systems can give us insight into how BCP self-assembly behavior evolves from AB diBCP, ABA triBCP to AB-type multiBCP. Experimentally, symmetric pentaBCPs are considered easier to synthesize than tetraBCPs because symmetric pentaBCPs grown via difunctional initiators require fewer chain extensions than tetraBCPs grown from monofunctional initiators.³⁸ As a result, study of pentaBCPs offers a more practical approach to understanding BCP self-assembly behavior.

In this paper, we investigate the equilibrium phase behavior of symmetric linear pentaBCPs with alternating sequences of A and B blocks – A_xB_yA_zB_yA_x and B_xA_yB_zA_yB_x – using a combined approach of SCFT calculation and MD simulation. We first present the SCFT calculations that lead to phase diagrams of the A_xB_yA_zB_yA_x and B_xA_yB_zA_yB_x pentaBCPs at different χN and various end vs. middle block volume fractions. Using a coarse-grained model for the pentaBCPs in MD simulations with chains of finite lengths and analogous design parameters as used in SCFT, we quantify distributions of chain conformations and spatial organization of the blocks within the A and B domains to explain the observed phase behavior in the SCFT phase diagram.

II. Approach

A. Pentablock copolymer (pentaBCP) design

We consider melts of linear A_xB_yA_zB_yA_x and B_xA_yB_zA_yB_x pentaBCPs with the design space shown in Fig. 1. For the remainder of this paper, we will denote these two cases of pentaBCPs as A₁B₁A₂B₂A₃ and B₁A₁B₂A₂B₃ pentaBCPs, identifying the order of the blocks in the subscript rather than the degree of polymerization of the block.

Fig. 1 Schematics of (a) symmetric A₁B₁A₂B₂A₃ pentaBCP and (b) symmetric B₁A₁B₂A₂B₃ pentaBCP, where N_γ is degree of polymerization of γ type block in the polymer chain ((a) γ = {A₁, B₁, A₂, B₂, A₃}, and (b) γ = {B₁, A₁, B₂, A₂, B₃}), and the A and B monomers have the same monomer volume. The polymer design parameters for each pentaBCP are defined below the schematics. The f_A is the overall volume fraction of A blocks in the polymer chain, and τ_A2 is the volume fraction of the middle A₂ block to the total A blocks for the A₁B₁A₂B₂A₃ pentaBCP, and τ_B2 is the volume fraction of the middle B₂ block to the total B blocks for the B₁A₁B₂A₂B₃ pentaBCP.

We only consider the case where the degree of polymerization of the two end blocks are equal and the degree of polymerization of the second and forth blocks are equal. In other words, N_A1 = N_A3 and N_B1 = N_B2 for the A₁B₁A₂B₂A₃ pentaBCPs, and N_B1 = N_B3 and N_A1 = N_A2 for the B₁A₁B₂A₂B₃ pentaBCPs, where N_γ is defined as degree of polymerization of γ type block in the polymer chain. We also assume that the A and B monomers have the same statistical segment length, b_A = b_B; the analogous cases where b_A ≠ b_B will be the focus of a future publication.

Our chosen designs reduce the number of parameters for the block compositions to two independent parameters, the overall volume fraction of A block, f_A, and the fraction of the middle A₂ (or B₂) block to the total A (or B) blocks for the A₁B₁A₂B₂A₃ block (for the B₁A₁B₂A₂B₃) copolymers, denoted as τ_A2 and τ_B2 as shown in Fig. 1.

In short, we study the phase behavior of linear A₁B₁A₂B₂A₃ and B₁A₁B₂A₂B₃ pentaBCP melts by varying the polymer design parameters, f_A, τ_A2, τ_B2, and segregation strength between A and B monomers, χN.

B. SCFT formalism for pentaBCP melt

We use SCFT formalism³⁹ to predict equilibrium phase behavior of pentaBCP melts in a system volume V. The polymers are modeled as flexible Gaussian chains with a total degree of polymerization N and A and B monomers with equal statistical segment length b and segment volume ν. With this model, the polymer design parameter τ_A2 (or τ_B2) is interpreted as the number ratio of statistical segments in the middle bridging blocks, A₂ (or B₂), to the total A (or B) blocks for A₁B₁A₂B₂A₃ (or B₁A₁B₂A₂B₃) pentaBCPs.

As we are considering single-component melts, we use canonical ensemble SCFT for calculating Helmholtz free energies of the pentaBCP melts. In the SCFT formalism, the partial partition function q(r, s) of the chain fragment is calculated by solving the following modified diffusion equation (eqn (1)). The parameter s is a continuous parameter s ∈ [0, N] indicating the specific segment along the chain contour with the free first segment denoted as s = 0 and the last segment as s at position r in q(r, s). The modified diffusion equation is written as follows:

with an initial condition q(r, 0) = 1. Here, ω_γ(r) is the chemical potential field acting on a segment of type γ = {A, B} positioned at r, where the type γ depends on the variable s. The SCFT formalism provides self-consistent mean field equations to determine the mean fields ω_A(r) and ω_B(r), which we explain later. The conjugate partial partition function q^†(r, s) for the remaining chain fragment, which has the free end (s = N) and the other end fixed at position r, is obtained by solving the following modified equationwith initial condition q^†(r, N) = 1.

From the computed partial partition functions, the density of monomers belonging to type γ blocks is determined by

where γ = {A₁, B₁, A₂, B₂, A₃} for the A₁B₁A₂B₂A₃ pentaBCPs, and γ = {B₁, A₁, B₂, A₂, B₃} for the B₁A₁B₂A₂B₃ pentaBCPs. Here, Q is the total partition function of the polymer chain, which is expressed as . The incompressibility condition requires the total density of monomers at position r satisfies .

In SCFT, the mean field solutions for the densities and chemical potential fields satisfy the following equations,

ω_A(r) = χϕ_B(r) + ξ(r) (4)

ω_B(r) = χϕ_A(r) + ξ(r) (5)

where χ is the Flory–Huggins interaction parameter between A and B type segment and ξ(r) is a Lagrangian multiplier pressure field to enforce the incompressibility of the polymer melts expressed as ϕ_A(r) + ϕ_B(r) = 1, where ϕ_A(r) and ϕ_B(r) is the total density of segment type A and B, respectively. After finding the mean field solutions, the Helmholtz free energy of the system per segment is calculated from the following expression,where k_B is the Boltzmann factor and T is the temperature of the system. Since the SCFT polymer chain model is invariant to the absolute values of N and b, we choose N = 1 and b = 1, which is the standard convention in SCFT formalism.³⁹ With this choice, the total number of segments in the system, V/ν, is equal to the total number of chains of length N, and therefore, the Helmholtz free energy per segment is interpreted as Helmholtz free energy per chain of length N.

C. SCFT calculations

The calculations of the above SCFT formalism for pentaBCPs are performed using the open-source polymer self-consistent field (PSCF) software.³⁹ The diffusion equations (eqn (1) and (2)) are numerically solved using the pseudo-spectral method⁴⁰ with periodic boundary conditions. For the integration, we use the integration step of ds ≈ 0.005 and a large number of spatial grid points (Table 1) to keep accuracy of numerical solutions as high as possible while keeping computational costs reasonable. To find the self-consistent mean field solutions, the diffusion equations are iteratively solved with chemical potential fields ω_A(r) and ω_B(r), which are updated by the self-consistent mean field equations (eqn (4) and (5)) for the next iteration cycle. The Anderson mixing iteration scheme is implemented when updating the fields iteratively to accelerate the convergence to the self-consistent mean field solutions.⁴¹ The iterative computation stops when the errors in the self-consistent mean field equations (eqn (4) and (5)) are less than the specified error tolerance. The SCFT simulation details implemented in PSCF software are provided in ref. 39.

Table 1 The candidate phases considered in the SCFT calculations, structure details, and spatial discretization. Body-centered cubic sphere (BCC), hexagonal-packed cylinders (C₆), double gyroid (DG), Fddd orthorhombic (O⁷⁰), and lamellae (L) morphologies are presented in the table

Name	Structure	Space group	Crystal system	SCFT grid size
BCC		Im3̄m	Cubic	64 × 64 × 64
C6		p6mm		64 × 64
DG		Ia3̄d	Cubic	64 × 64 × 64
O^70		Fddd	Orthorhombic	32 × 64 × 128
L				128

---

To determine the equilibrium phase at given polymer design parameters, the Helmholtz free energies for the ordered phases, lamellae (L), double gyroid (DG), Fddd orthorhombic network (O⁷⁰), hexagonal-packed cylinders (C₆), and body-centered cubic spheres (BCC), and homogeneous disordered phase (Dis), are compared (ESI† Fig. S2). The morphologies of the candidate phases we consider in SCFT calculations are provided in Table 1. We did not consider non-classical phases such as single gyroid and square array of cylinders; those phases have been predicted by SCFT to be stable in AB-type pentaBCP melts with highly asymmetric end blocks.^33,34 For the pentablock copolymer designs we consider in this work, the two end blocks have the same length where such non-classical phases are not favored.

D. Coarse-grained (CG) model used in molecular dynamics (MD) simulations

To complement SCFT morphological predictions with quantitative information of chain conformations in the observed equilibrium morphologies, we use CG MD simulations.

We model pentaBCPs as CG bead-spring chains composed of two types of beads, A and B, where each bead represents a Kuhn segment of the copolymer.⁴² Such models have been used extensively in past studies of MD simulations of copolymers as well as polymer blends to predict melt morphology as a function of polymer composition, architecture, and total chain length.^42–50 To match the SCFT design space, we maintain both bead types to be of the same statistical segment size 1d; all other simulation lengths are normalized to this 1d.

Pairs of bonded beads in each chain have a harmonic bond potential of the form

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

where k_bond is the force constant = 50k_BT/d² in all cases, r is the distance between the centers of the two bonded beads, and r₀ is the equilibrium length equal to 1d.

Pairs of non-bonded beads interact with each other via the cut-and-shifted Lennard-Jones (LJ) potential⁵¹ which has the form

in which the i and j subscripts refer to the type of non-bonded beads in the pair, ε_ij is the depth of the potential well, σ is the bead diameter taken to be 1d for both bead types, and r_cut is the potential cutoff distance of 2d. The depths of the self- and cross-interaction potentials between the monomers (ε_AA, ε_AB, and ε_BB) are related to χ viawhere z is the average bead coordination number within the first coordination shell defined as beads' center-to-center distance of ≤1.5d. The value of the z parameter is taken to be 6, guided by reported values in literature.^42,44,49,52 To model BCPs with positive values of χ (for mimicking segregation between A and B monomers), we choose the values of ε_AB ≤ ε_AA = ε_BB.

E. MD simulation

System sizes. It is a known challenge that the simulated morphology of BCPs can be altered by the simulation box size and its relationship with the expected periodicity in the ordered morphology.^53–56 To ensure the morphology we observe is the true equilibrium morphology and not one obtained due to our selected system size, we utilize three system sizes for each pentaBCP design.

In the first system, we simulate 400 chains of each polymer design at χN = 60 within a cubic simulation box with side length L^MD_target calculated as

where n_chains is the number of chains in simulation (n_chains = 400 for this first simulation set), N is the degree of polymerization (corresponding to the number of beads per chain in simulation), and η is the desired volume fraction of beads within the simulation box volume, taken to be 0.45 to represent melt conditions.^42,44,49,52

In the second system, the chain end-to-end distance R^MD_ee sampled at χN = 60 from the first system is used to specify the target cube side length L^MD_target. We take inspiration from the SCFT calculations in which the length of the unit cell L^SCFT_uc is expressed as

L^SCFT_uc = C^SCFT_ucR^SCFT_ee (11)

with R^SCFT_ee = bN^1/2 being the ideal chain end-to-end distance and C^SCFT_uc a dimensionless unit cell parameter used in SCFT calculations. The latter is used as a proportionality constant to specify the unit cell in MD simulations as L^MD_uc = C^SCFT_ucR^MD_ee in analogy to eqn (11) and using R^MD_ee sampled for each chain design from the first system. The side length L^MD_target used in the second system is an integer multiple of L^MD_uc. This L^MD_target is then used to specify n_chains according to eqn (10). We set the integer ratio between L^MD_target and L^MD_uc such that n_chains ≥ 400. We simulate each melt at χN from 40 to 120 in increments of 10, with the final simulation state at each χN used as the initial state of the next value of χN. We use this simulated annealing procedure to prevent kinetic trapping of the morphology.

The third system size is a cube of size L^MD_target which is a multiple of L_q* = 2π/q*, corresponding to the position of the major peak of the A domain structure factor S_AA(q*) sampled at χN = 120 for each chain design in the second system. We found that the L^MD_uc from the second system was often mismatched with L_q*, indicating that the simulation box size was incommensurate with the morphology periodicity. We quantify this discrepancy between the unit cell size found in SCFT and MD by comparing C^SCFT_uc with the MD unit cell parameter C^MD_uc, found as the quotient of L_q*, from the second system, and R^MD_ee from the first system. We found that C^MD_uc corresponded to the C^SCFT_uc for chain designs SCFT predicted would form the lamella (L) morphology, but often disagreed for those predicted to form C₆ and DG (cylinders and double gyroid) morphologies; this prompted us to use the third system with L^MD_uc = for chain designs which are predicted by SCFT calculations to form DG at χN = 60, and with L^MD_uc = L_q* if any other morphology is predicted. We set the integer ratio between L^MD_target and L^MD_uc such that n_chains ≥ 1000 to reduce the impact of finite-size effects. We simulate increasing χN from 40 to 120 in increments of 10 to sample the pentaBCP morphology and chain conformations.

For readers interested in the values of C^SCFT_uc used for the second system size, the R^MD_ee and L_q* sampled from the first and second system sizes, and C^MD_uc calculated from the results of the first and second system sizes, we provide this information in Table S1 in the ESI.†

Simulation details. We use the LAMMPS software package⁵⁷ to perform these MD simulations of melts of pentaBCPs using the CG model described earlier in each of the above system sizes. We use a timestep size of 0.01τ (where τ is the dimensionless LJ unit of time) for all stages of simulation. We use the MolTemplate software program to initialize the chains and specify an initial simulation volume;⁵⁸ we then generate each initial configuration by first shrinking each dimension of the simulation box to the L^MD_target under the NVT ensemble at T* = 1 (where T* is the reduced LJ temperature) using the Nosé–Hoover thermostat^59,60 over 5 × 10⁶ timesteps. At this stage, we set the depth of the potential well for each LJ interaction (ε_AA, ε_AB, ε_BB) to 0.25k_BT to permit the chains to mix (as segregation strength is 0). Once the target density is reached, we simulate each melt at constant volume in a cubic simulation volume of side length L^MD_target for another 5 × 10⁶ timesteps. This is followed by increasing ε_AA, ε_AB, and ε_BB to 0.5k_BT and continuing the simulation for another 5 × 10⁶ timesteps. These steps are simply meant to create a well-mixed melt configuration of the pentaBCPs.

We then simulate the melts at each segregation strength χN in the NPT ensemble (with T* = 1, P* = 1) using the Nosé–Hoover thermostat and barostat. We maintain a cubic simulation box by coupling all three side lengths. We set ε_AA and ε_BB equal to 1k_BT and select ε_AB so as to achieve our target χ using eqn (9). We first equilibrate the melt for 1.5 × 10⁷ timesteps, which we verify to be sufficient by ensuring the convergence of the profile of mean-squared internal distances:

where x₀ and x_i are the positions of the first and the i-th bead of each chain, respectively, and n_bonds is the number of bonds separating the first and the i-th bead. The profile of internal distances averaged across all chains in the melt will fluctuate around a steady shape once the melt has equilibrated. Once equilibrated, we sample the configuration of the melt every 5 × 10⁵ timesteps over the next 5 × 10⁶ timesteps.

Workflow. For every pentaBCP system, we perform three simulation trials for each of the three system sizes using different random number seeds for the initial velocities of all beads. All quantities sampled from each simulation are averaged across the three trials, including R^MD_ee from the first system size, L_q* from the second system size, and the morphologies and chain conformations from the third system size. We note that the averaged results from three separate runs of the third system size are reported in Section III of this manuscript and the ESI.†

Lastly, we also conduct tests of stability of the MD simulation and SCFT predicted morphologies using a new computational approach – rapid analysis of polymer structure and inverse design strategy (RAPSIDY) – that we describe in a separate publication.⁶¹ Briefly, in this new approach we begin MD simulations using initial configurations where we pre-place chains directly into a desired morphologies (e.g., lamellae or double gyroid or others) to rapidly test the stability or metastability of that morphology for the specific polymer design. We describe this method in detail in the ESI.†

Analysis. The melt morphologies we access in our simulations include disordered (Dis), disordered microphase-separated (DM), lamella (L), hexagonal-packed cylinders (C₆), body-centered cubic spheres (BCC), double gyroid (DG), and Fddd (O⁷⁰). We use the A–B bead pair radial distribution function g_AB(r) to determine if a melt is in a disordered morphology or disordered with microphase-separation. As the pentaBCP chains de-mix, the height of the contact peak of g_AB(r ≅ 1d) decreases. We define a melt to be in a disordered morphology if g_AB(r ≅ 1d) = 1.76 ± 0.06 (i.e., the average value ±95% confidence interval), which is the height of the contact peak sampled from a completely mixed pentaBCP melt with χ = 0 (see ESI† Fig. S3). We define a melt to be in a disordered microphase-separated (DM) morphology if its contact peak height lies g_AB(r ≅ 1d) < 1.76 and outside of the 95% confidence interval of the above value without it forming an ordered morphology (e.g. L, DG, C₆). The DM morphology is analogous to the “disordered micelle” state studied in sphere-forming BCP melts⁶² in that it features significant microphase separation without the formation of an ordered morphology. However, some DM melts contain continuous and asymmetric (“defective”) structures rather than a micelle population. We did not evaluate if the DM morphology shares any of the unique characteristics with the disordered micelle state (e.g., “memory” of previous nonequilibrium structures upon crossing (χN)_ODT and then returning to order).

We identify ordered microphase-separated morphologies of pentaBCP melts visually using isosurface snapshots and via the emergence of the microphase peak (q* peak) in the A domain structure factor S_AA(q). We employ a procedure derived from that described by Brisard and Levitz⁶³ to efficiently calculate the structure factor while accounting for finite simulation volume. We treat each bead as a point scatterer and calculate the scattering amplitude as

where N_beads is the total number of simulated beads, ρ_k is the scattering length density of the bead k, and x_k is the position of bead k. To find the structure factor of the A domains, we take ρ_k = 1 if bead k is of type A and ρ_k = 0 if bead k is of type B. We account for finite-size effects to find the corrected scattering amplitude A′(q) by subtracting the form factor of the simulation volume:where q_i and L_i are the magnitude of the q vector and the length of the simulation volume along axis i, respectively. The corrected scattering amplitude A′(q) is calculated for 300 vectors oriented in all directions using the Fibonacci sphere algorithm.⁶⁴ Finally, the A domain structure factor profile is calculated by averaging the square of the corrected scattering amplitudes across all 300 vectors such that S_AA(q) is only a function of the scattering vector length q:The structure factor is particularly important for identifying triply periodic morphologies, including DG and O⁷⁰, as they can be difficult to visually distinguish from the DM morphology. To test whether or not a specific melt has either of these morphologies, we quantify the ratio of the simulation box length to the length corresponding to the q* peak, L_q*, which should be an integer multiple of  : 1 and  : 1 for the DG and O⁷⁰ morphologies, respectively.^65–67

To analyze the conformations of simulated chains, we identify the A and B domains of the simulated melts using isosurfaces drawn with a Gaussian density mesh method⁶⁸ implemented in the OVITO software.⁶⁹ We identify A domains by generating isosurfaces around beads of type A with a resolution of 30, radius scaling of 100%, and iso value of 10; this combination of parameters correctly yields the fraction of the simulation volume occupied by A domains to be equal to f_A. We then quantify the fraction of each type of block (e.g., the end A blocks, B blocks, and middle A block for the A₁B₁A₂B₂A₃ pentaBCP) in A domains. For comparison, we also quantify the distribution of each block type in each domain using a particle mesh Ewald (PME) method developed by Essmann et al.,⁷⁰ and we present the results of both methods in the ESI.† This PME calculation is used to distribute bead charges onto a voxel grid for the efficient calculation of electrostatic interactions; as in previous simulation work on BCPs,⁷¹ we use this method to distribute the mass of each bead onto a voxel grid and then assign each voxel to the domain type (A or B) corresponding to the majority of the bead mass in that voxel. This allows us to confirm that our conclusions about the chain conformations are consistent and independent of the type of calculation. We also use these bead mass grids to create 1D and 2D distributions of each block type within A and B domains for comparison with SCFT data.

We further distinguish chain conformations by using the chain end-to-end distance R_ee to estimate the number of domains accessed by each chain. We analyze the R_ee distribution of simulated melts to determine the proportion of each conformation.

F. Simulated pentaBCP designs

The advantage of combining SCFT and MD simulations is that SCFT accelerates the design parameter sweep as compared to MD simulations. Thus, running MD simulations at every point in the design space defeats the purpose of this accelerated exploration with SCFT. Instead, using the results from SCFT (discussed in Section III.A), we select a few pentaBCP designs to simulate and analyze their chain conformations. The design parameters for the twenty simulated pentaBCP designs are listed in Table 2. Each simulated system is comprised of identical chains of specified N, fraction of A-type beads f_A, and the fraction τ_i represented by the chain's middle block of the total beads of its type. The composition of pentaBCPs of sequence A₁B₁A₂B₂A₃ and B₁A₁B₂A₂B₃ are specified by the parameters τ_A2 and τ_B2, respectively.

Table 2 Design parameters for the pentaBCPs used in CG MD simulations. N is the degree of polymerization of the simulated chains (i.e. the number of beads per chain), f_A is the fraction of A-type beads in each chain, and τ_A2 and τ_B2 are the fraction of the middle block of beads of A or B type, respectively

Block sequence	N	f A	τ A2 OR τB2
A1B1A2B2A3	51
A1B1A2B2A3	50	0.4
B1A1B2A2B3	48
B1A1B2A2B3	51

---

III. Results

A. SCFT phase diagrams of ABABA and BABAB pentaBCPs

In Fig. 2 we show the SCFT phase diagrams of the A₁B₁A₂B₂A₃ pentaBCPs at different χN values with respect to the design parameters, τ_A2 and f_A. We observe disordered phases (Dis) only in the phase diagrams calculated at χN = 35 and 40. Between χN = 40 and 60, all design parameters forming the disordered phases at χN = 40 undergo a transition from disorder to order and the phase diagram at χN = 60 only shows ordered morphologies.

Fig. 2 The SCFT phase diagrams with respect to the polymer design parameters τ_A2 and f_A (as defined in Fig. 1) for the A₁B₁A₂B₂A₃ pentaBCPs at (a) χN = 35, (b) χN = 40, and (c) χN = 60. The legend on top of the figure shows which symbols represent the different equilibrium phases identified by SCFT free energy calculations: homogeneous disorder (Dis), body-centered cubic spheres (BCC), hexagonal-packed cylinders (C₆), double gyroid (DG), Fddd (O⁷⁰), and lamellae (L) phases. The conventional presentation of SCFT phase diagrams with the continuous phase boundaries are presented in ESI† Fig. S4.

The values of τ_A2 = 0 and τ_A2 = 1 correspond to the symmetric ABA and BAB triBCPs, respectively. At these extreme values of τ_A2, for all three values of χN, as f_A increases, we observe this sequence of ordered phases: C₆ → DG → L. At intermediate values of τ_A2 and χN values of 35 and 40 (Fig. 2a and b), we observe the transition from Dis → BCC → C₆ → DG or O⁷⁰ → L. This order of sequences is similar to that observed with the AB diBCP phase diagram because as A/B block composition becomes more symmetric, the polymers assemble into morphologies with less curved A–B interfaces. While this order of sequences C₆ → DG → L with increasing f_A remains consistent at all values of τ_A2, the value of f_A where the OOT occurs depends on the value of τ_A2. At χN = 60 (Fig. 2c) we observe no disordered phase, and all the chain designs exhibit ordered phases (BCC, C₆, DG, and L).

The phase behavior A₁B₁A₂B₂A₃ pentaBCPs at the three values of χN varies the most between 0.3 ≤ τ_A2 ≤ 0.5. At χN = 35 (Fig. 2a), a wide window of disordered phase is observed around 0.3 ≤ τ_A2 ≤ 0.5. The ordered phases next to the disordered region of the phase diagram are weakly segregated with a broader A–B interface compared to the ordered phases in τ_A2 = 0 and τ_A2 = 1 at the same f_A (ESI† Fig. S6). In this weakly segregated regime, the DG phase stability is outcompeted by the O⁷⁰ phase stability, which is a characteristic feature of AB diBCP phase behavior near ODTs. At χN of 60 (Fig. 2c) and 0.3 ≤ τ_A2 ≤ 0.5, all phases are ordered over the f_A range we present, and the L stability region extends to low f_A values (f_A ≈ 0.38) with well segregated A and B domains (ESI† Fig. S7). In contrast, the width and position of the phase stability windows for the τ_A2 ≈ 0 and τ_A2 ≈ 1 are insensitive to χN, indicating that (χN)_ODT from SCFT calculations is well below 35. This suggests that among A₁B₁A₂B₂A₃ pentaBCP designs with the same amount A monomers, the pentaBCP designs with relatively evenly distributed A monomers among the three A blocks (0.3 ≤ τ_A2 ≤ 0.5) have lower ODT temperatures and higher (χN)_ODT than the designs where A monomers are predominantly distributed to either the end or the middle blocks (τ_A2 ≈ 0 and τ_A2 ≈ 1). Considering that τ_A2 = 0 and τ_A2 = 1 correspond to the ABA and BAB triBCPs, respectively, this finding indicates that the A₁B₁A₂B₂A₃ pentaBCPs have higher (χN)_ODT (or lower ODT temperature) than the triBCPs. The higher (χN)_ODT in A₁B₁A₂B₂A₃ pentaBCPs is likely because there are more available configurations in the pentaBCPs than the triBCPs when A and B blocks are microphase-separated by looping or bridging the intermediate blocks, which results in higher degree of A–B mixing for the pentaBCPs than the triBCPs at the same χN.

At all three values of χN for a fixed f_A value as we increase τ_A2, we observe a reentrant phase sequence; e.g., at f_A = 0.45, going from τ_A2 = 0 to 1 at χN = 60 (Fig. 2c), the equilibrium phase goes from L → DG → L → DG → L (see ESI† Fig. S8 for the f_A = 0.45 line overlaid on Fig. 2c; along this line the reentrant phase sequence is observed). Similar reentrant phase behavior has been observed in ABA triBCPs,^23,25 and ABAB tetraBCPs,^27,28,30 but not in simple AB diBCPs. In the case of ABA triBCPs,²³ the SCFT phase diagram of the ABA triBCPs is plotted with respect to τ and f_A, where the asymmetry parameter τ is the volume fraction of one A block to all A blocks; τ = 0 corresponds to the AB diBCP and τ = 0.5 corresponds to a symmetric ABA triBCP. For the ABA triBCPs, the reentrant phase transition (L → DG → C₆ → DG → L) and phase boundary shape with decreasing τ from 0.5 to 0 is similar to the shape of the phase boundaries and the ordered phase transitions observed in the half portion of our A₁B₁A₂B₂A₃ pentaBCP phase diagram (0 ≤ τ_A2 ≤ 0.5) in Fig. 2. For the ABA triBCP, the reentrant phase transition was explained by the tendency of asymmetric ABA triBCP to form curved A–B interfaces found in DG and C₆ while the AB diBCP (τ = 0) and the symmetric ABA triBCP (τ = 0.5) prefer the zero A–B interface curvature found in the L phase. The authors of that SCFT study²³ suggested that at the intermediate values of τ between 0 and 0.5, the A blocks of different lengths spatially organize both near and far from the A–B interfaces, increasing the tendency to curve the A–B interface toward the inner parts of the A domains to reduce the overall stretching penalty. This discussion of the triBCP is relevant because one could consider the A₁B₁A₂B₂A₃ pentaBCP design as two identical triBCPs connected in series. We assume that there is higher tendency to form curved geometry phases for the pentaBCP designs (0 < τ_A2 < 0.5 and 0.5 < τ_A2 < 1) equivalent to the two asymmetric triBCPs connected, caused by spatial organization of A blocks of different lengths, than for the symmetric ABA triBCP (τ_A2 = 0), symmetric BAB triBCP (τ_A2 = 1), and A₁B₁A₂B₂A₃ pentaBCP of τ_A2 = 0.5 which is equivalent to the two symmetric ABA triBCP connected in serial. For the case of A₁B₁A₂B₂A₃ pentaBCP of τ_A2 = 0.5, the middle A₂ block has twice the length of the end A blocks, but the looped or bridged A₂ block conformations make the stretched chain length similar between the middle A₂ block and the two end A blocks, making it less efficient to spatially organize for forming the curved A domains. The nonmonotonic trend in interfacial curvature due to different spatial organizations of blocks at different τ_A2 values likely cause the reentrant phase behavior. We discuss this effect of τ_A2 on the adopted chain conformations and spatial organization of the A and B blocks in Section III.B.

Fig. 3 shows the B₁A₁B₂A₂B₃ pentaBCP phase diagrams at three different χN. Similar to the A₁B₁A₂B₂A₃ phase diagrams in Fig. 2, B₁A₁B₂A₂B₃ pentaBCPs also exhibit reentrant phase behavior at fixed f_A values, and with decreasing χN, the ordered phases are destabilized exhibiting disordered phases at the intermediate τ_B2 values. In fact, if the identities of A and B blocks are interchanged, the B₁A₁B₂A₂B₃ pentaBCP phase diagrams correspond to those of the A₁B₁A₂B₂A₃ pentaBCP phase diagrams over the range 0.55 < f_A < 0.85. The combined phase diagrams of Fig. 2 and 3 are provided in ESI† Fig. S9. One key difference between Fig. 2 and 3 is that the ordered phase stability windows for the curved structures (BCC, C₆, DG, and O⁷⁰) are mostly located in f_A < 0.40 for the B₁A₁B₂A₂B₃ pentaBCPs while a larger stability window of those phases is observed in f_A > 0.40 for the A₁B₁A₂B₂A₃ pentaBCPs (Fig. 2). For example, the DG stability windows of B₁A₁B₂A₂B₃ pentaBCPs at χN = 60 (Fig. 3c) are located in f_A ≤ 0.38 while a large portion of DG stability windows of A₁B₁A₂B₂A₃ pentaBCPs at the same χN = 60 (Fig. 2c) is located over f_A = 0.38, especially with DG windows at the highest f_A for τ_A2 = 0.2 and τ_A2 = 0.8.

Fig. 3 Same figure caption as Fig. 2 but for B₁A₁B₂A₂B₃ pentaBCPs at (a) χN = 35, (b) χN = 40, and (c) χN = 60.

The differences in the position of the phase stability windows between A₁B₁A₂B₂A₃ and B₁A₁B₂A₂B₃ pentaBCPs implies that A₁B₁A₂B₂A₃ pentaBCPs have a higher tendency to form curved A–B interfaces than B₁A₁B₂A₂B₃ pentaBCPs at the same f_A and χN. For the A₁B₁A₂B₂A₃ pentaBCPs, the A domains are formed by the two free end blocks (A₁ and A₃ blocks) and the one looped or bridged middle block (A₂ block), while for the B₁A₁B₂A₂B₃ pentaBCPs the A domains are formed by looped or bridged A₁ and A₂ blocks of equal length. Local segregation of same component blocks within same domains has been well known as an important mechanism for accommodating interfacial curvatures of complex spherical phases.^28,72,73 We use MD simulations to calculate the spatial distributions of the A and B blocks in both A and B domains and present those results in Section III.C.

B. Density profiles of A and B blocks for ABABA pentaBCPs

Density profiles and hypothesized chain conformations. To understand the polymer chain conformations that could explain the reentrant phase behavior (Fig. 2c), we plot in Fig. 4a–e the distribution of block densities along one direction within the equilibrium phases of the A₁B₁A₂B₂A₃ pentaBCPs at χN = 60 with varying τ_A2 at fixed f_A = 0.4 (see ESI† Fig. S8 for the f_A = 0.4 line overlaid on Fig. 2c along which the reentrant phase sequence is observed). The density distributions at τ_A2 = 0.1 (Fig. 4a) show that the A domains in the DG phase are formed by the end A blocks (A₁ and A₃ blocks) and the shorter middle A₂ blocks are present in the B domains. Although the presence of A blocks in the B domains, which we call “solubilized A blocks” throughout this paper, results in unfavorable contacts of A and B monomers; the enthalpic energy penalty from the A and B mixing is likely compensated for by the configurational entropy gain from the relaxed B blocks (see hypothesized chain conformations in Fig. 4f). Contrarily, the conformations in which the short middle A₂ blocks reside in the A domains along with the A end blocks would cost higher configurational entropy of B blocks because the B blocks have to be looped or bridged (as in the following hypothesized majority chain conformations in Fig. 4g, i, and j).

Fig. 4 (a–e) The spatial SCFT density distributions for all B blocks (ϕ_B = ϕ_B1 + ϕ_B2, blue lines), the middle A blocks (ϕ_A2, black short-dashed lines), and the end A blocks (ϕ_A1 + ϕ_A3, red long-dashed lines) for the equilibrium phases at f_A = 0.4 and (a) τ_A2 = 0.1 (DG), (b) τ_A2 = 0.2 (C₆), (c) τ_A2 = 0.5 (L), (d) τ_A2 = 0.8 (C₆) and (e) τ_A2 = 0.9 (DG), identified on the A₁B₁A₂B₂A₃ pentaBCP phase diagram at χN = 60 (Fig. 2c). The A domain (ϕ_A > 0.5) and B domain (ϕ_B > 0.5) are represented with the shaded red and blue regions, respectively. The SCFT equilibrium morphologies (DG → C₆ → L → C₆ → DG) are presented next to the density profiles. The red lines on the DG and C₆ morphologies depict the directions along which the densities are plotted. For the DG phases, the densities are plotted along the [111] direction of the unit cells, which passes though the four three-fold connectors of DG network domains. For the C₆ phases, the densities are plotted along the center-to-center direction of the cylinders. (f–l) Our hypothesized majority and minority chain conformations forming the A/B domain interfaces, corresponding to from τ_A2 = 0.1 to 0.9. The A/B domain interfaces in the DG phases are depicted as flat to exaggerate the interfacial curvature change in the reentrant phase behavior.

The density distributions at τ_A2 = 0.2 (Fig. 4b) show the spatial segregation of A blocks within the A domains of the C₆ morphology. The long A end blocks form the core of the cylinder domains and the short middle A₂ blocks form the shell of the cylinder domains with their highest concentration being near the A–B interface. Some amount of middle A₂ block density (ϕ_A2 = 0.048) remains at the center of the B domains; this suggests that while majority of the chains prefer A and B blocks to be segregated in their respective domains by looping or bridging the B blocks (see hypothesized chain conformations in Fig. 4g), a small fraction of chains adopt configurations where middle A₂ blocks is solubilized within the B domains (see hypothesized chain conformations in Fig. 4h).

The density distributions at τ_A2 = 0.5 (Fig. 4c) show all end and middle A blocks exhibit similar spatial density distributions in the A lamellar domains; this is because the volume of two end A blocks is equal to that of the one middle A₂ block. We assume that at τ_A2 = 0.5, the degree of stretching in the end A blocks and middle A block are relatively comparable than the case where the length ratio between the two types of A blocks is highly asymmetric (e.g. τ_A2 = 0.1 or 0.2), as indicated in the similar spatial density distributions within the A domains. The spatial segregation of A blocks within the A domains of the C₆ morphology becomes stable again at τ_A2 = 0.8, where the long middle A₂ block and short end A blocks are preferentially localized into core and shell regions within the A cylindrical domains (Fig. 4d). We hypothesize that in a majority of cases, the A and B blocks segregate in the different domains (as shown in our schematic in Fig. 4j), and for a small fraction of chain conformations the short end A blocks are solubilized in the B domains (schematic in Fig. 4k). We also hypothesize that at extremely strong segregation regime (χN → ∞), the minority chain conformations (Fig. 4h and k) would disappear for both τ_A2 = 0.2 and τ_A2 = 0.8 due to significant enthalpic A–B contact penalty.

The density distributions at τ_A2 = 0.9 (Fig. 4e) show that the A domains in the DG phase are mostly formed by the middle A₂ blocks and the short end A blocks are solubilized in the B domains. As in the τ_A2 = 0.1 case, the enthalpic loss from the unfavorable A–B mixing is compensated by the configurational entropy gain from the relaxed B blocks (see our hypothesized conformation in Fig. 4l).

A–B interfacial curvature promoted by chain conformations. In the reentrant phase sequence presented in Fig. 4, the C₆ morphology have the highest interfacial curvature and becomes stable when the spatial segregation of A blocks within the A domains occurs. Both majority chain conformations (Fig. 4g and j) inducing the spatial segregations favor the curved A–B interfaces as indicated by higher interfacial curvatures in the C₆ phases (τ_A2 = 0.2 and 0.8) than in the DG phases (τ_A2 = 0.1 and 0.9) and in the L phase (τ_A2 = 0.5). It has been shown in prior literature^23,28,74 that spatial segregations of blocks of different length within the same domain reduces the overall stretching energy of the blocks and relieves the packing frustration in the domain, and in turn, induces a tendency towards higher A–B interfacial curvature than in domains formed by blocks with monodisperse length. We hypothesize that in the C₆ morphology the majority of the polymer chains adopt the chain conformations with the spatially segregated A blocks within the A domains (as in Fig. 4g and j) than with the solubilized A blocks (Fig. 4h and k) in order to minimize the A–B contact energetic penalty. The configurational entropy loss from looping or bridging B blocks is minimized by curving the A–B interface and favoring the C₆ phase with core–shell distribution of A blocks. The L phase becomes stable at τ_A2 = 0.5 (Fig. 4c) when the density distributions of the A blocks exhibit similar spatial density distributions. The spatial segregation of A blocks is suppressed by the comparable degree of stretching in the end A blocks and middle A block (shown schematically in Fig. 4i) and such similar spatial distributions of all three A blocks is likely to stabilize the L phase over the C₆ phase.

Motivation for MD simulation analysis. In the hypothesized majority chain conformation schematics in Fig. 4g, i, and j, we skip the possibility of the chains adopting a bridging configuration in more than one B block. When we analyze the density distributions in Fig. 4a–e, we assume that the population of chains adopting such a “bridging” conformation is small (negligible). In SCFT, it is possible to calculate bridging/looping fractions of block by solving partial partition functions with modified initial conditions which are defined by Voronoi unit cells.^20,75,76 However, the calculation is not straightforward for determining the fractions of multiple different chain conformations if the conformations are defined with more than one block adopting bridged or looped configurations, because it requires multiple evaluations of partial partition functions with a series of dependent initial conditions. Moreover, as the dimensionality of morphology increases, the number of different possible conformations defined by different Voronoi unit cells increases rapidly, and even determining the Voronoi unit cell is not trivial for the complex morphologies such as DG and O⁷⁰ with SCFT density solutions. The complexity in calculation of pentaBCP chain conformations by SCFT is discussed in ESI† Fig. S10–S13, where we presented the fractions of A₁B₁A₂B₂A₃ pentaBCP conformations calculated by SCFT in the L morphology. Even in the simplest case of one-dimensional L morphology, there are six possible chain conformations (ESI† Fig. S10), and for the C₆ morphology, evaluating the fraction of each chain conformation requires enumerating all possible initial conditions belonging to different Voronoi unit cells (ESI† Fig. S13). On the other hand, MD simulations is a straightforward approach for determining the population of different chain conformations, which will be presented in the following Section III.D. Therefore, to quantify the chain conformations and spatial distribution of the A and B blocks, we use CG MD simulations. Specifically, we use CG MD simulations to test our hypothesized conformations that explain the observed reentrant phase transition (DG → C₆ → L → C₆ → DG) for A₁B₁A₂B₂A₃ pentaBCPs with increasing τ_A2 from 0.1 to 0.9.

C. Chain conformations adopted by pentaBCPs in various morphologies

Ordered morphologies in simulation. Before analyzing chain conformations in MD simulation, we compare the morphologies predicted by SCFT calculations and MD simulations in ESI† Fig. S14. We present the morphologies from SCFT for each segregation strength studied and from MD starting at χN = 90, as this is the lowest χN at which an ordered morphology appeared in simulation. As seen in Fig. 2 and 3, the SCFT calculations predict that melts of each chain design will be ordered at χN = 60, whereas MD simulations suggest that each will remain in the disordered microphase separated (DM) morphology at this χN; this discrepancy in the disorder to order transition from SCFT and MD simulation has been seen in multiple previous studies.^77–81 We observe that each simulated melt adopts a DM morphology at 40 ≤ χN ≤ 80. Overall, we find good agreement in the morphologies predicted by SCFT and MD at the highest χN examined with either method (χN = 60 for SCFT and 120 for MD) for equivalent chain designs. In cases where there is disagreement and where a DM morphology was observed with MD at χN ≤ 90, we use RAPSIDY testing as described in the ESI† to verify the stability of a specific morphology for a particular chain design (see ESI† Fig. S15 and S16). Next, we use the chain-level detail available in MD simulations to estimate the proportions of the majority/minority chain conformations in Fig. 4 for select chain designs.

PentaBCP block distributions from MD simulation and SCFT. We plot the fractions of the end blocks (A₁ and A₃) and middle A₂ blocks and the B blocks (B₁ and B₂) in A domains of the A₁B₁A₂B₂A₃ pentaBCPs with f_A = 0.4 and several τ_A2 from MD simulation in Fig. 5. We provide the fraction of each type of block in A domains for every simulated pentaBCP design in ESI† Fig. S17. Additionally, to demonstrate that our results are consistent and independent of our chosen analysis method, we compare in ESI† Fig. S18 the A domain fractions from the isosurface and PME methods (described in Section II.E).

Fig. 5 The fractions of the end A blocks A₁ and A₃ (open red circle), the B blocks B₁ and B₂ (filled blue squares), and the middle A₂ block (filled black circles) within A domains, for A₁B₁A₂B₂A₃ pentaBCP copolymers with f_A = 0.4 at several τ_A2 at χN = 120 from MD simulation. The plotted fraction and error bars represent the average and standard deviation of the computed fractions from each sampled frame from all three trials using third system size (as described in Section II.E).

We find that smaller blocks are more likely to solubilize, leading to the lowest average A block solubilization when the end and middle A block sizes are balanced at τ_A2 = 0.5, as predicted in Fig. 4i. The end A blocks (A₁ and A₃) are increasingly solubilized in B domains as τ_A2 increases and, therefore, the end block size relative to the middle A₂ block shrinks. In contrast, middle A₂ blocks are less likely to reside in B domains as τ_A2 increases. The fraction of B blocks (B₁ and B₂) found in A domains is greatest when τ_A2 ∼ 0.5 and decreases at extreme values of τ_A2. We speculate that this is because the greater solubilization of the small A blocks at either extreme (i.e., τ_A2 = 0.2 and τ_A2 = 0.8) permits the B blocks to reside further away from the A–B interface on average and therefore decreases their solubilization in A domains. Conversely, the B blocks become more solubilized in A domains on average when the end and middle A blocks are both sufficiently large to remain in the A domains (at τ_A2 ∼ 0.5).

Using the preferred domain location of each block type, we can now evaluate the hypothesized proportions of each A₁B₁A₂B₂A₃ pentaBCP chain conformation in Fig. 4. For all τ_A2, the B blocks of the A₁B₁A₂B₂A₃ pentaBCPs mostly remain in B domains (see blue symbols in Fig. 5) which supports the hypothesized B block arrangement in chain conformations in Fig. 4. For all τ_A2 besides 0.2 and 0.8, more than 90% of the end and middle A blocks in A domains, which is also in accordance with the hypothesized chain conformations for τ_A2 = 0.5 (Fig. 4i). For τ_A2 of 0.2, the middle A₂ block for ∼30% of the chains (i.e., minority) resides in the B domain; this also supports the hypothesized “majority” and “minority” conformations in Fig. 4g and h. For τ_A2 of 0.8, the majority of end A blocks reside in A domains but chains have at least one “solubilized” end block; as each pentaBCP has two end blocks, it has a ∼ (70%)² = 49% chance of having at least one end A block in an A domain. This means that the most frequent conformation for chains with τ_A2 of 0.8 would be a combination of the hypothesized conformations in Fig. 4j and k, with one end A block in a B domain.

To further compare the distribution of each block type between SCFT and MD, we present one- and two-dimensional volume fraction distributions of the end A blocks, B blocks, and middle A₂ blocks of the A₁B₁A₂B₂A₃ pentaBCP with f_A = 0.4 and τ_A2 = 0.8 in Fig. 6; this polymer design assembles into the C₆ morphology at χN = 120 (Fig. 6a). We calculate one-dimensional distributions of each block, analogous to those derived from the SCFT calculations in Fig. 4, by sampling the volume fractions in a plane parallel to (one of) the periodic axis of the morphology; this periodic axis is indicated with a dashed line in Fig. 6a. We present the distributions derived from the SCFT calculations for the same chain designs, originally given in Fig. 4d, in Fig. 6b for comparison. We present the same data in ESI† Fig. S19 for a simulated melt of the A₁B₁A₂B₂A₃ pentaBCP with f_A = 0.4 and τ_A2 = 0.5, which forms the L phase.

Fig. 6 (a) Snapshot of 3D isosurface, (b and c) one-dimensional volume fraction profiles of end A blocks, B blocks, and middle A blocks from (b) SCFT and (c) MD, and (d–f) two-dimensional volume fraction distributions of (d) end A blocks, (e) B blocks, and (f) middle A blocks from MD for the A₁B₁A₂B₂A₃ pentaBCP melt with f_A = 0.4 and τ_A2 = 0.8 at χN = 120. The one-dimensional distributions are sampled within a center-to-center interval (solid line segment with end marks in Fig. 6a) of a diagonal plane parallel to the cylinder axes (dashed line) and spatially averaged across the plane in the direction of the solid arrow. They are plotted against the distance r from the origin along the dashed arrow. The two-dimensional distributions are averaged along the solid arrow; note the different scalebar ranges. All volume fraction distributions are averaged between all sampled frames across three trials.

We observe that the one-dimensional MD volume fraction distributions exhibit many of the same characteristics as the SCFT-derived distributions, including those which we took as evidence of the chain conformations hypothesized in Fig. 4. For the SCFT and MD distributions for the cylinder-forming melt (Fig. 6b and c respectively), the volume fraction of the end A blocks (ϕ_A1 + ϕ_A3) is highest near the A–B interface (ϕ_A1 + ϕ_A3 ∼ 0.16), declines towards the center of the B domains, and becomes nearly zero at the middle of the A domains. This suggests that the end A blocks segregate towards the outside of the A domains and can also be present in B domains in simulated melts of this chain design, as predicted by SCFT calculations and reflected in the chain conformations in Fig. 4j and k. The volume fraction of B blocks ϕ_B reaches a maximum in the middle of the B domains and a minimum in the A domains. Conversely, ϕ_A2 is highest at the middle of each A domain and lowest, but not zero, at the middle of the B domains. We note that ϕ_A2 ∼0.05 in the MD data and ∼0.01 in the SCFT data in the B domains.

The two-dimensional spatial distributions of the blocks provide further evidence that the spatial segregation of A blocks in A domains in simulation agrees with SCFT predictions. The two-dimensional distributions of end A blocks (ϕ_A1 + ϕ_A3) in Fig. 6d show that end A blocks (A₁ and A₃) preferentially segregate towards the outside of each cylindrical A domain near the A–B interface and have a sparse presence towards the center of each A domain and in the B domain between them. Fig. 6e and f show that the B blocks (B₁ and B₂) and middle A blocks (A₂) reside primarily in the B domains and the interior of A domains, respectively.

We find that the one-dimensional block volume fraction distributions for the L-forming melt in ESI† Fig. S19c are also consistent with SCFT results. The volume fractions of the end A blocks (ϕ_A1 + ϕ_A3) and middle A blocks (ϕ_A2) both reach a maximum of ∼0.45 at the middle of A domains and decline to nearly zero towards the center of B domains. These volume fraction distributions suggest that both end and middle A blocks distribute relatively evenly throughout A domains as reflected in the hypothesized conformations in Fig. 4i and do not exhibit the segregated distribution as observed in Fig. 6. Conversely, the volume fraction of B blocks (ϕ_B) is highest at the middle of B domains and declines to ∼0.08 at the middle of A domains. These distributions demonstrate that each block type mostly resides in domains of the same type in simulation, as reflected in Fig. 5, which also supports the conformations in Fig. 4i. The two-dimensional block distributions for the end and middle A blocks (ESI† Fig. S19d and f) further suggest that both block types distribute evenly throughout A domains.

The correspondence between the 1D block distributions from SCFT and MD and the local segregation present in the 2D distributions from MD both support the hypothesized chain conformations in Fig. 4. However, these data are insufficient to determine which A domains each A block resides in; the majority and minority chain conformations in Fig. 4 are further distinguished by whether or not the end blocks of a pentaBCP are in the same domain.

D. Determining the population of each chain conformation in MD simulations

We use the chain end-to-end distance R_ee to understand if both end blocks are located in the same domain or not. We define “looping” chains as chains where both second and fourth “connecting” blocks (e.g., B₁ and B₂ for A₁B₁A₂B₂A₃ and A₁ and A₂ for B₁A₁B₂A₂B₃) adopt looping configurations in domains of the same type, such that the end and middle A blocks all reside in A domains. Chains with the looping conformation often have a “bundled” appearance and, therefore, low R_ee. We define “bridging” chains as chains where both connecting blocks (e.g., B₁ and B₂ for A₁B₁A₂B₂A₃) adopt bridging configurations to stretch across distinct domains and therefore allow each of the chain's end and middle blocks to occupy separate domains as well (i.e., all five blocks located in their own domains), leading to relatively high R_ee. In contrast, the connecting blocks of chains with the “hairpin” conformation (described in the ESI†) occupy the same domain, leading to a much lower R_ee; we found this conformation was rare for all chain designs and we therefore do not systematically distinguish it. PentaBCPs may also have a connecting block in either configuration, which leads to a “hybrid” conformation with one looping and one bridging block and an intermediate chain R_ee.

In Fig. 7a, we present a snapshot of a L morphology formed by a melt of A₁B₁A₂B₂A₃ pentaBCP chains with f_A = 0.4 and τ_A2 = 0.5 with several example chains adopting looping, hybrid, and bridging conformations within the A and B domains of the L phase. We include the distribution of end-to-end distance R_ee in Fig. 7b for the same pentaBCP design (also examined in Fig. 7e) to show the low, medium, and high R_ee modes, separated at R_ee = 7d and 16d. The proportion of chains with R_ee < 7d, 7d ≤ R_ee < 16d, and 16d ≤ R_ee at several τ_A2 are shown in Fig. 7c–g, along with snapshots of representative chains and their R_ee in parenthesis from each population.

Fig. 7 (a) Representative snapshot of A₁B₁A₂B₂A₃ chains adopting looping, hybrid, and bridging conformations in a lamellar melt. (b) Example probability density histogram of chain end-to-end distance R_ee for the A₁B₁A₂B₂A₃ pentaBCP at χN = 120 with f_A = 0.4 and τ_A2 = 0.5, which is also represented in Fig. 7e. Black vertical lines mark the boundaries at R_ee = 7d and 16d between the low, medium, and high R_ee chain populations, and the gold vertical line marks the real-space distance L_q* sampled from S_AA(q). The R_ee distributions for each design are presented in ESI† Fig. S20. (c–g) The percentages of chains within a specific range of R_ee is denoted as [R_ee] for R_ee < 7d, 7d ≤ R_ee < 16d, and 16d ≤ R_ee for A₁B₁A₂B₂A₃ pentaBCP at χN = 120 with f_A = 0.4 and (c) τ_A2 = 0.2, (d) τ_A2 = 0.4, (e) τ_A2 = 0.5, (f) τ_A2 = 0.6, and (g) τ_A2 = 0.8. Representative chains are presented along with their R_ee for each population; the end A blocks of the chains (light red) are highlighted against the middle A blocks (dark red). The morphologies observed with MD simulation for these values of τ_A2 (0.2 to 0.8) are DG, L, L, L, and C₆, respectively (ESI† Fig. S14).

We use a gold line to mark the characteristic domain size of each melt on the distributions in Fig. 7b and S20,† which we identify as the real-space length corresponding to the major peak in the melt structure factor (L_q*). We observe that the low, medium, and high R_ee populations have peaks at ∼0.5, 1, and 2 times L_q*, which corresponds to the expected relationship between L_q* and R_ee for looping, hybrid, and bridging chains. Looping chains have small R_ee as their end blocks occupy the same domain, although it is possible for their R_ee ∼ L_q*. The ends of hybrid and bridging chains are separated by one and two morphology periods, respectively, which suggests that chains of these conformations form the medium and high R_ee populations. Inspection of the example chains also suggests that these low, medium, and high R_ee populations mostly consist of chains with the looping, hybrid, and bridging conformations, respectively.

The proportion of each population is broadly consistent for each given chain design: high R_ee chains are always the smallest fraction, ranging from ∼3–10% of all chains; small R_ee chains are the next largest population at ∼27–30%; and medium R_ee chains constitute ∼63–68% of all chains. The small proportion of high R_ee chains for each chain design support the assumption that the conformations with two B blocks bridging are rare and can be neglected in the hypothesized conformations in Fig. 4f–l. Further, the relative proportions of low and medium R_ee chains suggest that the hybrid conformation, the right-side conformation in Fig. 4g, i, and j, is most prevalent in melts of each chain design.

Although the size of each population is insensitive to τ_A2, we observe some trends that relate the domain fraction results presented in Fig. 5 to the prevalence of looping, hybrid, and bridging conformations. The bridging conformation, which has the highest average R_ee, relies on both end blocks remaining in separate domains and the middle block bridging to opposite sides of its domain. Accordingly, this conformation is suppressed in the melt of pentaBCP chains with τ_A2 = 0.2, for which the middle A block is the smallest among the simulated chain designs and therefore has a high stretching energy penalty in the bridging configuration; the high R_ee peak for the corresponding distribution in ESI† Fig. S20a is subtle. The distributions for the L-forming melts with τ_A2 = 0.4, 0.5, and 0.6 (Fig. S20b–d†) appear relatively similar; for all three systems, the low R_ee population decreases and the high R_ee population increases as τ_A2 increases, in tandem with the increase and decrease in solubilization of the end and middle A blocks, respectively. These data demonstrate that the bridging conformation is promoted as the size of the middle block increases and the looping conformation is suppressed as the size of the end blocks decreases, even within the same morphology. The proportions of each R_ee population do not continue these trends monotonically for the pentaBCP melt with τ_A2 = 0.8, as the low and high R_ee populations are instead larger and smaller, respectively, as compared to the τ_A2 = 0.6 melt.

IV. Discussions

A. Design rules for stabilizing DG and observing L morphologies

If one is targeting transport-favorable morphologies such as DG or L, the phase diagrams on Fig. 2 and 3 give valuable insights to experimentalists for the choice of polymer designs that gives high propensity for forming DG and L phases. Here, we discuss the design parameters that give the widest DG window and the stable L phase formation at the low f_A based on the phase diagrams (Fig. 2 and 3).

The maximum DG window width (Δf_A ≈ 0.06) is found around τ_A2 = 0.8 for the A₁B₁A₂B₂A₃ pentaBCPs for all three χN, which is slightly wider than that of diblock copolymers (Δf_A = 0.03–0.04). The higher stability of DG in the A₁B₁A₂B₂A₃ pentaBCPs can be explained as follows. In AB diBCPs, only one block per polymer chain forms the gyroid structure (e.g., A block if A is the minority composition block) while in the A₁B₁A₂B₂A₃ pentaBCP, the three A blocks form the gyroid structure. In the case of τ_A2 = 0.8, where the DG window width is the widest, the middle A block is eight times longer than the end A block. It has been shown in prior literature that dispersity in block lengths relieves the packing frustration in the complex network phases.^12,82 Packing frustration is high in the three-fold connectors of DG which have high variations in domain thickness and curvature compared to the classical L, C₆, and BCC phases. We speculate that the differences in the block length of the end and middle A blocks helps to efficiently organize the A blocks within A domains to relieve packing frustration and widen the DG stability window.

Moreover, compared to the AB diBCP design, the A₁B₁A₂B₂A₃ pentaBCP design introduces one additional free end of A blocks at the expense of one free end loss of B blocks. Due to the two free A ends in A₁B₁A₂B₂A₃ pentaBCP design, A blocks are more flexible in adopting various configurations than the B blocks tethered to the A blocks. The configurational entropy of A blocks helps to enhance the stability of curved phases (BCC, C₆, and DG) versus L phase because the total stretching energy of pentaBCPs is minimized by curving the A–B interface toward A domains at expense of stretching A blocks in order to relax the B blocks. The stability of curved phases over L phase is supported by our observation that in the A₁B₁A₂B₂A₃ pentaBCP phase diagram (Fig. 2), the DG windows are located around 0.45 < f_A < 0.5, which are relatively high values compared to those of the diBCPs and symmetric triBCPs (0.3 < f_A < 0.35).^20,21 We speculate that the along with the dispersity in the A blocks, the enhanced tendency to curve the A–B interface contributes to the stabilization of DG phase, and thus widening DG stability windows.

In contrast, among the polymer designs in the B₁A₁B₂A₂B₃ pentaBCPs, a comparably wide DG window width (Δf_A ≈ 0.04) is only found for χN = 60 around τ_A2 = 0.5 at low f_A values (0.34 < f_A < 0.38). The results suggest that for stabilizing DG phases in pentaBCPs, A₁B₁A₂B₂A₃ pentaBCPs around τ_A2 = 0.8 is the best design choice, but the B₁A₁B₂A₂B₃ pentaBCP with τ_B2 = 0.5 at higher χN is better if lower volume fraction DG structures are preferred.

In diBCP melt systems, the L phase has wide stability windows spanning over symmetric A/B compositional fractions (e.g., f_A ∼ 0.5). Generally, using a symmetric composition of A and B monomers is a reliable strategy for achieving L phase in AB-type BCP melts. However, in case one of the A and B type monomers have higher costs or are environmentally less favorable, it may be desired to find polymer design rules which provides stable L phases at highly asymmetric compositional fractions, where BCC or C₆ phases are usually preferred. The phase diagrams in Fig. 3 show that the L phase windows for the B₁A₁B₂A₂B₃ pentaBCPs at τ_B2 = 0.1 and 0.9 extends to the lowest f_A values among all the pentaBCP and triBCP designs we investigated. We note that the L phase is stable in B₁A₁B₂A₂B₃ pentaBCPs with highly asymmetric compositional fraction (f_A = 0.25) at τ_B2 = 0.1 and χN = 60 (Fig. 3c). Therefore, among all the pentaBCP designs we investigated, including the ABA/BAB triblocks, A₁B₁A₂B₂A₃ pentaBCPs, and B₁A₁B₂A₂B₃ pentaBCPs, the B₁A₁B₂A₂B₃ pentaBCPs at τ_B2 = 0.1 and 0.9 are the best design choice to stabilize L with the lowest content of the A monomer.

B. Underlying mechanism for the reentrant phase behavior

In Section III, we hypothesize that the spatial segregation of A blocks in A domains promotes the A–B interfacial curvature at intermediate τ_A2 values (0 < τ_A2 < 0.5 and 0.5 < τ_A2 < 1), which leads to the reentrant phase behavior. The MD simulation results for the spatial distributions of each block in Fig. 6c and S19c† confirm that the partitioning of A blocks within A domains is governed by their relative size and qualitatively consistent with SCFT density profiles in Fig. 4. It is notable that even with the different chain models (Gaussian thread model in SCFT and bead-spring model in MD simulation) and the field-based versus particle-based approach, we observe the similar reentrant phase behavior and the consistent spatial distributions of blocks within the various domains in both SCFT calculations and MD simulations.

The consistent distribution of R_ee across τ_A2 in ESI† Fig. S20 and the corresponding chain conformations in Fig. 7 suggest that R_ee is mainly determined by the configuration of B blocks (bridge, loop, or hybrid) but less correlated with τ_A2, which quantifies the relative size of the end and middle A blocks. The similar shape of probability density histograms of R_ee across τ_A2 in ESI† Fig. S20 indicates that the population of chain conformations is robust to the change of τ_A2 and morphology. Thus, it is likely that the spatial segregations of end A block and middle A blocks within A domains is mostly responsible for the reentrant phase behavior. Assuming that the configurations of B blocks are insensitive to the change of τ_A2, the overall stretching energy of polymer chains are affected by spatial organization of A blocks. At intermediate τ_A2 values (0 < τ_A2 < 0.5 and 0.5 < τ_A2 < 1), the spatial segregation of end and middle A blocks of different length within the same domains makes it easier to reduce the overall stretching energy of the blocks and relieves the packing frustration in the curved domains than the similar spatial distributions of A blocks at τ_A2 = 0.5, which we think is the driving factor for the reentrant phase behavior.

V. Conclusions

We employed a combined approach of SCFT calculations and MD simulations to study the equilibrium phase behavior of symmetric pentablock copolymer (pentaBCP) melts, with A₁B₁A₂B₂A₃ and B₁A₁B₂A₂B₃ sequence, where the first and last block have same volume fraction and second and fourth block have same volume fraction. We construct the SCFT phase diagrams for the self-assembled morphologies to identify the ideal pentaBCP designs for targeting the widest DG windows and the L phase stabilization at the lowest f_A. We also performed CG MD simulations with a set of pentaBCP designs selected from the SCFT phase diagrams and analyzed the SCFT and MD chain conformations in several morphologies. Our work leverages the chain-level detail provided by MD simulations to evaluate the underlying mechanisms of complex reentrant phase behavior observed in pentaBCP morphology predictions made with SCFT. This study establishes several design rules that may aid researchers in synthesizing pentaBCPs with transport-favorable morphologies and demonstrates the utility of complementary computational methods in the study of multiBCPs.

Data availability

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

Author contributions

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

Conflicts of interest

The authors declare that there is no known conflict of interests.

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 Caviness supercomputing cluster at University of Delaware (UD). Authors also thank Brent Sumerlin and Kaden Stevens at University of Florida for useful discussions about the limitations of various multiblock polymer synthesis techniques.

References

1. F. S. Bates and G. H. Fredrickson, Phys. Today, 1999, 52, 32–38.
2. H. Feng, X. Lu, W. Wang, N.-G. Kang and J. W. Mays, Polymers, 2017, 9, 494.
3. C. Tang, E. M. Lennon, G. H. Fredrickson, E. J. Kramer and C. J. Hawker, Science, 2008, 322, 429–432.
4. E. J. W. Crossland, M. Kamperman, M. Nedelcu, C. Ducati, U. Wiesner, D. M. Smilgies, G. E. S. Toombes, M. A. Hillmyer, S. Ludwigs, U. Steiner and H. J. Snaith, Nano Lett., 2009, 9, 2807–2812.
5. M. Stefik, S. Guldin, S. Vignolini, U. Wiesner and U. Steiner, Chem. Soc. Rev., 2015, 44, 5076–5091.
6. A. C. Chang, T. Mulia, Y. S. Wu, Y. H. Weng, Y. C. Lin and W. C. Chen, J. Polym. Sci., 2023, 61, 2633–2654.
7. L. Li, L. Schulte, L. D. Clausen, K. M. Hansen, G. E. Jonsson and S. Ndoni, ACS Nano, 2011, 5, 7754–7766.
8. H. S. Kang, S. W. Han, C. Park, S. W. Lee, H. Eoh, J. Baek, D. G. Shin, T. H. Park, J. Huh, H. Lee, D. E. Kim, D. Y. Ryu, E. L. Thomas, W. G. Koh and C. Park, Sci. Adv., 2020, 6, eabb5769.
9. M. W. Matsen and F. S. Bates, J. Polym. Sci., Part B: Polym. Phys., 1997, 35, 945–952.
10. L. Ludwik, Macromolecules, 1980, 13, 1602–1617.
11. M. W. Matsen, J. Phys.: Condens. Matter, 2002, 14, R21–R47.
12. M. W. Matsen and F. S. Bates, Macromolecules, 1996, 29, 7641–7644.
13. M. Iijima, D. Ulkoski, S. Sakuma, D. Matsukuma, N. Nishiyama, H. Otsuka and C. Scholz, Polym. Int., 2016, 65, 1132–1141.
14. N. Singh, M. v. S. Tureau and T. H. Epps, Soft Matter, 2009, 5, 4757–4762.
15. M. W. Bates, S. M. Barbon, A. E. Levi, R. M. Lewis, H. K. Beech, K. M. Vonk, C. Zhang, G. H. Fredrickson, C. J. Hawker and C. M. Bates, ACS Macro Lett., 2020, 9, 396–403.
16. S. J. Park, F. S. Bates and K. D. Dorfman, Phys. Rev. Mater., 2023, 7, 105601.
17. S. J. Park, G. K. Cheong, F. S. Bates and K. D. Dorfman, Macromolecules, 2021, 54, 9063–9070.
18. S. Mekcham and K. Nomura, J. Am. Chem. Soc., 2023, 145, 17001–17006.
19. M. B. Koo, S. W. Lee, J. M. Lee and K. T. Kim, J. Am. Chem. Soc., 2020, 142, 14028–14032.
20. M. W. Matsen and R. B. Thompson, J. Chem. Phys., 1999, 111, 7139–7146.
21. M. W. Matsen, Macromolecules, 2012, 45, 2161–2165.
22. A. M. Mayes and M. Olvera De La Cruz, J. Chem. Phys., 1989, 91, 7228–7235.
23. M. W. Matsen, J. Chem. Phys., 2000, 113, 5539–5544.
24. M. W. Hamersky, S. D. Smith, A. O. Gozen and R. J. Spontak, Phys. Rev. Lett., 2005, 95, 168306.
25. S. Sakurai, K. Shirouchi, S. Munakata, H. Kurimura, S. Suzuki, J. Watanabe, T. Oda, N. Shimizu, K. Tanida and K. Yamamoto, Macromolecules, 2017, 50, 8647–8657.
26. B. R. Chapman, M. W. Hamersky, J. M. Milhaupt, C. Kostelecky, T. P. Lodge, E. D. Von Meerwall and S. D. Smith, Macromolecules, 1998, 31, 4562–4573.
27. B. Zhao, W. Jiang, L. Chen, W. Li, F. Qiu and A. C. Shi, ACS Macro Lett., 2018, 7, 95–99.
28. B. Zhao, C. Wang, Y. Chen and M. Liu, Langmuir, 2021, 37, 5642–5650.
29. S. Li, W. Tao, K. Gao, N. Athir, F. Li, Y. Chen, J. Liu, L. Zhang and M. Tsige, RSC Adv., 2019, 9, 42029–42042.
30. S. Ahn, J. K. Kim, B. Zhao, C. Duan and W. Li, Macromolecules, 2018, 51, 4415–4421.
31. S. Ahn, Y. Seo, J. K. Kim, C. Duan, L. Zhang and W. Li, Macromolecules, 2019, 52, 9039–9044.
32. B. Zhao, Q. Dong, Y. Wei and Y. Xu, Molecules, 2023, 28, 3536.
33. Q. Xie, Y. Qiang and W. Li, ACS Macro Lett., 2022, 11, 205–209.
34. K. Yuan, Z. Xu, X. Huang and W. Li, Chem. – Eur. J., 2023, 29, e202301043.
35. Q. Dong and W. Li, Macromolecules, 2021, 54, 203–213.
36. L. Li and W. Li, Giant, 2021, 7, 100065.
37. S. J. Park, F. S. Bates and K. D. Dorfman, ACS Macro Lett., 2022, 11, 643–650.
38. M. E. Lott, L. Trachsel, E. Schué, C. L. G. Davidson, R. A. Olson S, D. I. Pedro, F. Chang, Y. Hong, W. G. Sawyer and B. S. Sumerlin, Macromolecules, 2024, 57, 4007–4015.
39. 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.
40. A. Ranjan, J. Qin and D. C. Morse, Macromolecules, 2008, 41, 942–954.
41. M. W. Matsen, Eur. Phys. J. E: Soft Matter Biol. Phys., 2009, 30, 361–369.
42. G. S. Grest, M.-D. Lacasse, K. Kremer and A. M. Gupta, J. Chem. Phys., 1996, 105, 10583–10594.
43. G. Fredrickson, The Equilibrium Theory of Inhomogeneous Polymers, 2005.
44. A. Kulshreshtha, R. C. Hayward and A. Jayaraman, Macromolecules, 2022, 55, 2675–2690.
45. J. H. Ryu, H. S. Wee and W. B. Lee, Phys. Rev. E, 2016, 94, 032501.
46. B. Steinmüller, M. Müller, K. R. Hambrecht, G. D. Smith and D. Bedrov, Macromolecules, 2011, 45, 1107–1117.
47. M. A. Horsch, Z. Zhang, C. R. Iacovella and S. C. Glotzer, J. Chem. Phys., 2004, 121, 11455–11462.
48. J. R. Brown, Y. Seo, S. W. Sides and L. M. Hall, Macromolecules, 2017, 50, 5619–5626.
49. H. Guo and K. Kremer, J. Chem. Phys., 2003, 119, 9308–9320.
50. L. A. Moreira, G. Zhang, F. Müller, T. Stuehn and K. Kremer, Macromol. Theory Simul., 2015, 24, 419–431.
51. J. E. Jones, Proc. R. Soc. London, Ser. A, 1924, 106, 463–477.
52. M. Murat, G. S. Grest and K. Kremer, Macromolecules, 1999, 32, 595–609.
53. A. Arora, D. C. Morse, F. S. Bates and K. D. Dorfman, Soft Matter, 2015, 11, 4862–4867.
54. J. Huh, C. Park and Y. K. Kwon, J. Chem. Phys., 2010, 133, 114903.
55. F. J. Martinez-Veracoechea and F. A. Escobedo, J. Chem. Phys., 2006, 125, 104907.
56. J. Škvor and Z. Posel, Macromol. Theory Simul., 2014, 24, 141–151.
57. 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.
58. A. I. Jewett, D. Stelter, J. Lambert, S. M. Saladi, O. M. Roscioni, M. Ricci, L. Autin, M. Maritan, S. M. Bashusqeh, T. Keyes, R. T. Dame, J. E. Shea, G. J. Jensen and D. S. Goodsell, J. Mol. Biol., 2021, 433, 166841.
59. S. Nosé, Mol. Phys., 2006, 52, 255–268.
60. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
61. V. Liao, T. Myers and A. Jayaraman, Soft Matter, 2024 10.1039/D4SM01037J.
62. K. Kim, A. Arora, R. M. Lewis, 3rd, M. Liu, W. Li, A. C. Shi, K. D. Dorfman and F. S. Bates, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 847–854.
63. S. Brisard and P. Levitz, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 013305.
64. Á. González, Math. Geosci., 2009, 42, 49–64.
65. D. A. Hajduk, P. E. Harper, S. M. Gruner, C. C. Honeker, G. Kim, E. L. Thomas and L. J. Fetters, Macromolecules, 2002, 27, 4063–4075.
66. T. H. Epps, E. W. Cochran, C. M. Hardy, T. S. Bailey, R. S. Waletzko and F. S. Bates, Macromolecules, 2004, 37, 7085–7088.
67. H. Lee, S. Kwon, J. Min, S. M. Jin, J. H. Hwang, E. Lee, W. B. Lee and M. J. Park, Science, 2024, 383, 70–76.
68. M. Krone, J. E. Stone, T. Ertl and K. Schulten, Eurographics Conference on Visualization, 2012.
69. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
70. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
71. F. A. Detcheverry, H. Kang, K. C. Daoulas, M. Müller, P. F. Nealey and J. J. de Pablo, Macromolecules, 2008, 41, 4989–5001.
72. M. Liu, Y. Qiang, W. Li, F. Qiu and A. C. Shi, ACS Macro Lett., 2016, 5, 1167–1171.
73. A. J. Mueller, A. P. Lindsay, A. Jayaraman, T. P. Lodge, M. K. Mahanthappa and F. S. Bates, ACS Macro Lett., 2020, 9, 576–582.
74. Z. Xu and W. Li, Chin. J. Chem., 2022, 40, 1083–1090.
75. F. Drolet and G. H. Fredrickson, Macromolecules, 2001, 34, 5317–5324.
76. R. K. W. Spencer and M. W. Matsen, Macromolecules, 2017, 50, 1681–1687.
77. F. S. Bates, M. F. Schulz, A. K. Khandpur, S. Förster, J. H. Rosedale, K. Almdal and K. Mortensen, Faraday Discuss., 1994, 98, 7–18.
78. T. M. Beardsley and M. W. Matsen, J. Chem. Phys., 2021, 154, 124902.
79. M. W. Matsen, T. M. Beardsley and J. D. Willis, Phys. Rev. Lett., 2023, 130, 248101.
80. M. W. Matsen, T. M. Beardsley and J. D. Willis, Phys. Rev. Mater., 2023, 7, 105605.
81. O. N. Vassiliev and M. W. Matsen, J. Chem. Phys., 2003, 118, 7700–7713.
82. C. To Lai and A.-C. Shi, Macromol. Theory Simul., 2021, 30, 2100019.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4me00138a

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

---