Morphological Diversity in Diblock Copolymer Solutions: A Molecular Dynamics Study

: Coarse-grained molecular dynamics simulations that incorporate explicit water-mediated hydrophilic/hydrophobic interactions are employed to track spatiotemporal evolution of diblock copolymer aggregation in initially homogeneous solutions. A phase portrait of the observed morphologies and their quantitative geometric features such as aggregation numbers, packing parameters, and radial distribution functions of solvent/monomers are presented. Energetic and entropic measures relevant to self-assembly such as speciﬁc solvent accessible surface area (SASA) and probability distribution functions (pdfs) of segmental stretch of copolymer chains are analyzed. The simulations qualitatively capture experimentally observed morphological diversity in diblock copolymer solutions. Topologically simpler structures predicted include spherical micelles, vesicles (polymersomes), lamellae (bilayers), linear wormlike micelles, and tori. More complex morphologies observed for larger chain lengths and nearly symmetric copolymer compositions include branched wormlike micelles with Y-shaped junctions and cylindrical micelle networks. For larger concentrations, vesicle strands, held together by hydrogen bonds, and “giant” composite aggregates that consist of lamellar, mixed hydrophobic/hydrophilic regions and percolating water cores are predicted. All structures are dynamic and exhibit diffuse domain boundaries. Morphology transitions across topologically simpler structures can be rationalized based on specific SASA measurements. PDFs of segmental stretch within vesicular assemblies appear to follow a log-normal distribution conducive for maximizing configuration entropy.

Chemical heterogeneity of BCPs finds applications in precision manufacturing of advanced thermoplastic elastomers [5]. Further, they have been extensively used in industry as emulsifiers and compatibilizers. In such roles, they act as traditional surfactants by reducing the interfacial tension between organic and aqueous phases. For instance, assemblies of diblock copolymers comprising hydrophilic and hydrophobic portions could be used to separate and encapsulate hydrocarbons from subterranean porous rock formations to aid enhanced oil recovery processes. Further, BCPs are typically non-ionic. This makes them environmentally less-intrusive candidates for isolating and scavenging spilled oil or organic contaminants from the environment. A distinct advantage of BCPs over traditional surfactants is the ability to vary their chain lengths and comparative chemical affinities of their hydrophilic and hydrophobic segments. Consequently, BCPs have been employed as molecular bridges in the synthesis of mesoporous materials. For instance, BCP molecules have been used to control the overlap distance between the hydrophobic hydration spheres during the growth of inorganic-organic composite structures in zeolite synthesis to manipulate the structure and catalytic function of the final product [21].
Renewed interest in BCP self-assembly in solution was catalyzed in the 1990s by the pioneering discoveries of vesicular morphologies, referred to as polymersomes [8][9][10][11][12][13][14]: see reviews by Bleul et al. [22], Mohammadi et al. [23], Zhang and Zhang [24], Discher et al. [25], and Antonietti and Förster [15]. The nomenclature (polymersomes) is borrowed from the literature on biological vesicles, or liposomes, which are made up of amphiphilic phospholipid molecules. However, as summarized in the review by Bleul et al. [22], BCPs offer robust means of size control of the vesicles by varying the chemical composition and/or molecular weight of the copolymer chain, changing the polymer concentration, the addition of secondary molecules capable of inserting themselves into the assembly and act as spacers, as well as microfluidics-based flow manipulation to produce droplets of controllable sizes. As a result, polymersomes of sizes that range between several nanometers to tens of micrometers can be routinely manufactured with low degrees of polydispersity. Such morphologies form the building blocks of complex nanostructures that offer significant promise in nanomedicine for applications as drug carriers for targeted delivery of therapeutic agents to treat diseases including many cancers, in gene therapy, for breaching biological barriers such as the blood-brain barrier and intestinal epithelium, and as diagnostic agents: see for instanc Table 1 in Karayianni and Pispas [3]. For example, polymersomes and rodlike micelles of PB-PEO copolymers loaded with the anti-cancer drug paclitaxel have experimentally displayed the controlled release profile and long-term stability necessary for inhibiting the growth of MCF-7 human breast cancer cells [26]. Diblock copolymer micelles are also being widely used in targeted drug delivery for lipophilic drugs such as FK506, administered in organ transplant recipients for immunosuppression [27,28]. The versatility and controllability of physical properties as well as biochemical activity of biodegradable polymersome-based nano-constructs are being actively explored for cancer theranostics, an approach that integrates diagnosis (e.g., by imaging) and therapeutics [23].
Despite the insights gained from extensive experimental studies on the phase behavior of BCPs in solution, several questions regarding the thermodynamic and kinetic pathways underlying their self-assembly leading to diverse morphologies remain unanswered. First, real time visualization of spontaneous formation of molecular aggregates in such systems is beyond the spatiotemporal resolution of advanced imaging technologies. Consequently, morphology characterization is achieved often by using cryo-TEM experiments although progress in NMR imaging has made it possible to track the trajectories of individually tagged molecules in native solution as they navigate within and across molecular aggregates [29,30]. Second, experimental observations of shape transitions are often rationalized based on elegant, yet simplified geometric arguments, such as the one based on the packing parameter, which was developed to explain spherical, cylindrical, and vesicular aggregates in amphiphilic surfactants [31]. However, experiments show that self-assembly of BCPs in solution leads to extraordinarily diverse morphologies. Further, a clear demarcation of hydrophilic and hydrophobic domains does not exist in BCP aggregates: e.g., coarse-grained molecular dynamics (CGMD) simulations show significant degrees of chain entanglement and interdigitation of the hydrophobic segments of the upper and lower leaflets of BCP bilayers [32,33].
Questions regarding entropic (e.g., polymer chain stretch and orientation) and energetic (e.g., solvent accessible surface area) changes as well as solvent transport mechanisms underlying morphology transitions, especially rod to lamella as well as lamella to vesicle transitions, remain largely unresolved. For instance, the proposed mechanisms of vesicle formation have been at best speculative. One proposed pathway in which a sufficiently large lamellar structure develops curvature, bends, and folds into a vesicle is based on energetic considerations. Specifically, it is envisioned that during such morphology change, the lamella discards the energy associated with the unfavorable hydrophobic interactions of its edge groups with the surrounding solvent molecules at the expense of gaining curvature energy [15,22,34]. Another proposed pathway entails solvent diffusion into large spherical micelles creating a polymer-free core region [22]. On the other hand, CGMD simulations conducted by Markvoort and coworkers [35] on phospholipid systems predict an increase in the overall potential energy associated with bilayer to vesicle transition, suggesting the existence of entropic driving forces. Configurational entropy associated with the distributions of chain stretch and orientation of hydrophobic/hydrophilic segments within BCP aggregates could influence their macroscopic mechanical properties, e.g., bending rigidity of bilayers. This in turn could have pronounced implications on shape transitions. To date, the relationships among internal molecular architecture of BCP aggregates, their mechanical properties, and shape dynamics remain poorly understood.
The overarching goal of the present work is to develop faithful molecular dynamics (MD) simulations to study spontaneous self-assembly in a model (PB-PEO) diblock copolymer solution leading to the formation of diverse morphologies. The initial conditions of our simulations correspond to a homogeneous BCP solution at a prescribed concentration in which polymer molecules of a given chain length and composition are randomly distributed. No pre-assembly is employed, and spatiotemporal evolution of structures is tracked in real time. In the past, both CGMD simulations with pre-assembled initial conditions [32,33] and dissipative particle dynamics techniques [36,37] have been employed to study self-assembly, bilayer to vesicle transition and kinetics of copolymer exchange between self-assembled aggregates in BCP solutions. In the present study, we employ CGMD simulations which are adapted from our previous investigations of self-assembly and shape transitions in as well as shear/elongational flow dynamics and rheology of cationic surfactant solutions [38][39][40][41][42][43][44][45]. Our simulations utilize the coarse graining approach and force fields provided by the MARTINI framework [46]. Utilizing coarse-grained molecular representations allow for conducting simulations that span sufficiently large length and time scales to capture spontaneous self-assembly processes in the presence of explicit solvent-polymer interactions.
Simulation details and data analytics methods are discussed in Section 2. In Section 3, we summarize the results of a set of over 200 (excluding replicates) CGMD simulations conducted by systematically varying the chain length, hydrophobic to hydrophilic segmental composition, and polymer concentration. Overall, the simulations not only capture experimentally observed morphologies but also render significant geometric, topological, and configurational details to our understanding of the structure of BCP aggregates. Various morphologies predicted by the simulations are organized into a phase diagram. Insights and intuitions provided by this work, together with those gained from experiments, enhance our understanding of morphological diversity in BCP solutions and could lead to the development of more accurate theories of BCP aggregation in the future.

Coarse-Grained (CG) Polymer Model
Even with the computational prowess available today, it is impractical to access, using all atom descriptions of the polymer and solvent, the length and time scales over which spontaneous copolymer self-assembly takes place in solution. On the other hand, mesoscopic models such as the ones employed in Brownian [47] or dissipative particle dynamics simulations [48] neglect explicit inter-molecular and solvent-mediated physicochemical interactions that govern the self-assembly process. For instance, in Brownian dynamics, the influence of the solvent on polymer motion is modeled by a random force in accordance with the fluctuation-dissipation theorem. In this work, our goal is to simulate structure evolution in the presence of explicit solvent-mediated forces, starting from an initially homogeneous copolymer solution, i.e., without resorting to pre-assembly of molecules as done in the past [32,33]. Indeed, atomistic molecular dynamics (MD) simulations in which the positions and momenta of solvent and polymer atoms are evolved in time subject to Newton's second law of motion would provide the most detailed portrait of the self-assembly process. However, the number of degrees of freedom and ultrasmall time steps (typically much smaller than the time scale of inter-atomic bond fluctuations) required in atomistic MD make it computationally intractable for conducting simulations that span length and time scales sufficient to capture spontaneous self-assembly of polymeric molecules. Further, interatomic interactions within monomeric units do not have an appreciable influence on thermodynamic self-assembly induced primarily by hydrophobic interactions. Hence, to realize an ensemble of simulations which collectively can provide statistically meaningful results on the effect of concentration, molecular weight, and composition of the copolymer on morphology transitions in monodisperse diblock BCP solutions, we employ a coarsegrained (CG) MD framework. Such CGMD models have been used successfully to probe the role of solvent-and salt-mediated interactions on the aqueous-phase self-assembly of surfactants [44,45], study the topology and energetics of wormlike micelles and micelle networks [41,43], understand flow-induced deformation and rupture of micelles [40,42], and explore the structure and rheology of micelle-nanoparticle complexes [38,39]. A similar approach is employed in the present study.
CG molecular modeling approaches originate from the pioneering work of Levitt and Warshel in the 1970s [49,50]. In the present study, we employ the MARTINI force field, developed by Marrink and coworkers at the University of Groningen, initially used in 2004 for molecular dynamics simulation of lipids [51] and later (2007) extended to various other molecules [46]. In this representation, several atoms are represented by a single CG bead. Based on their interaction sites, the beads are classified into four main types, namely Q, P, N, and C representing beads with charged, polar, nonpolar, and apolar characteristics, respectively. Each bead is further classified into many subtypes, which are used to represent the chemical nature of the underlying atomic structure more accurately [46]. Within one of the four main types, subtypes are either distinguished by their hydrogen-bonding capabilities (d: donor, a: acceptor, da: both donor and acceptor, 0: none), or by a number indicating the degree of polarity, which ranges from 1 (low polarity) to 5 (high polarity). Ring structures are denoted by the prefix 'S' [46].
In the present study, water (four molecule cluster, W), PB (one monomer), and PEO (one monomer) are represented by bead types P4, C4, and SNda, respectively [46,52]. The solvent consists of 90.9% of P4 water beads along with 9.1% of "antifreeze" beads (WF), which are of type BP4 and slightly bigger than the P4 beads. The WF beads are employed to disturb the lattice packing of the uniformly sized solvent beads and decrease the solvent freezing point to less than 250 K [46,53]. Figure 1 shows an example of the CG representation of a single copolymer chain that is composed of 6 butadiene and 6 ethylene oxide monomers: see Table S1 for additional details. In the simulations, the number of monomers is varied systematically to study the effect of chain length and chain composition on morphology evolution. Similar representation of PEO was used by Yang [54] in CGMD simulations to study dynamics of PEO chains in uniform shear flow.
Colloids Interfaces 2023, 7, x 5 of 23 morphology evolution. Similar representation of PEO was used by Yang [54] in CGMD simulations to study dynamics of PEO chains in uniform shear flow. The stretching and bending interactions between chemically bonded beads are modeled by weak harmonic potentials Vb and Vθ, respectively, by Equations (1) and (2) as where b0 and θ0 represent the equilibrium bond distance and equilibrium bond angle at the minimum energy configuration, respectively, and Kb and Kθ represent constants associated with bond stretching and bending energies, respectively. The non-bonded interactions between beads i and j are described by a 6-12 Lennard-Jones (LJ) potential given by 12 − (6) 6 . ( Equation (3) can also be written as where (12) = 4 12 and (6) = 4 6 .
In the above equations, r represents the distance between beads i and j, is the depth (minimum) of the potential well (energy function), and is the distance at which is zero. and are also referred to as the interaction strength and cutoff length, respectively.
MARTINI force field requires the LJ potential to be shifted smoothly to 0 between 0.9 and 1.2 nm, rather than using a hard cutoff to avoid singularities in the computation of the inter-particle forces. This is achieved by using a potential-switch function SV given by where = 6 12. The stretching and bending interactions between chemically bonded beads are modeled by weak harmonic potentials V b and V θ , respectively, by Equations (1) and (2) as where b 0 and θ 0 represent the equilibrium bond distance and equilibrium bond angle at the minimum energy configuration, respectively, and K b and K θ represent constants associated with bond stretching and bending energies, respectively. The non-bonded interactions between beads i and j are described by a 6-12 Lennard-Jones (LJ) potential given by Equation (3) can also be written as where C In the above equations, r represents the distance between beads i and j, ε ij is the depth (minimum) of the potential well (energy function), and σ ij is the distance at which V LJ is zero. ε ij and σ ij are also referred to as the interaction strength and cutoff length, respectively.
MARTINI force field requires the LJ potential to be shifted smoothly to 0 between 0.9 and 1.2 nm, rather than using a hard cutoff to avoid singularities in the computation of the inter-particle forces. This is achieved by using a potential-switch function S V given by where α = 6 or 12.
In Equation (7), the parameters A, B, and C are defined by and where values of r shi f t = 0.9 nm and r cut = 1.2 nm are used. The force field parameters used in the present simulations are listed in Tables S2-S4 [46,52,54].

Simulation Details
The simulations are performed by using the Gromacs 2020.2 MD software. First, simulations are conducted to minimize the energy of the initial system. This is followed by a short (a few ns in duration) NVT simulation to equilibrate the system at a desired temperature to ensure algorithmic stability. Subsequently, a sufficiently long NPT simulation (production run) is carried out. Velocity generation at the outset of a simulation, which is carried out by sampling from a Gaussian distribution that yields the mean temperature, is imperfect. When coupled with a barostat, initial velocity distributions thus generated can frequently lead to numerical instabilities. Hence, equilibration is better performed for an NVT ensemble for a short period of time to get the correct velocity distribution. Proceeding NPT simulations lead to the establishment of the appropriate system density.
The initial simulation box is cubic with a linear dimension of 40 nm. In comparison, the contour length of the longest copolymer chain is 6.58 nm. The number of water beads (W+WF) is held constant at 440,000, and the number of copolymer beads (chains) is varied to achieve a desired polymer concentration. During the NPT simulation, the box size changes to between 39 and 41 nm to match the correct system density.
To generate a phase diagram of self-assembled morphologies in the chain length (N)-mole fraction (x) space, the length of the copolymer chains is varied from 6 to 15 beads. For each chain length, we carry out simulations for different compositions. The composition sequence is denoted by 2PB(N-2)PEO, 3PB(N-3)PEO, . . . . . . , and (N-2)PB2PEO. For example, for a chain length N = 8, we perform simulations for copolymers with compositions denoted by 2PB6PEO, 3PB5PEO, 4PB4PEO, 5PB3PEO, and 6PB2PEO. Additionally, we conduct simulations for 0PB11PEO, 1PB10PEO, 10PB1PEO, and 11PB0PEO. For each composition, simulations are replicated with three different initial conditions. Hence, the total number of simulations used to construct the phase diagram in the present study is 237. The polymer concentration is held approximately constant at 10 mol% in these simulations.
The above choice of polymer concentration was made based on (i) considerations of computational expenses which increases rapidly with increasing number of degrees of freedom and (ii) experimental observations that a variety of self-assembled morphologies manifests at comparable polymer concentrations in PB-PEO and similar systems [11][12][13][14]. Under these simulation conditions, the number of copolymer beads ranges from 7500 to 4000, thus making the total number of monomers between 42,000 and 60,000. The chain lengths used in the simulations translate to polymer molecular weights (MWs) that range between 264 Da (6PEO) and 810 Da (15PB). Hence, these are relatively short polymers with MWs comparable to surfactants with long hydrocarbon chains, e.g., cetyl-trimethyl ammonium chloride, a C-16 cationic surfactant studied extensively in experiments and our previous CGMD simulations [38][39][40][41][42][43][44][45] has a MW of 320 Da. However, as shown in Section 3, the ability to manipulate the relative lengths of the hydrophilic (PEO) and hydrophobic (PB) segments within a polymer chain allows for the generation of diverse morphologies in copolymer solutions as compared to the limited set of structures observed in surfactant systems. We note that the persistence lengths of the PEO and PB chains are comparable to their bond lengths. This renders significant chain flexibility, a property that makes the copolymers amenable to a broad range of packing protocols.
To study concentration effects, we use the symmetric 6PB6PEO chain and vary the copolymer molar concentrations from 5.8% to 38.0%, which correspond to mass concentra-tions ranging from 4.0% to 29.5%. The reference pressure and temperature are 1 bar and 300 K, respectively. A v-rescale thermostat [55] with a time constant of 1.0 ps is used for temperature coupling. Berendsen barostat [56] with a time constant of 4.0 ps is used for isotropic pressure coupling during the NPT simulations. The simulations employ periodic boundary conditions along all three spatial coordinates. The time step for equilibration runs is 20 fs, and the run time is 2 ns. The time step for production runs is 40 fs, and the run time is up to 600 ns.

Data Analytics: Methods and Tools
Definition of a cluster: any two beads with center-to-center separation less than the cutoff distance are regarded as part of a cluster (aggregate). Further, within any bead in each cluster, there exists at least one bead at a distance less than the cutoff value. In this study, we use 0.5 nm as the cutoff distance. This is based on the cutoff distance of 0.47 nm for non-bonded LJ interactions between PB and PEO monomers.
Radius of gyration (R g ): the radius of gyration is used to quantify the size of selfassembled clusters. Radius of gyration is the root mean square distance of the beads in the cluster from its center of mass and is given by where r i is the distance from the ith bead to the center of mass. Radial distribution function, RDF (g(r)): radial distribution function describes how the density of B, EO, or W beads varies as a function of distance from the center of mass of a cluster. A larger value of g(r) signifies tighter packing of the beads. RDF is given by where < ρ(r) > is the average density of beads at a distance r from the center of mass, ρ = 3n m /(4πr 3 max ) is the average cluster density, n m is the total number of monomer beads contained in the cluster, r max is the distance of the farthest bead from the center of mass of the cluster, V(r) = (4/3)π (r + ∆r) 3 − r 3 and In the summation in Equation (12), beads with centers at the inner boundary are included while those at the outer boundary are excluded and ∆r = 0.1 nm.
Solvent accessible surface area (SASA): The algorithm for computing SASA is illustrated in Figure S1. Specifically, a probe is slid along the periphery of the beads that constitute a cluster, and the area circumscribed by the motion of the center of the probe is equal to the SASA: see Eisenhaber et al. [57] for details. The radii of the probe, PEO monomer bead, and PB monomer bead are 0.25 nm, 0.215 nm, and 0.235 nm, respectively. Note that the probe radius is slightly greater than one half of the cutoff distance of 0.47 nm for non-bonded interaction between water and PB/PEO. Linear dimension, surface area, volume, and packing parameter of the aggregates: The data generated from MD simulations include a set of spatial coordinates of the various beads (i.e., the so-called points cloud). The coordinates of its inner and outer bounding surfaces are obtained by using the alpha-shape algorithm, illustrated in Figure 2 [58]. The volume V and surface area A of molecular aggregates are computed from a 3-dimensional mesh constructed from the boundary coordinates. Matlab TM 2023a software is employed for this purpose. The cluster radius is calculated as the average distance between the surface points and the center of mass. Sphericity S is calculated as S = (6π 1/2 V) 2/3 /A. For a sphere, S = 1.
For monolayer aggregates, the packing parameter, P, is defined as P ≡ v/al, where v, a and l denote the volume of the hydrophobic segment, surface area of the hydrophobic-hydrophilic interface, the hydrophobic chain length, respectively [3,31]. In this work, we employ the same definition for monolayer structures. For bilayers (micelles and vesicles) and mixed structures, a used in packing parameter calculations is the area of the outer PEO-PB interface and the area of the solvent-aggregate interface, respectively. These definitions yield p values of 1/3 for monolayers and mixed structures of spherical shape. For a perfectly spherical bilayer consisting of PB and PEO segments of equal length, the above definition yields p ≈ 0.48.
Colloids Interfaces 2023, 7, x 8 of 23 surface points and the center of mass. Sphericity S is calculated as = (6 1/2 ) 2/3 / . For a sphere, S = 1. For monolayer aggregates, the packing parameter, P, is defined as P ≡ v/al, where v, a and l denote the volume of the hydrophobic segment, surface area of the hydrophobic-hydrophilic interface, the hydrophobic chain length, respectively [3,31]. In this work, we employ the same definition for monolayer structures. For bilayers (micelles and vesicles) and mixed structures, a used in packing parameter calculations is the area of the outer PEO-PB interface and the area of the solvent-aggregate interface, respectively. These definitions yield p values of 1/3 for monolayers and mixed structures of spherical shape. For a perfectly spherical bilayer consisting of PB and PEO segments of equal length, the above definition yields p ≈ 0.48.

Software Used
The initial configurations are built by using PACKMOL 20.14.0 [59]. The simulations are performed using Gromacs 2020.2 [60]. The molecular visualization is performed by using VMD 1.9.3 [61]. Data analysis and graphing are performed using Matlab TM 2023a, WPS office TM 2021, C++, and OriginPro TM 2017.

Vesicles (Polymersomes)
We first describe the various self-assembled morphologies observed in the simulations for different xPEO values for copolymer concentration of 10 mol%. Figure 3a shows the cross section (inset) of a representative vesicular structure and the corresponding RDFs of PB, PEO, and water, respectively, for a 5PB6PEO copolymer. RDFs of vesicles formed by copolymers of different segmental compositions are also shown in Figure 3bd. Table S5 summarizes the geometric features of the vesicles whose RDFs are shown in

Software Used
The initial configurations are built by using PACKMOL 20.14.0 [59]. The simulations are performed using Gromacs 2020.2 [60]. The molecular visualization is performed by using VMD 1.9.3 [61]. Data analysis and graphing are performed using Matlab TM 2023a, WPS office TM 2021, C++, and OriginPro TM 2017.

Vesicles (Polymersomes)
We first describe the various self-assembled morphologies observed in the simulations for different x PEO values for copolymer concentration of 10 mol%. Figure 3a shows the cross section (inset) of a representative vesicular structure and the corresponding RDFs of PB, PEO, and water, respectively, for a 5PB6PEO copolymer. RDFs of vesicles formed by copolymers of different segmental compositions are also shown in Figure 3b-d. Table S5 summarizes the geometric features of the vesicles whose RDFs are shown in Figure 3a-d.
Common features of the vesicular morphologies are the presence of a polymer-free core (center region), two maxima in the PEO RDF, and a single maximum in the PB RDF in between the two peaks in the PEO RDF. Water density in the core region exhibits a maximum close to the center and the inner PEO-water interface of the vesicle. This is attributed to the existence of hydrogen-bonded water-ethylene oxide pairs. The PB layer is almost devoid of water, as expected. Further, the shapes are nearly spherical with S ranging between 0.94 and 0.96 (Table S5). However, the RDFs show that the distribution of the hydrophobic and hydrophilic entities within the vesicle are far from those expected based on ideal geometric depictions often employed in the literature. For instance, the interface between PB and PEO segments is diffuse as inferred from the overlap between their RDFs. Second, the distribution of monomers depends on the chain composition (x PEO ) with chains with PB to PEO ratio of unity (cases 4d) exhibiting a nearly symmetric distribution of PB segments. Third, we observe the second (outer) peak in the PEO RDF to be smaller than the first. This is attributed to the more crowded packing of PEO segments around the water core. As seen from Table S5, the density of water within the vesicle core is smaller than that in the bulk. This is a consequence of hydrogen bonding interactions between PEO and water beads and the steric exclusion due to the confinement provided by the PEO-water interface, which prevents a homogeneous distribution of water within the nanometer-sized cores of the vesicles observed in our simulations.
Colloids Interfaces 2023, 7, x 9 of 23 Common features of the vesicular morphologies are the presence of a polymer-free core (center region), two maxima in the PEO RDF, and a single maximum in the PB RDF in between the two peaks in the PEO RDF. Water density in the core region exhibits a maximum close to the center and the inner PEO-water interface of the vesicle. This is attributed to the existence of hydrogen-bonded water-ethylene oxide pairs. The PB layer is almost devoid of water, as expected. Further, the shapes are nearly spherical with S ranging between 0.94 and 0.96 (Table S5). However, the RDFs show that the distribution of the hydrophobic and hydrophilic entities within the vesicle are far from those expected based on ideal geometric depictions often employed in the literature. For instance, the interface between PB and PEO segments is diffuse as inferred from the overlap between their RDFs. Second, the distribution of monomers depends on the chain composition (xPEO) with chains with PB to PEO ratio of unity (cases 4d) exhibiting a nearly symmetric distribution of PB segments. Third, we observe the second (outer) peak in the PEO RDF to be smaller than the first. This is attributed to the more crowded packing of PEO segments around the water core. As seen from Table S5, the density of water within the vesicle core is smaller than that in the bulk. This is a consequence of hydrogen bonding interactions between PEO and water beads and the steric exclusion due to the confinement provided by the PEO-water interface, which prevents a homogeneous distribution of water within the nanometer-sized cores of the vesicles observed in our simulations. Table S5 shows that the packing parameters of the vesicles are close to 1/2, which is the lower limit of p for polymersomes. This is expected since the sizes of the water core and the copolymer bilayer (corona) are comparable to each other. For instance, for a symmetric copolymer chain (5PB5PEO), the calculated value of p is in good agreement with that for a spherical bilayer structure (≈0.48). Further, the cluster density of all the vesicles is nearly constants at approximately 10 beads per nm 3 . However, we see significant variations in the aggregation numbers ranging from 500 to 900 depending on the chain length and ratio PB to PEO composition within the copolymer. To probe whether such variability in n might be rationalized based on thermodynamic considerations, we calculate the  Table S5 shows that the packing parameters of the vesicles are close to 1/2, which is the lower limit of p for polymersomes. This is expected since the sizes of the water core and the copolymer bilayer (corona) are comparable to each other. For instance, for a symmetric copolymer chain (5PB5PEO), the calculated value of p is in good agreement with that for a spherical bilayer structure (≈0.48). Further, the cluster density of all the vesicles is nearly constants at approximately 10 beads per nm 3 . However, we see significant variations in the aggregation numbers ranging from 500 to 900 depending on the chain length and ratio PB to PEO composition within the copolymer. To probe whether such variability in n might be rationalized based on thermodynamic considerations, we calculate the probability distribution function (pdf) of the end-to-end distance Q, p(Q), of the hydrophobic as well as the inner and outer hydrophilic polymer segments within the aggregates. The information entropy H A associated with the ensemble of molecules that makes up the aggregate can be inferred from such data. Specifically, H A ≡ −∑ p ln p, where p is the pdf associated with the distribution of particles in the ensemble [62]. Figure 4 shows the pdfs of s ≡ Q max − Q for a representative vesicle. The distributions are skewed towards lower s values and have long tails. The best fit of p(s) to a log-normal distribution is also shown. The qualitative agreement with a log-normal distribution of the chain-stretch parameter (s), which was also observed for vesicular structures with different N and x PEO in our simulations, is intriguing. It is well known from probability theory that for a random variable with given expectation and variance, log-normal distribution maximizes the information entropy [63]. Further investigations are needed to establish whether such entropy maximization principles constitute self-assembly motifs in BCP solutions.
lower s values and have long tails. The best fit of p(s) to a log-normal distribution is also shown. The qualitative agreement with a log-normal distribution of the chain-stretch parameter (s), which was also observed for vesicular structures with different N and xPEO in our simulations, is intriguing. It is well known from probability theory that for a random variable with given expectation and variance, log-normal distribution maximizes the information entropy [63]. Further investigations are needed to establish whether such entropy maximization principles constitute self-assembly motifs in BCP solutions.   Table S6. It is evident from Figure 5 that for a given chain length, compositional variations have a pronounced influence on the micelle structure. For the 2PB4PEO system (xPEO ≈ 0.67), a PEO-dominated micelle with a small PB core, surrounded by a wellmixed outer region of PB and PEO monomers, are observed. However, for 3PB3PEO chains (xPEO = 0.5), bilayer micelles with well-separated hydrophilic and hydrophobic domains are seen. Further, when xPEO is decreased to ≈0.33 (4PB2PEO system), PB-dominated multilayer micelles with a PEO-rich core, surrounded by a well-mixed intermediate domain and an external region consisting of a PB-rich shell and a PEO-solvent interface, are detected. As seen from Table S6, the aggregation numbers vary significantly across these composite micelle structures: asymmetric chains (4PB2PEO and 2PB4PEO) that form mixed multilayer structures have significantly higher n values as compared to that of the bilayer formed by symmetric chains (3PB3PEO). The packing parameter values are less than 1/3, as expected for nearly spherical assemblies. However, the cluster density is approximately 10 beads/nm 3 , as it is the case for vesicles observed in this study. This parameter is determined by the bond lengths and non-bonded LJ interaction forces. Figure 6 shows representative cross sections and RDFs of nearly spherical micelles made up of relatively longer BCP chains (N ≥ 11). Corresponding geometric parameters are reported in Table S7. The cluster formed by chains with a large hydrophobic fraction (xPEO ≈ 0.15), has a PB-rich core devoid of water and a relatively sharp PEO-water interface. In contrast, for chains with larger hydrophilic fraction (xPEO ≈ 0.42 and 0.55), a bilayer structure with a well-defined PEO core is observed. These structures exist adjacent to vesicular and lamellar morphologies in the N-xPEO phase diagram shown in Figure 11. Structures   Table S6. It is evident from Figure 5 that for a given chain length, compositional variations have a pronounced influence on the micelle structure. For the 2PB4PEO system (x PEO ≈ 0.67), a PEO-dominated micelle with a small PB core, surrounded by a well-mixed outer region of PB and PEO monomers, are observed. However, for 3PB3PEO chains (x PEO = 0.5), bilayer micelles with well-separated hydrophilic and hydrophobic domains are seen. Further, when x PEO is decreased to ≈0.33 (4PB2PEO system), PB-dominated multilayer micelles with a PEO-rich core, surrounded by a well-mixed intermediate domain and an external region consisting of a PB-rich shell and a PEO-solvent interface, are detected. As seen from Table S6, the aggregation numbers vary significantly across these composite micelle structures: asymmetric chains (4PB2PEO and 2PB4PEO) that form mixed multilayer structures have significantly higher n values as compared to that of the bilayer formed by symmetric chains (3PB3PEO). The packing parameter values are less than 1/3, as expected for nearly spherical assemblies. However, the cluster density is approximately 10 beads/nm 3 , as it is the case for vesicles observed in this study. This parameter is determined by the bond lengths and non-bonded LJ interaction forces.     Table S7. The cluster formed by chains with a large hydrophobic fraction (x PEO ≈ 0.15), has a PB-rich core devoid of water and a relatively sharp PEO-water interface. In contrast, for chains with larger hydrophilic fraction (x PEO ≈ 0.42 and 0.55), a bilayer structure with a well-defined PEO core is observed. These structures exist adjacent to vesicular and lamellar morphologies in the N-x PEO phase diagram shown in Figure 11. Structures similar to those shown in Figure 6a,c or Figure 6e are referred to as PB-dominated micelles (PB-S) and spherical micelles, respectively, in the phase diagram.   panels (a,b)) and spherical micelles (7PB5PEO, panels (c,d); 5PB6PEO, panels (e,f)).

Irregular-Shaped PEO-Rich Micelles
For longer chains, increasing the hydrophilic fraction leads to the formation of irregular-shaped, single-layer structures with a hydrophobic core and a hydrophilic solvent interface, as shown in Figure 7a,b (3PB8PEO, xPEO = 0.73). The aggregation numbers for these structures range between ≈ 150 and 450 (n = 213 for the cluster shown in Figure 7a). In Figure 7c,d, we show the compact PB-dominated spherical micelle observed for 8PB3PEO (xPEO = 0.27) to contrast it with the irregular-shaped aggregate realized for the complementary chain composition xPEO = 0.73. Overall, BCP chains with larger hydrophilic fractions tend to form water interfaces with greater PEO stretch and mobility. This reduces the sphericity and packing parameter of the BCP assembly as seen in Table S8.  panels (a,b)) and spherical micelles (7PB5PEO, panels (c,d); 5PB6PEO, panels (e,f)).

Irregular-Shaped PEO-Rich Micelles
For longer chains, increasing the hydrophilic fraction leads to the formation of irregularshaped, single-layer structures with a hydrophobic core and a hydrophilic solvent interface, as shown in Figure 7a,b (3PB8PEO, x PEO = 0.73). The aggregation numbers for these structures range between ≈ 150 and 450 (n = 213 for the cluster shown in Figure 7a). In Figure 7c,d, we show the compact PB-dominated spherical micelle observed for 8PB3PEO (x PEO = 0.27) to contrast it with the irregular-shaped aggregate realized for the complementary chain composition x PEO = 0.73. Overall, BCP chains with larger hydrophilic fractions tend to form water interfaces with greater PEO stretch and mobility. This reduces the sphericity and packing parameter of the BCP assembly as seen in Table S8.

Linear and Branched Wormlike Micelles, and Micelle Networks
For chain lengths N > 9 and hydrophilic fraction xPEO > 0.4, cylindrical aggregates are observed. These aggregates can grow into flexible, wormlike micelles as shown in Figure  8a. Flexibility is inferred from the large ratio of the contour length of the micelle to its persistence length. To calculate the latter, the micelle is partitioned into cylindrical segments with linear dimension of 2 nm and reconstructed as a chain that connects the centers of mass of these segments. The persistence length is obtained by fitting the average of the cosine of the angle subtended by the vectors connecting the center of masses of adjacent cylindrical segments to a decaying exponential function [64]. A small persistence length implies a fast loss of orientational correlation along the chain contour and consequently a high degree of chain flexibility. RDFs of PB and PEO monomers along the (circular) cross section of the wormlike structure show a PB-dominated core and a PEO corona. However, there is significant overlap between the radial distribution functions of the hydrophobic

Linear and Branched Wormlike Micelles, and Micelle Networks
For chain lengths N > 9 and hydrophilic fraction x PEO > 0.4, cylindrical aggregates are observed. These aggregates can grow into flexible, wormlike micelles as shown in Figure 8a. Flexibility is inferred from the large ratio of the contour length of the micelle to its persistence length. To calculate the latter, the micelle is partitioned into cylindrical segments with linear dimension of 2 nm and reconstructed as a chain that connects the centers of mass of these segments. The persistence length is obtained by fitting the average of the cosine of the angle subtended by the vectors connecting the center of masses of adjacent cylindrical segments to a decaying exponential function [64]. A small persistence length implies a fast loss of orientational correlation along the chain contour and consequently a high degree of chain flexibility. RDFs of PB and PEO monomers along the (circular) cross section of the wormlike structure show a PB-dominated core and a PEO corona. However, there is significant overlap between the radial distribution functions of the hydrophobic (PB) and hydrophilic (PEO) segments: the PB-PEO interface is diffuse and the overlap region spans ≈ 2 nm, which is greater than four times the length of the bond connecting the PEO and PB monomers within the copolymer. Further, the RDFs show that the micelle core is free of water.
observed. These aggregates can grow into flexible, wormlike micelles as shown in Figure  8a. Flexibility is inferred from the large ratio of the contour length of the micelle to its persistence length. To calculate the latter, the micelle is partitioned into cylindrical segments with linear dimension of 2 nm and reconstructed as a chain that connects the centers of mass of these segments. The persistence length is obtained by fitting the average of the cosine of the angle subtended by the vectors connecting the center of masses of adjacent cylindrical segments to a decaying exponential function [64]. A small persistence length implies a fast loss of orientational correlation along the chain contour and consequently a high degree of chain flexibility. RDFs of PB and PEO monomers along the (circular) cross section of the wormlike structure show a PB-dominated core and a PEO corona. However, there is significant overlap between the radial distribution functions of the hydrophobic (PB) and hydrophilic (PEO) segments: the PB-PEO interface is diffuse and the overlap region spans ≈ 2 nm, which is greater than four times the length of the bond connecting the PEO and PB monomers within the copolymer. Further, the RDFs show that the micelle core is free of water.  Cylindrical micelles of amphiphilic surfactants have been known to form branched structures, especially in the presence of added salt. Branch formation has been hypothesized as a mechanism of stress relaxation in micellar fluids [65]. Specifically, the micelle branches, unlike polymeric ones which are held together by strong covalent bonds, can slide along longer backbones of a network structure, thereby healing internal stresses. This mechanism has been confirmed by CGMD simulations [41]. The present study reveals the existence of branched wormlike micelles as well as micelle networks in BCP solutions as shown in Figure 8c,d, respectively. The former structures consist of junctions with Y-shaped branches as seen in the experiments of Jain and Bates [12] which are illustrated in Figure 8e. In fact, these authors suggested that lamellar, rodlike, and Y-shaped structures form the building blocks of complex assemblies seen in BCP solutions. The micelle network closely resembles the bi-continuous rodlike structures seen experimentally in PS-PAA diblock BCP solutions: see Figure 4c in Mai and Eisenberg [11]. The micelle structures shown in Figure 8 are dynamic. To illustrate this, we track the end-to-end distance and the orientation, as quantified by the angle θ subtended by the projection of the end-to-vector on the x-y plane and the z axis, of a wormlike micelle as function of time. During the time interval shown in Figure 9, the micelle can be seen to grow by merging with smaller aggregates. Further, the orientation angle undergoes temporal fluctuations with a relaxation time of approximately 40 ns. The relaxation time is calculated as the inverse of the decay rate of the best exponential fit of the autocorrelation function ACF of θ vs. time during the time interval over which Q remains practically constant. This procedure follows the algorithm employed by Sambasivam et al. [40] for rodlike surfactant micelles.
PAA diblock BCP solutions: see Figure 4c in Mai and Eisenberg [11]. The micelle structures shown in Figure 8 are dynamic. To illustrate this, we track the end-to-end distance and the orientation, as quantified by the angle θ subtended by the projection of the endto-vector on the x-y plane and the z axis, of a wormlike micelle as function of time. During the time interval shown in Figure 9, the micelle can be seen to grow by merging with smaller aggregates. Further, the orientation angle undergoes temporal fluctuations with a relaxation time of approximately 40 ns. The relaxation time is calculated as the inverse of the decay rate of the best exponential fit of the autocorrelation function ACF of θ vs. time during the time interval over which Q remains practically constant. This procedure follows the algorithm employed by Sambasivam et al. [40] for rodlike surfactant micelles.

Tori and Lamellae (Bilayers)
Figure 10a-c show a torus structure, PB organization within the torus, and the RDFs of PB and PEO monomers. This structure consists of a hydrophilic core and a corona. The intermediate region is rich in hydrophobic moieties. There is significant overlap between the PB and PEO RDFs suggesting the presence of diffuse hydrophobic-hydrophilic interfaces. Such diffuse interfaces are also observed in lamellar structures as illustrated in Figure 10d,e. Absence of sharp interfaces in bilayer structures was also reported in previous CGMD studies [32,33].

Phase Diagram
The BCP morphologies observed in our simulations are organized into a phase portrait in the N vs. xPEO space in Figure 11. For xPEO < 0.15, PB-dominated spherical micelles (Figure 6a) are observed for all chain lengths. In the limit as xPEO → 1, no micellization is observed, as expected. The phase portrait reveals that intermediate chain compositions

Phase Diagram
The BCP morphologies observed in our simulations are organized into a phase portrait in the N vs. x PEO space in Figure 11. For x PEO < 0.15, PB-dominated spherical micelles (Figure 6a) are observed for all chain lengths. In the limit as x PEO → 1, no micellization is observed, as expected. The phase portrait reveals that intermediate chain compositions (0.4 < x PEO < 0.7) lead to remarkable diversity in morphology. Specifically, consistent with experimental observations, simulations predict that rodlike, lamellar, and vesicular structures share phase boundaries [14]. Concentration, thermal, or mechanical (e.g., by flow shear) perturbations of solutions near the phase boundaries can induce topological transitions in BCP solutions. Moreover, the phase portrait shows that two or more morphologies coexist except for very low fractions of PEO or PB and very short chains. Furthermore, for slightly asymmetric hydrophilic chain compositions, i.e., around the neighborhood of x PEO = 0.6, the morphological complexity increases substantially with increasing chain length: see for examples regions Y and Z in the phase diagram where micelle networks ( Figure 8d) coexist with vesicles, spherical micelles, and/or tori. Pathways of morphology transitions among these structures need to be investigated further. We have provided representative real-time visualizations of structure evolution constructed from computed MD trajectories in Supplementary Materials (Videos S1-S4).

SASA and Morphology Transitions
In Figures 12 and 13, we report changes in the cluster size, as quantified by the radius of gyration Rg of the largest aggregate, and solvent accessible surface area (SASA) per monomer for the entire system, respectively. Figure 12 shows that Rg remains practically constant for xPEO values ≤ 0.6, although spherical micelle to vesicle transition is predicted for xPEO ≈ 0.3. However, Rg increases for xPEO > 0.6 as spherical structures merge to form rodlike micelles. Further increase in PEO composition leads to the breakage of rods into smaller PEO-dominated micelles of irregular shapes (Figure 7a). Figure 13 shows that specific SASA increases almost linearly with increasing xPEO in the region where spherical and vesicular structures are observed (xPEO < 0.6). This can be attributed to the enhanced hydrophilicity of the BCP chains in the case of spherical micelles and the availability of an

SASA and Morphology Transitions
In Figures 12 and 13, we report changes in the cluster size, as quantified by the radius of gyration R g of the largest aggregate, and solvent accessible surface area (SASA) per monomer for the entire system, respectively. Figure 12 shows that R g remains practically constant for x PEO values ≤ 0.6, although spherical micelle to vesicle transition is predicted for x PEO ≈ 0.3. However, R g increases for x PEO > 0.6 as spherical structures merge to form rodlike micelles. Further increase in PEO composition leads to the breakage of rods into smaller PEO-dominated micelles of irregular shapes (Figure 7a). Figure 13 shows that specific SASA increases almost linearly with increasing x PEO in the region where spherical and vesicular structures are observed (x PEO < 0.6). This can be attributed to the enhanced hydrophilicity of the BCP chains in the case of spherical micelles and the availability of an additional water-PEO interface for vesicles. As cylindrical structures form, the solvent-PEO interface also grows rapidly proportional to the lateral surface area of the flexible wormlike micelles (x PEO > 0.6).

Effects of Copolymer Concentration on Morphology
To explore the effect of copolymer concentration C, expressed in wt%, on the shape and internal structure of self-assembled structures, we perform simulations for C ranging from 4 to 29.5 for N = 12 and xPEO = 0.5, around which pronounced morphological diversity is observed. This point corresponds to the red circle in Figure 11. Figure 14 shows Rg of the largest cluster observed for each copolymer concentration and the number of clusters NC as functions of C. When C reaches 29.5, polymers self-assemble into "giant structures", that span the entire (periodic) simulation box. Hence Rg is not reported for this value of C

Effects of Copolymer Concentration on Morphology
To explore the effect of copolymer concentration C, expressed in wt%, on the shape and internal structure of self-assembled structures, we perform simulations for C ranging from 4 to 29.5 for N = 12 and xPEO = 0.5, around which pronounced morphological diversity is observed. This point corresponds to the red circle in Figure 11. Figure 14 shows Rg of the largest cluster observed for each copolymer concentration and the number of clusters NC as functions of C. When C reaches 29.5, polymers self-assemble into "giant structures", that span the entire (periodic) simulation box. Hence Rg is not reported for this value of C Figure 13. Specific SASA depends on composition and increases rapidly as sphere to rod transition occurs (N: inset).

Effects of Copolymer Concentration on Morphology
To explore the effect of copolymer concentration C, expressed in wt%, on the shape and internal structure of self-assembled structures, we perform simulations for C ranging from 4 to 29.5 for N = 12 and x PEO = 0.5, around which pronounced morphological diversity is observed. This point corresponds to the red circle in Figure 11. Figure 14 shows R g of the largest cluster observed for each copolymer concentration and the number of clusters N C as functions of C. When C reaches 29.5, polymers self-assemble into "giant structures", that span the entire (periodic) simulation box. Hence R g is not reported for this value of C as it is not a meaningful measure. Structure coarsening via plausible coalescence of smaller aggregates can be inferred from the data: N C is seen to decrease from 14 to 2 as C is increased from 4 to 29.5. This is also evident from the images of the cross sections of representative structures shown in Figure 15. Geometric features of the dominant structures observed for different C values are reported in Table S9. For the lowest C (=4), irregular-shaped micelles (S ≈ 0.85) are formed. Increasing C to 7.7 results in the formation of vesicles ((n ≈ 900, R ≈ 6.3 nm), which get larger with further increase in C to 11.1 (n ≈ 2700, R ≈ 9.4 nm) with a core thrice the size of the vesicle formed at C = 7.7. The examination of the pdfs of the magnitude of the segmental end-to-end vector Q for the two vesicles, shown in Figure 16, do not reveal appreciable differences between their means and variances. Hence, the increase in vesicle size with increasing C is a consequence of the enhancement in the aggregation number. as it is not a meaningful measure. Structure coarsening via plausible coalescence of smaller aggregates can be inferred from the data: NC is seen to decrease from 14 to 2 as C is increased from 4 to 29.5. This is also evident from the images of the cross sections of representative structures shown in Figure 15. Geometric features of the dominant structures observed for different C values are reported in Table S9. For the lowest C (=4), irregular-shaped micelles (S ≈ 0.85) are formed. Increasing C to 7.7 results in the formation of vesicles ((n ≈ 900, R ≈ 6.3 nm), which get larger with further increase in C to 11.1 (n ≈ 2700, R ≈ 9.4 nm) with a core thrice the size of the vesicle formed at C = 7.7. The examination of the pdfs of the magnitude of the segmental end-to-end vector Q for the two vesicles, shown in Figure 16, do not reveal appreciable differences between their means and variances. Hence, the increase in vesicle size with increasing C is a consequence of the enhancement in the aggregation number.   as it is not a meaningful measure. Structure coarsening via plausible coalescence of smaller aggregates can be inferred from the data: NC is seen to decrease from 14 to 2 as C is increased from 4 to 29.5. This is also evident from the images of the cross sections of representative structures shown in Figure 15. Geometric features of the dominant structures observed for different C values are reported in Table S9. For the lowest C (=4), irregular-shaped micelles (S ≈ 0.85) are formed. Increasing C to 7.7 results in the formation of vesicles ((n ≈ 900, R ≈ 6.3 nm), which get larger with further increase in C to 11.1 (n ≈ 2700, R ≈ 9.4 nm) with a core thrice the size of the vesicle formed at C = 7.7. The examination of the pdfs of the magnitude of the segmental end-to-end vector Q for the two vesicles, shown in Figure 16, do not reveal appreciable differences between their means and variances. Hence, the increase in vesicle size with increasing C is a consequence of the enhancement in the aggregation number.    For larger C values (=14. 3,20), we observe concatenated vesicles, or vesicle strands, as shown in Figure 17. Each substructure consists of a water core and PEO-PB-PEO shell. The RDFs (Figure 17d) of PB and PEO monomers clearly show that the molecular bridge that attaches constituent vesicles in the strand to each other is compositionally richer in the hydrophilic (PEO) component. We compare the segmental stretches of PEO and PB within the bridge with those farther from it on the vesicle surface. As shown in Figure 18, the segmental Q statistics are practically identical for the two groups. This suggests that the molecular bridges are held together by hydrogen bonding interactions between the PEO segments. Vesicle strands predicted by our simulations have been observed in experiments as well: see, for instance, Figures 7 and 8 in Antonietti and Förster [15]. These authors have rationalized, based on considerations of bending energy minimization, the genesis of vesicular structures of non-spherical shapes and complex topologies/symmetries via mechanisms including swelling/deswelling, budding, and fusion in copolymer systems with bilayer symmetry. A more comprehensive set of simulations and data analytics are necessary to elucidate and quantify the physico-chemical interactions that lead to the formation of such composite structures. For larger C values (=14.3, 20), we observe concatenated vesicles, or vesicle strands, as shown in Figure 17. Each substructure consists of a water core and PEO-PB-PEO shell. The RDFs (Figure 17d) of PB and PEO monomers clearly show that the molecular bridge that attaches constituent vesicles in the strand to each other is compositionally richer in the hydrophilic (PEO) component. We compare the segmental stretches of PEO and PB within the bridge with those farther from it on the vesicle surface. As shown in Figure 18, the segmental Q statistics are practically identical for the two groups. This suggests that the molecular bridges are held together by hydrogen bonding interactions between the PEO segments. Vesicle strands predicted by our simulations have been observed in experiments as well: see, for instance, Figures 7 and 8 in Antonietti and Förster [15]. These authors have rationalized, based on considerations of bending energy minimization, the genesis of vesicular structures of non-spherical shapes and complex topologies/symmetries via mechanisms including swelling/deswelling, budding, and fusion in copolymer systems with bilayer symmetry. A more comprehensive set of simulations and data analytics are necessary to elucidate and quantify the physico-chemical interactions that lead to the formation of such composite structures.
For the largest C value (=29.5) accessible to our present simulations, we observe "giant" structures that consists of lamellar regions, domains with mixed compositions of PEO and PB and water cores: see Figure S2 for representative visualizations of 2D cross sections and RDFs of PB and PEO monomers in the lamellar region. The 3D architectures of such composite structures are quite complex: the loci of the lamella and water core change as a function of depth within the interior of the structure: see Figure S3 and 3D visualizations provided in the Supplementary Materials (Videos S5-S7). Further investigations are needed to quantify their topological and percolation characteristics. We provide representative real-time visualizations of structure evolution of such "giant" structures from computed MD trajectories in Supplementary Materials (Videos S8-S10).  For the largest C value (=29.5) accessible to our present simulations, we observe "giant" structures that consists of lamellar regions, domains with mixed compositions of PEO and PB and water cores: see Figure S2 for representative visualizations of 2D cross sections and RDFs of PB and PEO monomers in the lamellar region. The 3D architectures of such composite structures are quite complex: the loci of the lamella and water core change as a function of depth within the interior of the structure: see Figure S3 and 3D visualizations provided in the Supplementary Materials (Videos S5-S7). Further investigations are needed to quantify their topological and percolation characteristics. We provide representative real-time visualizations of structure evolution of such "giant" structures from computed MD trajectories in Supplementary Materials (Videos S8-S10).

Conclusions and Outlook
Coarse-grained molecular dynamics simulations that incorporate explicit water-mediated hydrophilic/hydrophobic interactions are employed to systematically investigate  For the largest C value (=29.5) accessible to our present simulations, we observe "giant" structures that consists of lamellar regions, domains with mixed compositions of PEO and PB and water cores: see Figure S2 for representative visualizations of 2D cross sections and RDFs of PB and PEO monomers in the lamellar region. The 3D architectures of such composite structures are quite complex: the loci of the lamella and water core change as a function of depth within the interior of the structure: see Figure S3 and 3D visualizations provided in the Supplementary Materials (Videos S5-S7). Further investigations are needed to quantify their topological and percolation characteristics. We provide representative real-time visualizations of structure evolution of such "giant" structures from computed MD trajectories in Supplementary Materials (Videos S8-S10).

Conclusions and Outlook
Coarse-grained molecular dynamics simulations that incorporate explicit water-mediated hydrophilic/hydrophobic interactions are employed to systematically investigate

Conclusions and Outlook
Coarse-grained molecular dynamics simulations that incorporate explicit water-mediated hydrophilic/hydrophobic interactions are employed to systematically investigate self-assembly in diblock (PEO-PB) copolymer solutions. The simulations do not employ pre-assembled structures as initial conditions. Instead, they track spatiotemporal evolution of copolymer aggregation starting from an initial state that corresponds to a homogeneous BCP solution at a prescribed concentration. The effects of polymer chain length (N), mole fraction of the hydrophilic (PEO) block (x PEO ), and polymer concentration (C) are investigated. A phase portrait of the observed morphologies and a detailed discussion of their quantitative geometric features such as aggregation numbers, packing parameters, and radial distribution functions of solvent/monomers are presented. Simulation data are analyzed to glean insights into the energetics and entropy measures relevant to self-assembly such as specific SASA and pdfs of segmental stretch of hydrophilic and hydrophobic segments within the aggregates.
Overall, the simulations qualitatively capture experimentally observed morphological diversity in diblock BCP solutions [9][10][11][12][13]15,27,28]. Topologically, simpler structures predicted include spherical micelles, vesicles (polymersomes), lamellae (bilayers), flexible linear (wormlike) micelles, and tori. More complex morphologies observed for relatively large chain lengths and hydrophilic to hydrophobic molar composition ratio slightly greater than unity include branched wormlike micelles with Y-shaped junctions [12] and networks of cylindrical micelles [11]. For polymer concentrations ranging from 14 to 20, wt% formation of vesicle strands is predicted. The molecular bridges that concatenate the vesicles are compositionally rich in the hydrophilic component and are held together by hydrogen bonds. For polymer concentrations ≈ 30 wt%, structure coarsening leads to the formation of "giant" composite aggregates that consist of lamellar, mixed PEO/PB regions and percolating water cores.
The copolymer morphologies are dynamic, i.e., structure growth or dissociation due to addition/removal of copolymers, reorganization of constituent molecules within the aggregates, and chain length/orientation fluctuations are observed in the simulations. Further, even within a given morphology classification, e.g., spherical micelles or vesicles, significant variability exists in terms of aggregation numbers, sphericity, and copolymer distribution/packing, depending on N and x PEO . However, the cluster density, dictated by the minimum in the potential energy of intermolecular interactions, is constant across these structures. Simulations also reveal that clear demarcation between hydrophobic and hydrophilic domains does not exist within spherical/wormlike micelles, vesicles, lamellae, or tori. Significant overlap is predicted between the RDFs of ethylene oxide and butadiene monomers within the molecular assemblies suggesting that the domain boundaries are diffuse. This is consistent with previous studies [25,32,33,36]. A similar feature is observed for water-polymer interfaces as well.
Morphology transitions across topologically simpler structures can be rationalized based on specific SASA measurements. Simulations also suggest that pdfs of segmental stretch of PEO and PB within vesicular assemblies follow a log-normal distribution, i.e., one that maximizes entropy for a random variable with given mean and variance. The universality of this observation needs to be tested by using larger sets of simulation data. Moreover, the mechanisms of topological transformations which include copolymer exchange, fragmentation or coalescence of aggregates, curvature fluctuations, bending/folding of lamellae, and water/copolymer diffusion within micelles need to be investigated. The present study provides a reliable platform to guide future research along such directions.