The self-assembly of vertically aligned block copolymer (BCP) cylinders becomes more difficult when their minority block possesses a lower surface tension than majority block. Because of this specific characteristic, these inverse cylinder systems are highly susceptible to kinetic trapping in thin film geometries. While coarse-grained (CG) molecular dynamics (MD) simulations utilizing hard-wall substrates — where substrate beads are frozen in their initial place but modeled to mimic a random copolymer brush that facilitates vertical orientation — successfully capture the formation of vertical inverse cylinders, their rigid boundaries do not readily accommodate localized relaxation, limiting the statistical likelihood of sampling the fully aligned morphology. To facilitate the sampling of vertically oriented inverse cylinders in MD, we introduced an extended CGMD model featuring a mechanically compliant, harmonically tethered soft substrate. Across an ensemble of independent replicas, the soft boundary allowed interfacial chains to relieve vertical compressive strain, evidenced by a significantly reduced radius of gyration. By enabling localized reorganization, this mechanical buffering shifted the morphological yield toward higher structural fidelity compared to the hard-wall baseline. Ultimately, by substantially increasing both the maximum vertical orientation metric and the peak topological cylinder count, this study demonstrates that incorporating substrate compliance is a critical computational strategy for promoting structural rearrangement and reliably sampling defect-free vertical arrays of inverse cylinders in MD simulations.

Keywords: Block Copolymer, Self-Assembly, Inverse Cylinders, Molecular Dynamics Simulation

The self-assembly of block copolymers (BCPs) is a cornerstone technology for nanoscale device fabrication due to its ability to spontaneously form regular periodic morphologies [1,2]. For linear diblock copolymers, the resulting morphology is dictated by the constituent volume fractions and phase separation thermodynamics. Vertically oriented cylinders are particularly desirable for nanomembranes and templates due to their high aspect ratios and uniform nanoscale pores [3-6].
However, achieving precise vertical orientation requires balancing the interfacial energies at both the free surface and the underlying substrate [7-9]. In conventional systems like poly(styrene)-block-poly(methyl methacrylate) (PS-b-PMMA), the constituent blocks possess inherent surface tension differences under standard conditions. To neutralize the top free surface, thermal annealing is often performed at specific temperatures (~225°C) where the surface tensions of the two blocks match and neutralize the boundary. Simultaneously, to prevent preferential wetting at the bottom boundary, the substrate is typically modified with a neutral random copolymer brush (e.g., PS-r-PMMA). This bottom interfacial selectivity — and ultimately the vertical orientation of the overlying BCP film — is experimentally governed by precisely tuning the chemical composition, specifically the fraction of styrene, within this random copolymer brush.
This precise chemical tuning reveals a stark thermodynamic contrast between standard and inverse morphologies. For regular cylindrical systems—where the higher surface tension PMMA forms the minority cylinder block—the system exhibits a relatively broad processing window for substrate neutrality, requiring a styrene fraction in the brush ranging from 0.59 to 0.72 [10,11]. Conversely, inverse cylinder systems, where the lower surface tension PS forms the minority block, present a severe thermodynamic challenge. Experimentally, the processing window for achieving neutral wetting in these inverse systems narrows drastically, restricting the required styrene fraction to merely 0.55 to 0.57 [10].
Molecular dynamics (MD) simulations are widely used to investigate BCP self-assembly [12-15]. Our previous MD studies utilizing hard-wall substrate models — where substrate beads stay in their initial place but exhibit tunable selectivity toward BCP beads to mimic a random copolymer brush — have successfully captured and reproduced this precise physical phenomenon [16]. In these MD simulations, generic B-b-A BCP are used as representative chains. This substrate selectivity is parameterized into a variable, Γ, which is modeled to equivalently capture the fraction of PMMA in PS-r-PMMA (1-F_styrene). Within this framework, Γ = 0 represents a substrate that exclusively wets the majority B blocks, whereas Γ = 1 exclusively wets the minority A blocks. Our prior hard-wall MD models accurately demonstrated that vertical inverse cylinders do successfully form, and confirmed that their theoretical processing window collapses to a narrow margin of just 0.35≤Γ≤0.4 compared to a wider window for regular cylinders (0.3≤Γ≤0.45), qualitatively mirroring experimental trends [16]. However, because the phase separation of inverse cylinders is highly susceptible to kinetic trapping, not every independent simulation run successfully reaches this ideal vertical state. While the hard-wall model accurately identifies the thermodynamic Γ window, the statistical likelihood of sampling perfectly vertical inverse cylinders across multiple runs remains limited by localized trapping in metastable states (such as parallel or island-and-hole configurations).
In physical experiments, the random copolymer brush layer is an organic, flexible medium possessing inherent mechanical buffering capacities. Therefore, the objective of this study is to introduce a movable soft substrate model to simulate this mechanical compliance within MD simulations. Using an extended coarse-grained framework, harmonic tethering potentials were incorporated to allow the substrate beads to fluctuate thermally and mechanically. Across an ensemble of independent replica simulations, the resulting morphological yields were systematically evaluated by analyzing interfacial chain packing dynamics, a normalized vertical orientation metric, and topological cylinder quantification. Ultimately, this study investigates how mechanical buffering promotes localized structural reorganization, highlighting the potential of substrate compliance as an effective computational strategy to significantly improve the statistical sampling rate of fully aligned inverse cylinders.

2.1 Block Copolymer Chain Model
In this study, MD simulations of BCP thin films were performed using an extended Kremer-Grest model [17]. The BCP chains were modeled as linear B-b-A diblock copolymers, with each chain consisting of a total of 20 coarse-grained beads. To represent a typical cylinder-forming composition, the minority block (A) was assigned 5 beads, while the majority block (B) was assigned 15 beads. This yields a minority block volume fraction of f_A, which naturally forms a cylindrical morphology in the bulk state.
Adjacent beads connected along the polymer backbone interact via a combination of the Finitely Extensible Nonlinear Elastic (FENE) potential and a purely repulsive Lennard-Jones (LJ) potential. The bonded interaction is mathematically expressed as:



The first term represents the FENE spring, which prevents unphysical stretching of the bonds, utilizing standard parameters of spring constant k = 30ε/σ² and maximum extent R₀ = 1.5σ. The second term represents the excluded volume effect, preventing the physical overlap of constituent beads. To ensure this interaction remains purely repulsive, it is truncated at a distance of r_c = 2^1/6 σ≈1.1225σ with a corresponding shift factor of C₁ = 0.25.
For non-bonded interactions between beads, as well as interactions beyond the repulsive cutoff distance, an attractive LJ potential was employed to simulate the cohesive and microphase separation behaviors of the polymer melt:



To accurately capture the long-range attractive interactions critical for BCP self-assembly, the cutoff distance for the non-bonded potential was set to r_c = 2.5σ. The potential was shifted smoothly to zero at the cutoff boundary using the constant C₂.
The interaction parameter ε_ij is the critical factor governing the thermodynamic driving force for phase separation and the distinct surface tensions of the respective blocks. To capture the unique interfacial properties of the inverse cylinder morphology, self-cohesion parameters are set to ε_AA = 0.99ε and ε_BB = 1.01ε. This specific energy configuration dictates that the minority block (A) possesses a lower surface tension compared to the majority block (B). Physically, this drives the minority block to preferentially segregate toward the free surface and substrate interface. Furthermore, the cross-interaction energy between the differing blocks was fixed at ε_AB = 0.5ε. This value primarily drives phase separation.
2.2 Block Copolymer Thin Film Model
To simulate the self-assembly behavior of the BCP thin films, a three-dimensional simulation box was constructed with periodic boundary conditions applied in the x and y directions, and a fixed boundary condition applied in the z direction. The dimensions of the simulation box were set to L_x=57σ, and L_y=90σ. Within this space, 1,800 BCP chains (totaling 36,000 polymer beads) were deposited onto the underlying substrate. The underlying substrate, designed to mimic a substrate modified with random copolymer brush, was modeled as a single layer of hexagonally packed substrate beads (S-type) positioned at the bottom of the simulation box. The lattice spacing between adjacent substrate beads was set to 1.0σ.
To evaluate the impact of substrate mechanical compliance on the self-assembly kinetics, two distinct thermodynamic boundaries were established. For the hard-wall benchmark model, the substrate beads were excluded from the NVT ensemble time integration. This establishes an idealized, immovable boundary that successfully isolates and identifies the fundamental thermodynamic processing windows of the BCP film. To explore the kinetic pathways and mechanical buffering effects present in physical experiments, a soft-substrate model was introduced, wherein the substrate beads were fully incorporated into the NVT ensemble. To prevent macroscopic drift while allowing localized mechanical buffering, each substrate bead was tethered to its initial equilibrium lattice position (r₀) via a harmonic spring potential:



A spring constant of k = 10ε/σ² was applied, providing the exact mechanical flexibility required to mimic an organic polymer brush layer capable of dynamically absorbing localized compressive strains.
To control the interfacial energy between the substrate and the BCP thin film, a substrate selectivity parameter, G, was introduced. The interaction strengths between the substrate beads and the respective polymer blocks (ε_SA and ε_SB) were defined through a linear interpolation function based on Γ:

According to this relationship, Γ = 0 yields a strongly B-preferential surface, while Γ = 1 yields a strongly A-preferential surface. Based on prior computational studies, Γ = 0.375 was adopted as the fixed neutral parameter for the inverse BCP cylinder system, serving as the thermodynamic baseline upon which the mechanical buffering effects were evaluated.
Finally, because the microphase separation of inverse cylinders within this narrow G window is highly stochastic and prone to kinetic trapping, a single trajectory is insufficient to evaluate the true morphological yield. To rigorously quantify the probability of sampling the perfectly vertical state, 10 independent replica simulations were initially executed for both the hard-wall and soft-substrate conditions. Each replica was initialized with a uniquely randomized velocity seed at the start of the thermal annealing process. To mitigate the influence of extreme stochastic outliers on the overall structural trends, the highest and lowest performing replicas from each ensemble were subsequently excluded. This statistical methodology yielded 8 independent, robust trajectories per condition, forming the final ensemble used to calculate all average morphological metrics and statistical variances.
All simulations were carried out using the Large-scale Atomistic/Molecular Massively Parallel Simulator (LAMMPS) and molecular snapshots obtained using the Visual Molecular Dynamics (VMD) software package. To define the excluded volume interactions across the system, the effective bead diameters for identical pairs (σ_ii) and the cross-interaction distances for distinct pairs (σ_ij) were uniformly set to 1.0σ. To drive the self-assembly process, annealing simulations were initially performed at a reduced temperature of T = 1.0ε/k_B. This was followed by a thermal quench, where the system was cooled to and subsequently equilibrated at T = 0.3ε/k_B.
2.3 Morphological Characterization and Orientation Metrics
To systematically evaluate both the molecular-level confor- mational relaxation and the macroscopic structural fidelity of the inverse cylinder thin films, this study employed three distinct analytical metrics: the interfacial radius of gyration, a normalized vertical orientation metric, and a topological cylinder count.
First, to quantify the degree of localized packing strain and molecular relaxation near the substrate boundary, the conformational size of individual BCP chains at the annealing temperature was evaluated by calculating their radius of gyration (R_g):



where N = 20 is the total number of beads per chain, r_i is the position vector of the i-th bead, and r_com is the center-of-mass position vector of that specific chain. To isolate the mechanical buffering effect of the substrate, the R_g values of the chains were computed as a function of their perpendicular distance from the substrate interface. By comparing the R_g profiles of chains immediately adjacent to the boundary against those situated in the bulk film, the degree of spatial confinement and the relief of vertical compressive strain afforded by the soft substrate could be directly quantified.
Second, to strictly quantify the macroscopic degree of vertical alignment across the self-assembled film, a normalized orientation metric (Φ) was calculated after quenching. If the inverse cylinders are perfectly perpendicular to the substrate, their projection onto the xy plane results in completely distinct, non-overlapping A and B block areas. However, if the cylinders tilt, deform, or become trapped in parallel con- figurations, the projected paths of the A and B beads will inevitably overlap. Based on this geometric principle, the ratio of the exclusive projected area — occupied solely by A beads without any B bead overlap — relative to the total cross-sectional area of the simulation box (Area_box = L_xL_y) was determined. This ratio was then normalized by the theoretical volume fraction of the minority block (f_A):



where A_proj and B_proj represent the projected areas of the respective blocks, and the intersection A_proj∩B^c_proj isolates the pure A domain projection. Theoretically, Φ approaches 1.0 for a mathematically perfect vertical cylinder array. However, due to natural molecular-level interfacial roughness and thermal fluctuations in MD simulations, slight perimeter overlaps are unavoidable. Consequently, the highest practical Φ values denoting exceptional vertical alignment of cylinders across the entire simulation box typically range between 0.7 and 0.8.
Finally, because a high Φ value does not inherently guarantee that the resulting morphology consists of discrete, well-separated cylinders (e.g., fused domains could yield similar projection scores). By applying a spatial clustering algorithm to the A-block density map, the total number of discrete, continuous A domains spanning the film thickness was explicitly counted. By simultaneously evaluating the relief of chain-level packing strain (R_g), the geometric straightness of the domains (Φ), and the successful microphase separation into discrete columns (cylinder count), a comprehensive statistical framework was established to compare the morphological yields of the hard-wall and soft-substrate ensembles.

3.1 Interfacial Chain Packing and Conformational Relaxation
The thermodynamic stability of vertically oriented inverse cylinders is heavily dependent on the localized packing behavior of the BCP chains at the substrate interface (Fig. 1 and 2). To evaluate the mechanical buffering effect of the substrate on molecular conformations, the radius of gyration (R_g) of the BCP chains was computed as a function of their perpendicular distance from the substrate boundary.
In the hard-wall baseline model, the immovable lattice establishes a strict, idealized geometric boundary (Fig. 3). Consequently, a well-defined exclusion zone was observed, with the closest BCP chains maintaining a separation distance of approximately 1.0 to 1.5s from the substrate. Because the rigid lattice cannot deform to accommodate localized pressure, the chains adjacent to the hard wall maintain a relatively uniform R_g profile that mirrors the chains in the bulk. This uniform profile reflects the exact, unyielding spatial constraints of the hard-wall system, which effectively isolates the fundamental thermodynamics but provides limited physical free volume for localized conformational relaxation.
Conversely, the soft-wall substrate model fundamentally altered the interfacial packing dynamics by introducing localized mechanical compliance. Because the tethered substrate beads dynamically respond to the pressure exerted by the overlying polymer melt, they allow for subtle geometric rearrangements. As a result, BCP chains were able to approach the substrate to distances below 1.0s. Furthermore, several replica simulations exhibited minor interpenetration of the BCP chains into the substrate layer. This interpenetration closely mimics the physical chain entanglement and wetting behaviors characteristically observed in experimental organic random copolymer brushes.
In addition, the radii of gyration of the BCP chains located immediately adjacent to the soft substrate were found to be smaller than those of the chains situated in the bulk film. This measurable reduction in R_g indicates that the mechanical compliance of the soft substrate enables the adjacent polymer chains to adopt more compact molecular conformations. By providing the degrees of freedom required for the polymer chains to relieve localized vertical compressive strain, the soft substrate effectively absorbs interfacial stress.
3.2 Statistical Analysis of Vertical Orientation and Cylinder Yield
Because the self-assembly of inverse cylinders is a highly stochastic process in MD simulations, analyzing a single trajectory is insufficient to capture the true morphological yield. Following the exclusion of the maximum and minimum stochastic outliers, the 8 robust replicas for both the hard-wall and soft-wall models were analyzed. To effectively visualize this statistical distribution, the average values, standard errors, and individual replica data points were plotted for both the vertical fraction (Φ) and the topological cylinder count.
The structural fidelity of the resulting morphologies revealed distinct statistical trends between the two substrate conditions. For the hard-wall substrates, the vertical orientation distribution was relatively narrow and constrained (Fig. 4). The ensemble yielded an average Φ of approximately 0.36, with the highest-performing individual replica reaching a maximum Φ of 0.46. This tight distribution indicates that while the rigid baseline model successfully captures the thermodynamic formation of inverse cylinders, the rigid boundary geometrically restricts the extent of morphological evolution, making it statistically difficult for the system to fully escape metastable intermediate states (such as localized parallel orientations or terracing).
In contrast, introducing mechanical compliance via the soft-wall substrate noticeably shifted the orientation probability toward higher alignment fidelity. The soft-wall ensemble achieved an elevated average Φ of approximately 0.4 and, more importantly, unlocked a maximum Φ of 0.54. The highest-performing replicas in the soft-wall system corresponded to comparatively well-aligned ordered vertical cylinder arrays. It is important to note that the soft-wall ensemble also exhibited a broader variance, including a highly frustrated trajectory with a Φ of just 0.15. This specific outlier is physically significant: it demonstrates that kinetic trapping is an inherent thermodynamic reality of inverse BCPs, rather than an artifact of the hard wall. However, the significantly elevated maximums confirm that the soft substrate effectively promotes structural reorganization. By relaxing localized packing strain (as evidenced by the R_g reduction), the compliant boundary greatly increases the statistical probability of the system successfully navigating complex energy landscapes to access high-fidelity vertical orientations.
This trend in vertical fraction was directly corroborated by the topological cylinder count (Fig. 5). The hard-wall model yielded a constrained cylinder distribution, peaking at 11 discrete domains for its most successful trajectory. This ceiling reflects the difficulty of fully separating and standing the domains upright under rigid confinement in the entire simulation box. The soft-wall substrate, however, facilitated much denser and more distinct domain formation, yielding up to 18 discrete cylinders in its peak trajectories. The correlated increases in both the maximum vertical fraction and the peak cylinder count demonstrate that incorporating an organic, flexible boundary is a critical factor for computationally reproducing and statistically sampling the optimal vertical alignment observed in physical BCP thin film experiments.

The self-assembly of vertically aligned inverse BCP cylinders in thin films is inherently challenging due to a drastically narrowed processing window for both substrate neutrality and surface tension mismatch. While physical experiments typically mitigate this through thermal annealing to neutralize the free surface, this computational study deliberately maintained a slight surface tension mismatch during the annealing phase. This stringent thermodynamic environment was designed to rigorously investigate the system by exacerbating the kinetic trapping that typically prevents fully aligned morphologies in MD models with hard-wall substrates.
To overcome these sampling limitations, we introduced a mechanically compliant soft-wall substrate. By allowing interfacial chains to relieve localized compressive strain, evidenced by a reduced radius of gyration, this flexible boundary enabled localized reorganization. Consequently, the soft-wall model substantially increased both the maximum vertical orientation fraction and peak topological cylinder count, highlighting the potential impact of substrate compliance as a computational strategy to improve the sampling of inverse cylinders.
To further validate this potential, ongoing studies are evaluating the soft boundary's impact across varying substrate selectivity (Γ) and interaction strength (ε_AB). Future work will determine if this mechanical buffering universally promotes the assembly of both regular and inverse BCP chains, regardless of block volume fraction (f_A).

This work was supported by a 2-Year Research Grant of Pusan National University.