*3.2. Descriptor*

As mentioned earlier, equally important to the simulation itself is the analysis of the results, which can be "on the fly" or in a post-simulation step. Monte Carlo simulations, such as the ones presented here, provide no dynamical information, so the emphasis is placed on the study of the local and global structure, organization, topology, and phase behavior. Over recent decades, conceptually different approaches have led to the development and application of descriptors and analyzers of local structure in computer-generated configurations or of digitally processed experimental samples [81–93].

Here, since we are particularly interested in studying entropy- and energy-driven crystallization of polymer-based systems under extreme conditions, we propose a modeling scheme where the MC-based simulator is connected to a descriptor of the local structure ("descriptor") in the form of the CCE norm [60,61]. The version adopted in Simu-D is very similar to the one we presented very recently, so the concept, methodology, and technical implementation, reported in detail in [60], are all also applicable to the present context. Thus, in the continuation, we will provide a brief description on the main aspects of the CCE norm descriptor and the new features, as implemented in Simu-D. Given an atomic or particulate system in two or three dimensions, the CCE norm proceeds by comparing the local environment around each site with the ideal ones of specific reference crystals.

The main concept behind the CCE norm descriptor is that each ideal crystal is uniquely identified by a set of symmetry operations (elements of its point group) [94–97]. The identification of the totality of these crystallographic operations, or of an equally discriminating subset of them, and their application to the nearest neighbors of an atom or particle is key in the implementation of the CCE algorithm.

As explained in detail in Ref. [60], the CCE norm is defined with respect to a specific crystal *X*. Once the reference crystal *X* is selected, the point group is identified along with the generating symmetry elements. Given a site (atom or particle), *i*, the *Nvor*(*i*) nearest neighbors are identified through Voronoi tessellation. The *Nvor*(*i*) population is then compared against the coordination number of the reference crystal *X*, *Ncoord*(*X*). The latter is, for example, equal to 12 for the face-centered cubic (FCC) and hexagonal close-packed (HCP) crystals in 3-D and 6 for the triangular (TRI) crystal in 2-D. If the coordination number is larger than the number of nearest neighbors, *Ncoord*(*X*) > *Nvor*(*i*), a penalty function is applied [60]. In the opposite case, *Ncoord*(*X*) < *Nvor*(*i*), only the *Ncoord*(*X*) closest neighbors are kept for the successive CCE-based analysis. The characteristic crystallographic element(s) is(are) identified and the corresponding actions for each one of them are applied to the coordinates of the neighbor atoms relative to the given site. For example the HCP crystal, with *Ncoord*(HCP) = 12, has one geometric symmetry element in the form of a sixfold rotoinversion axis, while the body centered cubic (BCC), with *Ncoord*(BCC) = 8, has five such elements: Four three-fold roto-inversion axes and one inversion center.

One important point in the CCE norm analysis in a 3-D system is that the orientation of each symmetry axis, or at least of a sub-set of them, is not known a priory. Accordingly, we scan the orientation space SO(3) around the given site with a mesh of discretization width *ϕstep*, which is the same for the azimuthal and polar angles.

For a given orientation, the actions of the symmetry element are executed. This procedure is then repeated over all symmetry elements. The goal of these crystallographic operations is to map the real coordinate system (given site *i* and *Ncoord*(*X*) neighbors) into the ideal one of the reference crystal *X*. Once this is completed, the algorithm proceeds to the next point of the discrete mesh until the whole orientation space is examined. This mapping allows to simultaneously quantify the orientational and radial similarity of the given local environment with respect to the ideal one of crystal *X*. This is realized through the calculation of a norm (see Equation (2) of [60]). The CCE-based norm for the given atom *i* with respect to reference crystal *X*, *ε<sup>X</sup> i* , is the one that corresponds to the global minimum of the norms as calculated over all possible orientations of all symmetry elements (axes). The same process is repeated over all particles or atoms of the systems and all reference crystals. Currently, the CCE descriptor, as implemented in Simu-D, includes the following crystals: Face-centered cubic (FCC), hexagonal close-packed (HCP), bodycentered cubic (BCC), and hexagonal (HEX) for 3-D systems and honeycomb (HON), square (SQU), and triangular (TRI) for 2-D systems. Additionally, the local structure can be quantified with respect to fivefold (FIV) and pentagonal (PEN) local symmetries, in 3-D and 2-D, respectively. The lower the value of the CCE norm, the higher the similarity of the local environment to the reference crystal. A site is labelled *X*-type when its minimum CCE norm is lower than a critical threshold, *<sup>ε</sup>thres*, i.e., *ε<sup>X</sup> i* ≤ *<sup>ε</sup>thres*. By construction, as the characteristic crystallographic elements and operations constitute a distinctive feature for a crystal, the CCE norm is highly discriminatory, so that when the CCE norm with respect to crystal *X* is low, the corresponding norm for other crystal types is high.

An extensive analysis of the underlying concept, the minimum distinguishable set of symmetry elements and corresponding actions for each reference crystal, the algorithmic implementation on the CCE-norm descriptor, the required computational time, and the optimal selection of parameters are all discussed in detail in Ref. [60]. The Simu-D version contains certain additional features. As an option, the "on-the-fly" implementation allows, during the scanning of the spherical space, for the CCE analysis to stop when the norm is found to be lower than the pre-set threshold and pass to the next atom so as to expedite the process and provide a preliminary structural identification. Additionally, the CCE descriptor further identifies the clusters of all atoms that bear the same similarity. For example, it detects clusters of ordered sites, calculates their size (in number of atoms), as well as their shape. The cluster-based analysis functions with the same proximity criterion as the cluster identification used for the moves of the simulator component. An additional condition for the cluster identification is that it should contain sites that have all the same similarity (with respect to a single crystal type *X* or to a pair of them ( *X* or *Y*)). Finally, the CCE descriptor provides information on the shape, size, and statistics of the Voronoi cells, as extracted from the Voronoi tessellation.

A table with a summary of the main variables used by the Simu-D suite along with a brief description can be found in Appendix A (Table A1).

#### **4. Simu-D: Applications**

In this section, we briefly present polymer-based systems that can be simulated and successively analyzed with the Simu-D suite. Emphasis is placed on the simulation of systems under extreme conditions: These can range from a very high concentration (packing density), extreme confinement, or presence of nanofillers with dimensions significantly larger than the monomer size or any combination of the above. In the case of entropy- or energy-driven phase transitions, the corresponding analysis takes place through the CCE descriptor on the frames and trajectories generated by the simulator component.

The main point to be highlighted is that the Simu-D suite is built in a modularbased approach with the goals of general applicability and simplicity. So, all examples in the continuation have been or can be simulated and successively analyzed without any modification of the code being required from the end user. Here, it is not our intention to exhaustively analyze the physical behavior of each reported case but rather to provide evidence that such systems can be modeled and then characterized by the Simu-D software.

#### *4.1. Packing Efficiency of Semi-Flexible Athermal Polymers (3-D)*

How atoms, particles, or macroscopic objects are packed in the most efficient way is a topic of paramount importance in various fields and applications. Ordered packings of nonoverlapping spheres in 3-D have an upper limit in the volume fraction, which corresponds to the one reached by the HCP or FCC crystals [98,99]. For disordered systems of the same entities, the corresponding packing density is globally accepted to be approximately 10–12% lower and corresponds to the Random Close-Packed (RCP) limit [100,101] or its equivalent Maximally Random Jammed (MRJ) state [102]. In the very first application of the MC-based code that served as the initial seed for the Simu-D suite, it was demonstrated that freely jointed chains of tangent hard spheres can be packed as efficiently as monomeric analogs [103]. However, the corresponding state of semi-flexible polymers or even of freely rotating chains is still a subject for investigation [104,105]. To this end, we used the simulator component to generate and successively equilibrate random packings of semi-flexible chains with a varied equilibrium angle, degree of stiffness, as quantified by the spring constant in Equation (2), average chain length, and volume fraction. Exploring the combined effect of the physical variables stated previously requires the conduction of numerous simulations starting from dilute systems all the way up to the RCP/MRJ limit. The range of the latter is expected to be a function of the rigidity of the chain and thus depend strongly on the bending constant and equilibrium angle [104].

Figure 4 shows bulk system configurations for semi-flexible hard-sphere chains with average length *N* = 100 ( *Nat* = 4800) with an equilibrium (supplement) angle of *θ*0 = 120◦ at different packing densities of *ϕ* = 0.001, 0.1, and 0.60.

Using the Simu-D generation-equilibration modules, structures of semi-flexible athermal polymers can be simulated at very high densities, which are comparable to the densest ones observed for fully flexible (freely jointed) polymers [66,103] or monomeric counterparts [102]. The acceptance rate of the employed local and chain-connectivity-altering move as a function of packing density for the 48-chain *N* = 100 system with *θ*0 = 108◦ is shown in Figure 5 and is reminiscent of the one obtained for freely jointed chains [59]. As expected, the acceptance rate of local moves is significantly reduced as the system reaches progressively higher concentrations. Towards this, the configurational bias scheme aids in reducing this drop. The reduction for semi-flexible chains is especially apparent for the two variants of the reptation move. In sharp contrast, the acceptance rate of chainconnectivity-altering moves shows opposite trends: The higher the concentration, the higher the acceptance rate. Especially for the simplified End-Bridging at low volume fractions, acceptance is very small. This is expected as in such a dilute system there are very few or no pairs of atoms that can trigger the move. As the concentration increases, the population of such pairs also increases because chains start to feel each other and the contact network around each site becomes richer. In parallel, none of the CCAMs, as incorporated in Simu-D, requires the displacement of any atoms. Thus, their performance is enhanced at very high packing densities, and especially near the MRJ state.

According to the RCP/MRJ definition, the maximum-density state should correspond to the densest structures, which are characterized by the maximum randomness or, equivalently, the minimum order [102]. The concept of rattlers [102] and flippers [103] can be invoked to quantify the fraction of individual sites and groups of them, which are able to perform movements in their local vicinity for monomeric and polymeric packings, respectively. In both cases, it is well demonstrated that the flipper/rattler population diminishes as we approach the MRJ state. Alternatively, one could attempt to quantify the lack of order in the system through the proper definition of corresponding parameters. Towards this, we employ the CCE norm (descriptor module of Simu-D) to calculate the similarity of the local environment around each monomer site to the close-packed (HCP and FCC) crystals, which are the dominant ones in the crystallization of hard sphere packings at high volume fractions [72,106]. The absence of such crystals should correspond to a highly disordered but densely packed medium near or at the RCP/MRJ state.

**Figure 4.** Bulk system configurations of semi-flexible chains of tangent hard spheres of uniform size with average length of *N* = 100 and an equilibrium angle of *θ*0 = 120◦ at progressively higher volume fractions, *ϕ*: (**top left**) 0.001, (**top right**) 0.10 and (**bottom left**) 0.60. (**bottom right**): All three system configurations shown together allowing for a visual comparison of their dimensions. Monomers are colored according to their parent chain. Sphere monomers are shown with coordinates of their centers subjected to periodic boundary conditions. Image created with VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

Figure 6 shows the final configuration for the 48-chain *N* = 100 system with an equilibrium bending angle of *θ*0 = 120◦ at a density of approximately 0.64, which corresponds to the range of RCP/MRJ, as established for monomers and freely jointed chains. The left panel shows monomers colored according to the parent molecule, while the right one uses a coloring scheme according to the values of the CCE norm. More precisely, blue and red correspond to sites with HCP (*ε*HCP *i* ≤ *εthres* = 0.245) and FCC (*ε*FCC *i* ≤ *εthres* = 0.245) similarity, respectively, while green is used to represent FIV-like (*ε*FIV*i* ≤ *εthres* = 0.245) sites. All remaining amorphous (AMO) ones, which constitute most of the system, are shown in yellow with reduced dimensions in a 2:5 scale for clarity purposes. More accurately, amorphous (AMO) designates sites that show no similarity to any of the reference 3-D crystal (HCP, FCC, HEX, BCC) or local symmetry (FIV). This does not exclude the possibility that a specific site showing similarity to another "unknown" crystal not included in the reference list. Still, as mentioned earlier and given the very high concentration of the

generated athermal packings, the presence of non-compact crystals can be excluded. This is evident as no traces of BCC or HEX crystals are detected in any of the nearly jammed polymer configurations, such as the ones visualized in Figure 6.

**Figure 5.** Percentage of acceptance as a function of packing density for the local and chainconnectivity-altering moves employed in the MC simulation of 48 chains of *N* = 100, *θ*0 = 108◦ in the bulk. Vertical dash lines denote the end of the regime where *<sup>n</sup>*dis trial configurations are attempted per local MC move.

**Figure 6.** Jammed packing of semi-flexible chains of tangent hard spheres of uniform size with average length of *N* = 100 and an equilibrium angle of *θ*0 = 120◦at a packing density of *ϕ* = 0.637. (**Left panel**): Monomers are colored according to the parent chain; (**Right panel**): Monomers are colored according to the lowest value of the CCE norm. Blue, red, and green denote HCP, FCC, and FIV similarity, respectively. Amorphous (AMO) ones are colored yellow with reduced dimensions for clarity. Sphere monomers are shown with coordinates of their centers being subjected to periodic boundary conditions. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

Visual inspection of the jammed configuration in Figure 6 suggests a predominantly amorphous structure with few ordered HCP and FCC sites randomly distributed along the whole volume of the simulation cell. In fact, one can observe that the population of FIV-like sites is higher than that of the close packed crystal ones. Moving on to quantitative analysis based on the CCE order parameter [60] for the specific structure shown in Figure 6, the HCP, FCC, and FIV fractions are 0.022, 0.021, and 0.053, respectively, further demonstrating the predominance of disorder. The random character of the maximally jammed state for semi-flexible chains is in perfect qualitative agreemen<sup>t</sup> with the one exhibited by freely jointed analogs; the same can be stated for the growth of fivefold local symmetry with an increasing concentration as observed for monomeric counterparts [107,108] as well as for freely jointed chains [70].

#### *4.2. Entropy-Driven Crystallization of Semi-Flexible Athermal Polymers*

The presence of fivefolds in a random particulate packing acts as an inhibitor to crystallization [107,108], especially as the concentration approaches that of the jamming state. However, after a critical volume fraction (melting point) is exceeded, and if the observation (here simulation) time is sufficiently long, hard sphere packings crystallize. Similar phase behavior is observed for freely jointed chains, albeit with differences in the critical packing density and the morphology of the established crystals, both depending strongly on the gaps between bonded atoms [69]. Using the Simu-D suite we can extend the simulations to capture the effect of chain stiffness. As an example, Figure 7 shows the phase transition as first simulated and then identified by the CCE-based analysis for the 100-chain *N* = 12 system at *ϕ* = 0.58 with an equilibrium angle of *θ*0 = 90◦. In the left panel of Figure 7, the initial configuration is presented, as produced through the generation module, while in the right panel, the final one after the execution of 3 × 10<sup>11</sup> MC steps is shown. In both system states, monomers are colored according to the value of the CCE norm. It can be unmistakably concluded that the specific system shows crystallization, with the final stable configuration being of defect-ridden, fivefold-free, alternating HCP and FCC layers. Given that the hard-sphere chain system is athermal, such a phase transition is dictated solely by the increase in the total entropy of the system through a mechanism similar to the one observed in freely jointed chains where the local environment around each ordered site becomes more symmetric in the crystal phase [71,72,109].

#### *4.3. Phase Behavior of Athermal Blends (Polymers and Monomers)*

Through the incorporation of MC moves involving individual monomers and polymer chains (IdEx1, IdEx2, and IdEx3), Simu-D software can tackle blends of chains and monomers of varied relative fractions and different interactions between species. Example cases include a high-density athermal blend of polymers and monomers with a varied number of chains as can be seen in the panels of Figure 8. The system consists of 54,000 interacting sites at a packing density of *ϕ* = 0.56 and an average chain length of *N* = 1000. The polymer relative concentration ranges from 0 (0 chains and 54,000 monomers), 0.185 (10 chains and 44,000 monomers), 0.741 (40 chains and 14,000 monomers) to 1 (54 chains and 0 monomers). The objective here is to study how the relative concentration of the different entities (single versus chain monomers) could affect crystallization. This is motivated by the fact that the selected volume fraction is below and above the melting point of strictly tangent chains and individual monomers, respectively.

#### *4.4. Energy-Driven Cluster and Crystal Formation of Attractive Chains*

Up to this point, all systems studied were athermal with all interactions being of the hard sphere type. In the following case, we employ the square well (SW) attractive potential. Additionally, we carry out the simulations at a constant volume (NVT) for chain systems and at a constant (and high) pressure (NPT) for monomeric ones. In both cases, the starting configuration corresponds to a low-density (*ϕ* = 0.05) hard sphere system where we activate the SW interactions between all sites. Given the attraction felt between the monomers, clusters start to form, which, depending on the applied intensity and range of interactions, could further lead to ordered morphologies or amorphous glasses [63]. From the technical point of view in NVT simulations, one should activate collective, clusterrelated moves since, especially at a low concentration and a high attraction intensity, small and isolated clusters could be created, disallowing further inter-cluster aggregation and the eventual formation of a single entity. The phase diagram of SW chains can be surprisingly rich with different crystals and amorphous morphologies resulting as a function of the attraction range. As an example, Figure 9 hosts the final system configuration obtained from NVT simulations on chains (*<sup>ε</sup>SW* = 1.2, *σ*2 = 1.15, *N* = 12) and NPT simulations on monomers (*<sup>ε</sup>SW* = 2.1, *σ*2 = 1.15), both having *Nat* = 1200 interacting sites. For the given pairs of intensity and range of attraction, the established morphologies consist of HCP and FCC crystallites with random stacking directions. In the case of the chain cluster (left panel of Figure 9), the presence of fivefold sites in the form of twin defects is particularly evident in the meeting points of the HCP and FCC planes.

#### *4.5. Polymers under Confinement*

In a recent publication [75], we demonstrated the ability of the early version of the Simu-D suite to create polymer configurations under tube-like and plate-like confinement. Extreme conditions correspond to the case where the distance between the confining agents/surfaces is approximately equal to the diameter of the monomers. For example, for plates, this extreme condition corresponds effectively to a 2-D polymer system. The latest implementation of Simu-D allows for flexibility in the applied geometry of confinement departing potentially from orthogonal cells. Figures 10 and 11 show system configurations of freely jointed chains of tangent hard spheres (*N* = 12, *Nch* = 60) being confined in cylindrical (closed ends) and spherical geometries, respectively.

**Figure 7.** Snapshots of the semi-flexible *N* = 12 system (*θ*0 = 90◦) at *ϕ* = 0.58. (**Left panel**): Initial configuration as produced by the generator module of Simu-D. (**Right panel**): Final configuration of the simulation after the execution of 3 × 10<sup>11</sup> MC steps of the simulator module. Monomers are colored according to the lowest value of the CCE norm (descriptor module). Blue, red, and green colors denote HCP, FCC, and FIV similarity, respectively. Amorphous (AMO) ones are colored yellow with reduced dimensions for clarity. Spherical monomers are shown with coordinates of their centers subjected to periodic boundary conditions. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

**Figure 8.** Bulk systems of 54,000 interacting hard spheres with varied relative concentrations of polymer content with chains having an average length of *N* = 1000 at *ϕ* = 0.58. The polymer relative concentration ranges from 0, 0.185, 0.741 to 1 (from left to right). In the pure polymer configuration (rightmost panel), sites are colored according to the parent chain. In all other cases, single and chain monomers are shown in red and blue, respectively. For chain concentrations of 0.185 and 0.741, sphere dimensions of dominant species are shown in a 2:5 scale for clarity. Image created with the VMD visualization software [54].

**Figure 9.** Final configurations of systems whose sites interact with the attractive square well potential. (**Left panel**) NVT simulations on chains (*<sup>ε</sup>*SW = 1.2, *σ*2 = 1.15, *N* = 12, *N*ch = 100, *ϕ* = 0.05); (**Right panel**) NPT simulations on monomers (*<sup>ε</sup>*SW = 2.1, *σ*2 = 1.15, *N*at = 1200). Sites are colored according to the CCE norm: Blue, red, green, cyan, and purple correspond to sites with HCP, FCC, FIV, BCC, and HEX similarity, respectively. Amorphous (AMO) sites are colored yellow and shown with reduced dimensions for clarity. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

**Figure 10.** System snapshots of linear, fully flexible chains ( *N*ch = 60, *N* = 12, *ϕ* = 0.40) under cylindrical confinement with closed ends with a length-to-diameter ratio equal to 2 (**left panel**) and 10 (**right panel**). Monomers are colored according to the parent chain. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

**Figure 11.** System snapshots of linear, fully flexible chains (*N*ch = 60, *N* = 12) under spherical confinement at a packing density of (from left to right): *ϕ* = 0.30, 0.40 and 0.50. Monomers are colored according to the parent chain. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

In the cylindrical confinement, the cell (length to diameter) aspect ratio increases from 2 (left panel) to 10 (right panel) while the volume fraction remains the same (*ϕ* = 0.40). In the spherical one, the packing density changes from 0.30 (leftmost panel) up to 0.50 (rightmost panel). Reaching such densities allows the investigation of crystal nucleation and the growth of chain systems and eventually the comparison with monomeric analogs, as recently reported in [110,111].

#### *4.6. Polymer Nanocomposites*

The latest implementation of Simu-D allows for the simulation of polymer-based nanocomposites (PNCs). The nanofillers can exist as compact objects of cylindrical or spherical forms and at various concentrations and interactions with the chain monomers. Figure 12 shows examples of polymer nanocomposites where all entities interact through the hard-core potential. A single nanofiller is inserted, which is a sphere of size *dsph* (in units of monomer diameter, *σ*). The nanoparticle is positioned so that the coordinates of its center coincide with the center of the simulation cell. The left panel shows a PNC with *dsph* = 5 at *ϕ* = 9.9 × 10−<sup>3</sup> while *dsph* = 20 at *ϕ* = 0.29 is presented in the right panel. Taking into account the presence of the nanosphere, the effective packing densities are *ϕeff* = 0.01 and 0.55 for the systems in the left and right panels of Figure 12, respectively. The minimal difference for the former case is due to the small nanoparticle size (*dsph* = 5) compared to the large volume of the simulation cell, while in the latter case, the nanoparticle, due to its massive size (*dsph* = 20), has a profound effect on the reduction of the available volume.

Another example from MC simulations on PNCs is provided in Figure 13. This time the nanofiller takes the form of a single, infinitely long cylinder whose direction coincides with the direction of one of the axes of the simulation cell. The diameter of the cylinder is *dcyl* = 5 and is dispersed in a polymer matrix (*N* = 100, *Nch* = 48), which consists of (left panel) freely jointed and (right panel) semi-flexible, rod-like (*θ*0 = 0◦) chains.

#### *4.7. Comparison with Independent Algorithms*

A relevant task is to compare the results of any simulation suite with existing ones, preferably using different simulation methods. Here we should mention that with respect to jamming, our Simu-D produces very dense and nearly jammed random packings of hard spheres (chains or monomers) with volume fractions very close to the ones reported in the literature from independent sources on the RCP/MRJ state: *ϕ*MRJ ≈ 0.64–0.65, with the exact value and the salient characteristics being very dependent on the employed protocol [112–114].

**Figure 12.** System snapshots of polymer nanocomposite (*N* = 100, *N*ch = 48) at different effective packing density, *ϕ*eff. The nanofiller, shown in red, corresponds to a single, impenetrable sphere with diameter *<sup>d</sup>*sp<sup>h</sup> (in units of *σ*) whose center is located at the center of the simulation cell: (**left**) *ϕ*eff = 0.01, *<sup>d</sup>*sp<sup>h</sup> = 5 and (**right**) *ϕ*eff = 0.55, *<sup>d</sup>*sp<sup>h</sup> = 20. Monomers are colored according to the parent chain and are shown as semitransparent spheres for clarity. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

**Figure 13.** System snapshots of polymer nanocomposite (*N* = 100, *N*ch = 48, *ϕ*eff = 0.10). The nanofiller, shown in blue, corresponds to a single, impenetrable cylinder with diameter *<sup>d</sup>*cy<sup>l</sup> = 5 (in units of σ) and infinite length. The cylinder is oriented along the direction of one of the cell axes. (**left**) Freely jointed chains and (**right**) semi-flexible, rod-like chains (*θ*0 = 0◦). Monomers are colored according to the parent chain and are shown as semitransparent spheres for clarity. Image created with the VMD visualization software [54]. Figure panels are also available as interactive, 3-D images in Supplementary Materials.

In parallel, the melting point of monomeric hard spheres, as determined by simulations conducted with the present protocol, coincides with the well-established value available in the literature [115]. With respect to the crystallization of athermal polymers and the effect of bond gaps (or bond tangency), the results based on the application of the Simu-D suite [69] are in excellent agreemen<sup>t</sup> with independent studies using event-driven Molecular Dynamics (MD) simulations [116].

Furthermore, as an additional testbed, we use the crystallization of monomeric hard spheres. Towards this, we use the same amorphous system configuration, composed of 54,000 hard spheres at a volume fraction of *ϕ* = 0.56. Using this initial structure, we embark on two different kinds of simulation: (i) Stochastic MC using the Simu-D suite and (ii) event-driven, collision-based MD. Given that the reference system consists of hard spheres at an elevated concentration, total crystallinity can be considered as the sum of the fractions of sites with HCP and FCC-like similarity, with FIV local symmetry acting as a structural competitor to compact crystals. In both cases, the local structure is quantified through the CCE metric. Even if the two simulation methods are distinctly different, one (MD) based on collision-based dynamics the other (MC) being completely stochastic, the corresponding trends on the evolution of crystallinity, as seen in Figure 14, are strikingly similar, not only in qualitative but also in quantitative terms.

**Figure 14.** Evolution of crystallinity, *τ*c, and of the fraction of sites with fivefold (FIV) local symmetry, *S*FIV, as a function of MC steps (**left panel**) and MD collisions (**right panel**). Both the MC simulation, performed through the Simu-D suite, and the independent, event-driven MD simulation, are conducted on the same random initial configuration of 54,000 monomeric hard spheres of uniform size at a packing density of *ϕ* = 0.56. Total crystallinity is calculated here as the sum of fractions of sites with FCC and HCP character, as quantified by the CCE norm descriptor.

## **5. Conclusions**

We present the latest implementation of Simu-D, a simulator-descriptor suite used to model and successively analyze/describe polymer-based systems under extreme conditions of concentration (packing density), confinement, and nanofiller content. The simulator part is based on Monte Carlo (MC) algorithms, including localized, chain-connectivityaltering, identity-exchange, and cluster moves in various statistical ensembles. The descriptor is based on the characteristic crystallographic element (CCE) norm, which is a metric to gauge the local structure by comparing it with reference crystals in two and three dimensions. The suite has a modular approach, allowing the addition of features, and is built considering efficiency, general applicability, and ease of use. Monomers/atoms/particles are presented as spheres, which interact through standard bonded and pairwise, nonbonded terms.

We have provided examples ranging from applications on bulk, pure macromolecular systems, of blends with monomers, under various conditions of confinement to polymerbased nanocomposites. Through such simulations one can study general phase transitions,

packing ability, and local and global structure as a function of the aforementioned parameters. In the examples provided, emphasis is placed on the simplified hard-sphere and square-well models, but chemically realistic systems can be simulated as well.

Presently, Simu-D is further expanded to tackle more complex systems including the simulation of terminally grafted nanoparticles anchored on polymer chains as seen in Figure 15, and polymer adsorption on flat or nanostructured surfaces.

**Figure 15.** Terminally grafted nanoparticles on polymer chains at a volume fraction of *ϕ* = 0.50 as simulated through the Simu-D suite. Each nanoparticle, shown in red and in semitransparent format, has a size of *d*nano = 8 and is anchored to a single polymer chain. Macromolecules are represented as freely jointed chains of tangent hard spheres with an average length of *N* = 100.

Our simulator-descriptor suite is rather lacking with respect to the available potentials and interactions, especially compared to latest simulators [24]. However, as explained earlier, our intention is to simulate general, but still coarse-grained, representations of atomic and particulate systems with emphasis placed on extreme conditions, such as jamming, confinement, anchoring, presence of nanoparticles or all possible combinations of the above.

Chemical reactions can also be studied by assigning reactant and product types and implementing identity-change algorithms with the MC simulations being cast in the proper reactive ensemble [117].

In a coarse-to-fine approach, the present MC suite could be further benefited by algorithms that allow the simulation of chemically complex, all-atom systems through reversible, adaptive, or bijective mapping [34–37,118]. Furthermore, the suite should be compatible with independent and efficient MC algorithms such as the event-chain ones [40,41], parallel techniques [119,120], but also with different analyzers. The latter can be in the form of geometric [121,122], stochastic [123], or energy-based [124,125] codes for the topological analysis of the primitive path network of entanglements as abundantly encountered in densely packed systems of long polymer chains.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/ijms222212464/s1. The manuscript is available in interactive, 3-D format. With the exception of Figure 8, all panels corresponding to system configurations are also available as stand-alone, interactive 3-D pdf files.

**Author Contributions:** Conceptualization: N.C.K. and M.L.; methodology: K.F., N.C.K. and M.L.; software (Simu-D): M.H., D.M.-F., P.M.R., K.F., N.C.K. and M.L.; data curation: M.H., D.M.-F. and P.M.R.; visualization: M.H., D.M.-F. and N.C.K.; funding acquisition: K.F., N.C.K. and M.L.; writing— original draft preparation: M.H. and N.C.K.; writing—review and editing: D.M.-F., P.M.R., K.F. and M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by MINECO/FEDER (Ministerio de Economia y Competitividad, Fondo Europeo de Desarrollo Regional), gran<sup>t</sup> numbers "MAT2015-70478-P" and "RTI2018-097338-B-I00" and by UPM and Santander Bank, "Programa Propio UPM Santander".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Presented simulation trajectories and corresponding data are fully available upon request.

**Acknowledgments:** Very fruitful discussions with Doros N. Theodorou, Martin Kröger, Mukta Tripathy, and Robert S. Hoy are deeply appreciated. Authors deeply thank Clara Pedrosa, Jorge Rey, Javier Benito, Olia Bouzid, and Manuel Santiago for helpful interactions. The authors acknowledge support through projects "MAT2015-70478-P" and "RTI2018-097338-B-I00" of MINECO/FEDER (Ministerio de Economia y Competitividad, Fondo Europeo de Desarrollo Regional). M.H. and D.M.F. acknowledge financial support through "Programa Propio UPM Santander" of UPM and Santander Bank. The authors gratefully acknowledge the Universidad Politécnica de Madrid for providing computing resources on the Magerit Supercomputer through projects "p208", "r553", "s341", and "t736".

**Conflicts of Interest:** The authors declare no conflict of interest.
