Stabilizing a Double Gyroid Network Phase with 2 nm Feature Size by Blending of Lamellar and Cylindrical Forming Block Oligomers.

Molecular dynamics simulations are used to study binary blends of an AB-type diblock and an AB2-type miktoarm triblock amphiphiles (also known as high-χ block oligomers) consisting of sugar-based (A) and hydrocarbon (B) blocks. In their pure form, the AB diblock and AB2 triblock amphiphiles self-assemble into ordered lamellar (LAM) and cylindrical (CYL) structures, respectively. At intermediate compositions, however, the AB2-rich blend (0.2 ≤ x AB ≤ 0.4) forms a double gyroid (DG) network, whereas perforated lamellae (PL) are observed in the AB-rich blend (0.5 ≤ x AB ≤ 0.8). All of the ordered mesophases present domain pitches under 3 nm, with 1 nm feature sizes for the polar domains. Structural analyses reveal that the nonuniform interfacial curvatures of DG and PL structures are supported by local composition variations of the LAM- and CYL-forming amphiphiles. Self-consistent mean field theory calculations for blends of related AB and AB2 block polymers also show the DG network at intermediate compositions, when A is the minority block, but PL is not stable. This work provides molecular-level insights into how blending of shape-filling molecular architectures enables network phase formation with extremely small feature sizes over a wide composition range.

Introduction

Exploiting molecular self-assembly is an attractive means for the production of nanostructured functional materials and is essential for miniaturization, reaching sub-2 nm feature sizes. A variety of morphologies, such as lamellae (LAM), hexagonally-packed cylinders (CYL), body-centered cubic micellar (BCC), and three-dimensional network (NET) structures, can be formed by self-assembly of amphiphilic block oligomers and related block polymers. Among these possible structures, NET structures with their bicontinuous and interconnected domains are promising for applications as nanoporous membranes,[1−3] ion transport media,[4,5] drug delivery devices,[6,7] and in many other specialty applications.[8−11] However, packing frustration often leads to a narrow composition window for stable NET phases, since the negative Gaussian curvatures of the NET connectors differ profoundly from those found within the struts, whereas the simpler LAM and CYL phases exhibit a constant mean curvature.[12] Various strategies have been introduced to drive the stable NET structures over wider composition and/or temperature ranges, including the design of amphiphiles with shape-filling architectures among ABC triblock polymers,[13] coil–brush/bottlebrush polymers,[14−18] and glycolipids,[19−21] or by addition of solvent to relieve packing frustration.[22,23]

There is also strong interest in miniaturization of the self-assembled NET feature sizes down to sub-2 nm for applications such as nanofiltration and nanopatterning.[24−27] In diblock polymer systems, microphase segregation in the mean-field limit requires the product of the chain length (N) and Flory–Huggins interaction parameter (χ) to be greater than 10.5 (for the LAM morphology), and successful efforts have been made to prepare “high χ–low N” materials with reduced domain sizes.[28−37] For example, NET phases with sub-5 nm domain spacing are observed in glycolipids and T-shaped liquid crystals.[38−42] In our previous simulation and experimental studies, a series of block oligomers containing sugar-based (A) and hydrocarbon-based (B) blocks were found to self-assemble into ordered thermotropic phases including LAM, perforated lamellae (PL), and CYL with domain periods as small as 1.2 nm,[43−45] and good agreement is observed for the domain spacing with experimental data.[44,46]

Blending or mixing block polymers has been extensively studied to create new morphologies that are absent in neat systems.[47−52] In the present study, molecular dynamics (MD) simulations are utilized to investigate blends of two block oligomers containing a hydrophilic tetraol headgroup (with four CHOH repeat units, abbreviated here as H4) covalently connected to one or two hydrophobic alkyl tails (with CH, CH2, and CH3 units, abbreviated here as T): H4T9 and H4T(T8)2 (Figure ). Individually, these two molecules can self-assemble into thermotropic liquid crystalline LAM and CYL structures, respectively.[43,45] We also perform self-consistent mean field theory (SCFT) calculations on blends of AB diblock and AB2 miktoarm triblock polymers (with NA = NB for both architectures) to explore the difference between the entropy-driven self-assembly of flexible block polymer blends and the enthalpy-driven self-assembly of stiff high-χ block oligomer blends.

Figure 1

Chemical structures of the diblock AB (n-tridecan-1,2,3,4-tetraol) and triblock AB2 (5-octyl-tridecan-1,2,3,4-tetraol) amphiphiles.

Results and Discussion

Phase Behavior of AB/AB2 Blends

The H4T9 and H4T(T8)2 blends are probed by MD simulations at TSIM = 460 K, a temperature that is near the middle of the stability windows for the LAM and CYL morphologies formed by the neat compounds, respectively. A sequential simulation workflow is employed to deduce the appropriate system size (i.e., numbers of AB and AB2 amphiphiles) and to set up guiding fields so that the stability of different network phases can be evaluated (see Methods section). Conceptually, the H4T9 and H4T(T8)2 amphiphiles are comprised of either a ditopic or a tritopic central linker bead that is connected to a hydrophilic block comprised of four CHOH repeat units and either one or two lipophilic blocks comprised of four C2H4 repeat units. With identical characteristic dimensions of A and individual B blocks, the volume fractions of the polar block (fA) are 0.213 and 0.338 for neat H4T(T8)2 and H4T9 amphiphiles, respectively. As can be observed from the static structure factors (see Figure ; details of the calculation of the structure are provided in the Supporting Information), the neat H4T(T8)2 system (xAB = 0.0, the mole fraction of the AB amphiphile in AB/AB2 blends) yields an ordered CYL phase with the characteristic structure factor peaks with q/q* ratios of , , and (where q and q* are the magnitude of the scattering wave vector and the location of the first peak, respectively), and the d10-spacing is found to be 2.04 nm. The neat H4T9 system (xAB = 1.0) self-assembles into an ordered LAM phase with q/q* = 1, 2, 3, ... and d = 2.59 nm from the computed structure factor. These morphologies are also evident in snapshots of the simulated systems showing the dividing surfaces between polar and nonpolar regions (see Figure ) and the molecular configurations in ball-and-stick representations (see Figure S1 in the Supporting Information).

Figure 2

(a) Static structure factors S(q) and (b) snapshots (only the minority block volume is shown as a surface mesh) for the equilibrium morphologies observed at TSIM = 460 K for the different compositions (with two orientations shown for the neat AB2 system, xAB = 0.0). The labels denote the composition, ordered morphology, and d-spacing (d10, d211, d10, and d10 for CYL, DG, PL, and LAM, respectively).

When only a small fraction of the other component is added to the predominantly H4T(T8)2 or H4T9 systems (i.e., xAB = 0.1 and 0.9), the morphologies remain unaltered relative to the neat systems (see Figure ), but a small degree of disorder is evident in the snapshot for xAB = 0.1 where “bridges” are present between some of the cylinders. At intermediate compositions, however, stable double gyroid (DG) structures are observed for xAB = 0.2, 0.3, and 0.4 with the d211-spacing increasing from 2.18 to 2.35 nm as more of the AB diblock is added. The characteristic peak position ratios at , , , , , , , and , consistent with the Ia3̅d (Q230) space group symmetry,[13] are found in the structure factors obtained at these three compositions (see Figure a). Furthermore, this morphology assignment is supported by representative slices from the structures at xAB = 0.4 (see Figure e and also Figures S2e and S3e for xAB = 0.2 and 0.3) demonstrating that the internal DG structure is well-preserved. The propensity for these intermediate compositions to self-assemble into a stable DG morphology is also reflected from the disorganized configurations achieved rapidly from random initial structures (see Figures a, S2a, and S3a), where the prevalence of 3-fold connectors rules out other NET structures with 4-fold connectors (e.g., double diamond (DD), Pn3̅m, Q224) and 6-fold connectors (e.g., Plumber’s nightmare (P), also referred to as double primitive, Im3̅m, Q229). Furthermore, single gyroid and double diamond structures generated through system-size tuning and guiding fields for the mixtures at xAB = 0.2 to 0.4 are found to be rather unstable and turn into disordered structures within 10 ns after switching off the guiding fields (see Figure S4).

Figure 3

Snapshots of the minority-region surface meshes at xAB = 0.4 and TSIM = 460 K for (a) a nonequilibrium, disordered bicontinuous structure with system size fine-tuned for a DG structure with 8 unit cells. (b) Interaction sites residing in the subvolume of a DG network that are used as the guiding field to aid the self-assembly process. (c) Equilibrium structure under the applied guiding field. (d) Equilibrium structure reached after the guiding field is removed. (e) Slices with thickness of 15 Å in (111), (110), and (211) directions of the equilibrated DG structure without guiding field.

For xAB = 0.5 to 0.8, self-assembly into perforated lamellae (PL) is observed. The extra peaks in the structure factors that are very close to q* (see Figure a for xAB = 0.5 to 0.7) reflect the average lateral domain periods between neighboring perforations. Notably, the extent of perforations increases with increasing AB2 content (see Figure b) to allow for higher overall interfacial curvature, as the shape-filling architecture of AB2 by itself favors the formation of a CYL phase. When the density of perforations is high, they can arrange into hexagonally ordered patterns. Together with the order transverse to the lamellae, the hexagonal order in the perforations imparts three-dimensional periodicity, just like NET and BCC structures, and requires system size tuning to achieve ordered abab... (HPLab) and abcabc... (HPLabc) stackings. Because of the prohibitively large unit cell dimensions, fine-tuning of system sizes to achieve the HPL morphology[45] is not attempted, and these morphologies are referred to as PL throughout this work. Prior work on diblock polymer melts showed two additional peaks for hexagonally perforated layers (HPL),[53] but our previous simulations for linear BAB triblock oligomers also yield a single additional peak for a well-ordered HPL morphology.[45] Thus, it is possible that the different packing constraints for the stiff oligomers lead to a close match of lamellar spacing and perforation distance resulting in one of the minor peaks being hidden within the main peak.

To test that the PL structure for the mixture at xAB = 0.5 is stable, system-size tuning and guiding field are applied to generate a DG structure. However, we observe (see Figure S5) that, upon removal of the guiding field, the system rapidly evolves over only 10 ns into a disordered bicontinuous structure with flatter (i.e., less spherical) cross-sections of the struts than found for the DG structure at xAB = 0.4. After a much longer trajectory (≈800 ns), a PL-like morphology is recovered for this simulation at xAB = 0.5.

A convolutional neural network for image recognition in Fourier space (FTCNN, see Methods) designed to distinguish various ordered morphologies and trained on ordered and disordered structures of diblock oligomers and of star triblock oligomers (but not including any blends and only synthetic data for network morphologies)[54] also distinguishes among the PL structures (see Figure ). It should be noted that the PL morphology with disordered perforations was not included as a class during the FTCNN training. For the mixtures at xAB = 0.5 and 0.6 with their high density of perforations, the FTCNN classifies individual configurations encountered through the trajectory as HPL (see Figure ), whereas the configurations for xAB = 0.7 and 0.8 with low densities of perforations are classified as LAM. A small fraction of the configurations for the xAB = 0.7 mixture are misclassified as CYL. The configurations from the mixtures at xAB = 0.0 and 0.1 are classified overwhelmingly as CYL with the occasional misclassification as body-centered cubic micelles. Every configuration for the mixtures at xAB = 0.2, 0.3, and 0.4 is classified as DG, and zero probabilities are assigned to the other three types of canonical network phases in the training set, thereby providing further support for the stability of the DG morphology. Similarly, the configurations for the mixtures at xAB = 0.2, 0.3, and 0.4 are with a very high degree of certainty classified as LAM.

Figure 4

Stack plots of the softmax classification probabilities obtained from the Fourier transform convolutional neural network (FTCNN) model for all compositions. For each xAB, 50 simulation frames separated by 2 ns are selected for inference: body-centered cubic micelles (BCC), double diamond (DD), double gyroid (DG), disordered (DIS), hexagonally packed cylinders (CYL), hexagonally perforated lamellar (HPL), lamellar (LAM), plumber’s nightmare (P), and single gyroid (SG).

Domain Spacing and Amphiphile Packing

Although the lengths of the A and B blocks in the H4T9 and H4T(T8)2 amphiphiles are matched, the domain spacing (d = 2π/q*) varies considerably among the ordered structures. A plot of the domain spacing as a function of xAB (see Figure S6) indicates an approximately linear increase with (∂d/∂xAB) = 0.08 nm for xAB ≤ 0.4 covering the CYL and DG morphologies. Similarly, we find an approximately linear increase but with larger (∂d/∂xAB) = 0.11 nm for the PL region (0.5 ≤ xAB ≤ 0.8). The transition from DG to PL morphology (xAB = 0.4 and 0.5, respectively) leads to a discontinuous drop in d by ≈0.1 nm. Interestingly, d = 2.65 nm is found for both the PL structure at xAB = 0.8 and the LAM structure at xAB = 0.9, and d actually appears to decrease as the AB2 molecules are removed to reach the neat AB system. For the CYL morphology, the d10 spacing corresponds to the distance between rows of cylinders and not to the distance between the cylinders. The radial distribution function (RDF) for the terminal oxygen atoms (i.e., primary hydroxyl groups) of the H4T9 and H4T(T8)2 amphiphiles is another metric to probe the spacing between the minority blocks (see Figure S7). The broad fourth peak represents the head-to-head spacing between two domains. There is a clear trend of the fourth peak position shifting toward larger values with increasing xAB. The positions of the fourth peak are approximately 20% larger than the d-spacings, indicating that the primary hydroxyl group is not anchored to the center of the polar domain. As an aside, the lower height of the fourth peak for for xAB = 0.1 reflects the disorder introduced by bridges between the cylinders (see Figures b and S1).

The contour length of an H4T9 amphiphile is ≈1.9 nm, which is considerably larger than d/2 = 1.3 nm and the 1.6 nm deduced from the fourth peak position of the LAM phases. The simulation snapshots (see Figure S1) illustrate a significant degree of interdigitation of the lipophilic tails for the CYL morphology (where even the center of the triangular space between three neighboring cylinders shows relatively high occupancy by tail segments. Only the amphiphiles in the simulation box are shown, and thus, the central part of the tilted box gives the best impression of the packing, whereas periodic images would be needed toward the edges) and also for the PL phases with xAB ≤ 0.7 (where the central region between two polar sheets also shows high occupancy). In contrast, the snapshots for the PL morphology at xAB = 0.8 and for the LAM morphology clearly indicate a region of low density at the midpoint between the leaflets; i.e., there is less interdigitation. The cross-section of the “branched” headgroup with its CHOH repeat units is larger than that of a single nonpolar tail but smaller than that occupied by two tails. Thus, a higher degree of tilt of the nonpolar tails away from the direction normal to the lamellar plane is needed as the fraction of AB2 amphiphiles is reduced, thereby explaining the small decrease in d10 observed for xAB = 0.9 and 1.0.

Since the AB and AB2 amphiphiles share a common tetraol (H4) headgroup, we focus on the packing of the tail groups that occupy a volume fraction of about 75%. Specifically, we surmise that differences in the shape of the T9 and T(T8)2 tail groups reflected in the orientational distributions of neighboring alkyl tails lead to different preferred local interfacial curvature between the A-rich and B-rich regions. Here, we define this interface through the positions of “junction” beads, , placed at the center of the C–C bond connecting the H and T segments. Analysis of the – RDF (see Figure S8) shows remarkably consistent features for all compositions (despite their different morphologies) with the first peak and minimum located at 5.5 and 7 Å, respectively; 7 Å is used as the distance cutoff to define the first coordination sphere of a bead (see Figure b). The orientation of an alkyl tail is then determined by the vector R⃗AB or that points from to each terminal methyl group (see Figure a). The angle θ measures the relative alignment of two alkyl tail vectors belonging to two amphiphiles with neighboring beads; AB–AB, AB–AB2, and AB2–AB2 pairs yield one, two, and four distinct θ angles, respectively.

Figure 5

(a) Illustration of the junction sites J and R⃗AB and vectors for the diblock and triblock amphiphiles. (b) Zoomed-in snapshot of a perforation for the PL structure at xAB = 0.6. The minority (polar) region is shown as red surface mesh, and the nonpolar segments of AB and AB2 amphiphiles are displayed in cyan and yellow, respectively. A spherical region with radius of 7 Å centered at one junction is indicated by the dashed line. (c) Angular distribution functions for the tail group end-to-end vectors from neighboring sites for all compositions (xAB) at TSIM = 460 K.

Angular distribution functions (ADF) for tail group end-to-end vectors belonging to neighboring beads are shown in Figure c (see also Figure S9 for the corresponding heatmap). The ADFs fall into four groups with distinct features. For the LAM phases formed at xAB = 1.0 and 0.9, parallel orientations (θ = 0°) are strongly preferred, angles between 90 and 150° are least favorable, and a small fraction of neighboring vectors are found with antiparallel orientation. Given the tight cutoff distance of 7 Å, the antiparallel orientations arise mostly from AB molecules with their A segments in the same sheet and their tails pointing into different nonpolar regions of the LAM phase. At xAB = 0.8 and 0.7, the PL phases with a low density of perforations show a weaker preference for parallel orientations, but a minimum and a secondary maximum in the ADF are now observed for θ ≈ 40 and 50°, respectively. Marked changes are evident for the transition to PL phases with a high density of perforations (xAB = 0.6 and 0.5); here, the ADFs are approximately symmetric with a preference for angles between 50 and 130°.

The ADFs for the systems with DG morphology (xAB = 0.4 to 0.2) fall in between those of the PL phases with low and high densities of perforations. That is, angles with θ < 90° are clearly preferred over those with θ > 90° (i.e., the ADFs are far from symmetric), angles near 50° are now most preferred, and there is only a weak secondary maximum for parallel orientations. This secondary maximum gradually becomes weaker for the DG phases as xAB decreases and is no longer present for the CYL phase containing only AB2 amphiphiles. Interestingly, as xAB decreases, the population for angles up to 130° increases, whereas the fraction of antiparallel orientations decreases. These subtle changes are likely caused by the cross-section of the polar cores of the cylinders and circular struts possessing a diameter that is slightly larger than the thickness of the polar regions with low curvature forming the nodes and lamellar sheets.[55]Figures S10–S12 show 2-D radial-angular distribution functions (RADFs). When accounting for the increase in volume for radial shells (i.e., distance-normalized), then the strongest orientational preferences are observed for – distances, r, near 5.5 Å, the location of the maximum in the RDFs. The RADFs also show that the antiparallel orientations found for the LAM phases are associated with larger r values near the distance cutoff. The fairly sharp boundaries for | cos θ | > 0.8 observed for the CYL phases can be attributed to the double tails of the AB2 amphiphiles where a pair of neighboring beads results in the calculations of four angles that are constrained by the preference for certain intramolecular conformations.

The angular distribution for the intramolecular CH3––CH3 angle of the AB2 amphiphiles is also analyzed to assess whether changes in self-assembly morphology are accompanied by changes in the molecular shape. Despite the short length of the alkyl tails, the CH3––CH3 distribution is broad with a weak preference for angles near 100° and a shoulder near 20° for xAB ≤ 0.6 (see Figure S13). As the morphology transitions to PL with low density of perforations and LAM (xAB ≥ 0.7), the distribution shifts toward angles between 20° and 80° that allow for better packing in these low-curvature morphologies.

Local Composition Enhancement

For a blend, an important question is whether a network morphology containing nodes and struts with nonconstant Gaussian interfacial curvature can be stabilized by local segregation of the different shape-filling amphiphiles. In the DG morphology, the minority domain consists of two opposite-handed, interpenetrating node-strut networks formed by the polar A blocks, while the majority B blocks fill the remaining space. For the simulation box containing 23 unit cells, we define spherical subvolumes (Figure a) with a radius of rcut around the centers of all 128 nodes (each DG unit cell contains 16 equivalent node positions according to the Ia3̅d symmetry). The local mole fraction within a given subvolume is then calculated as , where NO(rcut) is the number of oxygen atoms within the subvolume belonging to a given type of amphiphile. Some caution is required for the choice of the subvolume. When rcut is too small, it does not cover all of the node region and small NO values lead to large statistical uncertainties. When rcut is set to half of the distance between two neighboring nodes, dnode = a/23/2 (where a is the length of the cubic unit cell), then the spherical subvolumes of adjancent nodes will touch and higher rcut values are not appropriate. For the AB/AB2 blends investigated here, dnode ranges from 1.89 nm (xAB = 0.2) to 2.03 nm (xAB = 0.4). As shown in Figure b, the local composition in the vicinity of the node locations is enhanced in polar groups from the AB amphiphiles compared to the global mole fraction, xAB. The enhancement is larger for the blend with xAB = 0.4 than for xAB = 0.2 and diminishes with increasing rcut. It should be noted here that only a small subregion at the node center is characterized by a flat interface with zero normal curvature, whereas the node edges between the struts possess high normal curvature and negative Gaussian curvature. Nevertheless, there is clear indication that the DG morphology is stabilized by a significant degree of local segregation. Local segregation has also been observed in SCFT calculations for blends of a gyroid-forming diblock polymer and a homopolymer-like diblock polymer, where the latter is enriched at the node.[56]

Figure 6

(a) Snapshot of the DG structure at xAB = 0.3. The red and orange dots highlight oxygen atoms from AB and AB2 amphiphiles, respectively, that are found within spherical subvolumes centered at all the node locations with a radius equivalent to a quarter of the distance between adjacent nodes (rcut = dnode/4). The remaining oxygen atoms are displayed as green dots. Blue lines represent the DG skeleton. The zoomed-in area shows the surface mesh of a node and connecting strut. The simulation box contains eight unit cells. (b) Local composition enhancement for oxygen atoms of AB amphiphiles in nodal subvolumes (xO,ABnode/xAB) as a function of cutoff distance. The shaded areas represent the 95% confidence intervals.

Phase Behavior of AB and AB Block Polymer Blends from SCFT Calculations

A comparison of the phase diagrams for the AB/AB2 amphiphile (from MD simulations) and block polymer (from SCFT calculations) blends is shown in Figure . The most striking difference between these two systems is the absence of a PL phase in the block polymers (with NA = NB for each block; i.e., fA = 0.5 and 0.333 for the diblock and star polymers, respectively). For the block polymer mixture, it is not surprising that the HPL phase is absent when the perforated lamellar structures are comprised of the A blocks. For pure diblock polymer melts, HPL is metastable due to the high packing frustration in the majority component domain;[12] the addition of majority-component homopolymers is predicted to reduce this packing frustration and stabilize the HPL phase.[57,58] For AB2 star polymer melts, SCFT predicts a narrow composition window for which the HPL morphology is stable at high segregation when B is the minority block.[59] When A is the minority block, however, the HPL phase becomes unfavorable for the AB2 miktoarm architecture due to the high entropic penalty of stretching the B blocks to fill the majority-component layers, which makes the HPL phase less stable than for linear diblock polymer melts. Considering the phase behavior of these two neat block polymers, we would not anticipate that mixing of AB and AB2 would give rise to an HPL phase when A is the minority block. In contrast, when B is the minority block, HPL is observed at high segregation in SCFT for the AB and AB2 blend system with a very small composition of the AB block polymers (see Figure S14). We thus conclude that the presence of PL in the block oligomer system is enthalpically driven since the entropic chain stretching penalties, which drive the frustration in the block polymer system and suppress PL, are less important for the stiff oligomers.

Figure 7

Phase diagrams as a function of the mole fraction for (top) AB and AB2 amphiphile mixtures from MD simulations at TSIM = 460 K and (bottom) AB diblock and AB2 star polymer mixtures with the same statistical segment lengths from SCFT calculations at different χN values, where N is the degree of polymerization of the AB chains. CYL, DG, PL, and LAM morphologies are represented by red up triangles, blue down triangles, green squares, and magenta diamonds, respectively. The A-block volume fractions are indicated at the tops of the graphs.

We also observe that the gyroid region for the block polymers is shifted to higher fA values compared to the block oligomer blends. In the block polymer system, the geometry of the microstructure is determined by a competition between the stretching penalties of the A and B blocks, balanced by the unfavorable energy associated with the interfacial area between the A-rich and the B-rich domains. For example, as the block composition of a diblock polymer becomes more asymmetric, a more curved geometry is preferred to minimize the total stretching energy.[60] Moreover, architecturally asymmetric miktoarm star block polymers spontaneously curve toward the minority domains.[61,62] Since the chain stretching entropy is an essential thermodynamic factor for the high-molecular-weight block polymer systems, the gyroid region positioned at relatively higher xAB (and fA) values can be rationalized by the spontaneous curvatures induced by the entropic driving force in the asymmetric polymer architectures, but the directional character of the strong hydrogen-bond interactions and the slight bottlebrush character present for the minority component of the amphiphiles certainly also plays a role. In the mean-field limit, the relative strength of the A–A and B–B interactions (εAA – εBB) does not directly impact the selection of the ordered state because the free energy term reflecting the difference in the self-interaction energies alters the free energies of all morphologies in the same way at a given volume fraction of A-monomers. However, the Flory–Huggins parameter χ expressed in terms of εAA, εBB, and εAB increases with increasing (εAA – εBB)/kBT, where ε is defined as the depth of the potential energy well, so the effect of increasing relative attractive A–A interactions on selection of the ordered state is the same as the effect of increasing χ.

Comparison to Blends of AB and A(B) Diblock Oligomers and Polymers

To assess the importance of the architectural difference (i.e., diblock and miktoarm triblock) of the oligomers and polymers on the phase behavior, we also carried out MD simulations for which the branched H4T(T8)2 amphiphile is replaced by its linear analogue H4T17 (n-heneicosan-1,2,3,4-tetraol abbreviated as A(B2)). For the stiff oligomers, we observe that the LAM morphology persists for linear diblock architectures (see Figure ) even when the volume fraction of the minority block is lowered. As should be expected, increasing xAB results in a decrease of the domain spacing as a longer linear amphiphile is replaced with a shorter one.

Figure 8

(a) Static structure factors S(q) and (b) snapshots (only the minority block volume is shown as a surface mesh) for the equilibrium morphologies observed at TSIM = 460 K for AB/A(B2) diblock blends for xAB ≤ 0.4. The labels denote the composition, ordered morphology, and d10-spacing.

Data for the analogous SCFT calculations considering blends of an AB diblock polymer (fA = 0.5) and a 1.5 times longer A(B2) diblock polymer (fA = 0.333) are shown in Figure S15. Compared to the AB/AB2 blend, the stability window for the LAM morphology for the diblock AB/A(B2) blend is significantly extended down to xAB ≈ 0.2, and the DG window is shifted to the A(B2)-rich region. For the neat A(B2) diblock polymer, the CYL morphology is stable for low χN values, whereas DG is found for high χN values, and the results are consistent with the phase behavior predicted for a neat diblock polymer with fA = 0.333.[56] Different binary blends of linear AB diblock copolymers have been examined previously by SCFT calculations, and Lai and Shi[56] predicted the formation of a double diamond phase for the binary mixture where one species is a homopolymer-like AB diblock. Their calculations illustrate that for diblock polymer mixture systems, the relative block composition plays an essential role in manipulating the stabilization of network morphologies.

Conclusion

MD simulations are utilized to probe the phase behavior for blends of low-molecular-weight AB diblock and AB2 miktoarm triblock amphiphiles containing polar, sugar-based (A) and nonpolar, hydrocarbon (B) segments. In their pure form, the AB and AB2 amphiphiles are capable of self-assembling into ordered LAM and CYL phases, respectively. Blending these two amphiphiles gives rise to large composition windows for stable DG networks and PL phases. Sub-3 nm d-spacings (or the included sphere for the B-region in DG) are achieved for all ordered morphologies. The different ordered morphologies are enabled due to the propensities of AB and AB2 amphiphiles for different B block orientational alignments and local demixing that, for example, yields a slightly enhanced fraction of AB amphiphiles in the nodes of the DG network. SCFT calculations for AB and AB2 block polymer blends also yield a window for a stable DG phase, yet a stable PL region is not observed when A is the minority block. This work illustrates the differences and similarities between diblock/triblock polymer mixtures and stiff/H-bonding oligomer mixtures and points to the need for adjustments of the design principles for tuning shape-filling oligomer/polymer architectures to achieve wide stability windows for network phases that provide domain sizes ranging from a few to many tens of nanometers.

Methods

Molecular Models and Simulation Details

The transferable potentials for phase equilibria united-atom (TraPPE-UA) force field is used to model the H4T9 and H4T(T8)2 amphiphiles (see Supporting Information and Tables S1 and S2 for force field details). Molecular dynamics (MD) simulations are performed in the isobaric–isothermal (NpT) ensemble using the GROMACS 2021.3 software.[63,64] The system sizes range from 1000 to 2000 molecules. The Nosé–Hoover thermostat[65,66] with a time constant τ = 0.4 ps and the Parrinello–Rahman barostat[67] with a time constant τ = 2 ps are used to control temperature and pressure, respectively. The particle-mesh Ewald method[68] is used to compute the electrostatic interactions. The p-LINCS algorithm is used to constrain the O–H bond length,[69] which enables a 2 fs time steps. The initial disordered configurations used as input for the MD simulations are generated from short Monte Carlo (MC) simulations in the canonical ensemble (TSIM = 3000 K). The simulated systems contain random mixtures of stereoisomers for each amphiphile because the coupled-decoupled configurational-bias MC moves[70] applied to a united-atom model do not preserve the tacticity (whereas specific synthesis routes may yield a preference for specific stereoisomers[44]). MD trajectories consisting of at least 300 ns are used to equilibrate the systems at TSIM = 460 K. Thereafter, trajectories of an additional 100 ns for LAM, PL, and CYL phases and of 600 ns for DG phases are performed and used for the structural analysis. Logarithmic plots of the mean-squared displacements obtained during the analysis period indicate that the AB and AB2 molecules reach the diffusive regime beyond 200 ns, and the amphiphiles diffuse on average by more than twice the domain spacing over 600 ns (see Figure S16). The simulation snapshots showing surface meshes and ball-and-stick representations are produced with the Ovito and VMD visualization packages.[71,72]

System Size Tuning and Workflow

For all the xAB values, initial simulations for 1000-molecule systems utilize orthorhombic simulation cells to allow for independent fluctuations in x-, y-, and z- dimensions to ameliorate incommensurability effects.[73,74] However, this strategy is only suitable for morphologies with periodicity in one or two dimensions, such as the LAM and CYL morphologies. In this case, the segregated domains can reach the preferred domain spacing through rotations, lateral expansions, or contractions in the perpendicular direction to the planes of lamellae or the axes of cylinders.

For the 1000-molecule systems at xAB values of 0.2, 0.3, and 0.4, locally segregated but globally disorganized networks are observed. Such defective NET structures are expected, because for three-dimensionally periodic structures, such as BCC and NET, the commensurability issue persists unless the simulation box contains exactly integer multiples of the unit cell which requires knowledge of the lattice parameter (a) and the number of molecules per unit cell (NUC). In addition, the structural orientation must match, or the system is likely to end up with distorted and metastable structures induced by the simulation box, which are unstable in the thermodynamic limit.[75−77] Thus, the system size has to be optimized for these structures. To fine-tune the simulation box so that it contains integer multiples of unit cells with accurate a and NUC that corresponds to a specific periodic structure, the structure factor, S(q), for the disorganized NET-like structures is first calculated, and the appropriate unit cell dimensions for a specific morphology (assuming cubic) is estimated using the position of the broad peak:[78]a = 2πm/q* where m corresponds to the first observable reflection spacing ratio of that morphology (e.g., for DG). The NUC can now be estimated using the average volume per molecule from the initial simulation. For each xAB forming a disorganized NET structure, the simulations were reinitiated with 8NUC molecules estimated for NET phases including double gyroid (DG), double diamond (DD), and single gyroid (SG) structures, within a cubic simulation box of length Lbox = 2a. However, even with a reasonable system size, equilibrating the system into one of the ordered NET structures is still challenging due to the small free energy difference between the disorganized and ordered NET phases.[76] Therefore, we utilize a MD workflow for 3D NET simulation,[79] which employs a guiding field generated by Gaussian interaction sites of different strengths for H and T beads that uniformly reside in the minority and majority domains of the candidate NET structures. For example, Figure b shows the guiding sites for the minority domain of a DG structure containing 8 unit cells. With the presence of the guiding field, the A and B blocks are aligned into two different domains of the ordered gyroid matrices (e.g., Figure c). Then the guiding field is removed, and the stability of each candidate NET morphology is assessed over a prolonged trajectory of ≈1 μs (e.g., see Figure d,e). The values of NUC and other details of the DG unit cells are reported in Tables S3 and S4.

Fourier Transform Convolutional Neural Network

The FTCNN is a three-dimensional convolutional neural network (CNN) that encodes representations of the discrete volumes to extract more discriminative features of each sample. For the training, we use the same point cloud data as in our previous work of a PointNet,[54] a deep neural network for simulation morphology detection using point clouds. These point clouds contain the coordinates representing repeat units of the minority component for nine distinct morphologies: body-centered cubic micelles, double diamond, double gyroid, disordered, hexagonally packed cylinders, hexagonally perforated lamellar, lamellar, plumber’s nightmare, and single gyroid. To ensure the rotational invariance for the classification results, a spatially uniform random rotation was applied on each point cloud, and the resulting point cloud was wrapped into the original box using the periodic boundary condition. All rotated point clouds were then voxelized into 16 × 16 × 16 volumetric occupancy grids (see Figure a). The value in each voxel grid represents the number of points inside the corresponding cubic subregion. Then, the discrete Fourier transform is applied on voxelized point clouds as described by the following equation:where f(x, y, z) denotes the values of voxel grids as a function of the grid indices; N, N, and N are the number of grids along x, y, and z directions; and u, v, and w are the indices in the reciprocal lattice. The fast Fourier transform algorithm[80] was used for efficient computation of F(u, v, w). For consistency, the zero-frequency component is shifted to the center of each spectrum. Examples of the processed voxels and their projections in x, y, and z directions or, equivalently, y – z, x – z, and x – y planes can be seen from Figure b. By such periodic transformation, the translational invariance of the structures under the periodic boundary conditions is ensured. Therefore, no extra preprocessing for the point clouds was required before being voxelized.

Figure 9

(a) Process of generating the voxel training set data (16 × 16 × 16) from a sample point cloud taken from a synthetic double gyroid structure. (b) Projections of the Fourier transformed 3-D occupancy grids (16 × 16 × 16) in x-, y-, and z-directions represented as heatmaps. (c) Schematics of the FTCNN network, including the data preprocessing stage. The 3-D voxels obtained from discrete Fourier transform of the occupancy grids or the corresponding point clouds are fed into two 3-D convolution layers accompanied by max-pooling layers that are followed by a fully connected layer and a softmax layer for the classification task.

The FTCNN is built with six hierarchically stacked layers including 3D convolutional layers, 3D max-pooling layers, and fully connected layers (see Figure c). The input voxels are first fed into two consecutive composite layers, each consisting of one 3D convolutional layer with a kernel size of 3 × 3 × 3 and one max-pooling layer with a 2 × 2 × 2 voxel window to down-sample the feature volumes. These two composite layers contain 16 and 64 filters, respectively. In contrast to the PointNet layers,[54,81,82] which make the prediction results invariant to the order of the input points, the convolution operation used in a CNN is index-ordered. This is because, for the Fourier transformed occupancy grids, the relative positions of the high-intensity grids are crucial to distinguish among different original structures. Subsequently, two fully connected layers with 128 and k neurons are used to output the classification labels. In this case, k = 9, which refers to nine distinct equilibrium morphologies. Leaky rectifed linear units (LReLU) with a slope of 0.01 for negative values were applied as activation functions after each convolutional and fully connected layer. Batch normalization and dropout with a 0.85 keep ratio are applied before the last fully connected layer. For each morphology, 3000 voxels are obtained by preprocessing the point clouds from the MD simulation trajectories and the synthetic network structures[54] and are split into training and test data in a 4:1 ratio. The model is trained over 1000 epochs with a batch size of 128 voxels. The Adam optimizer[83] is used with a learning rate of 0.001. The exponential decay rates for the first and the second moments are chosen as 0.9 and 0.999, respectively. The model is implemented in PyTorch.

SCFT Calculations

Canonical SCFT calculations for incompressible melts of AB and AB2 block polymer blends (NA = NB, i.e., fA = 0.5 for xAB = 1 and fA = 0.33 for xAB = 0) are performed using the CPU version of the open-source C++ PSCF software.[84] In the SCFT calculations, both AB and AB2 block polymers are modeled as flexible Gaussian chains with equal degrees of polymerization for each A and B block, NA = NB, and equal statistical segment lengths. The modified diffusion equation in SCFT[85] is numerically solved using the pseudospectral method[86] with periodic boundary conditions, integration steps of ds = 0.01, and spatial grid sizes of 403. Initial guesses for LAM, DG, and CYL are all readily available from prior work.[85] For rapid convergence of the SCFT algorithm, we generate the initial guess chemical potential fields for the perforated lamellae by a level set method,[85] which calculates the initial input structure by the level surface of the first nonzero symmetry-adapted basis function of the specific space group. The fields are then iteratively updated using an Anderson-mixing scheme, optimizing the variable unit-cell dimensions with stress relaxation.[87] The iterations are stopped when the errors in the self-consistent field equations drop below the specified tolerance of ε = 10–5. After obtaining the self-consistent mean field solutions, the free energies of the LAM, DG, CYL, and HPL phase (with ab and abc stacking) are calculated and compared to determine the equilibrium states.